(git:374b731)
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-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! **************************************************************************************************
18MODULE 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
53CONTAINS
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
497END 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
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:49
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