(git:374b731)
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-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! **************************************************************************************************
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
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 - 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
723END 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: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
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
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