Line data Source code
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 : ! **************************************************************************************************
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_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 :
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 39 : 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 39 : SELECT CASE (method)
79 : CASE DEFAULT
80 0 : CPABORT("Unsupported parametrization")
81 : CASE (pz_orig)
82 39 : p_string = 'ORIG'
83 : CASE (pz_dmc)
84 0 : p_string = 'DMC'
85 : CASE (pz_vmc)
86 0 : p_string = 'VMC'
87 : END SELECT
88 :
89 39 : IF (PRESENT(reference)) THEN
90 : reference = "J. P. Perdew and Alex Zunger," &
91 : //" Phys. Rev. B 23, 5048 (1981)" &
92 1 : //"["//TRIM(p_string)//"]"
93 1 : IF (.NOT. lsd) THEN
94 1 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
95 1 : reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
96 : END IF
97 : END IF
98 : END IF
99 39 : IF (PRESENT(shortform)) THEN
100 : shortform = "J. P. Perdew et al., PRB 23, 5048 (1981)" &
101 1 : //"["//TRIM(p_string)//"]"
102 1 : IF (.NOT. lsd) THEN
103 1 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
104 1 : shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
105 : END IF
106 : END IF
107 : END IF
108 39 : IF (PRESENT(needs)) THEN
109 38 : IF (lsd) THEN
110 6 : needs%rho_spin = .TRUE.
111 : ELSE
112 32 : needs%rho = .TRUE.
113 : END IF
114 : END IF
115 39 : IF (PRESENT(max_deriv)) max_deriv = 3
116 :
117 39 : 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 654 : 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 218 : POINTER :: dummy, e_0, e_rho, e_rho_rho, &
149 218 : e_rho_rho_rho, rho
150 : TYPE(xc_derivative_type), POINTER :: deriv
151 :
152 218 : CALL timeset(routineN, timer_handle)
153 218 : NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
154 :
155 218 : CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
156 :
157 : CALL xc_rho_set_get(rho_set, rho=rho, &
158 218 : local_bounds=bo, rho_cutoff=rho_cutoff)
159 218 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
160 :
161 218 : CALL pz_init(method, rho_cutoff)
162 :
163 218 : dummy => rho
164 :
165 218 : e_0 => dummy
166 218 : e_rho => dummy
167 218 : e_rho_rho => dummy
168 218 : e_rho_rho_rho => dummy
169 :
170 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
171 218 : allocate_deriv=.TRUE.)
172 218 : CALL xc_derivative_get(deriv, deriv_data=e_0)
173 218 : IF (order >= 1) THEN
174 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
175 212 : allocate_deriv=.TRUE.)
176 212 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
177 : END IF
178 218 : IF (order >= 2) THEN
179 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
180 0 : allocate_deriv=.TRUE.)
181 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
182 : END IF
183 218 : IF (order >= 3) THEN
184 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
185 0 : allocate_deriv=.TRUE.)
186 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
187 : END IF
188 218 : IF (order > 3) THEN
189 0 : CPABORT("derivatives bigger than 3 not implemented")
190 : END IF
191 :
192 218 : CALL pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
193 :
194 218 : CALL timestop(timer_handle)
195 :
196 218 : 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 218 : 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 218 : !$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 218 : 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 36 : 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 12 : POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
274 12 : eab, eabb, ebb, ebbb
275 : TYPE(xc_derivative_type), POINTER :: deriv
276 :
277 12 : CALL timeset(routineN, timer_handle)
278 12 : NULLIFY (a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
279 :
280 12 : 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 12 : local_bounds=bo, rho_cutoff=rho_cutoff)
284 12 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
285 :
286 12 : CALL pz_init(method, rho_cutoff)
287 :
288 12 : dummy => a
289 :
290 12 : e_0 => dummy
291 12 : ea => dummy
292 12 : eaa => dummy; eab => dummy; ebb => dummy
293 12 : eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
294 :
295 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
296 12 : allocate_deriv=.TRUE.)
297 12 : CALL xc_derivative_get(deriv, deriv_data=e_0)
298 12 : IF (order >= 1) THEN
299 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
300 6 : allocate_deriv=.TRUE.)
301 6 : CALL xc_derivative_get(deriv, deriv_data=ea)
302 : END IF
303 12 : IF (order >= 2) THEN
304 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
305 0 : allocate_deriv=.TRUE.)
306 0 : CALL xc_derivative_get(deriv, deriv_data=eaa)
307 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
308 0 : allocate_deriv=.TRUE.)
309 0 : CALL xc_derivative_get(deriv, deriv_data=eab)
310 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
311 0 : allocate_deriv=.TRUE.)
312 0 : CALL xc_derivative_get(deriv, deriv_data=ebb)
313 : END IF
314 12 : IF (order >= 3) THEN
315 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
316 0 : allocate_deriv=.TRUE.)
317 0 : CALL xc_derivative_get(deriv, deriv_data=eaaa)
318 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
319 0 : allocate_deriv=.TRUE.)
320 0 : CALL xc_derivative_get(deriv, deriv_data=eaab)
321 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
322 0 : allocate_deriv=.TRUE.)
323 0 : CALL xc_derivative_get(deriv, deriv_data=eabb)
324 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
325 0 : allocate_deriv=.TRUE.)
326 0 : CALL xc_derivative_get(deriv, deriv_data=ebbb)
327 : END IF
328 12 : IF (order > 3) THEN
329 0 : 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 12 : ebbb, npoints, order, sc)
334 :
335 12 : CALL timestop(timer_handle)
336 :
337 12 : 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 12 : 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 12 : order_ = ABS(order)
370 :
371 : !$OMP PARALLEL DO PRIVATE ( k, rho, ed ) DEFAULT(NONE)&
372 12 : !$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 12 : 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 230 : SUBROUTINE pz_init(method, cutoff)
417 :
418 : INTEGER, INTENT(IN) :: method
419 : REAL(KIND=dp), INTENT(IN) :: cutoff
420 :
421 230 : CALL set_util(cutoff)
422 :
423 230 : eps_rho = cutoff
424 :
425 230 : initialized = .FALSE.
426 :
427 230 : SELECT CASE (method)
428 :
429 : CASE DEFAULT
430 0 : CPABORT("Unknown method")
431 :
432 : CASE (pz_orig)
433 230 : CALL cite_reference(Perdew1981)
434 230 : ga(0) = -0.1423_dp; ga(1) = -0.0843_dp
435 230 : b1(0) = 1.0529_dp; b1(1) = 1.3981_dp
436 230 : b2(0) = 0.3334_dp; b2(1) = 0.2611_dp
437 230 : A(0) = 0.0311_dp; A(1) = 0.01555_dp
438 230 : B(0) = -0.048_dp; B(1) = -0.0269_dp
439 230 : C(0) = +0.0020_dp; C(1) = +0.0007_dp
440 230 : D(0) = -0.0116_dp; D(1) = -0.0048_dp
441 :
442 : CASE (pz_dmc)
443 0 : CALL cite_reference(Ortiz1994)
444 0 : ga(0) = -0.103756_dp; ga(1) = -0.065951_dp
445 0 : b1(0) = 0.56371_dp; b1(1) = 1.11846_dp
446 0 : b2(0) = 0.27358_dp; b2(1) = 0.18797_dp
447 0 : A(0) = 0.031091_dp; A(1) = 0.015545_dp
448 0 : B(0) = -0.046644_dp; B(1) = -0.025599_dp
449 0 : C(0) = -0.00419_dp; C(1) = -0.00329_dp
450 0 : D(0) = -0.00983_dp; D(1) = -0.00300_dp
451 :
452 : CASE (pz_vmc)
453 0 : CALL cite_reference(Ortiz1994)
454 0 : ga(0) = -0.093662_dp; ga(1) = -0.055331_dp
455 0 : b1(0) = 0.49453_dp; b1(1) = 0.93766_dp
456 0 : b2(0) = 0.25534_dp; b2(1) = 0.14829_dp
457 0 : A(0) = 0.031091_dp; A(1) = 0.015545_dp
458 0 : B(0) = -0.046644_dp; B(1) = -0.025599_dp
459 0 : C(0) = -0.00884_dp; C(1) = -0.00677_dp
460 230 : D(0) = -0.00688_dp; D(1) = -0.00093_dp
461 :
462 : END SELECT
463 :
464 230 : initialized = .TRUE.
465 :
466 230 : END SUBROUTINE pz_init
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param r ...
471 : !> \param z ...
472 : !> \param g ...
473 : !> \param order ...
474 : ! **************************************************************************************************
475 37362837 : 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 37362837 : IF (r >= 1.0_dp) THEN
485 :
486 37265424 : sr = SQRT(r)
487 : ! order 0 must always be calculated
488 37265424 : g(0) = ga(z)/(1.0_dp + b1(z)*sr + b2(z)*r)
489 37265424 : IF (order >= 1) THEN
490 : g(1) = -ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))/ &
491 36869492 : (1.0_dp + b1(z)*sr + b2(z)*r)**2
492 : END IF
493 37265424 : IF (order >= 2) THEN
494 0 : 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 0 : /(4.0_dp*(1.0_dp + b1(z)*sr + b2(z)*r)**2*rsr)
500 : END IF
501 37265424 : 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 0 : ((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 97413 : g(0) = A(z)*LOG(r) + B(z) + C(z)*r*LOG(r) + D(z)*r
515 97413 : IF (order >= 1) THEN
516 85783 : g(1) = A(z)/r + C(z)*LOG(r) + C(z) + D(z)
517 : END IF
518 97413 : IF (order >= 2) THEN
519 0 : rr = r*r
520 0 : g(2) = -A(z)/rr + C(z)/r
521 : END IF
522 97413 : IF (order >= 3) THEN
523 0 : g(3) = 2.0_dp*A(z)/(rr*r) - C(z)/rr
524 : END IF
525 :
526 : END IF
527 :
528 37362837 : END SUBROUTINE calc_g
529 :
530 : ! **************************************************************************************************
531 : !> \brief ...
532 : !> \param rho ...
533 : !> \param ed ...
534 : !> \param order ...
535 : !> \param sc ...
536 : ! **************************************************************************************************
537 36664213 : 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 36664213 : calc = .FALSE.
549 :
550 36664213 : order_ = order
551 36664213 : IF (order_ >= 0) THEN
552 109934389 : calc(0:order_) = .TRUE.
553 : ELSE
554 0 : order_ = -1*order_
555 0 : calc(order_) = .TRUE.
556 : END IF
557 :
558 36664213 : CALL calc_rs(rho, r(0))
559 36664213 : CALL calc_g(r(0), 0, e0, order_)
560 :
561 36664213 : IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
562 36664213 : IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
563 36664213 : IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
564 :
565 36664213 : m = 0
566 36664213 : IF (calc(0)) THEN
567 36664213 : ed(m) = sc*e0(0)
568 36664213 : m = m + 1
569 : END IF
570 36664213 : IF (calc(1)) THEN
571 36605963 : ed(m) = sc*e0(1)*r(1)
572 36605963 : m = m + 1
573 : END IF
574 36664213 : IF (calc(2)) THEN
575 0 : ed(m) = sc*e0(2)*r(1)**2 + sc*e0(1)*r(2)
576 0 : m = m + 1
577 : END IF
578 36664213 : IF (calc(3)) THEN
579 0 : 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 36664213 : 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 349312 : 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 349312 : calc = .FALSE.
607 :
608 349312 : order_ = 0
609 349312 : IF (PRESENT(order)) order_ = order
610 :
611 349312 : IF (order_ >= 0) THEN
612 873280 : calc(0:order_) = .TRUE.
613 : ELSE
614 0 : order_ = -1*order_
615 0 : calc(order_) = .TRUE.
616 : END IF
617 :
618 349312 : rho = a + b
619 :
620 349312 : CALL calc_fx(a, b, f(0:order_), order_)
621 349312 : CALL calc_rs(rho, r(0))
622 349312 : CALL calc_g(r(0), 0, e0(0:order_), order_)
623 349312 : CALL calc_g(r(0), 1, e1(0:order_), order_)
624 349312 : CALL calc_z(a, b, z(:, :), order_)
625 :
626 : !! calculate first partial derivatives
627 349312 : IF (order_ >= 1) THEN
628 174656 : r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
629 174656 : tr = e0(1) + (e1(1) - e0(1))*f(0)
630 174656 : tz = (e1(0) - e0(0))*f(1)
631 : END IF
632 :
633 : !! calculate second partial derivatives
634 349312 : IF (order_ >= 2) THEN
635 0 : r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
636 0 : trr = e0(2) + (e1(2) - e0(2))*f(0)
637 0 : trz = (e1(1) - e0(1))*f(1)
638 0 : tzz = (e1(0) - e0(0))*f(2)
639 : END IF
640 :
641 : !! calculate third derivatives
642 349312 : IF (order_ >= 3) THEN
643 0 : r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
644 0 : trrr = e0(3) + (e1(3) - e0(3))*f(0)
645 0 : trrz = (e1(2) - e0(2))*f(1)
646 0 : trzz = (e1(1) - e0(1))*f(2)
647 0 : tzzz = (e1(0) - e0(0))*f(3)
648 : END IF
649 :
650 349312 : m = 0
651 349312 : IF (calc(0)) THEN
652 349312 : ed(m) = e0(0) + (e1(0) - e0(0))*f(0)
653 349312 : ed(m) = ed(m)*sc
654 349312 : m = m + 1
655 : END IF
656 :
657 349312 : IF (calc(1)) THEN
658 174656 : ed(m) = tr*r(1) + tz*z(1, 0)
659 174656 : ed(m) = ed(m)*sc
660 174656 : ed(m + 1) = tr*r(1) + tz*z(0, 1)
661 174656 : ed(m + 1) = ed(m + 1)*sc
662 174656 : m = m + 2
663 : END IF
664 :
665 349312 : IF (calc(2)) THEN
666 : ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
667 0 : + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
668 0 : ed(m) = ed(m)*sc
669 : ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
670 0 : + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
671 0 : 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 0 : + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
674 0 : ed(m + 2) = ed(m + 2)*sc
675 0 : m = m + 3
676 : END IF
677 :
678 349312 : 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 0 : + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
684 0 : 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 0 : + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
692 0 : 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 0 : + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
700 0 : 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 0 : + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
706 0 : ed(m + 3) = ed(m + 3)*sc
707 : END IF
708 :
709 349312 : END SUBROUTINE pz_lsd_ed_loc
710 :
711 : END MODULE xc_perdew_zunger
|
