(git:db0abe2)
Loading...
Searching...
No Matches
xc_perdew_wang.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 Perdew-Wang correlation potential and
10!> energy density and ist derivatives with respect to
11!> the spin-up and spin-down densities up to 3rd order.
12!> \par History
13!> 18-MAR-2002, TCH, working version
14!> fawzi (04.2004) : adapted to the new xc interface
15!> \see functionals_utilities
16! **************************************************************************************************
17MODULE xc_perdew_wang
18
19 USE kinds, ONLY: dp
20 USE pw_types, ONLY: pw_r3d_rs_type
21 USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
22 xc_dset_get_derivative
23 USE xc_derivative_types, ONLY: xc_derivative_get, &
24 xc_derivative_type
25 USE xc_functionals_utilities, ONLY: calc_fx, &
26 calc_rs, &
27 calc_z, &
28 set_util
29 USE xc_input_constants, ONLY: pw_dmc, &
30 pw_orig, &
31 pw_vmc
32 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
33 USE xc_rho_set_types, ONLY: xc_rho_set_get, &
34 xc_rho_set_type
35 USE xc_derivative_desc, ONLY: deriv_rho, &
36 deriv_rhoa, &
37 deriv_rhob
38#include "../base/base_uses.f90"
39
40 IMPLICIT NONE
41 PRIVATE
42
43 REAL(KIND=dp), DIMENSION(-1:1) :: a, a1, b1, b2, b3, b4
44 REAL(KIND=dp), DIMENSION(-1:1) :: c0, c1, c2, c3
45 REAL(KIND=dp), DIMENSION(-1:1) :: d0, d1
46 REAL(KIND=dp) :: eps_rho
47 REAL(KIND=dp), PARAMETER :: &
48 epsilon = 5.e-13_dp, &
49 fpp = 0.584822362263464620726223866376013788782_dp ! d^2f(0)/dz^2
50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_perdew_wang'
51
52 PUBLIC :: perdew_wang_info, perdew_wang_lda_eval, perdew_wang_lsd_eval, perdew_wang_fxc_calc
53
54CONTAINS
55
56! **************************************************************************************************
57!> \brief Return some info on the functionals.
58!> \param method ...
59!> \param lsd ...
60!> \param reference full reference
61!> \param shortform short reference
62!> \param needs ...
63!> \param max_deriv ...
64!> \param scale ...
65! **************************************************************************************************
66 SUBROUTINE perdew_wang_info(method, lsd, reference, shortform, needs, &
67 max_deriv, scale)
68 INTEGER, INTENT(in) :: method
69 LOGICAL, INTENT(in) :: lsd
70 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
71 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
72 INTEGER, INTENT(out), OPTIONAL :: max_deriv
73 REAL(kind=dp), INTENT(in) :: scale
74
75 CHARACTER(len=3) :: p_string
76
77 SELECT CASE (method)
78 CASE DEFAULT
79 cpabort("Unsupported parametrization")
80 CASE (pw_orig)
81 p_string = 'PWO'
82 CASE (pw_dmc)
83 p_string = 'DMC'
84 CASE (pw_vmc)
85 p_string = 'VMC'
86 END SELECT
87
88 IF (PRESENT(reference)) THEN
89 reference = "J. P. Perdew and Yue Wang," &
90 //" Phys. Rev. B 45, 13244 (1992)" &
91 //"["//trim(p_string)//"]"
92 IF (scale /= 1._dp) THEN
93 WRITE (reference(len_trim(reference) + 1:len(reference)), "('s=',f5.3)") &
94 scale
95 END IF
96 IF (.NOT. lsd) THEN
97 IF (len_trim(reference) + 6 < len(reference)) THEN
98 reference(len_trim(reference) + 1:len_trim(reference) + 7) = ' {LDA}'
99 END IF
100 END IF
101 END IF
102 IF (PRESENT(shortform)) THEN
103 shortform = "J. P. Perdew et al., PRB 45, 13244 (1992)" &
104 //"["//trim(p_string)//"]"
105 IF (scale /= 1._dp) THEN
106 WRITE (shortform(len_trim(shortform) + 1:len(shortform)), "('s=',f5.3)") &
107 scale
108 END IF
109 IF (.NOT. lsd) THEN
110 IF (len_trim(shortform) + 6 < len(shortform)) THEN
111 shortform(len_trim(shortform) + 1:len_trim(shortform) + 7) = ' {LDA}'
112 END IF
113 END IF
114 END IF
115 IF (PRESENT(needs)) THEN
116 IF (lsd) THEN
117 needs%rho_spin = .true.
118 ELSE
119 needs%rho = .true.
120 END IF
121 END IF
122 IF (PRESENT(max_deriv)) max_deriv = 3
123
124 END SUBROUTINE perdew_wang_info
125
126! **************************************************************************************************
127!> \brief Initializes the functionals
128!> \param method name of the method used for parameters
129!> \param cutoff the cutoff density
130! **************************************************************************************************
131 SUBROUTINE perdew_wang_init(method, cutoff)
132
133 INTEGER, INTENT(IN) :: method
134 REAL(kind=dp), INTENT(IN) :: cutoff
135
136 INTEGER :: k
137
138 CALL set_util(cutoff)
139
140 eps_rho = cutoff
141
142 ! values for -ac are the same for all methods
143 a(-1) = 0.016887_dp
144 a1(-1) = 0.11125_dp
145 b1(-1) = 10.357_dp
146 b2(-1) = 3.6231_dp
147 b3(-1) = 0.88026_dp
148 b4(-1) = 0.49671_dp
149
150 SELECT CASE (method)
151
152 CASE DEFAULT
153 cpabort("Unknown method")
154
155 CASE (pw_orig)
156 a(0) = 0.031091_dp; a(1) = 0.015545_dp
157 a1(0) = 0.21370_dp; a1(1) = 0.20548_dp
158 b1(0) = 7.5957_dp; b1(1) = 14.1189_dp
159 b2(0) = 3.5876_dp; b2(1) = 6.1977_dp
160 b3(0) = 1.6382_dp; b3(1) = 3.3662_dp
161 b4(0) = 0.49294_dp; b4(1) = 0.62517_dp
162
163 CASE (pw_dmc)
164 a(0) = 0.031091_dp; a(1) = 0.015545_dp
165 a1(0) = 0.026481_dp; a1(1) = 0.022465_dp
166 b1(0) = 7.5957_dp; b1(1) = 14.1189_dp
167 b2(0) = 3.5876_dp; b2(1) = 6.1977_dp
168 b3(0) = -0.46647_dp; b3(1) = -0.56043_dp
169 b4(0) = 0.13354_dp; b4(1) = 0.11313_dp
170
171 CASE (pw_vmc)
172 a(0) = 0.031091_dp; a(1) = 0.015545_dp
173 a1(0) = -0.002257_dp; a1(1) = -0.009797_dp
174 b1(0) = 7.5957_dp; b1(1) = 14.1189_dp
175 b2(0) = 3.5876_dp; b2(1) = 6.1977_dp
176 b3(0) = -0.52669_dp; b3(1) = -0.91381_dp
177 b4(0) = 0.03755_dp; b4(1) = 0.01538_dp
178
179 END SELECT
180
181 DO k = -1, 1, 1
182 c0(k) = a(k)
183 c1(k) = -2.0_dp*c0(k)*log(2.0_dp*a(k)*b1(k))
184 c2(k) = a(k)*a1(k)
185 c3(k) = -2.0_dp*a(k)*(a1(k)*log(2.0_dp*a(k)*b1(k)) &
186 - (b2(k)/b1(k))**2 + (b3(k)/b1(k)))
187 d0(k) = a1(k)/b4(k)
188 d1(k) = a1(k)*b3(k)/(b4(k)**2)
189 END DO
190
191 END SUBROUTINE perdew_wang_init
192
193! **************************************************************************************************
194!> \brief Calculate the correlation energy and its derivatives
195!> wrt to rho (the electron density) up to 3rd order. This
196!> is the LDA version of the Perdew-Wang correlation energy
197!> If no order argument is given, then the routine calculates
198!> just the energy.
199!> \param method ...
200!> \param rho_set ...
201!> \param deriv_set ...
202!> \param order maximum order of derivatives to calculate
203!> order must lie between 0 and 3.
204!> \param scale ...
205! **************************************************************************************************
206 SUBROUTINE perdew_wang_lda_eval(method, rho_set, deriv_set, order, scale)
207
208 INTEGER, INTENT(in) :: method
209 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
210 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
211 INTEGER, INTENT(in) :: order
212 REAL(kind=dp), INTENT(in) :: scale
213
214 CHARACTER(len=*), PARAMETER :: routinen = 'perdew_wang_lda_eval'
215
216 INTEGER :: npoints, timer_handle
217 INTEGER, DIMENSION(2, 3) :: bo
218 REAL(kind=dp) :: rho_cutoff
219 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_rho, e_rho_rho, &
220 e_rho_rho_rho, rho
221 TYPE(xc_derivative_type), POINTER :: deriv
222
223 CALL timeset(routinen, timer_handle)
224 NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
225 CALL xc_rho_set_get(rho_set, rho=rho, &
226 local_bounds=bo, rho_cutoff=rho_cutoff)
227 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
228
229 CALL perdew_wang_init(method, rho_cutoff)
230
231 dummy => rho
232
233 e_0 => dummy
234 e_rho => dummy
235 e_rho_rho => dummy
236 e_rho_rho_rho => dummy
237
238 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
239 allocate_deriv=.true.)
240 CALL xc_derivative_get(deriv, deriv_data=e_0)
241 IF (order >= 1) THEN
242 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
243 allocate_deriv=.true.)
244 CALL xc_derivative_get(deriv, deriv_data=e_rho)
245 END IF
246 IF (order >= 2) THEN
247 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
248 allocate_deriv=.true.)
249 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
250 END IF
251 IF (order >= 3) THEN
252 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
253 allocate_deriv=.true.)
254 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
255 END IF
256 IF (order > 3) THEN
257 cpabort("derivatives bigger than 3 not implemented")
258 END IF
259
260 CALL perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
261 npoints, order, scale)
262
263 CALL timestop(timer_handle)
264
265 END SUBROUTINE perdew_wang_lda_eval
266
267! **************************************************************************************************
268!> \brief ...
269!> \param rho ...
270!> \param e_0 ...
271!> \param e_rho ...
272!> \param e_rho_rho ...
273!> \param e_rho_rho_rho ...
274!> \param npoints ...
275!> \param order ...
276!> \param scale ...
277! **************************************************************************************************
278 SUBROUTINE perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, scale)
279 !FM low level calc routine
280 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho
281 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
282 INTEGER, INTENT(in) :: npoints, order
283 REAL(kind=dp), INTENT(in) :: scale
284
285 INTEGER :: k
286 REAL(kind=dp), DIMENSION(0:3) :: ed
287
288!$OMP PARALLEL DO PRIVATE (k, ed) DEFAULT(NONE)&
289!$OMP SHARED(npoints,rho,eps_rho,scale,e_0,e_rho,e_rho_rho,e_rho_rho_rho,order)
290 DO k = 1, npoints
291
292 IF (rho(k) > eps_rho) THEN
293!! order_ is positive as it must be in this case:
294!! ec(:,2) needs ed(:,1) for example
295 CALL pw_lda_ed_loc(rho(k), ed, order)
296 ed(0:order) = scale*ed(0:order)
297
298 e_0(k) = e_0(k) + rho(k)*ed(0)
299 IF (order >= 1) THEN
300 e_rho(k) = e_rho(k) + ed(0) + rho(k)*ed(1)
301 END IF
302 IF (order >= 2) THEN
303 e_rho_rho(k) = e_rho_rho(k) + 2.0_dp*ed(1) + rho(k)*ed(2)
304 END IF
305 IF (order >= 3) THEN
306 e_rho_rho_rho(k) = e_rho_rho_rho(k) + 3.0_dp*ed(2) + rho(k)*ed(3)
307 END IF
308
309 END IF
310
311 END DO
312!$OMP END PARALLEL DO
313
314 END SUBROUTINE perdew_wang_lda_calc
315
316! **************************************************************************************************
317!> \brief Calculate the correlation energy and its derivatives
318!> wrt to rho (the electron density) up to 3rd order. This
319!> is the LSD version of the Perdew-Wang correlation energy
320!> If no order argument is given, then the routine calculates
321!> just the energy.
322!> \param method ...
323!> \param rho_set ...
324!> \param deriv_set ...
325!> \param order maximum order of derivatives to calculate
326!> order must lie between 0 and 3.
327!> \param scale ...
328! **************************************************************************************************
329 SUBROUTINE perdew_wang_lsd_eval(method, rho_set, deriv_set, order, scale)
330 INTEGER, INTENT(in) :: method
331 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
332 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
333 INTEGER, INTENT(IN), OPTIONAL :: order
334 REAL(kind=dp), INTENT(in) :: scale
335
336 CHARACTER(len=*), PARAMETER :: routinen = 'perdew_wang_lsd_eval'
337
338 INTEGER :: npoints, timer_handle
339 INTEGER, DIMENSION(2, 3) :: bo
340 REAL(kind=dp) :: rho_cutoff
341 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
342 eab, eabb, eb, ebb, ebbb
343 TYPE(xc_derivative_type), POINTER :: deriv
344
345 CALL timeset(routinen, timer_handle)
346 NULLIFY (a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
347 CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, &
348 local_bounds=bo, rho_cutoff=rho_cutoff)
349 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
350
351 CALL perdew_wang_init(method, rho_cutoff)
352
353 ! meaningful default for the arrays we don't need: let us make compiler
354 ! and debugger happy...
355 dummy => a
356
357 e_0 => dummy
358 ea => dummy; eb => dummy
359 eaa => dummy; eab => dummy; ebb => dummy
360 eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
361
362 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
363 allocate_deriv=.true.)
364 CALL xc_derivative_get(deriv, deriv_data=e_0)
365 IF (order >= 1) THEN
366 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
367 allocate_deriv=.true.)
368 CALL xc_derivative_get(deriv, deriv_data=ea)
369 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
370 allocate_deriv=.true.)
371 CALL xc_derivative_get(deriv, deriv_data=eb)
372 END IF
373 IF (order >= 2) THEN
374 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
375 allocate_deriv=.true.)
376 CALL xc_derivative_get(deriv, deriv_data=eaa)
377 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
378 allocate_deriv=.true.)
379 CALL xc_derivative_get(deriv, deriv_data=eab)
380 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
381 allocate_deriv=.true.)
382 CALL xc_derivative_get(deriv, deriv_data=ebb)
383 END IF
384 IF (order >= 3) THEN
385 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
386 allocate_deriv=.true.)
387 CALL xc_derivative_get(deriv, deriv_data=eaaa)
388 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
389 allocate_deriv=.true.)
390 CALL xc_derivative_get(deriv, deriv_data=eaab)
391 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
392 allocate_deriv=.true.)
393 CALL xc_derivative_get(deriv, deriv_data=eabb)
394 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
395 allocate_deriv=.true.)
396 CALL xc_derivative_get(deriv, deriv_data=ebbb)
397 END IF
398 IF (order > 3) THEN
399 cpabort("derivatives bigger than 3 not implemented")
400 END IF
401
402 CALL perdew_wang_lsd_calc(a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
403 ebbb, npoints, order, scale)
404
405 CALL timestop(timer_handle)
406
407 END SUBROUTINE perdew_wang_lsd_eval
408
409! **************************************************************************************************
410!> \brief ...
411!> \param rhoa ...
412!> \param rhob ...
413!> \param e_0 ...
414!> \param ea ...
415!> \param eb ...
416!> \param eaa ...
417!> \param eab ...
418!> \param ebb ...
419!> \param eaaa ...
420!> \param eaab ...
421!> \param eabb ...
422!> \param ebbb ...
423!> \param npoints ...
424!> \param order ...
425!> \param scale ...
426! **************************************************************************************************
427 SUBROUTINE perdew_wang_lsd_calc(rhoa, rhob, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
428 ebbb, npoints, order, scale)
429 !FM low-level computation routine
430 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob
431 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, &
432 eabb, ebbb
433 INTEGER, INTENT(in) :: npoints, order
434 REAL(kind=dp), INTENT(in) :: scale
435
436 INTEGER :: k
437 REAL(kind=dp) :: rho
438 REAL(kind=dp), DIMENSION(0:9) :: ed
439
440!$OMP PARALLEL DO PRIVATE (k, rho, ed) DEFAULT(NONE)&
441!$OMP SHARED(npoints,rhoa,rhob,eps_rho,order,e_0,ea,eb,eaa,eab,ebb,eaaa,eaab,eabb,ebbb,scale)
442 DO k = 1, npoints
443
444 rho = rhoa(k) + rhob(k)
445 IF (rho > eps_rho) THEN
446
447 ed = 0.0_dp
448 CALL pw_lsd_ed_loc(rhoa(k), rhob(k), ed, order)
449 ed = ed*scale
450
451 e_0(k) = e_0(k) + rho*ed(0)
452 IF (order >= 1) THEN
453 ea(k) = ea(k) + ed(0) + rho*ed(1)
454 eb(k) = eb(k) + ed(0) + rho*ed(2)
455 END IF
456 IF (order >= 2) THEN
457 eaa(k) = eaa(k) + 2.0_dp*ed(1) + rho*ed(3)
458 eab(k) = eab(k) + ed(1) + ed(2) + rho*ed(4)
459 ebb(k) = ebb(k) + 2.0_dp*ed(2) + rho*ed(5)
460 END IF
461 IF (order >= 3) THEN
462 eaaa(k) = eaaa(k) + 3.0_dp*ed(3) + rho*ed(6)
463 eaab(k) = eaab(k) + 2.0_dp*ed(4) + ed(3) + rho*ed(7)
464 eabb(k) = eabb(k) + 2.0_dp*ed(4) + ed(5) + rho*ed(8)
465 ebbb(k) = ebbb(k) + 3.0_dp*ed(5) + rho*ed(9)
466 END IF
467
468 END IF
469
470 END DO
471
472 END SUBROUTINE perdew_wang_lsd_calc
473
474! **************************************************************************************************
475!> \brief ...
476!> \param rho_a ...
477!> \param rho_b ...
478!> \param fxc_aa ...
479!> \param fxc_ab ...
480!> \param fxc_bb ...
481!> \param scalec ...
482!> \param eps_rho ...
483! **************************************************************************************************
484 SUBROUTINE perdew_wang_fxc_calc(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb, scalec, eps_rho)
485 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
486 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_ab, fxc_bb
487 REAL(kind=dp), INTENT(in) :: scalec, eps_rho
488
489 INTEGER :: i, j, k
490 INTEGER, DIMENSION(2, 3) :: bo
491 REAL(kind=dp) :: rho, rhoa, rhob, eaa, eab, ebb
492 REAL(kind=dp), DIMENSION(0:9) :: ed
493
494 CALL perdew_wang_init(pw_orig, eps_rho)
495 bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
496!$OMP PARALLEL DO PRIVATE(i,j,k,rho,rhoa,rhob,ed,eaa,eab,ebb) DEFAULT(NONE)&
497!$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_ab,fxc_bb,scalec,eps_rho)
498 DO k = bo(1, 3), bo(2, 3)
499 DO j = bo(1, 2), bo(2, 2)
500 DO i = bo(1, 1), bo(2, 1)
501 rhoa = rho_a%array(i, j, k)
502 rhob = rho_b%array(i, j, k)
503 rho = rhoa + rhob
504 IF (rho > eps_rho) THEN
505 ed = 0.0_dp
506 CALL pw_lsd_ed_loc(rhoa, rhob, ed, 2)
507 ed = ed*scalec
508 eaa = 2.0_dp*ed(1) + rho*ed(3)
509 eab = ed(1) + ed(2) + rho*ed(4)
510 ebb = 2.0_dp*ed(2) + rho*ed(5)
511 fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + eaa
512 fxc_ab%array(i, j, k) = fxc_ab%array(i, j, k) + eab
513 fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + ebb
514 END IF
515 END DO
516 END DO
517 END DO
518
519 END SUBROUTINE perdew_wang_fxc_calc
520
521! **************************************************************************************************
522!> \brief ...
523!> \param r ...
524!> \param z ...
525!> \param g ...
526!> \param order ...
527! **************************************************************************************************
528 PURE SUBROUTINE calc_g(r, z, g, order)
529
530! Calculates g and its derivatives wrt r up to 3rd order, where:
531!
532! g = .... for r < 1
533! g = .... for r > 100 and everywhere else
534! g = 2A(1+a1*r)ln(1+1/(2A(b1*r^1/2 + b2*r + b3*r^(3/2) + b4*r^2))).
535
536 REAL(kind=dp), INTENT(IN) :: r
537 INTEGER, INTENT(IN) :: z
538 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: g
539 INTEGER, INTENT(IN) :: order
540
541 REAL(kind=dp) :: a1_, a_, b1_, b2_, b3_, b4_, rr, rsr, &
542 sr, t11, t12, t14, t15, t16, t20, t22, &
543 t3, t40, t44, t45, t47, t48, t55, t56
544
545 a_ = a(z); a1_ = a1(z)
546 b1_ = b1(z); b2_ = b2(z); b3_ = b3(z); b4_ = b4(z)
547
548 sr = sqrt(r)
549 rsr = r*sr
550 rr = r*r
551
552 IF (r < 0.5_dp) THEN
553
554 ! order 0 must always be calculated
555 g(0) = c0(z)*log(r) - c1(z) + c2(z)*r*log(r) - c3(z)*r
556 IF (order >= 1) g(1) = c0(z)/r + c2(z)*log(r) + c2(z) - c3(z)
557 IF (order >= 2) g(2) = -c0(z)/rr + c2(z)/r
558 IF (order >= 3) g(3) = 2.0_dp*c0(z)/(rr*r) - c2(z)/rr
559
560 ELSE IF (r <= 100.0_dp) THEN
561
562 t3 = 1.0_dp + a1_*r
563 t11 = b1_*sr + b2_*r + b3_*rsr + b4_*rr
564 t12 = t11**2
565 t15 = 1.0_dp + 0.5_dp/a_/t11
566 t16 = log(t15)
567 t20 = 0.5_dp*b1_/sr + b2_ + 1.5_dp*b3_*sr + 2.0_dp*b4_*r
568
569 ! order 0 must always be calculated
570 g(0) = -2.0_dp*a_*t3*t16
571
572 IF (order >= 1) THEN
573
574 g(1) = -2.0_dp*a_*a1_*t16 + t3*t20/(t12*t15)
575
576 END IF
577
578 IF (order >= 2) THEN
579
580 t40 = -0.25_dp*b1_/rsr + 0.75_dp*b3_/sr + 2.0_dp*b4_
581
582 g(2) = 2.0_dp*a1_*t20/(t12*t15) &
583 - 2.0_dp*(t20**2)*t3/(t12*t11*t15) &
584 + t3*t40/(t12*t15) &
585 + 0.5_dp*t3*(t20**2)/(a_*(t12**2)*(t15**2))
586
587 END IF
588
589 IF (order >= 3) THEN
590
591 t14 = 1.0_dp/t12/t11
592 t22 = t20**2
593 t56 = t22*t20
594 t47 = t15**2
595 t48 = 1.0_dp/t47
596
597 t44 = t12**2
598 t45 = 1.0_dp/t44
599 t55 = t3*t45
600
601 g(3) = &
602 -6.0_dp*a1_*t14*t22/t15 &
603 + 3.0_dp*a1_*t40/(t15*t12) &
604 + 1.5_dp*a1_*t45*t22*t48/a_ &
605 + 6.0_dp*t55*t56/t15 &
606 - 6.0_dp*t3*t14*t20*t40/t15 &
607 - 3.0_dp*t3*t56*t48/(a_*t44*t11) &
608 + 0.375_dp*t3*(b1_/(rr*sr) - b3_/rsr)/(t12*t15) &
609 + 1.5_dp*t55*t40*t48*t20/a_ &
610 + 0.5_dp*t3*t56/((a_**2)*t44*t12*t47*t15)
611
612 END IF
613
614 ELSE
615
616 ! order 0 must always be calculated
617 g(0) = -d0(z)/r + d1(z)/rsr
618 IF (order >= 1) g(1) = d0(z)/rr - 1.5_dp*d1(z)/(rsr*r)
619 IF (order >= 2) g(2) = -2.0_dp*d0(z)/(rr*r) + 3.75_dp*d1(z)/(rsr*rr)
620 IF (order >= 3) g(3) = 6.0_dp*d0(z)/(rr*rr) - 13.125_dp*d1(z)/(rsr*rr*r)
621
622 END IF
623
624 END SUBROUTINE calc_g
625
626! **************************************************************************************************
627!> \brief ...
628!> \param rho ...
629!> \param ed ...
630!> \param order ...
631! **************************************************************************************************
632 SUBROUTINE pw_lda_ed_loc(rho, ed, order)
633
634 REAL(kind=dp), INTENT(IN) :: rho
635 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
636 INTEGER, INTENT(IN) :: order
637
638 INTEGER :: m, order_
639 LOGICAL, DIMENSION(0:3) :: calc
640 REAL(kind=dp), DIMENSION(0:3) :: e0, r
641
642 order_ = order
643 ed = 0
644 calc = .false.
645
646 IF (order_ >= 0) THEN
647 calc(0:order_) = .true.
648 ELSE
649 order_ = -1*order_
650 calc(order_) = .true.
651 END IF
652
653 CALL calc_rs(rho, r(0))
654 CALL calc_g(r(0), 0, e0, order_)
655
656 IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
657 IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
658 IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
659
660 m = 0
661 IF (calc(0)) THEN
662 ed(m) = e0(0)
663 m = m + 1
664 END IF
665 IF (calc(1)) THEN
666 ed(m) = e0(1)*r(1)
667 m = m + 1
668 END IF
669 IF (calc(2)) THEN
670 ed(m) = e0(2)*r(1)**2 + e0(1)*r(2)
671 m = m + 1
672 END IF
673 IF (calc(3)) THEN
674 ed(m) = e0(3)*r(1)**3 + e0(2)*3.0_dp*r(1)*r(2) + e0(1)*r(3)
675 END IF
676
677 END SUBROUTINE pw_lda_ed_loc
678
679! **************************************************************************************************
680!> \brief ...
681!> \param a ...
682!> \param b ...
683!> \param ed ...
684!> \param order ...
685! **************************************************************************************************
686 SUBROUTINE pw_lsd_ed_loc(a, b, ed, order)
687
688 REAL(kind=dp), INTENT(IN) :: a, b
689 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
690 INTEGER, INTENT(IN) :: order
691
692 INTEGER :: m, order_
693 LOGICAL, DIMENSION(0:3) :: calc
694 REAL(kind=dp) :: rho, tr, trr, trrr, trrz, trz, trzz, tz, &
695 tzz, tzzz
696 REAL(kind=dp), DIMENSION(0:3) :: ac, e0, e1, f, r
697 REAL(kind=dp), DIMENSION(0:3, 0:3) :: z
698
699 order_ = order
700 calc = .false.
701
702 IF (order_ > 0) THEN
703 calc(0:order_) = .true.
704 ELSE
705 order_ = -1*order_
706 calc(order_) = .true.
707 END IF
708
709 rho = a + b
710
711 CALL calc_fx(a, b, f(0:order_), order_)
712 CALL calc_rs(rho, r(0))
713 CALL calc_g(r(0), -1, ac(0:order_), order_)
714 CALL calc_g(r(0), 0, e0(0:order_), order_)
715 CALL calc_g(r(0), 1, e1(0:order_), order_)
716 CALL calc_z(a, b, z(0:order_, 0:order_), order_)
717
718!! calculate first partial derivatives
719 IF (order_ >= 1) THEN
720 r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
721 tr = e0(1) &
722 + fpp*ac(1)*f(0) &
723 - fpp*ac(1)*f(0)*z(0, 0)**4 &
724 + (e1(1) - e0(1))*f(0)*z(0, 0)**4
725 tz = fpp*ac(0)*f(1) &
726 - fpp*ac(0)*f(1)*z(0, 0)**4 &
727 - fpp*ac(0)*f(0)*4.0_dp*z(0, 0)**3 &
728 + (e1(0) - e0(0))*f(1)*z(0, 0)**4 &
729 + (e1(0) - e0(0))*f(0)*4.0_dp*z(0, 0)**3
730 END IF
731
732!! calculate second partial derivatives
733 IF (order_ >= 2) THEN
734 r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
735 trr = e0(2) &
736 + fpp*ac(2)*f(0) &
737 - fpp*ac(2)*f(0)*z(0, 0)**4 &
738 + (e1(2) - e0(2))*f(0)*z(0, 0)**4
739 trz = fpp*ac(1)*f(1) &
740 - fpp*ac(1)*f(1)*z(0, 0)**4 &
741 - fpp*ac(1)*f(0)*4.0_dp*z(0, 0)**3 &
742 + (e1(1) - e0(1))*f(1)*z(0, 0)**4 &
743 + (e1(1) - e0(1))*f(0)*4.0_dp*z(0, 0)**3
744 tzz = fpp*ac(0)*f(2) &
745 - fpp*ac(0)*f(2)*z(0, 0)**4 &
746 - fpp*ac(0)*f(1)*8.0_dp*z(0, 0)**3 &
747 - fpp*ac(0)*f(0)*12.0_dp*z(0, 0)**2 &
748 + (e1(0) - e0(0))*f(2)*z(0, 0)**4 &
749 + (e1(0) - e0(0))*f(1)*8.0_dp*z(0, 0)**3 &
750 + (e1(0) - e0(0))*f(0)*12.0_dp*z(0, 0)**2
751 END IF
752
753!! calculate third derivatives
754 IF (order_ >= 3) THEN
755
756 r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
757
758 trrr = e0(3) &
759 + fpp*ac(3)*f(0) &
760 - fpp*ac(3)*f(0)*z(0, 0)**4 &
761 + (e1(3) - e0(3))*f(0)*z(0, 0)**4
762
763 trrz = fpp*ac(2)*f(1) &
764 - fpp*ac(2)*f(1)*z(0, 0)**4 &
765 - fpp*ac(2)*f(0)*4.0_dp*z(0, 0)**3 &
766 + (e1(2) - e0(2))*f(1)*z(0, 0)**4 &
767 + (e1(2) - e0(2))*f(0)*4.0_dp*z(0, 0)**3
768
769 trzz = fpp*ac(1)*f(2) &
770 - fpp*ac(1)*f(2)*z(0, 0)**4 &
771 - fpp*ac(1)*f(1)*8.0_dp*z(0, 0)**3 &
772 - fpp*ac(1)*f(0)*12.0_dp*z(0, 0)**2 &
773 + (e1(1) - e0(1))*f(2)*z(0, 0)**4 &
774 + (e1(1) - e0(1))*f(1)*8.0_dp*z(0, 0)**3 &
775 + (e1(1) - e0(1))*f(0)*12.0_dp*z(0, 0)**2
776
777 tzzz = fpp*ac(0)*f(3) &
778 - fpp*ac(0)*f(3)*z(0, 0)**4 &
779 - fpp*ac(0)*f(2)*12.0_dp*z(0, 0)**3 &
780 - fpp*ac(0)*f(1)*36.0_dp*z(0, 0)**2 &
781 - fpp*ac(0)*f(0)*24.0_dp*z(0, 0) &
782 + (e1(0) - e0(0))*f(3)*z(0, 0)**4 &
783 + (e1(0) - e0(0))*f(2)*12.0_dp*z(0, 0)**3 &
784 + (e1(0) - e0(0))*f(1)*36.0_dp*z(0, 0)**2 &
785 + (e1(0) - e0(0))*f(0)*24.0_dp*z(0, 0)
786 END IF
787
788 m = 0
789 IF (calc(0)) THEN
790 ed(m) = e0(0) &
791 + fpp*ac(0)*f(0)*(1.0_dp - z(0, 0)**4) &
792 + (e1(0) - e0(0))*f(0)*z(0, 0)**4
793 m = m + 1
794 END IF
795 IF (calc(1)) THEN
796 ed(m) = tr*r(1) + tz*z(1, 0)
797 ed(m + 1) = tr*r(1) + tz*z(0, 1)
798 m = m + 2
799 END IF
800 IF (calc(2)) THEN
801 ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
802 + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
803 ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
804 + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
805 ed(m + 2) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(0, 1) &
806 + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
807 m = m + 3
808 END IF
809 IF (calc(3)) THEN
810 ed(m) = &
811 trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(1, 0) &
812 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(1, 0) + tr*r(3) &
813 + 3.0_dp*trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**3 &
814 + 3.0_dp*trz*r(1)*z(2, 0) &
815 + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
816 ed(m + 1) = &
817 trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(1, 0) + z(0, 1)) &
818 + 2.0_dp*trzz*r(1)*z(1, 0)*z(0, 1) &
819 + 2.0_dp*trz*(r(2)*z(1, 0) + r(1)*z(1, 1)) &
820 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(0, 1) + tr*r(3) &
821 + trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**2*z(0, 1) &
822 + 2.0_dp*tzz*z(1, 0)*z(1, 1) &
823 + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
824 ed(m + 2) = &
825 trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(0, 1) + z(1, 0)) &
826 + 2.0_dp*trzz*r(1)*z(0, 1)*z(1, 0) &
827 + 2.0_dp*trz*(r(2)*z(0, 1) + r(1)*z(1, 1)) &
828 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(1, 0) + tr*r(3) &
829 + trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**2*z(1, 0) &
830 + 2.0_dp*tzz*z(0, 1)*z(1, 1) &
831 + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
832 ed(m + 3) = &
833 trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(0, 1) &
834 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(0, 1) + tr*r(3) &
835 + 3.0_dp*trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**3 &
836 + 3.0_dp*trz*r(1)*z(0, 2) &
837 + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
838 END IF
839
840 END SUBROUTINE pw_lsd_ed_loc
841
842END MODULE xc_perdew_wang
xc_functionals_utilities::calc_fx
Definition xc_functionals_utilities.F:37
xc_functionals_utilities::calc_rs
Definition xc_functionals_utilities.F:41
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
pw_types
Definition pw_types.F:24
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_functionals_utilities::calc_z
pure subroutine, public calc_z(a, b, z, order)
...
Definition xc_functionals_utilities.F:348
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::pw_orig
integer, parameter, public pw_orig
Definition xc_input_constants.F:88
xc_input_constants::pw_dmc
integer, parameter, public pw_dmc
Definition xc_input_constants.F:88
xc_input_constants::pw_vmc
integer, parameter, public pw_vmc
Definition xc_input_constants.F:88
xc_perdew_wang
Calculate the Perdew-Wang correlation potential and energy density and ist derivatives with respect t...
Definition xc_perdew_wang.F:17
xc_perdew_wang::perdew_wang_fxc_calc
subroutine, public perdew_wang_fxc_calc(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb, scalec, eps_rho)
...
Definition xc_perdew_wang.F:485
xc_perdew_wang::perdew_wang_lda_eval
subroutine, public perdew_wang_lda_eval(method, rho_set, deriv_set, order, scale)
Calculate the correlation energy and its derivatives wrt to rho (the electron density) up to 3rd orde...
Definition xc_perdew_wang.F:207
xc_perdew_wang::perdew_wang_info
subroutine, public perdew_wang_info(method, lsd, reference, shortform, needs, max_deriv, scale)
Return some info on the functionals.
Definition xc_perdew_wang.F:68
xc_perdew_wang::perdew_wang_lsd_eval
subroutine, public perdew_wang_lsd_eval(method, rho_set, deriv_set, order, scale)
Calculate the correlation energy and its derivatives wrt to rho (the electron density) up to 3rd orde...
Definition xc_perdew_wang.F:330
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
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
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