(git:c641773)
Loading...
Searching...
No Matches
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-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the Perdew-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! **************************************************************************************************
17MODULE 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_dp
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
56CONTAINS
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 - maximum order of derivatives to calculate
129!> order must lie between 0 and 3.
130!> \param pz_params input parameter (scaling)
131!> \par History
132!> 01.2007 added scaling [Manuel Guidon]
133! **************************************************************************************************
134 SUBROUTINE pz_lda_eval(method, rho_set, deriv_set, order, pz_params)
135
136 INTEGER, INTENT(in) :: method
137 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
138 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
139 INTEGER, INTENT(in) :: order
140 TYPE(section_vals_type), POINTER :: pz_params
141
142 CHARACTER(len=*), PARAMETER :: routinen = 'pz_lda_eval'
143
144 INTEGER :: npoints, timer_handle
145 INTEGER, DIMENSION(2, 3) :: bo
146 REAL(kind=dp) :: rho_cutoff, sc
147 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
148 POINTER :: dummy, e_0, e_rho, e_rho_rho, &
149 e_rho_rho_rho, rho
150 TYPE(xc_derivative_type), POINTER :: deriv
151
152 CALL timeset(routinen, timer_handle)
153 NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
154
155 CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
156
157 CALL xc_rho_set_get(rho_set, rho=rho, &
158 local_bounds=bo, rho_cutoff=rho_cutoff)
159 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
160
161 CALL pz_init(method, rho_cutoff)
162
163 dummy => rho
164
165 e_0 => dummy
166 e_rho => dummy
167 e_rho_rho => dummy
168 e_rho_rho_rho => dummy
169
170 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
171 allocate_deriv=.true.)
172 CALL xc_derivative_get(deriv, deriv_data=e_0)
173 IF (order >= 1) THEN
174 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
175 allocate_deriv=.true.)
176 CALL xc_derivative_get(deriv, deriv_data=e_rho)
177 END IF
178 IF (order >= 2) THEN
179 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
180 allocate_deriv=.true.)
181 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
182 END IF
183 IF (order >= 3) THEN
184 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
185 allocate_deriv=.true.)
186 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
187 END IF
188 IF (order > 3) THEN
189 cpabort("derivatives bigger than 3 not implemented")
190 END IF
191
192 CALL pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
193
194 CALL timestop(timer_handle)
195
196 END SUBROUTINE pz_lda_eval
197
198! **************************************************************************************************
199!> \brief ...
200!> \param rho ...
201!> \param e_0 ...
202!> \param e_rho ...
203!> \param e_rho_rho ...
204!> \param e_rho_rho_rho ...
205!> \param npoints ...
206!> \param order ...
207!> \param sc ...
208! **************************************************************************************************
209 SUBROUTINE pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
210 !FM low level calc routine
211 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho
212 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
213 INTEGER, INTENT(in) :: npoints, order
214 REAL(kind=dp) :: sc
215
216 INTEGER :: k
217 REAL(kind=dp), DIMENSION(0:3) :: ed
218
219!$OMP PARALLEL DO PRIVATE ( k, ed ) DEFAULT(NONE)&
220!$OMP SHARED(npoints,rho,eps_rho,order,e_0,e_rho,e_rho_rho,e_rho_rho_rho,sc)
221 DO k = 1, npoints
222
223 IF (rho(k) > eps_rho) THEN
224
225 CALL pz_lda_ed_loc(rho(k), ed, abs(order), sc)
226
227 e_0(k) = e_0(k) + rho(k)*ed(0)
228 IF (order >= 1) THEN
229 e_rho(k) = e_rho(k) + ed(0) + rho(k)*ed(1)
230 END IF
231 IF (order >= 2) THEN
232 e_rho_rho(k) = e_rho_rho(k) + 2.0_dp*ed(1) + rho(k)*ed(2)
233 END IF
234 IF (order >= 3) THEN
235 e_rho_rho_rho(k) = e_rho_rho_rho(k) + 3.0_dp*ed(2) + rho(k)*ed(3)
236 END IF
237
238 END IF
239
240 END DO
241!$OMP END PARALLEL DO
242
243 END SUBROUTINE pz_lda_calc
244
245! **************************************************************************************************
246!> \brief Calculate the correlation energy and its derivatives
247!> wrt to rho (the electron density) up to 3rd order. This
248!> is the LSD version of the Perdew-Zunger correlation energy
249!> If no order argument is given, then the routine calculates
250!> just the energy.
251!> \param method ...
252!> \param rho_set ...
253!> \param deriv_set ...
254!> \param order INTEGER, OPTIONAL - maximum order of derivatives to calculate
255!> order must lie between 0 and 3.
256!> \param pz_params input parameter (scaling)
257!> \par History
258!> 01.2007 added scaling [Manuel Guidon]
259! **************************************************************************************************
260 SUBROUTINE pz_lsd_eval(method, rho_set, deriv_set, order, pz_params)
261 INTEGER, INTENT(in) :: method
262 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
263 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
264 INTEGER, INTENT(IN), OPTIONAL :: order
265 TYPE(section_vals_type), POINTER :: pz_params
266
267 CHARACTER(len=*), PARAMETER :: routinen = 'pz_lsd_eval'
268
269 INTEGER :: npoints, timer_handle
270 INTEGER, DIMENSION(2, 3) :: bo
271 REAL(kind=dp) :: rho_cutoff, sc
272 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
273 POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
274 eab, eabb, ebb, ebbb
275 TYPE(xc_derivative_type), POINTER :: deriv
276
277 CALL timeset(routinen, timer_handle)
278 NULLIFY (a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
279
280 CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
281
282 CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, &
283 local_bounds=bo, rho_cutoff=rho_cutoff)
284 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
285
286 CALL pz_init(method, rho_cutoff)
287
288 dummy => a
289
290 e_0 => dummy
291 ea => dummy
292 eaa => dummy; eab => dummy; ebb => dummy
293 eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
294
295 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
296 allocate_deriv=.true.)
297 CALL xc_derivative_get(deriv, deriv_data=e_0)
298 IF (order >= 1) THEN
299 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
300 allocate_deriv=.true.)
301 CALL xc_derivative_get(deriv, deriv_data=ea)
302 END IF
303 IF (order >= 2) THEN
304 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
305 allocate_deriv=.true.)
306 CALL xc_derivative_get(deriv, deriv_data=eaa)
307 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
308 allocate_deriv=.true.)
309 CALL xc_derivative_get(deriv, deriv_data=eab)
310 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
311 allocate_deriv=.true.)
312 CALL xc_derivative_get(deriv, deriv_data=ebb)
313 END IF
314 IF (order >= 3) THEN
315 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
316 allocate_deriv=.true.)
317 CALL xc_derivative_get(deriv, deriv_data=eaaa)
318 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
319 allocate_deriv=.true.)
320 CALL xc_derivative_get(deriv, deriv_data=eaab)
321 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
322 allocate_deriv=.true.)
323 CALL xc_derivative_get(deriv, deriv_data=eabb)
324 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
325 allocate_deriv=.true.)
326 CALL xc_derivative_get(deriv, deriv_data=ebbb)
327 END IF
328 IF (order > 3) THEN
329 cpabort("derivatives bigger than 3 not implemented")
330 END IF
331
332 CALL pz_lsd_calc(a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, &
333 ebbb, npoints, order, sc)
334
335 CALL timestop(timer_handle)
336
337 END SUBROUTINE pz_lsd_eval
338
339! **************************************************************************************************
340!> \brief ...
341!> \param a ...
342!> \param b ...
343!> \param e_0 ...
344!> \param ea ...
345!> \param eaa ...
346!> \param eab ...
347!> \param ebb ...
348!> \param eaaa ...
349!> \param eaab ...
350!> \param eabb ...
351!> \param ebbb ...
352!> \param npoints ...
353!> \param order ...
354!> \param sc ...
355! **************************************************************************************************
356 SUBROUTINE pz_lsd_calc(a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, &
357 ebbb, npoints, order, sc)
358 !FM low-level computation routine
359 REAL(kind=dp), DIMENSION(*), INTENT(in) :: a, b
360 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, ea, eaa, eab, ebb, eaaa, eaab, &
361 eabb, ebbb
362 INTEGER, INTENT(in) :: npoints, order
363 REAL(kind=dp), INTENT(IN) :: sc
364
365 INTEGER :: k, order_
366 REAL(kind=dp) :: rho
367 REAL(kind=dp), DIMENSION(0:9) :: ed
368
369 order_ = abs(order)
370
371!$OMP PARALLEL DO PRIVATE ( k, rho, ed ) DEFAULT(NONE)&
372!$OMP SHARED(order_,order,npoints,eps_rho,A,b,sc,e_0,ea,eaa,eab,ebb,eaaa,eaab,eabb,ebbb)
373 DO k = 1, npoints
374
375 rho = a(k) + b(k)
376
377 IF (rho > eps_rho) THEN
378
379 CALL pz_lsd_ed_loc(a(k), b(k), ed, order_, sc)
380
381 e_0(k) = e_0(k) + rho*ed(0)
382
383 IF (order >= 1) THEN
384 ea(k) = ea(k) + ed(0) + rho*ed(1)
385 ea(k) = ea(k) + ed(0) + rho*ed(2)
386 END IF
387
388 IF (order >= 2) THEN
389 eaa(k) = eaa(k) + 2.0_dp*ed(1) + rho*ed(3)
390 eab(k) = eab(k) + ed(1) + ed(2) + rho*ed(4)
391 ebb(k) = ebb(k) + 2.0_dp*ed(2) + rho*ed(5)
392 END IF
393
394 IF (order >= 3) THEN
395 eaaa(k) = eaaa(k) + 3.0_dp*ed(3) + rho*ed(6)
396 eaab(k) = eaab(k) + 2.0_dp*ed(4) + ed(3) + rho*ed(7)
397 eabb(k) = eabb(k) + 2.0_dp*ed(4) + ed(5) + rho*ed(8)
398 ebbb(k) = ebbb(k) + 3.0_dp*ed(5) + rho*ed(9)
399 END IF
400
401 END IF
402
403 END DO
404!$OMP END PARALLEL DO
405
406 END SUBROUTINE pz_lsd_calc
407
408! **************************************************************************************************
409!> \brief Initializes the functionals
410!> \param method CHARACTER(3) - name of the method used for parameters
411!> \param cutoff REAL(KIND=dp) - the cutoff density
412!> \par History
413!> 18-MAR-2002, TCH, working version
414!> \note see functionals_utilities
415! **************************************************************************************************
416 SUBROUTINE pz_init(method, cutoff)
417
418 INTEGER, INTENT(IN) :: method
419 REAL(kind=dp), INTENT(IN) :: cutoff
420
421 CALL set_util(cutoff)
422
423 eps_rho = cutoff
424
425 initialized = .false.
426
427 SELECT CASE (method)
428
429 CASE DEFAULT
430 cpabort("Unknown method")
431
432 CASE (pz_orig)
433 CALL cite_reference(perdew1981)
434 ga(0) = -0.1423_dp; ga(1) = -0.0843_dp
435 b1(0) = 1.0529_dp; b1(1) = 1.3981_dp
436 b2(0) = 0.3334_dp; b2(1) = 0.2611_dp
437 a(0) = 0.0311_dp; a(1) = 0.01555_dp
438 b(0) = -0.048_dp; b(1) = -0.0269_dp
439 c(0) = +0.0020_dp; c(1) = +0.0007_dp
440 d(0) = -0.0116_dp; d(1) = -0.0048_dp
441
442 CASE (pz_dmc)
443 CALL cite_reference(ortiz1994)
444 ga(0) = -0.103756_dp; ga(1) = -0.065951_dp
445 b1(0) = 0.56371_dp; b1(1) = 1.11846_dp
446 b2(0) = 0.27358_dp; b2(1) = 0.18797_dp
447 a(0) = 0.031091_dp; a(1) = 0.015545_dp
448 b(0) = -0.046644_dp; b(1) = -0.025599_dp
449 c(0) = -0.00419_dp; c(1) = -0.00329_dp
450 d(0) = -0.00983_dp; d(1) = -0.00300_dp
451
452 CASE (pz_vmc)
453 CALL cite_reference(ortiz1994)
454 ga(0) = -0.093662_dp; ga(1) = -0.055331_dp
455 b1(0) = 0.49453_dp; b1(1) = 0.93766_dp
456 b2(0) = 0.25534_dp; b2(1) = 0.14829_dp
457 a(0) = 0.031091_dp; a(1) = 0.015545_dp
458 b(0) = -0.046644_dp; b(1) = -0.025599_dp
459 c(0) = -0.00884_dp; c(1) = -0.00677_dp
460 d(0) = -0.00688_dp; d(1) = -0.00093_dp
461
462 END SELECT
463
464 initialized = .true.
465
466 END SUBROUTINE pz_init
467
468! **************************************************************************************************
469!> \brief ...
470!> \param r ...
471!> \param z ...
472!> \param g ...
473!> \param order ...
474! **************************************************************************************************
475 SUBROUTINE calc_g(r, z, g, order)
476
477 REAL(kind=dp), INTENT(IN) :: r
478 INTEGER, INTENT(IN) :: z
479 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: g
480 INTEGER, INTENT(IN) :: order
481
482 REAL(kind=dp) :: rr, rsr, sr
483
484 IF (r >= 1.0_dp) THEN
485
486 sr = sqrt(r)
487 ! order 0 must always be calculated
488 g(0) = ga(z)/(1.0_dp + b1(z)*sr + b2(z)*r)
489 IF (order >= 1) THEN
490 g(1) = -ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))/ &
491 (1.0_dp + b1(z)*sr + b2(z)*r)**2
492 END IF
493 IF (order >= 2) THEN
494 rsr = r*sr
495 g(2) = &
496 2.0_dp*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))**2 &
497 /(1.0_dp + b1(z)*sr + b2(z)*r)**3 &
498 + ga(z)*b1(z) &
499 /(4.0_dp*(1.0_dp + b1(z)*sr + b2(z)*r)**2*rsr)
500 END IF
501 IF (order >= 3) THEN
502 g(3) = &
503 -6.0_dp*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))**3/ &
504 (1.0_dp + b1(z)*sr + b2(z)*r)**4 &
505 - (3.0_dp/2.0_dp)*ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))*b1(z)/ &
506 ((1.0_dp + b1(z)*sr + b2(z)*r)**3*rsr) &
507 - (3.0_dp/8.0_dp)*ga(z)*b1(z)/ &
508 ((1.0_dp + b1(z)*sr + b2(z)*r)**2*r*rsr)
509 END IF
510
511 ELSE
512
513 ! order 0 must always be calculated
514 g(0) = a(z)*log(r) + b(z) + c(z)*r*log(r) + d(z)*r
515 IF (order >= 1) THEN
516 g(1) = a(z)/r + c(z)*log(r) + c(z) + d(z)
517 END IF
518 IF (order >= 2) THEN
519 rr = r*r
520 g(2) = -a(z)/rr + c(z)/r
521 END IF
522 IF (order >= 3) THEN
523 g(3) = 2.0_dp*a(z)/(rr*r) - c(z)/rr
524 END IF
525
526 END IF
527
528 END SUBROUTINE calc_g
529
530! **************************************************************************************************
531!> \brief ...
532!> \param rho ...
533!> \param ed ...
534!> \param order ...
535!> \param sc ...
536! **************************************************************************************************
537 SUBROUTINE pz_lda_ed_loc(rho, ed, order, sc)
538
539 REAL(kind=dp), INTENT(IN) :: rho
540 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
541 INTEGER, INTENT(IN) :: order
542 REAL(dp), INTENT(IN) :: sc
543
544 INTEGER :: m, order_
545 LOGICAL, DIMENSION(0:3) :: calc
546 REAL(kind=dp), DIMENSION(0:3) :: e0, r
547
548 calc = .false.
549
550 order_ = order
551 IF (order_ >= 0) THEN
552 calc(0:order_) = .true.
553 ELSE
554 order_ = -1*order_
555 calc(order_) = .true.
556 END IF
557
558 CALL calc_rs(rho, r(0))
559 CALL calc_g(r(0), 0, e0, order_)
560
561 IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
562 IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
563 IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
564
565 m = 0
566 IF (calc(0)) THEN
567 ed(m) = sc*e0(0)
568 m = m + 1
569 END IF
570 IF (calc(1)) THEN
571 ed(m) = sc*e0(1)*r(1)
572 m = m + 1
573 END IF
574 IF (calc(2)) THEN
575 ed(m) = sc*e0(2)*r(1)**2 + sc*e0(1)*r(2)
576 m = m + 1
577 END IF
578 IF (calc(3)) THEN
579 ed(m) = sc*e0(3)*r(1)**3 + sc*e0(2)*3.0_dp*r(1)*r(2) + sc*e0(1)*r(3)
580 END IF
581
582 END SUBROUTINE pz_lda_ed_loc
583
584! **************************************************************************************************
585!> \brief ...
586!> \param a ...
587!> \param b ...
588!> \param ed ...
589!> \param order ...
590!> \param sc ...
591! **************************************************************************************************
592 SUBROUTINE pz_lsd_ed_loc(a, b, ed, order, sc)
593
594 REAL(kind=dp), INTENT(IN) :: a, b
595 REAL(kind=dp), DIMENSION(0:), INTENT(OUT) :: ed
596 INTEGER, INTENT(IN), OPTIONAL :: order
597 REAL(kind=dp), INTENT(IN) :: sc
598
599 INTEGER :: m, order_
600 LOGICAL, DIMENSION(0:3) :: calc
601 REAL(kind=dp) :: rho, tr, trr, trrr, trrz, trz, trzz, tz, &
602 tzz, tzzz
603 REAL(kind=dp), DIMENSION(0:3) :: e0, e1, f, r
604 REAL(kind=dp), DIMENSION(0:3, 0:3) :: z
605
606 calc = .false.
607
608 order_ = 0
609 IF (PRESENT(order)) order_ = order
610
611 IF (order_ >= 0) THEN
612 calc(0:order_) = .true.
613 ELSE
614 order_ = -1*order_
615 calc(order_) = .true.
616 END IF
617
618 rho = a + b
619
620 CALL calc_fx(a, b, f(0:order_), order_)
621 CALL calc_rs(rho, r(0))
622 CALL calc_g(r(0), 0, e0(0:order_), order_)
623 CALL calc_g(r(0), 1, e1(0:order_), order_)
624 CALL calc_z(a, b, z(:, :), order_)
625
626!! calculate first partial derivatives
627 IF (order_ >= 1) THEN
628 r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
629 tr = e0(1) + (e1(1) - e0(1))*f(0)
630 tz = (e1(0) - e0(0))*f(1)
631 END IF
632
633!! calculate second partial derivatives
634 IF (order_ >= 2) THEN
635 r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
636 trr = e0(2) + (e1(2) - e0(2))*f(0)
637 trz = (e1(1) - e0(1))*f(1)
638 tzz = (e1(0) - e0(0))*f(2)
639 END IF
640
641!! calculate third derivatives
642 IF (order_ >= 3) THEN
643 r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
644 trrr = e0(3) + (e1(3) - e0(3))*f(0)
645 trrz = (e1(2) - e0(2))*f(1)
646 trzz = (e1(1) - e0(1))*f(2)
647 tzzz = (e1(0) - e0(0))*f(3)
648 END IF
649
650 m = 0
651 IF (calc(0)) THEN
652 ed(m) = e0(0) + (e1(0) - e0(0))*f(0)
653 ed(m) = ed(m)*sc
654 m = m + 1
655 END IF
656
657 IF (calc(1)) THEN
658 ed(m) = tr*r(1) + tz*z(1, 0)
659 ed(m) = ed(m)*sc
660 ed(m + 1) = tr*r(1) + tz*z(0, 1)
661 ed(m + 1) = ed(m + 1)*sc
662 m = m + 2
663 END IF
664
665 IF (calc(2)) THEN
666 ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
667 + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
668 ed(m) = ed(m)*sc
669 ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
670 + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
671 ed(m + 1) = ed(m + 1)*sc
672 ed(m + 2) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(0, 1) &
673 + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
674 ed(m + 2) = ed(m + 2)*sc
675 m = m + 3
676 END IF
677
678 IF (calc(3)) THEN
679 ed(m) = trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(1, 0) &
680 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(1, 0) &
681 + tr*r(3) + 3.0_dp*trzz*r(1)*z(1, 0)**2 &
682 + tzzz*z(1, 0)**3 + 3.0_dp*trz*r(1)*z(2, 0) &
683 + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
684 ed(m) = ed(m)*sc
685 ed(m + 1) = trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(1, 0) + z(0, 1)) &
686 + 2.0_dp*trzz*r(1)*z(1, 0)*z(0, 1) &
687 + 2.0_dp*trz*(r(2)*z(1, 0) + r(1)*z(1, 1)) &
688 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(0, 1) + tr*r(3) &
689 + trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**2*z(0, 1) &
690 + 2.0_dp*tzz*z(1, 0)*z(1, 1) &
691 + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
692 ed(m + 1) = ed(m + 1)*sc
693 ed(m + 2) = trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(0, 1) + z(1, 0)) &
694 + 2.0_dp*trzz*r(1)*z(0, 1)*z(1, 0) &
695 + 2.0_dp*trz*(r(2)*z(0, 1) + r(1)*z(1, 1)) &
696 + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(1, 0) + tr*r(3) &
697 + trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**2*z(1, 0) &
698 + 2.0_dp*tzz*z(0, 1)*z(1, 1) &
699 + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
700 ed(m + 2) = ed(m + 2)*sc
701 ed(m + 3) = trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(0, 1) &
702 + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(0, 1) + tr*r(3) &
703 + 3.0_dp*trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**3 &
704 + 3.0_dp*trz*r(1)*z(0, 2) &
705 + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
706 ed(m + 3) = ed(m + 3)*sc
707 END IF
708
709 END SUBROUTINE pz_lsd_ed_loc
710
711END MODULE xc_perdew_zunger
xc_functionals_utilities::calc_fx
Definition xc_functionals_utilities.F:37
xc_functionals_utilities::calc_rs
Definition xc_functionals_utilities.F:41
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::perdew1981
integer, save, public perdew1981
Definition bibliography.F:36
bibliography::ortiz1994
integer, save, public ortiz1994
Definition bibliography.F:36
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:1052
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:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition xc_functionals_utilities.F:15
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition xc_functionals_utilities.F:54
xc_functionals_utilities::calc_z
pure subroutine, public calc_z(a, b, z, order)
...
Definition xc_functionals_utilities.F:348
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::pz_dmc
integer, parameter, public pz_dmc
Definition xc_input_constants.F:92
xc_input_constants::pz_vmc
integer, parameter, public pz_vmc
Definition xc_input_constants.F:92
xc_input_constants::pz_orig
integer, parameter, public pz_orig
Definition xc_input_constants.F:92
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:135
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:261
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
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78