(git:98357aa)
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! **************************************************************************************************
20 USE kinds, ONLY: dp
21 USE mathconstants, ONLY: pi
22 USE xc_derivative_desc, ONLY: deriv_rho,&
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
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 IF (order >= 0) THEN
143 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
144 allocate_deriv=.true.)
145 CALL xc_derivative_get(deriv, deriv_data=e_0)
146
147 CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
148 END IF
149 IF (order >= 1 .OR. order == -1) THEN
150 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
151 allocate_deriv=.true.)
152 CALL xc_derivative_get(deriv, deriv_data=e_rho)
153
154 CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
155 END IF
156 IF (order >= 2 .OR. order == -2) THEN
157 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
158 allocate_deriv=.true.)
159 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
160
161 CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
162 END IF
163 IF (order >= 3 .OR. order == -3) THEN
164 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
165 allocate_deriv=.true.)
166 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
167
168 CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
169 END IF
170 IF (order > 3 .OR. order < -3) THEN
171 cpabort("derivatives bigger than 3 not implemented")
172 END IF
173 CALL timestop(handle)
174 END SUBROUTINE thomas_fermi_lda_eval
175
176! **************************************************************************************************
177!> \brief ...
178!> \param rho_set ...
179!> \param deriv_set ...
180!> \param order ...
181! **************************************************************************************************
182 SUBROUTINE thomas_fermi_lsd_eval(rho_set, deriv_set, order)
183 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
184 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
185 INTEGER, INTENT(in) :: order
186
187 CHARACTER(len=*), PARAMETER :: routinen = 'thomas_fermi_lsd_eval'
188 INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
189
190 INTEGER :: handle, i, ispin, npoints
191 INTEGER, DIMENSION(2, 3) :: bo
192 REAL(kind=dp) :: epsilon_rho
193 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
194 POINTER :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
195 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: rho, rho_1_3
196 TYPE(xc_derivative_type), POINTER :: deriv
197
198 CALL timeset(routinen, handle)
199 NULLIFY (deriv)
200 DO i = 1, 2
201 NULLIFY (rho(i)%array, rho_1_3(i)%array)
202 END DO
203
204 CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
205 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
206 rhob=rho(2)%array, &
207 rho_cutoff=epsilon_rho, &
208 local_bounds=bo)
209 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
210 CALL thomas_fermi_init(epsilon_rho)
211
212 DO ispin = 1, 2
213 IF (order >= 0) THEN
214 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
215 allocate_deriv=.true.)
216 CALL xc_derivative_get(deriv, deriv_data=e_0)
217
218 CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
219 e_0, npoints)
220 END IF
221 IF (order >= 1 .OR. order == -1) THEN
222 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
223 allocate_deriv=.true.)
224 CALL xc_derivative_get(deriv, deriv_data=e_rho)
225
226 CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
227 e_rho, npoints)
228 END IF
229 IF (order >= 2 .OR. order == -2) THEN
230 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
231 rho_spin_name(ispin)], allocate_deriv=.true.)
232 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
233
234 CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
235 e_rho_rho, npoints)
236 END IF
237 IF (order >= 3 .OR. order == -3) THEN
238 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
239 rho_spin_name(ispin), rho_spin_name(ispin)], &
240 allocate_deriv=.true.)
241 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
242
243 CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
244 e_rho_rho_rho, npoints)
245 END IF
246 IF (order > 3 .OR. order < -3) THEN
247 cpabort("derivatives bigger than 3 not implemented")
248 END IF
249 END DO
250 CALL timestop(handle)
251 END SUBROUTINE thomas_fermi_lsd_eval
252
253! **************************************************************************************************
254!> \brief ...
255!> \param rho ...
256!> \param r13 ...
257!> \param e_0 ...
258!> \param npoints ...
259! **************************************************************************************************
260 SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
261
262 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
263 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
264 INTEGER, INTENT(in) :: npoints
265
266 INTEGER :: ip
267
268!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
269!$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13)
270 DO ip = 1, npoints
271
272 IF (rho(ip) > eps_rho) THEN
273
274 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
275
276 END IF
277
278 END DO
279
280 END SUBROUTINE thomas_fermi_lda_0
281
282! **************************************************************************************************
283!> \brief ...
284!> \param rho ...
285!> \param r13 ...
286!> \param e_rho ...
287!> \param npoints ...
288! **************************************************************************************************
289 SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
290
291 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
292 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
293 INTEGER, INTENT(in) :: npoints
294
295 INTEGER :: ip
296 REAL(kind=dp) :: f
297
298 f = f53*flda
299
300!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
301!$OMP SHARED(npoints,rho,eps_rho,e_rho,f,r13)
302 DO ip = 1, npoints
303
304 IF (rho(ip) > eps_rho) THEN
305
306 e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
307
308 END IF
309
310 END DO
311
312 END SUBROUTINE thomas_fermi_lda_1
313
314! **************************************************************************************************
315!> \brief ...
316!> \param rho ...
317!> \param r13 ...
318!> \param e_rho_rho ...
319!> \param npoints ...
320! **************************************************************************************************
321 SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
322
323 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
324 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
325 INTEGER, INTENT(in) :: npoints
326
327 INTEGER :: ip
328 REAL(kind=dp) :: f
329
330 f = f23*f53*flda
331
332!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
333!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,f,r13)
334 DO ip = 1, npoints
335
336 IF (rho(ip) > eps_rho) THEN
337
338 e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
339
340 END IF
341
342 END DO
343
344 END SUBROUTINE thomas_fermi_lda_2
345
346! **************************************************************************************************
347!> \brief ...
348!> \param rho ...
349!> \param r13 ...
350!> \param e_rho_rho_rho ...
351!> \param npoints ...
352! **************************************************************************************************
353 SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
354
355 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13
356 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
357 INTEGER, INTENT(in) :: npoints
358
359 INTEGER :: ip
360 REAL(kind=dp) :: f
361
362 f = -f13*f23*f53*flda
363
364!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
365!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,f,r13)
366 DO ip = 1, npoints
367
368 IF (rho(ip) > eps_rho) THEN
369
370 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
371
372 END IF
373
374 END DO
375
376 END SUBROUTINE thomas_fermi_lda_3
377
378! **************************************************************************************************
379!> \brief ...
380!> \param rhoa ...
381!> \param r13a ...
382!> \param e_0 ...
383!> \param npoints ...
384! **************************************************************************************************
385 SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
386
387 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
388 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
389 INTEGER, INTENT(in) :: npoints
390
391 INTEGER :: ip
392
393!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
394!$OMP SHARED(npoints,rhoa,eps_rho,e_0,flsd,r13a)
395 DO ip = 1, npoints
396
397 IF (rhoa(ip) > eps_rho) THEN
398 e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
399 END IF
400
401 END DO
402
403 END SUBROUTINE thomas_fermi_lsd_0
404
405! **************************************************************************************************
406!> \brief ...
407!> \param rhoa ...
408!> \param r13a ...
409!> \param e_rho ...
410!> \param npoints ...
411! **************************************************************************************************
412 SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
413
414 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
415 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho
416 INTEGER, INTENT(in) :: npoints
417
418 INTEGER :: ip
419 REAL(kind=dp) :: f
420
421 f = f53*flsd
422
423!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
424!$OMP SHARED(npoints,rhoa,eps_rho,e_rho,f,r13a)
425 DO ip = 1, npoints
426
427 IF (rhoa(ip) > eps_rho) THEN
428 e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
429 END IF
430
431 END DO
432
433 END SUBROUTINE thomas_fermi_lsd_1
434
435! **************************************************************************************************
436!> \brief ...
437!> \param rhoa ...
438!> \param r13a ...
439!> \param e_rho_rho ...
440!> \param npoints ...
441! **************************************************************************************************
442 SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
443
444 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
445 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho
446 INTEGER, INTENT(in) :: npoints
447
448 INTEGER :: ip
449 REAL(kind=dp) :: f
450
451 f = f23*f53*flsd
452
453!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
454!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,r13a)
455
456 DO ip = 1, npoints
457
458 IF (rhoa(ip) > eps_rho) THEN
459 e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
460 END IF
461
462 END DO
463
464 END SUBROUTINE thomas_fermi_lsd_2
465
466! **************************************************************************************************
467!> \brief ...
468!> \param rhoa ...
469!> \param r13a ...
470!> \param e_rho_rho_rho ...
471!> \param npoints ...
472! **************************************************************************************************
473 SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
474
475 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a
476 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho
477 INTEGER, INTENT(in) :: npoints
478
479 INTEGER :: ip
480 REAL(kind=dp) :: f
481
482 f = -f13*f23*f53*flsd
483
484!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
485!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,f,r13a)
486 DO ip = 1, npoints
487
488 IF (rhoa(ip) > eps_rho) THEN
489 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
490 END IF
491
492 END DO
493
494 END SUBROUTINE thomas_fermi_lsd_3
495
496END MODULE xc_thomas_fermi
497
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Utility routines for the functional calculations.
subroutine, public set_util(cutoff)
...
contains the structure
contains the structure
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
Calculate the Thomas-Fermi kinetic energy functional.
subroutine, public thomas_fermi_lsd_eval(rho_set, deriv_set, order)
...
subroutine, public thomas_fermi_lda_eval(rho_set, deriv_set, order)
...
subroutine, public thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
...
represent a pointer to a contiguous 3d array
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation