Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 41 : 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 41 : SELECT CASE (method)
79 : CASE DEFAULT
80 0 : CPABORT("Unsupported parametrization")
81 : CASE (pz_orig)
82 41 : 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 41 : 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 41 : 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 41 : IF (PRESENT(needs)) THEN
109 40 : IF (lsd) THEN
110 6 : needs%rho_spin = .TRUE.
111 : ELSE
112 34 : needs%rho = .TRUE.
113 : END IF
114 : END IF
115 41 : IF (PRESENT(max_deriv)) max_deriv = 3
116 :
117 41 : 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 138 : 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 46 : POINTER :: dummy, e_0, e_rho, e_rho_rho, &
151 46 : e_rho_rho_rho, rho
152 : TYPE(xc_derivative_type), POINTER :: deriv
153 :
154 46 : CALL timeset(routineN, timer_handle)
155 46 : NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
156 :
157 46 : CALL section_vals_val_get(pz_params, "scale_c", r_val=sc)
158 :
159 : CALL xc_rho_set_get(rho_set, rho=rho, &
160 46 : local_bounds=bo, rho_cutoff=rho_cutoff)
161 46 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
162 :
163 46 : CALL pz_init(method, rho_cutoff)
164 :
165 46 : dummy => rho
166 :
167 46 : e_0 => dummy
168 46 : e_rho => dummy
169 46 : e_rho_rho => dummy
170 46 : e_rho_rho_rho => dummy
171 :
172 46 : IF (order >= 0) THEN
173 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
174 46 : allocate_deriv=.TRUE.)
175 46 : CALL xc_derivative_get(deriv, deriv_data=e_0)
176 : END IF
177 46 : IF (order >= 1 .OR. order == -1) THEN
178 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
179 38 : allocate_deriv=.TRUE.)
180 38 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
181 : END IF
182 46 : IF (order >= 2 .OR. order == -2) THEN
183 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
184 0 : allocate_deriv=.TRUE.)
185 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
186 : END IF
187 46 : IF (order >= 3 .OR. order == -3) THEN
188 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
189 0 : allocate_deriv=.TRUE.)
190 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
191 : END IF
192 46 : IF (order > 3 .OR. order < -3) THEN
193 0 : CPABORT("derivatives bigger than 3 not implemented")
194 : END IF
195 :
196 46 : CALL pz_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, sc)
197 :
198 46 : CALL timestop(timer_handle)
199 :
200 46 : 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 46 : 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 46 : !$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 46 : 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 36 : 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 12 : POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
282 12 : eab, eabb, ebb, ebbb
283 : TYPE(xc_derivative_type), POINTER :: deriv
284 :
285 12 : CALL timeset(routineN, timer_handle)
286 12 : NULLIFY (a, b, e_0, ea, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
287 :
288 12 : 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 12 : local_bounds=bo, rho_cutoff=rho_cutoff)
292 12 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
293 :
294 12 : CALL pz_init(method, rho_cutoff)
295 :
296 12 : dummy => a
297 :
298 12 : e_0 => dummy
299 12 : ea => dummy;
300 12 : eaa => dummy; eab => dummy; ebb => dummy
301 12 : eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
302 :
303 12 : IF (order >= 0) THEN
304 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
305 12 : allocate_deriv=.TRUE.)
306 12 : CALL xc_derivative_get(deriv, deriv_data=e_0)
307 : END IF
308 12 : IF (order >= 1 .OR. order == -1) THEN
309 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
310 6 : allocate_deriv=.TRUE.)
311 6 : CALL xc_derivative_get(deriv, deriv_data=ea)
312 : END IF
313 12 : IF (order >= 2 .OR. order == -2) THEN
314 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
315 0 : allocate_deriv=.TRUE.)
316 0 : CALL xc_derivative_get(deriv, deriv_data=eaa)
317 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
318 0 : allocate_deriv=.TRUE.)
319 0 : CALL xc_derivative_get(deriv, deriv_data=eab)
320 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
321 0 : allocate_deriv=.TRUE.)
322 0 : CALL xc_derivative_get(deriv, deriv_data=ebb)
323 : END IF
324 12 : IF (order >= 3 .OR. order == -3) THEN
325 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
326 0 : allocate_deriv=.TRUE.)
327 0 : CALL xc_derivative_get(deriv, deriv_data=eaaa)
328 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
329 0 : allocate_deriv=.TRUE.)
330 0 : CALL xc_derivative_get(deriv, deriv_data=eaab)
331 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
332 0 : allocate_deriv=.TRUE.)
333 0 : CALL xc_derivative_get(deriv, deriv_data=eabb)
334 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
335 0 : allocate_deriv=.TRUE.)
336 0 : CALL xc_derivative_get(deriv, deriv_data=ebbb)
337 : END IF
338 12 : IF (order > 3 .OR. order < -3) THEN
339 0 : 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 12 : ebbb, npoints, order, sc)
344 :
345 12 : CALL timestop(timer_handle)
346 :
347 12 : 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 12 : 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 12 : order_ = ABS(order)
380 :
381 : !$OMP PARALLEL DO PRIVATE ( k, rho, ed ) DEFAULT(NONE)&
382 12 : !$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 12 : 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 58 : SUBROUTINE pz_init(method, cutoff)
429 :
430 : INTEGER, INTENT(IN) :: method
431 : REAL(KIND=dp), INTENT(IN) :: cutoff
432 :
433 58 : CALL set_util(cutoff)
434 :
435 58 : eps_rho = cutoff
436 :
437 58 : initialized = .FALSE.
438 :
439 58 : SELECT CASE (method)
440 :
441 : CASE DEFAULT
442 0 : CPABORT("Unknown method")
443 :
444 : CASE (pz_orig)
445 58 : CALL cite_reference(Perdew1981)
446 58 : ga(0) = -0.1423_dp; ga(1) = -0.0843_dp
447 58 : b1(0) = 1.0529_dp; b1(1) = 1.3981_dp
448 58 : b2(0) = 0.3334_dp; b2(1) = 0.2611_dp
449 58 : A(0) = 0.0311_dp; A(1) = 0.01555_dp
450 58 : B(0) = -0.048_dp; B(1) = -0.0269_dp
451 58 : C(0) = +0.0020_dp; C(1) = +0.0007_dp
452 58 : D(0) = -0.0116_dp; D(1) = -0.0048_dp
453 :
454 : CASE (pz_dmc)
455 0 : CALL cite_reference(Ortiz1994)
456 0 : ga(0) = -0.103756_dp; ga(1) = -0.065951_dp
457 0 : b1(0) = 0.56371_dp; b1(1) = 1.11846_dp
458 0 : b2(0) = 0.27358_dp; b2(1) = 0.18797_dp
459 0 : A(0) = 0.031091_dp; A(1) = 0.015545_dp
460 0 : B(0) = -0.046644_dp; B(1) = -0.025599_dp
461 0 : C(0) = -0.00419_dp; C(1) = -0.00329_dp
462 0 : D(0) = -0.00983_dp; D(1) = -0.00300_dp
463 :
464 : CASE (pz_vmc)
465 0 : CALL cite_reference(Ortiz1994)
466 0 : ga(0) = -0.093662_dp; ga(1) = -0.055331_dp
467 0 : b1(0) = 0.49453_dp; b1(1) = 0.93766_dp
468 0 : b2(0) = 0.25534_dp; b2(1) = 0.14829_dp
469 0 : A(0) = 0.031091_dp; A(1) = 0.015545_dp
470 0 : B(0) = -0.046644_dp; B(1) = -0.025599_dp
471 0 : C(0) = -0.00884_dp; C(1) = -0.00677_dp
472 58 : D(0) = -0.00688_dp; D(1) = -0.00093_dp
473 :
474 : END SELECT
475 :
476 58 : initialized = .TRUE.
477 :
478 58 : END SUBROUTINE pz_init
479 :
480 : ! **************************************************************************************************
481 : !> \brief ...
482 : !> \param r ...
483 : !> \param z ...
484 : !> \param g ...
485 : !> \param order ...
486 : ! **************************************************************************************************
487 917427 : 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 917427 : IF (r >= 1.0_dp) THEN
497 :
498 885016 : sr = SQRT(r)
499 : ! order 0 must always be calculated
500 885016 : g(0) = ga(z)/(1.0_dp + b1(z)*sr + b2(z)*r)
501 885016 : IF (order >= 1) THEN
502 : g(1) = -ga(z)*(b1(z)/(2.0_dp*sr) + b2(z))/ &
503 473769 : (1.0_dp + b1(z)*sr + b2(z)*r)**2
504 : END IF
505 885016 : IF (order >= 2) THEN
506 0 : 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 0 : /(4.0_dp*(1.0_dp + b1(z)*sr + b2(z)*r)**2*rsr)
512 : END IF
513 885016 : 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 0 : ((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 32411 : g(0) = A(z)*LOG(r) + B(z) + C(z)*r*LOG(r) + D(z)*r
527 32411 : IF (order >= 1) THEN
528 20471 : g(1) = A(z)/r + C(z)*LOG(r) + C(z) + D(z)
529 : END IF
530 32411 : IF (order >= 2) THEN
531 0 : rr = r*r
532 0 : g(2) = -A(z)/rr + C(z)/r
533 : END IF
534 32411 : IF (order >= 3) THEN
535 0 : g(3) = 2.0_dp*A(z)/(rr*r) - C(z)/rr
536 : END IF
537 :
538 : END IF
539 :
540 917427 : END SUBROUTINE calc_g
541 :
542 : ! **************************************************************************************************
543 : !> \brief ...
544 : !> \param rho ...
545 : !> \param ed ...
546 : !> \param order ...
547 : !> \param sc ...
548 : ! **************************************************************************************************
549 218803 : 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 218803 : calc = .FALSE.
561 :
562 218803 : order_ = order
563 218803 : IF (order_ >= 0) THEN
564 582534 : calc(0:order_) = .TRUE.
565 : ELSE
566 0 : order_ = -1*order_
567 0 : calc(order_) = .TRUE.
568 : END IF
569 :
570 218803 : CALL calc_rs(rho, r(0))
571 218803 : CALL calc_g(r(0), 0, e0, order_)
572 :
573 218803 : IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
574 218803 : IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
575 218803 : IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
576 :
577 218803 : m = 0
578 218803 : IF (calc(0)) THEN
579 218803 : ed(m) = sc*e0(0)
580 218803 : m = m + 1
581 : END IF
582 218803 : IF (calc(1)) THEN
583 144928 : ed(m) = sc*e0(1)*r(1)
584 144928 : m = m + 1
585 : END IF
586 218803 : IF (calc(2)) THEN
587 0 : ed(m) = sc*e0(2)*r(1)**2 + sc*e0(1)*r(2)
588 0 : m = m + 1
589 : END IF
590 218803 : IF (calc(3)) THEN
591 0 : 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 218803 : 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 349312 : 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 349312 : calc = .FALSE.
619 :
620 349312 : order_ = 0
621 349312 : IF (PRESENT(order)) order_ = order
622 :
623 349312 : IF (order_ >= 0) THEN
624 873280 : calc(0:order_) = .TRUE.
625 : ELSE
626 0 : order_ = -1*order_
627 0 : calc(order_) = .TRUE.
628 : END IF
629 :
630 349312 : rho = a + b
631 :
632 349312 : CALL calc_fx(a, b, f(0:order_), order_)
633 349312 : CALL calc_rs(rho, r(0))
634 349312 : CALL calc_g(r(0), 0, e0(0:order_), order_)
635 349312 : CALL calc_g(r(0), 1, e1(0:order_), order_)
636 349312 : CALL calc_z(a, b, z(:, :), order_)
637 :
638 : !! calculate first partial derivatives
639 349312 : IF (order_ >= 1) THEN
640 174656 : r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
641 174656 : tr = e0(1) + (e1(1) - e0(1))*f(0)
642 174656 : tz = (e1(0) - e0(0))*f(1)
643 : END IF
644 :
645 : !! calculate second partial derivatives
646 349312 : IF (order_ >= 2) THEN
647 0 : r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
648 0 : trr = e0(2) + (e1(2) - e0(2))*f(0)
649 0 : trz = (e1(1) - e0(1))*f(1)
650 0 : tzz = (e1(0) - e0(0))*f(2)
651 : END IF
652 :
653 : !! calculate third derivatives
654 349312 : IF (order_ >= 3) THEN
655 0 : r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
656 0 : trrr = e0(3) + (e1(3) - e0(3))*f(0)
657 0 : trrz = (e1(2) - e0(2))*f(1)
658 0 : trzz = (e1(1) - e0(1))*f(2)
659 0 : tzzz = (e1(0) - e0(0))*f(3)
660 : END IF
661 :
662 349312 : m = 0
663 349312 : IF (calc(0)) THEN
664 349312 : ed(m) = e0(0) + (e1(0) - e0(0))*f(0)
665 349312 : ed(m) = ed(m)*sc
666 349312 : m = m + 1
667 : END IF
668 :
669 349312 : IF (calc(1)) THEN
670 174656 : ed(m) = tr*r(1) + tz*z(1, 0)
671 174656 : ed(m) = ed(m)*sc
672 174656 : ed(m + 1) = tr*r(1) + tz*z(0, 1)
673 174656 : ed(m + 1) = ed(m + 1)*sc
674 174656 : m = m + 2
675 : END IF
676 :
677 349312 : IF (calc(2)) THEN
678 : ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
679 0 : + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
680 0 : ed(m) = ed(m)*sc
681 : ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
682 0 : + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
683 0 : 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 0 : + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
686 0 : ed(m + 2) = ed(m + 2)*sc
687 0 : m = m + 3
688 : END IF
689 :
690 349312 : 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 0 : + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
696 0 : 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 0 : + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
704 0 : 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 0 : + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
712 0 : 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 0 : + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
718 0 : ed(m + 3) = ed(m + 3)*sc
719 : END IF
720 :
721 349312 : END SUBROUTINE pz_lsd_ed_loc
722 :
723 : END MODULE xc_perdew_zunger
|
