(git:34ef472)
xc_thomas_fermi.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculate the Thomas-Fermi kinetic energy functional
10 !> \note
11 !> Order of derivatives is: LDA 0; 1; 2; 3;
12 !> LSD 0; a b; aa bb; aaa bbb;
13 !> \par History
14 !> JGH (26.02.2003) : OpenMP enabled
15 !> fawzi (04.2004) : adapted to the new xc interface
16 !> \author JGH (18.02.2002)
17 ! **************************************************************************************************
18 MODULE xc_thomas_fermi
19  USE cp_array_utils, ONLY: cp_3d_r_cp_type
20  USE kinds, ONLY: dp
21  USE xc_derivative_desc, ONLY: deriv_rho,&
22  deriv_rhoa,&
23  deriv_rhob
24  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
25  xc_dset_get_derivative
26  USE xc_derivative_types, ONLY: xc_derivative_get,&
27  xc_derivative_type
28  USE xc_functionals_utilities, ONLY: set_util
29  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
30  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
31  xc_rho_set_type
32 #include "../base/base_uses.f90"
33 
34  IMPLICIT NONE
35 
36  PRIVATE
37 
38 ! *** Global parameters ***
39 
40  REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
41  REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
42  f23 = 2.0_dp*f13, &
43  f43 = 4.0_dp*f13, &
44  f53 = 5.0_dp*f13
45 
46  PUBLIC :: thomas_fermi_info, thomas_fermi_lda_eval, thomas_fermi_lsd_eval
47 
48  REAL(KIND=dp) :: cf, flda, flsd
49  REAL(KIND=dp) :: eps_rho
50 
51  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_thomas_fermi'
52 
53 CONTAINS
54 
55 ! **************************************************************************************************
56 !> \brief ...
57 !> \param cutoff ...
58 ! **************************************************************************************************
59  SUBROUTINE thomas_fermi_init(cutoff)
60 
61  REAL(KIND=dp), INTENT(IN) :: cutoff
62 
63  eps_rho = cutoff
64  CALL set_util(cutoff)
65 
66  cf = 0.3_dp*(3.0_dp*pi*pi)**f23
67  flda = cf
68  flsd = flda*2.0_dp**f23
69 
70  END SUBROUTINE thomas_fermi_init
71 
72 ! **************************************************************************************************
73 !> \brief ...
74 !> \param lsd ...
75 !> \param reference ...
76 !> \param shortform ...
77 !> \param needs ...
78 !> \param max_deriv ...
79 ! **************************************************************************************************
80  SUBROUTINE thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
81  LOGICAL, INTENT(in) :: lsd
82  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
83  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
84  INTEGER, INTENT(out), OPTIONAL :: max_deriv
85 
86  IF (PRESENT(reference)) THEN
87  reference = "Thomas-Fermi kinetic energy functional: see Parr and Yang"
88  IF (.NOT. lsd) THEN
89  IF (len_trim(reference) + 6 < len(reference)) THEN
90  reference(len_trim(reference):len_trim(reference) + 6) = ' {LDA}'
91  END IF
92  END IF
93  END IF
94  IF (PRESENT(shortform)) THEN
95  shortform = "Thomas-Fermi kinetic energy functional"
96  IF (.NOT. lsd) THEN
97  IF (len_trim(shortform) + 6 < len(shortform)) THEN
98  shortform(len_trim(shortform):len_trim(shortform) + 6) = ' {LDA}'
99  END IF
100  END IF
101  END IF
102  IF (PRESENT(needs)) THEN
103  IF (lsd) THEN
104  needs%rho_spin = .true.
105  needs%rho_spin_1_3 = .true.
106  ELSE
107  needs%rho = .true.
108  needs%rho_1_3 = .true.
109  END IF
110  END IF
111  IF (PRESENT(max_deriv)) max_deriv = 3
112 
113  END SUBROUTINE thomas_fermi_info
114 
115 ! **************************************************************************************************
116 !> \brief ...
117 !> \param rho_set ...
118 !> \param deriv_set ...
119 !> \param order ...
120 ! **************************************************************************************************
121  SUBROUTINE thomas_fermi_lda_eval(rho_set, deriv_set, order)
122  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
123  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
124  INTEGER, INTENT(in) :: order
125 
126  CHARACTER(len=*), PARAMETER :: routinen = 'thomas_fermi_lda_eval'
127 
128  INTEGER :: handle, npoints
129  INTEGER, DIMENSION(2, 3) :: bo
130  REAL(kind=dp) :: epsilon_rho
131  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
132  POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
133  r13, rho
134  TYPE(xc_derivative_type), POINTER :: deriv
135 
136  CALL timeset(routinen, handle)
137 
138  CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
139  local_bounds=bo, rho_cutoff=epsilon_rho)
140  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
141  CALL thomas_fermi_init(epsilon_rho)
142 
143  IF (order >= 0) THEN
144  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
145  allocate_deriv=.true.)
146  CALL xc_derivative_get(deriv, deriv_data=e_0)
147 
148  CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
149  END IF
150  IF (order >= 1 .OR. order == -1) THEN
151  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
152  allocate_deriv=.true.)
153  CALL xc_derivative_get(deriv, deriv_data=e_rho)
154 
155  CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
156  END IF
157  IF (order >= 2 .OR. order == -2) THEN
158  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
159  allocate_deriv=.true.)
160  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
161 
162  CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
163  END IF
164  IF (order >= 3 .OR. order == -3) THEN
165  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
166  allocate_deriv=.true.)
167  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
168 
169  CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
170  END IF
171  IF (order > 3 .OR. order < -3) THEN
172  cpabort("derivatives bigger than 3 not implemented")
173  END IF
174  CALL timestop(handle)
175  END SUBROUTINE thomas_fermi_lda_eval
176 
177 ! **************************************************************************************************
178 !> \brief ...
179 !> \param rho_set ...
180 !> \param deriv_set ...
181 !> \param order ...
182 ! **************************************************************************************************
183  SUBROUTINE thomas_fermi_lsd_eval(rho_set, deriv_set, order)
184  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
185  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
186  INTEGER, INTENT(in) :: order
187 
188  CHARACTER(len=*), PARAMETER :: routinen = 'thomas_fermi_lsd_eval'
189  INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
190 
191  INTEGER :: handle, i, ispin, npoints
192  INTEGER, DIMENSION(2, 3) :: bo
193  REAL(kind=dp) :: epsilon_rho
194  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
195  POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
196  TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
197  TYPE(xc_derivative_type), POINTER :: deriv
198 
199  CALL timeset(routinen, handle)
200  NULLIFY (deriv)
201  DO i = 1, 2
202  NULLIFY (rho(i)%array, rho_1_3(i)%array)
203  END DO
204 
205  CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
206  rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
207  rhob=rho(2)%array, &
208  rho_cutoff=epsilon_rho, &
209  local_bounds=bo)
210  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
211  CALL thomas_fermi_init(epsilon_rho)
212 
213  DO ispin = 1, 2
214  IF (order >= 0) THEN
215  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
216  allocate_deriv=.true.)
217  CALL xc_derivative_get(deriv, deriv_data=e_0)
218 
219  CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
220  e_0, npoints)
221  END IF
222  IF (order >= 1 .OR. order == -1) THEN
223  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
224  allocate_deriv=.true.)
225  CALL xc_derivative_get(deriv, deriv_data=e_rho)
226 
227  CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
228  e_rho, npoints)
229  END IF
230  IF (order >= 2 .OR. order == -2) THEN
231  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
232  rho_spin_name(ispin)], allocate_deriv=.true.)
233  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
234 
235  CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
236  e_rho_rho, npoints)
237  END IF
238  IF (order >= 3 .OR. order == -3) THEN
239  deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
240  rho_spin_name(ispin), rho_spin_name(ispin)], &
241  allocate_deriv=.true.)
242  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
243 
244  CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
245  e_rho_rho_rho, npoints)
246  END IF
247  IF (order > 3 .OR. order < -3) THEN
248  cpabort("derivatives bigger than 3 not implemented")
249  END IF
250  END DO
251  CALL timestop(handle)
252  END SUBROUTINE thomas_fermi_lsd_eval
253 
254 ! **************************************************************************************************
255 !> \brief ...
256 !> \param rho ...
257 !> \param r13 ...
258 !> \param e_0 ...
259 !> \param npoints ...
260 ! **************************************************************************************************
261  SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
262 
263  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
264  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
265  INTEGER, INTENT(in) :: npoints
266 
267  INTEGER :: ip
268 
269 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
270 !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13)
271  DO ip = 1, npoints
272 
273  IF (rho(ip) > eps_rho) THEN
274 
275  e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
276 
277  END IF
278 
279  END DO
280 
281  END SUBROUTINE thomas_fermi_lda_0
282 
283 ! **************************************************************************************************
284 !> \brief ...
285 !> \param rho ...
286 !> \param r13 ...
287 !> \param e_rho ...
288 !> \param npoints ...
289 ! **************************************************************************************************
290  SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
291 
292  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
293  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
294  INTEGER, INTENT(in) :: npoints
295 
296  INTEGER :: ip
297  REAL(kind=dp) :: f
298 
299  f = f53*flda
300 
301 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
302 !$OMP SHARED(npoints,rho,eps_rho,e_rho,f,r13)
303  DO ip = 1, npoints
304 
305  IF (rho(ip) > eps_rho) THEN
306 
307  e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
308 
309  END IF
310 
311  END DO
312 
313  END SUBROUTINE thomas_fermi_lda_1
314 
315 ! **************************************************************************************************
316 !> \brief ...
317 !> \param rho ...
318 !> \param r13 ...
319 !> \param e_rho_rho ...
320 !> \param npoints ...
321 ! **************************************************************************************************
322  SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
323 
324  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
325  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
326  INTEGER, INTENT(in) :: npoints
327 
328  INTEGER :: ip
329  REAL(kind=dp) :: f
330 
331  f = f23*f53*flda
332 
333 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
334 !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,f,r13)
335  DO ip = 1, npoints
336 
337  IF (rho(ip) > eps_rho) THEN
338 
339  e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
340 
341  END IF
342 
343  END DO
344 
345  END SUBROUTINE thomas_fermi_lda_2
346 
347 ! **************************************************************************************************
348 !> \brief ...
349 !> \param rho ...
350 !> \param r13 ...
351 !> \param e_rho_rho_rho ...
352 !> \param npoints ...
353 ! **************************************************************************************************
354  SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
355 
356  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
357  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
358  INTEGER, INTENT(in) :: npoints
359 
360  INTEGER :: ip
361  REAL(kind=dp) :: f
362 
363  f = -f13*f23*f53*flda
364 
365 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
366 !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,f,r13)
367  DO ip = 1, npoints
368 
369  IF (rho(ip) > eps_rho) THEN
370 
371  e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
372 
373  END IF
374 
375  END DO
376 
377  END SUBROUTINE thomas_fermi_lda_3
378 
379 ! **************************************************************************************************
380 !> \brief ...
381 !> \param rhoa ...
382 !> \param r13a ...
383 !> \param e_0 ...
384 !> \param npoints ...
385 ! **************************************************************************************************
386  SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
387 
388  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
389  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
390  INTEGER, INTENT(in) :: npoints
391 
392  INTEGER :: ip
393 
394 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
395 !$OMP SHARED(npoints,rhoa,eps_rho,e_0,flsd,r13a)
396  DO ip = 1, npoints
397 
398  IF (rhoa(ip) > eps_rho) THEN
399  e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
400  END IF
401 
402  END DO
403 
404  END SUBROUTINE thomas_fermi_lsd_0
405 
406 ! **************************************************************************************************
407 !> \brief ...
408 !> \param rhoa ...
409 !> \param r13a ...
410 !> \param e_rho ...
411 !> \param npoints ...
412 ! **************************************************************************************************
413  SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
414 
415  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
416  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
417  INTEGER, INTENT(in) :: npoints
418 
419  INTEGER :: ip
420  REAL(kind=dp) :: f
421 
422  f = f53*flsd
423 
424 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
425 !$OMP SHARED(npoints,rhoa,eps_rho,e_rho,f,r13a)
426  DO ip = 1, npoints
427 
428  IF (rhoa(ip) > eps_rho) THEN
429  e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
430  END IF
431 
432  END DO
433 
434  END SUBROUTINE thomas_fermi_lsd_1
435 
436 ! **************************************************************************************************
437 !> \brief ...
438 !> \param rhoa ...
439 !> \param r13a ...
440 !> \param e_rho_rho ...
441 !> \param npoints ...
442 ! **************************************************************************************************
443  SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
444 
445  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
446  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
447  INTEGER, INTENT(in) :: npoints
448 
449  INTEGER :: ip
450  REAL(kind=dp) :: f
451 
452  f = f23*f53*flsd
453 
454 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
455 !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,r13a)
456 
457  DO ip = 1, npoints
458 
459  IF (rhoa(ip) > eps_rho) THEN
460  e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
461  END IF
462 
463  END DO
464 
465  END SUBROUTINE thomas_fermi_lsd_2
466 
467 ! **************************************************************************************************
468 !> \brief ...
469 !> \param rhoa ...
470 !> \param r13a ...
471 !> \param e_rho_rho_rho ...
472 !> \param npoints ...
473 ! **************************************************************************************************
474  SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
475 
476  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
477  REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
478  INTEGER, INTENT(in) :: npoints
479 
480  INTEGER :: ip
481  REAL(kind=dp) :: f
482 
483  f = -f13*f23*f53*flsd
484 
485 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
486 !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,f,r13a)
487  DO ip = 1, npoints
488 
489  IF (rhoa(ip) > eps_rho) THEN
490  e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
491  END IF
492 
493  END DO
494 
495  END SUBROUTINE thomas_fermi_lsd_3
496 
497 END MODULE xc_thomas_fermi
498 
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition: xc_derivative_desc.F:26
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition: xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition: xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition: xc_derivative_set_types.F:66
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition: xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition: xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition: xc_functionals_utilities.F:15
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition: xc_functionals_utilities.F:54
xc_rho_cflags_types
contains the structure
Definition: xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition: xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition: xc_rho_set_types.F:281
xc_thomas_fermi
Calculate the Thomas-Fermi kinetic energy functional.
Definition: xc_thomas_fermi.F:18
xc_thomas_fermi::thomas_fermi_lsd_eval
subroutine, public thomas_fermi_lsd_eval(rho_set, deriv_set, order)
...
Definition: xc_thomas_fermi.F:184
xc_thomas_fermi::thomas_fermi_lda_eval
subroutine, public thomas_fermi_lda_eval(rho_set, deriv_set, order)
...
Definition: xc_thomas_fermi.F:122
xc_thomas_fermi::thomas_fermi_info
subroutine, public thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
...
Definition: xc_thomas_fermi.F:81