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