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