(git:1f285aa)
xc_perdew_zunger.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-Zunger 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_zunger
18  USE bibliography, ONLY: ortiz1994,&
19  perdew1981,&
20  cite_reference
21  USE input_section_types, ONLY: section_vals_type,&
22  section_vals_val_get
23  USE kinds, ONLY: dp
24  USE xc_derivative_desc, ONLY: deriv_rho,&
25  deriv_rhoa,&
26  deriv_rhob
27  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
28  xc_dset_get_derivative
29  USE xc_derivative_types, ONLY: xc_derivative_get,&
30  xc_derivative_type
31  USE xc_functionals_utilities, ONLY: calc_fx,&
32  calc_rs,&
33  calc_z,&
34  set_util
35  USE xc_input_constants, ONLY: pz_dmc,&
36  pz_orig,&
37  pz_vmc
38  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
39  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
40  xc_rho_set_type
41 #include "../base/base_uses.f90"
42 
43  IMPLICIT NONE
44  PRIVATE
45 
46  LOGICAL :: initialized = .false.
47  REAL(KIND=dp), DIMENSION(0:1) :: a = 0.0_dp, b = 0.0_dp, c = 0.0_dp, d = 0.0_dp, &
48  b1 = 0.0_dp, b2 = 0.0_dp, ga = 0.0_dp
49 
50  REAL(KIND=dp), PARAMETER :: epsilon = 5.e-13
51  REAL(KIND=dp) :: eps_rho
52 
53  PUBLIC :: pz_info, pz_lda_eval, pz_lsd_eval
54  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_perdew_zunger'
55 
56 CONTAINS
57 
58 ! **************************************************************************************************
59 !> \brief Return some info on the functionals.
60 !> \param method ...
61 !> \param lsd ...
62 !> \param reference CHARACTER(*), INTENT(OUT), OPTIONAL - full reference
63 !> \param shortform CHARACTER(*), INTENT(OUT), OPTIONAL - short reference
64 !> \param needs ...
65 !> \param max_deriv ...
66 !> \par History
67 !> 18-MAR-2002, TCH, working version
68 ! **************************************************************************************************
69  SUBROUTINE pz_info(method, lsd, reference, shortform, needs, max_deriv)
70  INTEGER, INTENT(in) :: method
71  LOGICAL, INTENT(in) :: lsd
72  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
73  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
74  INTEGER, INTENT(out), OPTIONAL :: max_deriv
75 
76  CHARACTER(len=4) :: p_string
77 
78  SELECT CASE (method)
79  CASE DEFAULT
80  cpabort("Unsupported parametrization")
81  CASE (pz_orig)
82  p_string = 'ORIG'
83  CASE (pz_dmc)
84  p_string = 'DMC'
85  CASE (pz_vmc)
86  p_string = 'VMC'
87  END SELECT
88 
89  IF (PRESENT(reference)) THEN
90  reference = "J. P. Perdew and Alex Zunger," &
91  //" Phys. Rev. B 23, 5048 (1981)" &
92  //"["//trim(p_string)//"]"
93  IF (.NOT. lsd) THEN
94  IF (len_trim(reference) + 6 < len(reference)) THEN
95  reference(len_trim(reference):len_trim(reference) + 6) = ' {LDA}'
96  END IF
97  END IF
98  END IF
99  IF (PRESENT(shortform)) THEN
100  shortform = "J. P. Perdew et al., PRB 23, 5048 (1981)" &
101  //"["//trim(p_string)//"]"
102  IF (.NOT. lsd) THEN
103  IF (len_trim(shortform) + 6 < len(shortform)) THEN
104  shortform(len_trim(shortform):len_trim(shortform) + 6) = ' {LDA}'
105  END IF
106  END IF
107  END IF
108  IF (PRESENT(needs)) THEN
109  IF (lsd) THEN
110  needs%rho_spin = .true.
111  ELSE
112  needs%rho = .true.
113  END IF
114  END IF
115  IF (PRESENT(max_deriv)) max_deriv = 3
116 
117  END SUBROUTINE pz_info
118 
119 ! **************************************************************************************************
120 !> \brief Calculate the correlation energy and its derivatives
121 !> wrt to rho (the electron density) up to 3rd order. This
122 !> is the LDA version of the Perdew-Zunger correlation energy
123 !> If no order argument is given, then the routine calculates
124 !> just the energy.
125 !> \param method ...
126 !> \param rho_set ...
127 !> \param deriv_set ...
128 !> \param order INTEGER, OPTIONAL - order of derivatives to calculate
129 !> order must lie between -3 and 3. If it is negative then only
130 !> that order will be calculated, otherwise all derivatives up to
131 !> that order will be calculated.
132 !> \param pz_params input parameter (scaling)
133 !> \par History
134 !> 01.2007 added scaling [Manuel Guidon]
135 ! **************************************************************************************************
136  SUBROUTINE pz_lda_eval(method, rho_set, deriv_set, order, pz_params)
137 
138  INTEGER, INTENT(in) :: method
139  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
140  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
141  INTEGER, INTENT(in) :: order
142  TYPE(section_vals_type), POINTER :: pz_params
143 
144  CHARACTER(len=*), PARAMETER :: routinen = 'pz_lda_eval'
145 
146  INTEGER :: npoints, timer_handle
147  INTEGER, DIMENSION(2, 3) :: bo
148  REAL(kind=dp) :: rho_cutoff, sc
149  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
150  POINTER :: dummy, e_0, e_rho, e_rho_rho, &
151  e_rho_rho_rho, rho
152  TYPE(xc_derivative_type), POINTER :: deriv
153 
154  CALL timeset(routinen, timer_handle)
155  NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
156 
157  CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
158 
159  CALL xc_rho_set_get(rho_set, rho=rho, &
160  local_bounds=bo, rho_cutoff=rho_cutoff)
161  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
162 
163  CALL pz_init(method, rho_cutoff)
164 
165  dummy => rho
166 
167  e_0 => dummy
168  e_rho => dummy
169  e_rho_rho => dummy
170  e_rho_rho_rho => dummy
171 
172  IF (order >= 0) THEN
173  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
174  allocate_deriv=.true.)
175  CALL xc_derivative_get(deriv, deriv_data=e_0)
176  END IF
177  IF (order >= 1 .OR. order == -1) THEN
178  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
179  allocate_deriv=.true.)
180  CALL xc_derivative_get(deriv, deriv_data=e_rho)
181  END IF
182  IF (order >= 2 .OR. order == -2) THEN
183  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
184  allocate_deriv=.true.)
185  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
186  END IF
187  IF (order >= 3 .OR. order == -3) THEN
188  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
189  allocate_deriv=.true.)
190  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
191  END IF
192  IF (order > 3 .OR. order < -3) THEN
193  cpabort("derivatives bigger than 3 not implemented")
194  END IF
195 
196  CALL pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
197 
198  CALL timestop(timer_handle)
199 
200  END SUBROUTINE pz_lda_eval
201 
202 ! **************************************************************************************************
203 !> \brief ...
204 !> \param rho ...
205 !> \param e_0 ...
206 !> \param e_rho ...
207 !> \param e_rho_rho ...
208 !> \param e_rho_rho_rho ...
209 !> \param npoints ...
210 !> \param order ...
211 !> \param sc ...
212 ! **************************************************************************************************
213  SUBROUTINE pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
214  !FM low level calc routine
215  REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho
216  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
217  INTEGER, INTENT(in) :: npoints, order
218  REAL(kind=dp) :: sc
219 
220  INTEGER :: k
221  REAL(kind=dp), DIMENSION(0:3) :: ed
222 
223 !$OMP PARALLEL DO PRIVATE ( k, ed ) DEFAULT(NONE)&
224 !$OMP SHARED(npoints,rho,eps_rho,order,e_0,e_rho,e_rho_rho,e_rho_rho_rho,sc)
225  DO k = 1, npoints
226 
227  IF (rho(k) > eps_rho) THEN
228 
229  CALL pz_lda_ed_loc(rho(k), ed, abs(order), sc)
230 
231  IF (order >= 0) THEN
232  e_0(k) = e_0(k) + rho(k)*ed(0)
233  END IF
234  IF (order >= 1 .OR. order == -1) THEN
235  e_rho(k) = e_rho(k) + ed(0) + rho(k)*ed(1)
236  END IF
237  IF (order >= 2 .OR. order == -2) THEN
238  e_rho_rho(k) = e_rho_rho(k) + 2.0_dp*ed(1) + rho(k)*ed(2)
239  END IF
240  IF (order >= 3 .OR. order == -3) THEN
241  e_rho_rho_rho(k) = e_rho_rho_rho(k) + 3.0_dp*ed(2) + rho(k)*ed(3)
242  END IF
243 
244  END IF
245 
246  END DO
247 !$OMP END PARALLEL DO
248 
249  END SUBROUTINE pz_lda_calc
250 
251 ! **************************************************************************************************
252 !> \brief Calculate the correlation energy and its derivatives
253 !> wrt to rho (the electron density) up to 3rd order. This
254 !> is the LSD version of the Perdew-Zunger correlation energy
255 !> If no order argument is given, then the routine calculates
256 !> just the energy.
257 !> \param method ...
258 !> \param rho_set ...
259 !> \param deriv_set ...
260 !> \param order INTEGER, OPTIONAL - order of derivatives to calculate
261 !> order must lie between -3 and 3. If it is negative then only
262 !> that order will be calculated, otherwise all derivatives up to
263 !> that order will be calculated.
264 !> \param pz_params input parameter (scaling)
265 !> \par History
266 !> 01.2007 added scaling [Manuel Guidon]
267 ! **************************************************************************************************
268  SUBROUTINE pz_lsd_eval(method, rho_set, deriv_set, order, pz_params)
269  INTEGER, INTENT(in) :: method
270  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
271  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
272  INTEGER, INTENT(IN), OPTIONAL :: order
273  TYPE(section_vals_type), POINTER :: pz_params
274 
275  CHARACTER(len=*), PARAMETER :: routinen = 'pz_lsd_eval'
276 
277  INTEGER :: npoints, timer_handle
278  INTEGER, DIMENSION(2, 3) :: bo
279  REAL(kind=dp) :: rho_cutoff, sc
280  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
281  POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
282  eab, eabb, ebb, ebbb
283  TYPE(xc_derivative_type), POINTER :: deriv
284 
285  CALL timeset(routinen, timer_handle)
286  NULLIFY (a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
287 
288  CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
289 
290  CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, &
291  local_bounds=bo, rho_cutoff=rho_cutoff)
292  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
293 
294  CALL pz_init(method, rho_cutoff)
295 
296  dummy => a
297 
298  e_0 => dummy
299  ea => dummy;
300  eaa => dummy; eab => dummy; ebb => dummy
301  eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
302 
303  IF (order >= 0) THEN
304  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
305  allocate_deriv=.true.)
306  CALL xc_derivative_get(deriv, deriv_data=e_0)
307  END IF
308  IF (order >= 1 .OR. order == -1) THEN
309  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
310  allocate_deriv=.true.)
311  CALL xc_derivative_get(deriv, deriv_data=ea)
312  END IF
313  IF (order >= 2 .OR. order == -2) THEN
314  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
315  allocate_deriv=.true.)
316  CALL xc_derivative_get(deriv, deriv_data=eaa)
317  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
318  allocate_deriv=.true.)
319  CALL xc_derivative_get(deriv, deriv_data=eab)
320  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
321  allocate_deriv=.true.)
322  CALL xc_derivative_get(deriv, deriv_data=ebb)
323  END IF
324  IF (order >= 3 .OR. order == -3) THEN
325  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
326  allocate_deriv=.true.)
327  CALL xc_derivative_get(deriv, deriv_data=eaaa)
328  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
329  allocate_deriv=.true.)
330  CALL xc_derivative_get(deriv, deriv_data=eaab)
331  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
332  allocate_deriv=.true.)
333  CALL xc_derivative_get(deriv, deriv_data=eabb)
334  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
335  allocate_deriv=.true.)
336  CALL xc_derivative_get(deriv, deriv_data=ebbb)
337  END IF
338  IF (order > 3 .OR. order < -3) THEN
339  cpabort("derivatives bigger than 3 not implemented")
340  END IF
341 
342  CALL pz_lsd_calc(a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, &
343  ebbb, npoints, order, sc)
344 
345  CALL timestop(timer_handle)
346 
347  END SUBROUTINE pz_lsd_eval
348 
349 ! **************************************************************************************************
350 !> \brief ...
351 !> \param a ...
352 !> \param b ...
353 !> \param e_0 ...
354 !> \param ea ...
355 !> \param eaa ...
356 !> \param eab ...
357 !> \param ebb ...
358 !> \param eaaa ...
359 !> \param eaab ...
360 !> \param eabb ...
361 !> \param ebbb ...
362 !> \param npoints ...
363 !> \param order ...
364 !> \param sc ...
365 ! **************************************************************************************************
366  SUBROUTINE pz_lsd_calc(a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, &
367  ebbb, npoints, order, sc)
368  !FM low-level computation routine
369  REAL(kind=dp), DIMENSION(*), INTENT(in) :: a, b
370  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, ea, eaa, eab, ebb, eaaa, eaab, &
371  eabb, ebbb
372  INTEGER, INTENT(in) :: npoints, order
373  REAL(kind=dp), INTENT(IN) :: sc
374 
375  INTEGER :: k, order_
376  REAL(kind=dp) :: rho
377  REAL(kind=dp), DIMENSION(0:9) :: ed
378 
379  order_ = abs(order)
380 
381 !$OMP PARALLEL DO PRIVATE ( k, rho, ed ) DEFAULT(NONE)&
382 !$OMP SHARED(order_,order,npoints,eps_rho,A,b,sc,e_0,ea,eaa,eab,ebb,eaaa,eaab,eabb,ebbb)
383  DO k = 1, npoints
384 
385  rho = a(k) + b(k)
386 
387  IF (rho > eps_rho) THEN
388 
389  CALL pz_lsd_ed_loc(a(k), b(k), ed, order_, sc)
390 
391  IF (order >= 0) THEN
392  e_0(k) = e_0(k) + rho*ed(0)
393  END IF
394 
395  IF (order >= 1 .OR. order == -1) THEN
396  ea(k) = ea(k) + ed(0) + rho*ed(1)
397  ea(k) = ea(k) + ed(0) + rho*ed(2)
398  END IF
399 
400  IF (order >= 2 .OR. order == -2) THEN
401  eaa(k) = eaa(k) + 2.0_dp*ed(1) + rho*ed(3)
402  eab(k) = eab(k) + ed(1) + ed(2) + rho*ed(4)
403  ebb(k) = ebb(k) + 2.0_dp*ed(2) + rho*ed(5)
404  END IF
405 
406  IF (order >= 3 .OR. order == -3) THEN
407  eaaa(k) = eaaa(k) + 3.0_dp*ed(3) + rho*ed(6)
408  eaab(k) = eaab(k) + 2.0_dp*ed(4) + ed(3) + rho*ed(7)
409  eabb(k) = eabb(k) + 2.0_dp*ed(4) + ed(5) + rho*ed(8)
410  ebbb(k) = ebbb(k) + 3.0_dp*ed(5) + rho*ed(9)
411  END IF
412 
413  END IF
414 
415  END DO
416 !$OMP END PARALLEL DO
417 
418  END SUBROUTINE pz_lsd_calc
419 
420 ! **************************************************************************************************
421 !> \brief Initializes the functionals
422 !> \param method CHARACTER(3) - name of the method used for parameters
423 !> \param cutoff REAL(KIND=dp) - the cutoff density
424 !> \par History
425 !> 18-MAR-2002, TCH, working version
426 !> \note see functionals_utilities
427 ! **************************************************************************************************
428  SUBROUTINE pz_init(method, cutoff)
429 
430  INTEGER, INTENT(IN) :: method
431  REAL(kind=dp), INTENT(IN) :: cutoff
432 
433  CALL set_util(cutoff)
434 
435  eps_rho = cutoff
436 
437  initialized = .false.
438 
439  SELECT CASE (method)
440 
441  CASE DEFAULT
442  cpabort("Unknown method")
443 
444  CASE (pz_orig)
445  CALL cite_reference(perdew1981)
446  ga(0) = -0.1423_dp; ga(1) = -0.0843_dp
447  b1(0) = 1.0529_dp; b1(1) = 1.3981_dp
448  b2(0) = 0.3334_dp; b2(1) = 0.2611_dp
449  a(0) = 0.0311_dp; a(1) = 0.01555_dp
450  b(0) = -0.048_dp; b(1) = -0.0269_dp
451  c(0) = +0.0020_dp; c(1) = +0.0007_dp
452  d(0) = -0.0116_dp; d(1) = -0.0048_dp
453 
454  CASE (pz_dmc)
455  CALL cite_reference(ortiz1994)
456  ga(0) = -0.103756_dp; ga(1) = -0.065951_dp
457  b1(0) = 0.56371_dp; b1(1) = 1.11846_dp
458  b2(0) = 0.27358_dp; b2(1) = 0.18797_dp
459  a(0) = 0.031091_dp; a(1) = 0.015545_dp
460  b(0) = -0.046644_dp; b(1) = -0.025599_dp
461  c(0) = -0.00419_dp; c(1) = -0.00329_dp
462  d(0) = -0.00983_dp; d(1) = -0.00300_dp
463 
464  CASE (pz_vmc)
465  CALL cite_reference(ortiz1994)
466  ga(0) = -0.093662_dp; ga(1) = -0.055331_dp
467  b1(0) = 0.49453_dp; b1(1) = 0.93766_dp
468  b2(0) = 0.25534_dp; b2(1) = 0.14829_dp
469  a(0) = 0.031091_dp; a(1) = 0.015545_dp
470  b(0) = -0.046644_dp; b(1) = -0.025599_dp
471  c(0) = -0.00884_dp; c(1) = -0.00677_dp
472  d(0) = -0.00688_dp; d(1) = -0.00093_dp
473 
474  END SELECT
475 
476  initialized = .true.
477 
478  END SUBROUTINE pz_init
479 
480 ! **************************************************************************************************
481 !> \brief ...
482 !> \param r ...
483 !> \param z ...
484 !> \param g ...
485 !> \param order ...
486 ! **************************************************************************************************
487  SUBROUTINE calc_g(r, z, g, order)
488 
489  REAL(kind=dp), INTENT(IN) :: r
490  INTEGER, INTENT(IN) :: z
491  REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: g
492  INTEGER, INTENT(IN) :: order
493 
494  REAL(kind=dp) :: rr, rsr, sr
495 
496  IF (r >= 1.0_dp) THEN
497 
498  sr = sqrt(r)
499  ! order 0 must always be calculated
500  g(0) = ga(z)/(1.0_dp + b1(z)*sr + b2(z)*r)
501  IF (order >= 1) THEN
502  g(1) = -ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))/ &
503  (1.0_dp + b1(z)*sr + b2(z)*r)**2
504  END IF
505  IF (order >= 2) THEN
506  rsr = r*sr
507  g(2) = &
508  2.0_dp*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))**2 &
509  /(1.0_dp + b1(z)*sr + b2(z)*r)**3 &
510  + ga(z)*b1(z) &
511  /(4.0_dp*(1.0_dp + b1(z)*sr + b2(z)*r)**2*rsr)
512  END IF
513  IF (order >= 3) THEN
514  g(3) = &
515  -6.0_dp*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))**3/ &
516  (1.0_dp + b1(z)*sr + b2(z)*r)**4 &
517  - (3.0_dp/2.0_dp)*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))*b1(z)/ &
518  ((1.0_dp + b1(z)*sr + b2(z)*r)**3*rsr) &
519  - (3.0_dp/8.0_dp)*ga(z)*b1(z)/ &
520  ((1.0_dp + b1(z)*sr + b2(z)*r)**2*r*rsr)
521  END IF
522 
523  ELSE
524 
525  ! order 0 must always be calculated
526  g(0) = a(z)*log(r) + b(z) + c(z)*r*log(r) + d(z)*r
527  IF (order >= 1) THEN
528  g(1) = a(z)/r + c(z)*log(r) + c(z) + d(z)
529  END IF
530  IF (order >= 2) THEN
531  rr = r*r
532  g(2) = -a(z)/rr + c(z)/r
533  END IF
534  IF (order >= 3) THEN
535  g(3) = 2.0_dp*a(z)/(rr*r) - c(z)/rr
536  END IF
537 
538  END IF
539 
540  END SUBROUTINE calc_g
541 
542 ! **************************************************************************************************
543 !> \brief ...
544 !> \param rho ...
545 !> \param ed ...
546 !> \param order ...
547 !> \param sc ...
548 ! **************************************************************************************************
549  SUBROUTINE pz_lda_ed_loc(rho, ed, order, sc)
550 
551  REAL(kind=dp), INTENT(IN) :: rho
552  REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
553  INTEGER, INTENT(IN) :: order
554  REAL(dp), INTENT(IN) :: sc
555 
556  INTEGER :: m, order_
557  LOGICAL, DIMENSION(0:3) :: calc
558  REAL(kind=dp), DIMENSION(0:3) :: e0, r
559 
560  calc = .false.
561 
562  order_ = order
563  IF (order_ >= 0) THEN
564  calc(0:order_) = .true.
565  ELSE
566  order_ = -1*order_
567  calc(order_) = .true.
568  END IF
569 
570  CALL calc_rs(rho, r(0))
571  CALL calc_g(r(0), 0, e0, order_)
572 
573  IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
574  IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
575  IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
576 
577  m = 0
578  IF (calc(0)) THEN
579  ed(m) = sc*e0(0)
580  m = m + 1
581  END IF
582  IF (calc(1)) THEN
583  ed(m) = sc*e0(1)*r(1)
584  m = m + 1
585  END IF
586  IF (calc(2)) THEN
587  ed(m) = sc*e0(2)*r(1)**2 + sc*e0(1)*r(2)
588  m = m + 1
589  END IF
590  IF (calc(3)) THEN
591  ed(m) = sc*e0(3)*r(1)**3 + sc*e0(2)*3.0_dp*r(1)*r(2) + sc*e0(1)*r(3)
592  END IF
593 
594  END SUBROUTINE pz_lda_ed_loc
595 
596 ! **************************************************************************************************
597 !> \brief ...
598 !> \param a ...
599 !> \param b ...
600 !> \param ed ...
601 !> \param order ...
602 !> \param sc ...
603 ! **************************************************************************************************
604  SUBROUTINE pz_lsd_ed_loc(a, b, ed, order, sc)
605 
606  REAL(kind=dp), INTENT(IN) :: a, b
607  REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
608  INTEGER, INTENT(IN), OPTIONAL :: order
609  REAL(kind=dp), INTENT(IN) :: sc
610 
611  INTEGER :: m, order_
612  LOGICAL, DIMENSION(0:3) :: calc
613  REAL(kind=dp) :: rho, tr, trr, trrr, trrz, trz, trzz, tz, &
614  tzz, tzzz
615  REAL(kind=dp), DIMENSION(0:3) :: e0, e1, f, r
616  REAL(kind=dp), DIMENSION(0:3, 0:3) :: z
617 
618  calc = .false.
619 
620  order_ = 0
621  IF (PRESENT(order)) order_ = order
622 
623  IF (order_ >= 0) THEN
624  calc(0:order_) = .true.
625  ELSE
626  order_ = -1*order_
627  calc(order_) = .true.
628  END IF
629 
630  rho = a + b
631 
632  CALL calc_fx(a, b, f(0:order_), order_)
633  CALL calc_rs(rho, r(0))
634  CALL calc_g(r(0), 0, e0(0:order_), order_)
635  CALL calc_g(r(0), 1, e1(0:order_), order_)
636  CALL calc_z(a, b, z(:, :), order_)
637 
638 !! calculate first partial derivatives
639  IF (order_ >= 1) THEN
640  r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
641  tr = e0(1) + (e1(1) - e0(1))*f(0)
642  tz = (e1(0) - e0(0))*f(1)
643  END IF
644 
645 !! calculate second partial derivatives
646  IF (order_ >= 2) THEN
647  r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
648  trr = e0(2) + (e1(2) - e0(2))*f(0)
649  trz = (e1(1) - e0(1))*f(1)
650  tzz = (e1(0) - e0(0))*f(2)
651  END IF
652 
653 !! calculate third derivatives
654  IF (order_ >= 3) THEN
655  r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
656  trrr = e0(3) + (e1(3) - e0(3))*f(0)
657  trrz = (e1(2) - e0(2))*f(1)
658  trzz = (e1(1) - e0(1))*f(2)
659  tzzz = (e1(0) - e0(0))*f(3)
660  END IF
661 
662  m = 0
663  IF (calc(0)) THEN
664  ed(m) = e0(0) + (e1(0) - e0(0))*f(0)
665  ed(m) = ed(m)*sc
666  m = m + 1
667  END IF
668 
669  IF (calc(1)) THEN
670  ed(m) = tr*r(1) + tz*z(1, 0)
671  ed(m) = ed(m)*sc
672  ed(m + 1) = tr*r(1) + tz*z(0, 1)
673  ed(m + 1) = ed(m + 1)*sc
674  m = m + 2
675  END IF
676 
677  IF (calc(2)) THEN
678  ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
679  + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
680  ed(m) = ed(m)*sc
681  ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
682  + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
683  ed(m + 1) = ed(m + 1)*sc
684  ed(m + 2) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(0, 1) &
685  + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
686  ed(m + 2) = ed(m + 2)*sc
687  m = m + 3
688  END IF
689 
690  IF (calc(3)) THEN
691  ed(m) = trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(1, 0) &
692  + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(1, 0) &
693  + tr*r(3) + 3.0_dp*trzz*r(1)*z(1, 0)**2 &
694  + tzzz*z(1, 0)**3 + 3.0_dp*trz*r(1)*z(2, 0) &
695  + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
696  ed(m) = ed(m)*sc
697  ed(m + 1) = trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(1, 0) + z(0, 1)) &
698  + 2.0_dp*trzz*r(1)*z(1, 0)*z(0, 1) &
699  + 2.0_dp*trz*(r(2)*z(1, 0) + r(1)*z(1, 1)) &
700  + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(0, 1) + tr*r(3) &
701  + trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**2*z(0, 1) &
702  + 2.0_dp*tzz*z(1, 0)*z(1, 1) &
703  + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
704  ed(m + 1) = ed(m + 1)*sc
705  ed(m + 2) = trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(0, 1) + z(1, 0)) &
706  + 2.0_dp*trzz*r(1)*z(0, 1)*z(1, 0) &
707  + 2.0_dp*trz*(r(2)*z(0, 1) + r(1)*z(1, 1)) &
708  + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(1, 0) + tr*r(3) &
709  + trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**2*z(1, 0) &
710  + 2.0_dp*tzz*z(0, 1)*z(1, 1) &
711  + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
712  ed(m + 2) = ed(m + 2)*sc
713  ed(m + 3) = trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(0, 1) &
714  + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(0, 1) + tr*r(3) &
715  + 3.0_dp*trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**3 &
716  + 3.0_dp*trz*r(1)*z(0, 2) &
717  + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
718  ed(m + 3) = ed(m + 3)*sc
719  END IF
720 
721  END SUBROUTINE pz_lsd_ed_loc
722 
723 END MODULE xc_perdew_zunger
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::perdew1981
integer, save, public perdew1981
Definition: bibliography.F:43
bibliography::ortiz1994
integer, save, public ortiz1994
Definition: bibliography.F:43
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition: xc_derivative_desc.F:26
xc_derivative_desc::deriv_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::pz_dmc
integer, parameter, public pz_dmc
Definition: xc_input_constants.F:88
xc_input_constants::pz_vmc
integer, parameter, public pz_vmc
Definition: xc_input_constants.F:88
xc_input_constants::pz_orig
integer, parameter, public pz_orig
Definition: xc_input_constants.F:88
xc_perdew_zunger
Calculate the Perdew-Zunger correlation potential and energy density and ist derivatives with respect...
Definition: xc_perdew_zunger.F:17
xc_perdew_zunger::pz_lda_eval
subroutine, public pz_lda_eval(method, rho_set, deriv_set, order, pz_params)
Calculate the correlation energy and its derivatives wrt to rho (the electron density) up to 3rd orde...
Definition: xc_perdew_zunger.F:137
xc_perdew_zunger::pz_lsd_eval
subroutine, public pz_lsd_eval(method, rho_set, deriv_set, order, pz_params)
Calculate the correlation energy and its derivatives wrt to rho (the electron density) up to 3rd orde...
Definition: xc_perdew_zunger.F:269
xc_perdew_zunger::pz_info
subroutine, public pz_info(method, lsd, reference, shortform, needs, max_deriv)
Return some info on the functionals.
Definition: xc_perdew_zunger.F:70
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