(git:98357aa)
Loading...
Searching...
No Matches
xc_tfw.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!> plus the von Weizsaecker term
11!> \par History
12!> JGH (26.02.2003) : OpenMP enabled
13!> fawzi (04.2004) : adapted to the new xc interface
14!> \author JGH (18.02.2002)
15! **************************************************************************************************
16MODULE xc_tfw
17 USE cp_array_utils, ONLY: cp_3d_r_cp_type
18 USE kinds, ONLY: dp
19 USE mathconstants, ONLY: pi
20 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
21 deriv_norm_drhoa,&
22 deriv_norm_drhob,&
23 deriv_rho,&
24 deriv_rhoa,&
25 deriv_rhob
26 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
27 xc_dset_get_derivative
28 USE xc_derivative_types, ONLY: xc_derivative_get,&
29 xc_derivative_type
30 USE xc_functionals_utilities, ONLY: set_util
31 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
32 USE xc_rho_set_types, ONLY: xc_rho_set_get,&
33 xc_rho_set_type
34#include "../base/base_uses.f90"
35
36 IMPLICIT NONE
37
38 PRIVATE
39
40! *** Global parameters ***
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 :: tfw_lda_info, tfw_lda_eval, tfw_lsd_info, tfw_lsd_eval
47
48 REAL(KIND=dp) :: cf, flda, flsd, fvw
49 REAL(KIND=dp) :: eps_rho
50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tfw'
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief ...
56!> \param cutoff ...
57! **************************************************************************************************
58 SUBROUTINE tfw_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 fvw = 1.0_dp/72.0_dp
69
70 END SUBROUTINE tfw_init
71
72! **************************************************************************************************
73!> \brief ...
74!> \param reference ...
75!> \param shortform ...
76!> \param needs ...
77!> \param max_deriv ...
78! **************************************************************************************************
79 SUBROUTINE tfw_lda_info(reference, shortform, needs, max_deriv)
80 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
81 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
82 INTEGER, INTENT(out), OPTIONAL :: max_deriv
83
84 IF (PRESENT(reference)) THEN
85 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}"
86 END IF
87 IF (PRESENT(shortform)) THEN
88 shortform = "TF+vW kinetic energy functional {LDA}"
89 END IF
90 IF (PRESENT(needs)) THEN
91 needs%rho = .true.
92 needs%rho_1_3 = .true.
93 needs%norm_drho = .true.
94 END IF
95 IF (PRESENT(max_deriv)) max_deriv = 3
96
97 END SUBROUTINE tfw_lda_info
98
99! **************************************************************************************************
100!> \brief ...
101!> \param reference ...
102!> \param shortform ...
103!> \param needs ...
104!> \param max_deriv ...
105! **************************************************************************************************
106 SUBROUTINE tfw_lsd_info(reference, shortform, needs, max_deriv)
107 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
108 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
109 INTEGER, INTENT(out), OPTIONAL :: max_deriv
110
111 IF (PRESENT(reference)) THEN
112 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional"
113 END IF
114 IF (PRESENT(shortform)) THEN
115 shortform = "TF+vW kinetic energy functional"
116 END IF
117 IF (PRESENT(needs)) THEN
118 needs%rho_spin = .true.
119 needs%rho_spin_1_3 = .true.
120 needs%norm_drho = .true.
121 END IF
122 IF (PRESENT(max_deriv)) max_deriv = 3
123
124 END SUBROUTINE tfw_lsd_info
125
126! **************************************************************************************************
127!> \brief ...
128!> \param rho_set ...
129!> \param deriv_set ...
130!> \param order ...
131! **************************************************************************************************
132 SUBROUTINE tfw_lda_eval(rho_set, deriv_set, order)
133 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
134 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
135 INTEGER, INTENT(in) :: order
136
137 CHARACTER(len=*), PARAMETER :: routinen = 'tfw_lda_eval'
138
139 INTEGER :: handle, npoints
140 INTEGER, DIMENSION(2, 3) :: bo
141 REAL(kind=dp) :: epsilon_rho
142 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
143 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
144 e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, e_rho_rho_rho, grho, &
145 r13, rho
146 TYPE(xc_derivative_type), POINTER :: deriv
147
148 CALL timeset(routinen, handle)
149
150 CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
151 norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
152 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
153 CALL tfw_init(epsilon_rho)
154
155 ALLOCATE (s(npoints))
156 CALL calc_s(rho, grho, s, npoints)
157
158 IF (order >= 0) THEN
159 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
160 allocate_deriv=.true.)
161 CALL xc_derivative_get(deriv, deriv_data=e_0)
162
163 CALL tfw_u_0(rho, r13, s, e_0, npoints)
164 END IF
165 IF (order >= 1 .OR. order == -1) THEN
166 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
167 allocate_deriv=.true.)
168 CALL xc_derivative_get(deriv, deriv_data=e_rho)
169 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
170 allocate_deriv=.true.)
171 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
172
173 CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
174 END IF
175 IF (order >= 2 .OR. order == -2) THEN
176 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
177 allocate_deriv=.true.)
178 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
179 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
180 allocate_deriv=.true.)
181 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
182 deriv => xc_dset_get_derivative(deriv_set, &
183 [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
184 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
185
186 CALL tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, &
187 e_ndrho_ndrho, npoints)
188 END IF
189 IF (order >= 3 .OR. order == -3) THEN
190 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
191 allocate_deriv=.true.)
192 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
193 deriv => xc_dset_get_derivative(deriv_set, &
194 [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.true.)
195 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
196 deriv => xc_dset_get_derivative(deriv_set, &
197 [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
198 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
199
200 CALL tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
201 e_rho_ndrho_ndrho, npoints)
202 END IF
203 IF (order > 3 .OR. order < -3) THEN
204 cpabort("derivatives bigger than 3 not implemented")
205 END IF
206
207 DEALLOCATE (s)
208 CALL timestop(handle)
209 END SUBROUTINE tfw_lda_eval
210
211! **************************************************************************************************
212!> \brief ...
213!> \param rho ...
214!> \param grho ...
215!> \param s ...
216!> \param npoints ...
217! **************************************************************************************************
218 SUBROUTINE calc_s(rho, grho, s, npoints)
219 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho, grho
220 REAL(kind=dp), DIMENSION(*), INTENT(out) :: s
221 INTEGER, INTENT(in) :: npoints
222
223 INTEGER :: ip
224
225!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
226!$OMP SHARED(npoints,rho,eps_rho,s,grho)
227 DO ip = 1, npoints
228 IF (rho(ip) < eps_rho) THEN
229 s(ip) = 0.0_dp
230 ELSE
231 s(ip) = grho(ip)*grho(ip)/rho(ip)
232 END IF
233 END DO
234 END SUBROUTINE calc_s
235
236! **************************************************************************************************
237!> \brief ...
238!> \param rho_set ...
239!> \param deriv_set ...
240!> \param order ...
241! **************************************************************************************************
242 SUBROUTINE tfw_lsd_eval(rho_set, deriv_set, order)
243 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
244 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
245 INTEGER, INTENT(in) :: order
246
247 CHARACTER(len=*), PARAMETER :: routinen = 'tfw_lsd_eval'
248 INTEGER, DIMENSION(2), PARAMETER :: &
249 norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
250 rho_spin_name = [deriv_rhoa, deriv_rhob]
251
252 INTEGER :: handle, i, ispin, npoints
253 INTEGER, DIMENSION(2, 3) :: bo
254 REAL(kind=dp) :: epsilon_rho
255 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s
256 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
257 POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, &
258 e_rho_ndrho, e_rho_ndrho_ndrho, &
259 e_rho_rho, e_rho_rho_ndrho, &
260 e_rho_rho_rho
261 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: norm_drho, rho, rho_1_3
262 TYPE(xc_derivative_type), POINTER :: deriv
263
264 CALL timeset(routinen, handle)
265 NULLIFY (deriv)
266 DO i = 1, 2
267 NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
268 END DO
269
270 CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
271 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
272 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
273 norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
274 local_bounds=bo)
275 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
276 CALL tfw_init(epsilon_rho)
277
278 ALLOCATE (s(npoints))
279
280 DO ispin = 1, 2
281 CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints)
282
283 IF (order >= 0) THEN
284 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
285 allocate_deriv=.true.)
286 CALL xc_derivative_get(deriv, deriv_data=e_0)
287
288 CALL tfw_p_0(rho(ispin)%array, &
289 rho_1_3(ispin)%array, s, e_0, npoints)
290 END IF
291 IF (order >= 1 .OR. order == -1) THEN
292 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
293 allocate_deriv=.true.)
294 CALL xc_derivative_get(deriv, deriv_data=e_rho)
295 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
296 allocate_deriv=.true.)
297 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
298
299 CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, &
300 rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints)
301 END IF
302 IF (order >= 2 .OR. order == -2) THEN
303 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
304 rho_spin_name(ispin)], allocate_deriv=.true.)
305 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
306 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
307 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
308 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
309 deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
310 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
311 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
312
313 CALL tfw_p_2(rho(ispin)%array, norm_drho(ispin)%array, &
314 rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho, &
315 e_ndrho_ndrho, npoints)
316 END IF
317 IF (order >= 3 .OR. order == -3) THEN
318 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
319 rho_spin_name(ispin), rho_spin_name(ispin)], &
320 allocate_deriv=.true.)
321 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
322 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
323 rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
324 allocate_deriv=.true.)
325 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
326 deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
327 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
328 allocate_deriv=.true.)
329 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
330
331 CALL tfw_p_3(rho(ispin)%array, norm_drho(ispin)%array, &
332 rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho, &
333 e_rho_ndrho_ndrho, npoints)
334 END IF
335 IF (order > 3 .OR. order < -3) THEN
336 cpabort("derivatives bigger than 3 not implemented")
337 END IF
338 END DO
339
340 DEALLOCATE (s)
341 CALL timestop(handle)
342 END SUBROUTINE tfw_lsd_eval
343
344! **************************************************************************************************
345!> \brief ...
346!> \param rho ...
347!> \param r13 ...
348!> \param s ...
349!> \param e_0 ...
350!> \param npoints ...
351! **************************************************************************************************
352 SUBROUTINE tfw_u_0(rho, r13, s, e_0, npoints)
353
354 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s
355 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
356 INTEGER, INTENT(in) :: npoints
357
358 INTEGER :: ip
359
360!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
361!$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13,s,fvw)
362 DO ip = 1, npoints
363
364 IF (rho(ip) > eps_rho) THEN
365
366 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw*s(ip)
367
368 END IF
369
370 END DO
371
372 END SUBROUTINE tfw_u_0
373
374! **************************************************************************************************
375!> \brief ...
376!> \param rho ...
377!> \param grho ...
378!> \param r13 ...
379!> \param s ...
380!> \param e_rho ...
381!> \param e_ndrho ...
382!> \param npoints ...
383! **************************************************************************************************
384 SUBROUTINE tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
385
386 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
387 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
388 INTEGER, INTENT(in) :: npoints
389
390 INTEGER :: ip
391 REAL(kind=dp) :: f
392
393 f = f53*flda
394
395!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
396!$OMP SHARED(npoints,rho,eps_rho,e_rho,e_ndrho,grho,s,r13,f,fvw)
397 DO ip = 1, npoints
398
399 IF (rho(ip) > eps_rho) THEN
400
401 e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip) - fvw*s(ip)/rho(ip)
402 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grho(ip)/rho(ip)
403
404 END IF
405
406 END DO
407
408 END SUBROUTINE tfw_u_1
409
410! **************************************************************************************************
411!> \brief ...
412!> \param rho ...
413!> \param grho ...
414!> \param r13 ...
415!> \param s ...
416!> \param e_rho_rho ...
417!> \param e_rho_ndrho ...
418!> \param e_ndrho_ndrho ...
419!> \param npoints ...
420! **************************************************************************************************
421 SUBROUTINE tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
422 npoints)
423
424 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
425 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
426 INTEGER, INTENT(in) :: npoints
427
428 INTEGER :: ip
429 REAL(kind=dp) :: f
430
431 f = f23*f53*flda
432
433!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
434!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho,grho,f,fvw)
435 DO ip = 1, npoints
436
437 IF (rho(ip) > eps_rho) THEN
438
439 e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip) + 2.0_dp*fvw*s(ip)/(rho(ip)*rho(ip))
440 e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip))
441 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rho(ip)
442
443 END IF
444
445 END DO
446
447 END SUBROUTINE tfw_u_2
448
449! **************************************************************************************************
450!> \brief ...
451!> \param rho ...
452!> \param grho ...
453!> \param r13 ...
454!> \param s ...
455!> \param e_rho_rho_rho ...
456!> \param e_rho_rho_ndrho ...
457!> \param e_rho_ndrho_ndrho ...
458!> \param npoints ...
459! **************************************************************************************************
460 SUBROUTINE tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
461 e_rho_ndrho_ndrho, npoints)
462
463 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s
464 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
465 e_rho_ndrho_ndrho
466 INTEGER, INTENT(in) :: npoints
467
468 INTEGER :: ip
469 REAL(kind=dp) :: f
470
471 f = -f13*f23*f53*flda
472
473!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
474!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,r13,s,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw)
475 DO ip = 1, npoints
476
477 IF (rho(ip) > eps_rho) THEN
478
479 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip)) &
480 - 6.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)*rho(ip))
481 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
482 + 4.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)*rho(ip))
483 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
484 - 2.0_dp*fvw/(rho(ip)*rho(ip))
485 END IF
486
487 END DO
488
489 END SUBROUTINE tfw_u_3
490
491! **************************************************************************************************
492!> \brief ...
493!> \param rhoa ...
494!> \param r13a ...
495!> \param sa ...
496!> \param e_0 ...
497!> \param npoints ...
498! **************************************************************************************************
499 SUBROUTINE tfw_p_0(rhoa, r13a, sa, e_0, npoints)
500
501 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a, sa
502 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
503 INTEGER, INTENT(in) :: npoints
504
505 INTEGER :: ip
506
507!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
508!$OMP SHARED(npoints, rhoa,eps_rho,e_0,r13a,sa,flsd,fvw)
509 DO ip = 1, npoints
510
511 IF (rhoa(ip) > eps_rho) THEN
512 e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip) + fvw*sa(ip)
513 END IF
514
515 END DO
516
517 END SUBROUTINE tfw_p_0
518
519! **************************************************************************************************
520!> \brief ...
521!> \param rhoa ...
522!> \param grhoa ...
523!> \param r13a ...
524!> \param sa ...
525!> \param e_rho ...
526!> \param e_ndrho ...
527!> \param npoints ...
528! **************************************************************************************************
529 SUBROUTINE tfw_p_1(rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints)
530
531 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
532 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
533 INTEGER, INTENT(in) :: npoints
534
535 INTEGER :: ip
536 REAL(kind=dp) :: f
537
538 f = f53*flsd
539
540!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
541!$OMP SHARED(npoints,rhoa,eps_rho,r13a,sa,fvw,grhoa,e_rho,e_ndrho,f)
542 DO ip = 1, npoints
543
544 IF (rhoa(ip) > eps_rho) THEN
545 e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip) - fvw*sa(ip)/rhoa(ip)
546 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grhoa(ip)/rhoa(ip)
547 END IF
548
549 END DO
550
551 END SUBROUTINE tfw_p_1
552
553! **************************************************************************************************
554!> \brief ...
555!> \param rhoa ...
556!> \param grhoa ...
557!> \param r13a ...
558!> \param sa ...
559!> \param e_rho_rho ...
560!> \param e_rho_ndrho ...
561!> \param e_ndrho_ndrho ...
562!> \param npoints ...
563! **************************************************************************************************
564 SUBROUTINE tfw_p_2(rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho, &
565 e_ndrho_ndrho, npoints)
566
567 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
568 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
569 INTEGER, INTENT(in) :: npoints
570
571 INTEGER :: ip
572 REAL(kind=dp) :: f
573
574 f = f23*f53*flsd
575
576!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
577!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,fvw,r13a,sa,e_rho_ndrho,e_ndrho_ndrho)
578 DO ip = 1, npoints
579
580 IF (rhoa(ip) > eps_rho) THEN
581 e_rho_rho(ip) = e_rho_rho(ip) &
582 + f/r13a(ip) + 2.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip))
583 e_rho_ndrho(ip) = e_rho_ndrho(ip) &
584 - 2.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip))
585 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rhoa(ip)
586 END IF
587
588 END DO
589
590 END SUBROUTINE tfw_p_2
591
592! **************************************************************************************************
593!> \brief ...
594!> \param rhoa ...
595!> \param grhoa ...
596!> \param r13a ...
597!> \param sa ...
598!> \param e_rho_rho_rho ...
599!> \param e_rho_rho_ndrho ...
600!> \param e_rho_ndrho_ndrho ...
601!> \param npoints ...
602! **************************************************************************************************
603 SUBROUTINE tfw_p_3(rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho, &
604 e_rho_ndrho_ndrho, npoints)
605
606 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa
607 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
608 e_rho_ndrho_ndrho
609 INTEGER, INTENT(in) :: npoints
610
611 INTEGER :: ip
612 REAL(kind=dp) :: f
613
614 f = -f13*f23*f53*flsd
615
616!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
617!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw,sa,grhoa)
618 DO ip = 1, npoints
619
620 IF (rhoa(ip) > eps_rho) THEN
621 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
622 + f/(r13a(ip)*rhoa(ip)) &
623 - 6.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
624 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
625 + 4.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
626 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
627 - 2.0_dp*fvw/(rhoa(ip)*rhoa(ip))
628 END IF
629
630 END DO
631
632 END SUBROUTINE tfw_p_3
633
634END MODULE xc_tfw
635
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_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
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_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
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_tfw
Calculate the Thomas-Fermi kinetic energy functional plus the von Weizsaecker term.
Definition xc_tfw.F:16
xc_tfw::tfw_lda_info
subroutine, public tfw_lda_info(reference, shortform, needs, max_deriv)
...
Definition xc_tfw.F:80
xc_tfw::tfw_lda_eval
subroutine, public tfw_lda_eval(rho_set, deriv_set, order)
...
Definition xc_tfw.F:133
xc_tfw::tfw_lsd_info
subroutine, public tfw_lsd_info(reference, shortform, needs, max_deriv)
...
Definition xc_tfw.F:107
xc_tfw::tfw_lsd_eval
subroutine, public tfw_lsd_eval(rho_set, deriv_set, order)
...
Definition xc_tfw.F:243
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