(git:374b731)
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-2024 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 order of derivatives to calculate
203!> order must lie between -3 and 3. If it is negative then only
204!> that order will be calculated, otherwise all derivatives up to
205!> that order will be calculated.
206!> \param scale ...
207! **************************************************************************************************
208 SUBROUTINE perdew_wang_lda_eval(method, rho_set, deriv_set, order, scale)
209
210 INTEGER, INTENT(in) :: method
211 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
212 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
213 INTEGER, INTENT(in) :: order
214 REAL(kind=dp), INTENT(in) :: scale
215
216 CHARACTER(len=*), PARAMETER :: routinen = 'perdew_wang_lda_eval'
217
218 INTEGER :: npoints, timer_handle
219 INTEGER, DIMENSION(2, 3) :: bo
220 REAL(kind=dp) :: rho_cutoff
221 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_rho, e_rho_rho, &
222 e_rho_rho_rho, rho
223 TYPE(xc_derivative_type), POINTER :: deriv
224
225 CALL timeset(routinen, timer_handle)
226 NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
227 CALL xc_rho_set_get(rho_set, rho=rho, &
228 local_bounds=bo, rho_cutoff=rho_cutoff)
229 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
230
231 CALL perdew_wang_init(method, rho_cutoff)
232
233 dummy => rho
234
235 e_0 => dummy
236 e_rho => dummy
237 e_rho_rho => dummy
238 e_rho_rho_rho => dummy
239
240 IF (order >= 0) THEN
241 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
242 allocate_deriv=.true.)
243 CALL xc_derivative_get(deriv, deriv_data=e_0)
244 END IF
245 IF (order >= 1 .OR. order == -1) THEN
246 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
247 allocate_deriv=.true.)
248 CALL xc_derivative_get(deriv, deriv_data=e_rho)
249 END IF
250 IF (order >= 2 .OR. order == -2) THEN
251 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
252 allocate_deriv=.true.)
253 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
254 END IF
255 IF (order >= 3 .OR. order == -3) THEN
256 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
257 allocate_deriv=.true.)
258 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
259 END IF
260 IF (order > 3 .OR. order < -3) THEN
261 cpabort("derivatives bigger than 3 not implemented")
262 END IF
263
264 CALL perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
265 npoints, order, scale)
266
267 CALL timestop(timer_handle)
268
269 END SUBROUTINE perdew_wang_lda_eval
270
271! **************************************************************************************************
272!> \brief ...
273!> \param rho ...
274!> \param e_0 ...
275!> \param e_rho ...
276!> \param e_rho_rho ...
277!> \param e_rho_rho_rho ...
278!> \param npoints ...
279!> \param order ...
280!> \param scale ...
281! **************************************************************************************************
282 SUBROUTINE perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, scale)
283 !FM low level calc routine
284 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho
285 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
286 INTEGER, INTENT(in) :: npoints, order
287 REAL(kind=dp), INTENT(in) :: scale
288
289 INTEGER :: abs_order, k
290 REAL(kind=dp), DIMENSION(0:3) :: ed
291
292 abs_order = abs(order)
293
294!$OMP PARALLEL DO PRIVATE (k, ed) DEFAULT(NONE)&
295!$OMP SHARED(npoints,rho,eps_rho,abs_order,scale,e_0,e_rho,e_rho_rho,e_rho_rho_rho,order)
296 DO k = 1, npoints
297
298 IF (rho(k) > eps_rho) THEN
299!! order_ is positive as it must be in this case:
300!! ec(:,2) needs ed(:,1) for example
301 CALL pw_lda_ed_loc(rho(k), ed, abs_order)
302 ed(0:abs_order) = scale*ed(0:abs_order)
303
304 IF (order >= 0) THEN
305 e_0(k) = e_0(k) + rho(k)*ed(0)
306 END IF
307 IF (order >= 1 .OR. order == -1) THEN
308 e_rho(k) = e_rho(k) + ed(0) + rho(k)*ed(1)
309 END IF
310 IF (order >= 2 .OR. order == -2) THEN
311 e_rho_rho(k) = e_rho_rho(k) + 2.0_dp*ed(1) + rho(k)*ed(2)
312 END IF
313 IF (order >= 3 .OR. order == -3) THEN
314 e_rho_rho_rho(k) = e_rho_rho_rho(k) + 3.0_dp*ed(2) + rho(k)*ed(3)
315 END IF
316
317 END IF
318
319 END DO
320!$OMP END PARALLEL DO
321
322 END SUBROUTINE perdew_wang_lda_calc
323
324! **************************************************************************************************
325!> \brief Calculate the correlation energy and its derivatives
326!> wrt to rho (the electron density) up to 3rd order. This
327!> is the LSD version of the Perdew-Wang correlation energy
328!> If no order argument is given, then the routine calculates
329!> just the energy.
330!> \param method ...
331!> \param rho_set ...
332!> \param deriv_set ...
333!> \param order order of derivatives to calculate
334!> order must lie between -3 and 3. If it is negative then only
335!> that order will be calculated, otherwise all derivatives up to
336!> that order will be calculated.
337!> \param scale ...
338! **************************************************************************************************
339 SUBROUTINE perdew_wang_lsd_eval(method, rho_set, deriv_set, order, scale)
340 INTEGER, INTENT(in) :: method
341 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
342 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
343 INTEGER, INTENT(IN), OPTIONAL :: order
344 REAL(kind=dp), INTENT(in) :: scale
345
346 CHARACTER(len=*), PARAMETER :: routinen = 'perdew_wang_lsd_eval'
347
348 INTEGER :: npoints, timer_handle
349 INTEGER, DIMENSION(2, 3) :: bo
350 REAL(kind=dp) :: rho_cutoff
351 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
352 eab, eabb, eb, ebb, ebbb
353 TYPE(xc_derivative_type), POINTER :: deriv
354
355 CALL timeset(routinen, timer_handle)
356 NULLIFY (a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
357 CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, &
358 local_bounds=bo, rho_cutoff=rho_cutoff)
359 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
360
361 CALL perdew_wang_init(method, rho_cutoff)
362
363 ! meaningful default for the arrays we don't need: let us make compiler
364 ! and debugger happy...
365 dummy => a
366
367 e_0 => dummy
368 ea => dummy; eb => dummy
369 eaa => dummy; eab => dummy; ebb => dummy
370 eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
371
372 IF (order >= 0) THEN
373 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
374 allocate_deriv=.true.)
375 CALL xc_derivative_get(deriv, deriv_data=e_0)
376 END IF
377 IF (order >= 1 .OR. order == -1) THEN
378 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
379 allocate_deriv=.true.)
380 CALL xc_derivative_get(deriv, deriv_data=ea)
381 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
382 allocate_deriv=.true.)
383 CALL xc_derivative_get(deriv, deriv_data=eb)
384 END IF
385 IF (order >= 2 .OR. order == -2) THEN
386 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
387 allocate_deriv=.true.)
388 CALL xc_derivative_get(deriv, deriv_data=eaa)
389 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
390 allocate_deriv=.true.)
391 CALL xc_derivative_get(deriv, deriv_data=eab)
392 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
393 allocate_deriv=.true.)
394 CALL xc_derivative_get(deriv, deriv_data=ebb)
395 END IF
396 IF (order >= 3 .OR. order == -3) THEN
397 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
398 allocate_deriv=.true.)
399 CALL xc_derivative_get(deriv, deriv_data=eaaa)
400 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
401 allocate_deriv=.true.)
402 CALL xc_derivative_get(deriv, deriv_data=eaab)
403 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
404 allocate_deriv=.true.)
405 CALL xc_derivative_get(deriv, deriv_data=eabb)
406 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
407 allocate_deriv=.true.)
408 CALL xc_derivative_get(deriv, deriv_data=ebbb)
409 END IF
410 IF (order > 3 .OR. order < -3) THEN
411 cpabort("derivatives bigger than 3 not implemented")
412 END IF
413
414 CALL perdew_wang_lsd_calc(a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
415 ebbb, npoints, order, scale)
416
417 CALL timestop(timer_handle)
418
419 END SUBROUTINE perdew_wang_lsd_eval
420
421! **************************************************************************************************
422!> \brief ...
423!> \param rhoa ...
424!> \param rhob ...
425!> \param e_0 ...
426!> \param ea ...
427!> \param eb ...
428!> \param eaa ...
429!> \param eab ...
430!> \param ebb ...
431!> \param eaaa ...
432!> \param eaab ...
433!> \param eabb ...
434!> \param ebbb ...
435!> \param npoints ...
436!> \param order ...
437!> \param scale ...
438! **************************************************************************************************
439 SUBROUTINE perdew_wang_lsd_calc(rhoa, rhob, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
440 ebbb, npoints, order, scale)
441 !FM low-level computation routine
442 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob
443 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, &
444 eabb, ebbb
445 INTEGER, INTENT(in) :: npoints, order
446 REAL(kind=dp), INTENT(in) :: scale
447
448 INTEGER :: abs_order, k
449 REAL(kind=dp) :: rho
450 REAL(kind=dp), DIMENSION(0:9) :: ed
451
452 abs_order = abs(order)
453
454!$OMP PARALLEL DO PRIVATE (k, rho, ed) DEFAULT(NONE)&
455!$OMP SHARED(npoints,rhoa,rhob,eps_rho,abs_order,order,e_0,ea,eb,eaa,eab,ebb,eaaa,eaab,eabb,ebbb,scale)
456 DO k = 1, npoints
457
458 rho = rhoa(k) + rhob(k)
459 IF (rho > eps_rho) THEN
460
461 ed = 0.0_dp
462 CALL pw_lsd_ed_loc(rhoa(k), rhob(k), ed, abs_order)
463 ed = ed*scale
464
465 IF (order >= 0) THEN
466 e_0(k) = e_0(k) + rho*ed(0)
467 END IF
468 IF (order >= 1 .OR. order == -1) THEN
469 ea(k) = ea(k) + ed(0) + rho*ed(1)
470 eb(k) = eb(k) + ed(0) + rho*ed(2)
471 END IF
472 IF (order >= 2 .OR. order == -2) THEN
473 eaa(k) = eaa(k) + 2.0_dp*ed(1) + rho*ed(3)
474 eab(k) = eab(k) + ed(1) + ed(2) + rho*ed(4)
475 ebb(k) = ebb(k) + 2.0_dp*ed(2) + rho*ed(5)
476 END IF
477 IF (order >= 3 .OR. order == -3) THEN
478 eaaa(k) = eaaa(k) + 3.0_dp*ed(3) + rho*ed(6)
479 eaab(k) = eaab(k) + 2.0_dp*ed(4) + ed(3) + rho*ed(7)
480 eabb(k) = eabb(k) + 2.0_dp*ed(4) + ed(5) + rho*ed(8)
481 ebbb(k) = ebbb(k) + 3.0_dp*ed(5) + rho*ed(9)
482 END IF
483
484 END IF
485
486 END DO
487
488 END SUBROUTINE perdew_wang_lsd_calc
489
490! **************************************************************************************************
491!> \brief ...
492!> \param rho_a ...
493!> \param rho_b ...
494!> \param fxc_aa ...
495!> \param fxc_ab ...
496!> \param fxc_bb ...
497!> \param scalec ...
498!> \param eps_rho ...
499! **************************************************************************************************
500 SUBROUTINE perdew_wang_fxc_calc(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb, scalec, eps_rho)
501 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
502 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_ab, fxc_bb
503 REAL(kind=dp), INTENT(in) :: scalec, eps_rho
504
505 INTEGER :: i, j, k
506 INTEGER, DIMENSION(2, 3) :: bo
507 REAL(kind=dp) :: rho, rhoa, rhob, eaa, eab, ebb
508 REAL(kind=dp), DIMENSION(0:9) :: ed
509
510 CALL perdew_wang_init(pw_orig, eps_rho)
511 bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
512!$OMP PARALLEL DO PRIVATE(i,j,k,rho,rhoa,rhob,ed,eaa,eab,ebb) DEFAULT(NONE)&
513!$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_ab,fxc_bb,scalec,eps_rho)
514 DO k = bo(1, 3), bo(2, 3)
515 DO j = bo(1, 2), bo(2, 2)
516 DO i = bo(1, 1), bo(2, 1)
517 rhoa = rho_a%array(i, j, k)
518 rhob = rho_b%array(i, j, k)
519 rho = rhoa + rhob
520 IF (rho > eps_rho) THEN
521 ed = 0.0_dp
522 CALL pw_lsd_ed_loc(rhoa, rhob, ed, 2)
523 ed = ed*scalec
524 eaa = 2.0_dp*ed(1) + rho*ed(3)
525 eab = ed(1) + ed(2) + rho*ed(4)
526 ebb = 2.0_dp*ed(2) + rho*ed(5)
527 fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + eaa
528 fxc_ab%array(i, j, k) = fxc_ab%array(i, j, k) + eab
529 fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + ebb
530 END IF
531 END DO
532 END DO
533 END DO
534
535 END SUBROUTINE perdew_wang_fxc_calc
536
537! **************************************************************************************************
538!> \brief ...
539!> \param r ...
540!> \param z ...
541!> \param g ...
542!> \param order ...
543! **************************************************************************************************
544 PURE SUBROUTINE calc_g(r, z, g, order)
545
546! Calculates g and its derivatives wrt r up to 3rd order, where:
547!
548! g = .... for r < 1
549! g = .... for r > 100 and everywhere else
550! g = 2A(1+a1*r)ln(1+1/(2A(b1*r^1/2 + b2*r + b3*r^(3/2) + b4*r^2))).
551
552 REAL(kind=dp), INTENT(IN) :: r
553 INTEGER, INTENT(IN) :: z
554 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: g
555 INTEGER, INTENT(IN) :: order
556
557 REAL(kind=dp) :: a1_, a_, b1_, b2_, b3_, b4_, rr, rsr, &
558 sr, t11, t12, t14, t15, t16, t20, t22, &
559 t3, t40, t44, t45, t47, t48, t55, t56
560
561 a_ = a(z); a1_ = a1(z)
562 b1_ = b1(z); b2_ = b2(z); b3_ = b3(z); b4_ = b4(z)
563
564 sr = sqrt(r)
565 rsr = r*sr
566 rr = r*r
567
568 IF (r < 0.5_dp) THEN
569
570 ! order 0 must always be calculated
571 g(0) = c0(z)*log(r) - c1(z) + c2(z)*r*log(r) - c3(z)*r
572 IF (order >= 1) g(1) = c0(z)/r + c2(z)*log(r) + c2(z) - c3(z)
573 IF (order >= 2) g(2) = -c0(z)/rr + c2(z)/r
574 IF (order >= 3) g(3) = 2.0_dp*c0(z)/(rr*r) - c2(z)/rr
575
576 ELSE IF (r <= 100.0_dp) THEN
577
578 t3 = 1.0_dp + a1_*r
579 t11 = b1_*sr + b2_*r + b3_*rsr + b4_*rr
580 t12 = t11**2
581 t15 = 1.0_dp + 0.5_dp/a_/t11
582 t16 = log(t15)
583 t20 = 0.5_dp*b1_/sr + b2_ + 1.5_dp*b3_*sr + 2.0_dp*b4_*r
584
585 ! order 0 must always be calculated
586 g(0) = -2.0_dp*a_*t3*t16
587
588 IF (order >= 1) THEN
589
590 g(1) = -2.0_dp*a_*a1_*t16 + t3*t20/(t12*t15)
591
592 END IF
593
594 IF (order >= 2) THEN
595
596 t40 = -0.25_dp*b1_/rsr + 0.75_dp*b3_/sr + 2.0_dp*b4_
597
598 g(2) = 2.0_dp*a1_*t20/(t12*t15) &
599 - 2.0_dp*(t20**2)*t3/(t12*t11*t15) &
600 + t3*t40/(t12*t15) &
601 + 0.5_dp*t3*(t20**2)/(a_*(t12**2)*(t15**2))
602
603 END IF
604
605 IF (order >= 3) THEN
606
607 t14 = 1.0_dp/t12/t11
608 t22 = t20**2
609 t56 = t22*t20
610 t47 = t15**2
611 t48 = 1.0_dp/t47
612
613 t44 = t12**2
614 t45 = 1.0_dp/t44
615 t55 = t3*t45
616
617 g(3) = &
618 -6.0_dp*a1_*t14*t22/t15 &
619 + 3.0_dp*a1_*t40/(t15*t12) &
620 + 1.5_dp*a1_*t45*t22*t48/a_ &
621 + 6.0_dp*t55*t56/t15 &
622 - 6.0_dp*t3*t14*t20*t40/t15 &
623 - 3.0_dp*t3*t56*t48/(a_*t44*t11) &
624 + 0.375_dp*t3*(b1_/(rr*sr) - b3_/rsr)/(t12*t15) &
625 + 1.5_dp*t55*t40*t48*t20/a_ &
626 + 0.5_dp*t3*t56/((a_**2)*t44*t12*t47*t15)
627
628 END IF
629
630 ELSE
631
632 ! order 0 must always be calculated
633 g(0) = -d0(z)/r + d1(z)/rsr
634 IF (order >= 1) g(1) = d0(z)/rr - 1.5_dp*d1(z)/(rsr*r)
635 IF (order >= 2) g(2) = -2.0_dp*d0(z)/(rr*r) + 3.75_dp*d1(z)/(rsr*rr)
636 IF (order >= 3) g(3) = 6.0_dp*d0(z)/(rr*rr) - 13.125_dp*d1(z)/(rsr*rr*r)
637
638 END IF
639
640 END SUBROUTINE calc_g
641
642! **************************************************************************************************
643!> \brief ...
644!> \param rho ...
645!> \param ed ...
646!> \param order ...
647! **************************************************************************************************
648 SUBROUTINE pw_lda_ed_loc(rho, ed, order)
649
650 REAL(kind=dp), INTENT(IN) :: rho
651 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
652 INTEGER, INTENT(IN) :: order
653
654 INTEGER :: m, order_
655 LOGICAL, DIMENSION(0:3) :: calc
656 REAL(kind=dp), DIMENSION(0:3) :: e0, r
657
658 order_ = order
659 ed = 0
660 calc = .false.
661
662 IF (order_ >= 0) THEN
663 calc(0:order_) = .true.
664 ELSE
665 order_ = -1*order_
666 calc(order_) = .true.
667 END IF
668
669 CALL calc_rs(rho, r(0))
670 CALL calc_g(r(0), 0, e0, order_)
671
672 IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
673 IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
674 IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
675
676 m = 0
677 IF (calc(0)) THEN
678 ed(m) = e0(0)
679 m = m + 1
680 END IF
681 IF (calc(1)) THEN
682 ed(m) = e0(1)*r(1)
683 m = m + 1
684 END IF
685 IF (calc(2)) THEN
686 ed(m) = e0(2)*r(1)**2 + e0(1)*r(2)
687 m = m + 1
688 END IF
689 IF (calc(3)) THEN
690 ed(m) = e0(3)*r(1)**3 + e0(2)*3.0_dp*r(1)*r(2) + e0(1)*r(3)
691 END IF
692
693 END SUBROUTINE pw_lda_ed_loc
694
695! **************************************************************************************************
696!> \brief ...
697!> \param a ...
698!> \param b ...
699!> \param ed ...
700!> \param order ...
701! **************************************************************************************************
702 SUBROUTINE pw_lsd_ed_loc(a, b, ed, order)
703
704 REAL(kind=dp), INTENT(IN) :: a, b
705 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
706 INTEGER, INTENT(IN) :: order
707
708 INTEGER :: m, order_
709 LOGICAL, DIMENSION(0:3) :: calc
710 REAL(kind=dp) :: rho, tr, trr, trrr, trrz, trz, trzz, tz, &
711 tzz, tzzz
712 REAL(kind=dp), DIMENSION(0:3) :: ac, e0, e1, f, r
713 REAL(kind=dp), DIMENSION(0:3, 0:3) :: z
714
715 order_ = order
716 calc = .false.
717
718 IF (order_ > 0) THEN
719 calc(0:order_) = .true.
720 ELSE
721 order_ = -1*order_
722 calc(order_) = .true.
723 END IF
724
725 rho = a + b
726
727 CALL calc_fx(a, b, f(0:order_), order_)
728 CALL calc_rs(rho, r(0))
729 CALL calc_g(r(0), -1, ac(0:order_), order_)
730 CALL calc_g(r(0), 0, e0(0:order_), order_)
731 CALL calc_g(r(0), 1, e1(0:order_), order_)
732 CALL calc_z(a, b, z(0:order_, 0:order_), order_)
733
734!! calculate first partial derivatives
735 IF (order_ >= 1) THEN
736 r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
737 tr = e0(1) &
738 + fpp*ac(1)*f(0) &
739 - fpp*ac(1)*f(0)*z(0, 0)**4 &
740 + (e1(1) - e0(1))*f(0)*z(0, 0)**4
741 tz = fpp*ac(0)*f(1) &
742 - fpp*ac(0)*f(1)*z(0, 0)**4 &
743 - fpp*ac(0)*f(0)*4.0_dp*z(0, 0)**3 &
744 + (e1(0) - e0(0))*f(1)*z(0, 0)**4 &
745 + (e1(0) - e0(0))*f(0)*4.0_dp*z(0, 0)**3
746 END IF
747
748!! calculate second partial derivatives
749 IF (order_ >= 2) THEN
750 r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
751 trr = e0(2) &
752 + fpp*ac(2)*f(0) &
753 - fpp*ac(2)*f(0)*z(0, 0)**4 &
754 + (e1(2) - e0(2))*f(0)*z(0, 0)**4
755 trz = fpp*ac(1)*f(1) &
756 - fpp*ac(1)*f(1)*z(0, 0)**4 &
757 - fpp*ac(1)*f(0)*4.0_dp*z(0, 0)**3 &
758 + (e1(1) - e0(1))*f(1)*z(0, 0)**4 &
759 + (e1(1) - e0(1))*f(0)*4.0_dp*z(0, 0)**3
760 tzz = fpp*ac(0)*f(2) &
761 - fpp*ac(0)*f(2)*z(0, 0)**4 &
762 - fpp*ac(0)*f(1)*8.0_dp*z(0, 0)**3 &
763 - fpp*ac(0)*f(0)*12.0_dp*z(0, 0)**2 &
764 + (e1(0) - e0(0))*f(2)*z(0, 0)**4 &
765 + (e1(0) - e0(0))*f(1)*8.0_dp*z(0, 0)**3 &
766 + (e1(0) - e0(0))*f(0)*12.0_dp*z(0, 0)**2
767 END IF
768
769!! calculate third derivatives
770 IF (order_ >= 3) THEN
771
772 r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
773
774 trrr = e0(3) &
775 + fpp*ac(3)*f(0) &
776 - fpp*ac(3)*f(0)*z(0, 0)**4 &
777 + (e1(3) - e0(3))*f(0)*z(0, 0)**4
778
779 trrz = fpp*ac(2)*f(1) &
780 - fpp*ac(2)*f(1)*z(0, 0)**4 &
781 - fpp*ac(2)*f(0)*4.0_dp*z(0, 0)**3 &
782 + (e1(2) - e0(2))*f(1)*z(0, 0)**4 &
783 + (e1(2) - e0(2))*f(0)*4.0_dp*z(0, 0)**3
784
785 trzz = fpp*ac(1)*f(2) &
786 - fpp*ac(1)*f(2)*z(0, 0)**4 &
787 - fpp*ac(1)*f(1)*8.0_dp*z(0, 0)**3 &
788 - fpp*ac(1)*f(0)*12.0_dp*z(0, 0)**2 &
789 + (e1(1) - e0(1))*f(2)*z(0, 0)**4 &
790 + (e1(1) - e0(1))*f(1)*8.0_dp*z(0, 0)**3 &
791 + (e1(1) - e0(1))*f(0)*12.0_dp*z(0, 0)**2
792
793 tzzz = fpp*ac(0)*f(3) &
794 - fpp*ac(0)*f(3)*z(0, 0)**4 &
795 - fpp*ac(0)*f(2)*12.0_dp*z(0, 0)**3 &
796 - fpp*ac(0)*f(1)*36.0_dp*z(0, 0)**2 &
797 - fpp*ac(0)*f(0)*24.0_dp*z(0, 0) &
798 + (e1(0) - e0(0))*f(3)*z(0, 0)**4 &
799 + (e1(0) - e0(0))*f(2)*12.0_dp*z(0, 0)**3 &
800 + (e1(0) - e0(0))*f(1)*36.0_dp*z(0, 0)**2 &
801 + (e1(0) - e0(0))*f(0)*24.0_dp*z(0, 0)
802 END IF
803
804 m = 0
805 IF (calc(0)) THEN
806 ed(m) = e0(0) &
807 + fpp*ac(0)*f(0)*(1.0_dp - z(0, 0)**4) &
808 + (e1(0) - e0(0))*f(0)*z(0, 0)**4
809 m = m + 1
810 END IF
811 IF (calc(1)) THEN
812 ed(m) = tr*r(1) + tz*z(1, 0)
813 ed(m + 1) = tr*r(1) + tz*z(0, 1)
814 m = m + 2
815 END IF
816 IF (calc(2)) THEN
817 ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
818 + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
819 ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
820 + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
821 ed(m + 2) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(0, 1) &
822 + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
823 m = m + 3
824 END IF
825 IF (calc(3)) THEN
826 ed(m) = &
827 trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(1, 0) &
828 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(1, 0) + tr*r(3) &
829 + 3.0_dp*trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**3 &
830 + 3.0_dp*trz*r(1)*z(2, 0) &
831 + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
832 ed(m + 1) = &
833 trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(1, 0) + z(0, 1)) &
834 + 2.0_dp*trzz*r(1)*z(1, 0)*z(0, 1) &
835 + 2.0_dp*trz*(r(2)*z(1, 0) + r(1)*z(1, 1)) &
836 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(0, 1) + tr*r(3) &
837 + trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**2*z(0, 1) &
838 + 2.0_dp*tzz*z(1, 0)*z(1, 1) &
839 + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
840 ed(m + 2) = &
841 trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(0, 1) + z(1, 0)) &
842 + 2.0_dp*trzz*r(1)*z(0, 1)*z(1, 0) &
843 + 2.0_dp*trz*(r(2)*z(0, 1) + r(1)*z(1, 1)) &
844 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(1, 0) + tr*r(3) &
845 + trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**2*z(1, 0) &
846 + 2.0_dp*tzz*z(0, 1)*z(1, 1) &
847 + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
848 ed(m + 3) = &
849 trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(0, 1) &
850 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(0, 1) + tr*r(3) &
851 + 3.0_dp*trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**3 &
852 + 3.0_dp*trz*r(1)*z(0, 2) &
853 + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
854 END IF
855
856 END SUBROUTINE pw_lsd_ed_loc
857
858END 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: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_functionals_utilities::calc_z
pure subroutine, public calc_z(a, b, z, order)
...
Definition xc_functionals_utilities.F:340
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:84
xc_input_constants::pw_dmc
integer, parameter, public pw_dmc
Definition xc_input_constants.F:84
xc_input_constants::pw_vmc
integer, parameter, public pw_vmc
Definition xc_input_constants.F:84
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:501
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:209
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:340
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: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