(git:e7e05ae)
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 ! **************************************************************************************************
17 MODULE 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 
54 CONTAINS
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 
858 END MODULE xc_perdew_wang
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