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-Wang 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_wang
18 : #:include "xc_perdew_wang.fypp"
19 :
20 : USE kinds, ONLY: dp
21 : USE pw_types, ONLY: pw_r3d_rs_type
22 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
23 : xc_dset_get_derivative
24 : USE xc_derivative_types, ONLY: xc_derivative_get, &
25 : xc_derivative_type
26 : USE xc_functionals_utilities, ONLY: calc_fx, &
27 : calc_rs, &
28 : calc_z, &
29 : set_util
30 : USE xc_input_constants, ONLY: pw_dmc, &
31 : pw_orig, &
32 : pw_vmc
33 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
34 : USE xc_rho_set_types, ONLY: xc_rho_set_get, &
35 : xc_rho_set_type
36 : USE xc_derivative_desc, ONLY: deriv_rho, &
37 : deriv_rhoa, &
38 : deriv_rhob
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 : PRIVATE
43 :
44 : @: global_var_pw92()
45 : REAL(KIND=dp), PARAMETER :: &
46 : epsilon = 5.E-13_dp, &
47 : fpp = 0.584822362263464620726223866376013788782_dp ! d^2f(0)/dz^2
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_perdew_wang'
49 :
50 : PUBLIC :: perdew_wang_info, perdew_wang_lda_eval, perdew_wang_lsd_eval, perdew_wang_fxc_calc
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Return some info on the functionals.
56 : !> \param method ...
57 : !> \param lsd ...
58 : !> \param reference full reference
59 : !> \param shortform short reference
60 : !> \param needs ...
61 : !> \param max_deriv ...
62 : !> \param scale ...
63 : ! **************************************************************************************************
64 337 : SUBROUTINE perdew_wang_info(method, lsd, reference, shortform, needs, &
65 : max_deriv, scale)
66 : INTEGER, INTENT(in) :: method
67 : LOGICAL, INTENT(in) :: lsd
68 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
69 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
70 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
71 : REAL(kind=dp), INTENT(in) :: scale
72 :
73 : CHARACTER(len=3) :: p_string
74 :
75 337 : SELECT CASE (method)
76 : CASE DEFAULT
77 0 : CPABORT("Unsupported parametrization")
78 : CASE (pw_orig)
79 337 : p_string = 'PWO'
80 : CASE (pw_dmc)
81 0 : p_string = 'DMC'
82 : CASE (pw_vmc)
83 0 : p_string = 'VMC'
84 : END SELECT
85 :
86 337 : IF (PRESENT(reference)) THEN
87 : reference = "J. P. Perdew and Yue Wang," &
88 : //" Phys. Rev. B 45, 13244 (1992)" &
89 11 : //"["//TRIM(p_string)//"]"
90 11 : IF (scale /= 1._dp) THEN
91 : WRITE (reference(LEN_TRIM(reference) + 1:LEN(reference)), "('s=',f5.3)") &
92 0 : scale
93 : END IF
94 11 : IF (.NOT. lsd) THEN
95 11 : IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
96 11 : reference(LEN_TRIM(reference) + 1:LEN_TRIM(reference) + 7) = ' {LDA}'
97 : END IF
98 : END IF
99 : END IF
100 337 : IF (PRESENT(shortform)) THEN
101 : shortform = "J. P. Perdew et al., PRB 45, 13244 (1992)" &
102 11 : //"["//TRIM(p_string)//"]"
103 11 : IF (scale /= 1._dp) THEN
104 : WRITE (shortform(LEN_TRIM(shortform) + 1:LEN(shortform)), "('s=',f5.3)") &
105 0 : scale
106 : END IF
107 11 : IF (.NOT. lsd) THEN
108 11 : IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
109 11 : shortform(LEN_TRIM(shortform) + 1:LEN_TRIM(shortform) + 7) = ' {LDA}'
110 : END IF
111 : END IF
112 : END IF
113 337 : IF (PRESENT(needs)) THEN
114 326 : IF (lsd) THEN
115 24 : needs%rho_spin = .TRUE.
116 : ELSE
117 302 : needs%rho = .TRUE.
118 : END IF
119 : END IF
120 337 : IF (PRESENT(max_deriv)) max_deriv = 3
121 :
122 337 : END SUBROUTINE perdew_wang_info
123 :
124 1250 : @: init_pw92()
125 :
126 : ! **************************************************************************************************
127 : !> \brief Calculate the correlation energy and its derivatives
128 : !> wrt to rho (the electron density) up to 3rd order. This
129 : !> is the LDA version of the Perdew-Wang correlation energy
130 : !> If no order argument is given, then the routine calculates
131 : !> just the energy.
132 : !> \param method ...
133 : !> \param rho_set ...
134 : !> \param deriv_set ...
135 : !> \param order maximum order of derivatives to calculate
136 : !> order must lie between 0 and 3.
137 : !> \param scale ...
138 : ! **************************************************************************************************
139 440 : SUBROUTINE perdew_wang_lda_eval(method, rho_set, deriv_set, order, scale)
140 :
141 : INTEGER, INTENT(in) :: method
142 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
143 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
144 : INTEGER, INTENT(in) :: order
145 : REAL(kind=dp), INTENT(in) :: scale
146 :
147 : CHARACTER(len=*), PARAMETER :: routineN = 'perdew_wang_lda_eval'
148 :
149 : INTEGER :: npoints, timer_handle
150 : INTEGER, DIMENSION(2, 3) :: bo
151 : REAL(KIND=dp) :: rho_cutoff
152 220 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_rho, e_rho_rho, &
153 220 : e_rho_rho_rho, rho
154 : TYPE(xc_derivative_type), POINTER :: deriv
155 :
156 220 : CALL timeset(routineN, timer_handle)
157 220 : NULLIFY (rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, dummy)
158 : CALL xc_rho_set_get(rho_set, rho=rho, &
159 220 : local_bounds=bo, rho_cutoff=rho_cutoff)
160 220 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
161 :
162 220 : CALL perdew_wang_init(method, rho_cutoff)
163 :
164 220 : dummy => rho
165 :
166 220 : e_0 => dummy
167 220 : e_rho => dummy
168 220 : e_rho_rho => dummy
169 220 : e_rho_rho_rho => dummy
170 :
171 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
172 220 : allocate_deriv=.TRUE.)
173 220 : CALL xc_derivative_get(deriv, deriv_data=e_0)
174 220 : IF (order >= 1) THEN
175 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
176 218 : allocate_deriv=.TRUE.)
177 218 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
178 : END IF
179 220 : IF (order >= 2) THEN
180 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
181 0 : allocate_deriv=.TRUE.)
182 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
183 : END IF
184 220 : IF (order >= 3) THEN
185 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
186 0 : allocate_deriv=.TRUE.)
187 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
188 : END IF
189 220 : IF (order > 3) THEN
190 0 : CPABORT("derivatives bigger than 3 not implemented")
191 : END IF
192 :
193 : CALL perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
194 220 : npoints, order, scale)
195 :
196 220 : CALL timestop(timer_handle)
197 :
198 220 : END SUBROUTINE perdew_wang_lda_eval
199 :
200 : ! **************************************************************************************************
201 : !> \brief ...
202 : !> \param rho ...
203 : !> \param e_0 ...
204 : !> \param e_rho ...
205 : !> \param e_rho_rho ...
206 : !> \param e_rho_rho_rho ...
207 : !> \param npoints ...
208 : !> \param order ...
209 : !> \param scale ...
210 : ! **************************************************************************************************
211 220 : SUBROUTINE perdew_wang_lda_calc(rho, e_0, e_rho, e_rho_rho, e_rho_rho_rho, npoints, order, scale)
212 : !FM low level calc routine
213 : REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rho
214 : REAL(KIND=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
215 : INTEGER, INTENT(in) :: npoints, order
216 : REAL(kind=dp), INTENT(in) :: scale
217 :
218 : INTEGER :: k
219 : REAL(KIND=dp), DIMENSION(0:3) :: ed
220 :
221 : !$OMP PARALLEL DO PRIVATE (k, ed) DEFAULT(NONE)&
222 220 : !$OMP SHARED(npoints,rho,eps_rho,scale,e_0,e_rho,e_rho_rho,e_rho_rho_rho,order)
223 : DO k = 1, npoints
224 :
225 : IF (rho(k) > eps_rho) THEN
226 : !! order_ is positive as it must be in this case:
227 : !! ec(:,2) needs ed(:,1) for example
228 : CALL pw_lda_ed_loc(rho(k), ed, order)
229 : ed(0:order) = scale*ed(0:order)
230 :
231 : e_0(k) = e_0(k) + rho(k)*ed(0)
232 : IF (order >= 1) THEN
233 : e_rho(k) = e_rho(k) + ed(0) + rho(k)*ed(1)
234 : END IF
235 : IF (order >= 2) THEN
236 : e_rho_rho(k) = e_rho_rho(k) + 2.0_dp*ed(1) + rho(k)*ed(2)
237 : END IF
238 : IF (order >= 3) THEN
239 : e_rho_rho_rho(k) = e_rho_rho_rho(k) + 3.0_dp*ed(2) + rho(k)*ed(3)
240 : END IF
241 :
242 : END IF
243 :
244 : END DO
245 : !$OMP END PARALLEL DO
246 :
247 220 : END SUBROUTINE perdew_wang_lda_calc
248 :
249 : ! **************************************************************************************************
250 : !> \brief Calculate the correlation energy and its derivatives
251 : !> wrt to rho (the electron density) up to 3rd order. This
252 : !> is the LSD version of the Perdew-Wang correlation energy
253 : !> If no order argument is given, then the routine calculates
254 : !> just the energy.
255 : !> \param method ...
256 : !> \param rho_set ...
257 : !> \param deriv_set ...
258 : !> \param order maximum order of derivatives to calculate
259 : !> order must lie between 0 and 3.
260 : !> \param scale ...
261 : ! **************************************************************************************************
262 40 : SUBROUTINE perdew_wang_lsd_eval(method, rho_set, deriv_set, order, scale)
263 : INTEGER, INTENT(in) :: method
264 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
265 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
266 : INTEGER, INTENT(IN), OPTIONAL :: order
267 : REAL(kind=dp), INTENT(in) :: scale
268 :
269 : CHARACTER(len=*), PARAMETER :: routineN = 'perdew_wang_lsd_eval'
270 :
271 : INTEGER :: npoints, timer_handle
272 : INTEGER, DIMENSION(2, 3) :: bo
273 : REAL(KIND=dp) :: rho_cutoff
274 20 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: a, b, dummy, e_0, ea, eaa, eaaa, eaab, &
275 20 : eab, eabb, eb, ebb, ebbb
276 : TYPE(xc_derivative_type), POINTER :: deriv
277 :
278 20 : CALL timeset(routineN, timer_handle)
279 20 : NULLIFY (a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, ebbb)
280 : CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, &
281 20 : local_bounds=bo, rho_cutoff=rho_cutoff)
282 20 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
283 :
284 20 : CALL perdew_wang_init(method, rho_cutoff)
285 :
286 : ! meaningful default for the arrays we don't need: let us make compiler
287 : ! and debugger happy...
288 20 : dummy => a
289 :
290 20 : e_0 => dummy
291 20 : ea => dummy; eb => dummy
292 20 : eaa => dummy; eab => dummy; ebb => dummy
293 20 : eaaa => dummy; eaab => dummy; eabb => dummy; ebbb => dummy
294 :
295 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
296 20 : allocate_deriv=.TRUE.)
297 20 : CALL xc_derivative_get(deriv, deriv_data=e_0)
298 20 : IF (order >= 1) THEN
299 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
300 20 : allocate_deriv=.TRUE.)
301 20 : CALL xc_derivative_get(deriv, deriv_data=ea)
302 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
303 20 : allocate_deriv=.TRUE.)
304 20 : CALL xc_derivative_get(deriv, deriv_data=eb)
305 : END IF
306 20 : IF (order >= 2) THEN
307 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
308 0 : allocate_deriv=.TRUE.)
309 0 : CALL xc_derivative_get(deriv, deriv_data=eaa)
310 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
311 0 : allocate_deriv=.TRUE.)
312 0 : CALL xc_derivative_get(deriv, deriv_data=eab)
313 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
314 0 : allocate_deriv=.TRUE.)
315 0 : CALL xc_derivative_get(deriv, deriv_data=ebb)
316 : END IF
317 20 : IF (order >= 3) THEN
318 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
319 0 : allocate_deriv=.TRUE.)
320 0 : CALL xc_derivative_get(deriv, deriv_data=eaaa)
321 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
322 0 : allocate_deriv=.TRUE.)
323 0 : CALL xc_derivative_get(deriv, deriv_data=eaab)
324 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
325 0 : allocate_deriv=.TRUE.)
326 0 : CALL xc_derivative_get(deriv, deriv_data=eabb)
327 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
328 0 : allocate_deriv=.TRUE.)
329 0 : CALL xc_derivative_get(deriv, deriv_data=ebbb)
330 : END IF
331 20 : IF (order > 3) THEN
332 0 : CPABORT("derivatives bigger than 3 not implemented")
333 : END IF
334 :
335 : CALL perdew_wang_lsd_calc(a, b, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
336 20 : ebbb, npoints, order, scale)
337 :
338 20 : CALL timestop(timer_handle)
339 :
340 20 : END SUBROUTINE perdew_wang_lsd_eval
341 :
342 : ! **************************************************************************************************
343 : !> \brief ...
344 : !> \param rhoa ...
345 : !> \param rhob ...
346 : !> \param e_0 ...
347 : !> \param ea ...
348 : !> \param eb ...
349 : !> \param eaa ...
350 : !> \param eab ...
351 : !> \param ebb ...
352 : !> \param eaaa ...
353 : !> \param eaab ...
354 : !> \param eabb ...
355 : !> \param ebbb ...
356 : !> \param npoints ...
357 : !> \param order ...
358 : !> \param scale ...
359 : ! **************************************************************************************************
360 20 : SUBROUTINE perdew_wang_lsd_calc(rhoa, rhob, e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, eabb, &
361 : ebbb, npoints, order, scale)
362 : !FM low-level computation routine
363 : REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob
364 : REAL(KIND=dp), DIMENSION(*), INTENT(inout) :: e_0, ea, eb, eaa, eab, ebb, eaaa, eaab, &
365 : eabb, ebbb
366 : INTEGER, INTENT(in) :: npoints, order
367 : REAL(kind=dp), INTENT(in) :: scale
368 :
369 : INTEGER :: k
370 : REAL(KIND=dp) :: rho
371 : REAL(KIND=dp), DIMENSION(0:9) :: ed
372 :
373 : !$OMP PARALLEL DO PRIVATE (k, rho, ed) DEFAULT(NONE)&
374 20 : !$OMP SHARED(npoints,rhoa,rhob,eps_rho,order,e_0,ea,eb,eaa,eab,ebb,eaaa,eaab,eabb,ebbb,scale)
375 : DO k = 1, npoints
376 :
377 : rho = rhoa(k) + rhob(k)
378 : IF (rho > eps_rho) THEN
379 :
380 : ed = 0.0_dp
381 : CALL pw_lsd_ed_loc(rhoa(k), rhob(k), ed, order)
382 : ed = ed*scale
383 :
384 : e_0(k) = e_0(k) + rho*ed(0)
385 : IF (order >= 1) THEN
386 : ea(k) = ea(k) + ed(0) + rho*ed(1)
387 : eb(k) = eb(k) + ed(0) + rho*ed(2)
388 : END IF
389 : IF (order >= 2) THEN
390 : eaa(k) = eaa(k) + 2.0_dp*ed(1) + rho*ed(3)
391 : eab(k) = eab(k) + ed(1) + ed(2) + rho*ed(4)
392 : ebb(k) = ebb(k) + 2.0_dp*ed(2) + rho*ed(5)
393 : END IF
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 :
405 20 : END SUBROUTINE perdew_wang_lsd_calc
406 :
407 : ! **************************************************************************************************
408 : !> \brief ...
409 : !> \param rho_a ...
410 : !> \param rho_b ...
411 : !> \param fxc_aa ...
412 : !> \param fxc_ab ...
413 : !> \param fxc_bb ...
414 : !> \param scalec ...
415 : !> \param eps_rho ...
416 : ! **************************************************************************************************
417 10 : SUBROUTINE perdew_wang_fxc_calc(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb, scalec, eps_rho)
418 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
419 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_ab, fxc_bb
420 : REAL(kind=dp), INTENT(in) :: scalec, eps_rho
421 :
422 : INTEGER :: i, j, k
423 : INTEGER, DIMENSION(2, 3) :: bo
424 : REAL(KIND=dp) :: rho, rhoa, rhob, eaa, eab, ebb
425 : REAL(KIND=dp), DIMENSION(0:9) :: ed
426 :
427 10 : CALL perdew_wang_init(pw_orig, eps_rho)
428 100 : bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
429 : !$OMP PARALLEL DO PRIVATE(i,j,k,rho,rhoa,rhob,ed,eaa,eab,ebb) DEFAULT(NONE)&
430 10 : !$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_ab,fxc_bb,scalec,eps_rho)
431 : DO k = bo(1, 3), bo(2, 3)
432 : DO j = bo(1, 2), bo(2, 2)
433 : DO i = bo(1, 1), bo(2, 1)
434 : rhoa = rho_a%array(i, j, k)
435 : rhob = rho_b%array(i, j, k)
436 : rho = rhoa + rhob
437 : IF (rho > eps_rho) THEN
438 : ed = 0.0_dp
439 : CALL pw_lsd_ed_loc(rhoa, rhob, ed, 2)
440 : ed = ed*scalec
441 : eaa = 2.0_dp*ed(1) + rho*ed(3)
442 : eab = ed(1) + ed(2) + rho*ed(4)
443 : ebb = 2.0_dp*ed(2) + rho*ed(5)
444 : fxc_aa%array(i, j, k) = fxc_aa%array(i, j, k) + eaa
445 : fxc_ab%array(i, j, k) = fxc_ab%array(i, j, k) + eab
446 : fxc_bb%array(i, j, k) = fxc_bb%array(i, j, k) + ebb
447 : END IF
448 : END DO
449 : END DO
450 : END DO
451 :
452 10 : END SUBROUTINE perdew_wang_fxc_calc
453 :
454 15954082 : @:calc_g()
455 :
456 : ! **************************************************************************************************
457 : !> \brief ...
458 : !> \param rho ...
459 : !> \param ed ...
460 : !> \param order ...
461 : ! **************************************************************************************************
462 11340400 : SUBROUTINE pw_lda_ed_loc(rho, ed, order)
463 :
464 : REAL(KIND=dp), INTENT(IN) :: rho
465 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: ed
466 : INTEGER, INTENT(IN) :: order
467 :
468 : INTEGER :: m, order_
469 : LOGICAL, DIMENSION(0:3) :: calc
470 : REAL(KIND=dp), DIMENSION(0:3) :: e0, r
471 :
472 11340400 : order_ = order
473 56702000 : ed = 0
474 11340400 : calc = .FALSE.
475 :
476 11340400 : IF (order_ >= 0) THEN
477 33930075 : calc(0:order_) = .TRUE.
478 : ELSE
479 0 : order_ = -1*order_
480 0 : calc(order_) = .TRUE.
481 : END IF
482 :
483 11340400 : CALL calc_rs(rho, r(0))
484 11340400 : CALL calc_g(r(0), 0, e0, order_)
485 :
486 11340400 : IF (order_ >= 1) r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
487 11340400 : IF (order_ >= 2) r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
488 11340400 : IF (order_ >= 3) r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
489 :
490 11340400 : m = 0
491 11340400 : IF (calc(0)) THEN
492 11340400 : ed(m) = e0(0)
493 11340400 : m = m + 1
494 : END IF
495 11340400 : IF (calc(1)) THEN
496 11249275 : ed(m) = e0(1)*r(1)
497 11249275 : m = m + 1
498 : END IF
499 11340400 : IF (calc(2)) THEN
500 0 : ed(m) = e0(2)*r(1)**2 + e0(1)*r(2)
501 0 : m = m + 1
502 : END IF
503 11340400 : IF (calc(3)) THEN
504 0 : ed(m) = e0(3)*r(1)**3 + e0(2)*3.0_dp*r(1)*r(2) + e0(1)*r(3)
505 : END IF
506 :
507 11340400 : END SUBROUTINE pw_lda_ed_loc
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param a ...
512 : !> \param b ...
513 : !> \param ed ...
514 : !> \param order ...
515 : ! **************************************************************************************************
516 1537894 : SUBROUTINE pw_lsd_ed_loc(a, b, ed, order)
517 :
518 : REAL(KIND=dp), INTENT(IN) :: a, b
519 : REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: ed
520 : INTEGER, INTENT(IN) :: order
521 :
522 : INTEGER :: m, order_
523 : LOGICAL, DIMENSION(0:3) :: calc
524 : REAL(KIND=dp) :: rho, tr, trr, trrr, trrz, trz, trzz, tz, &
525 : tzz, tzzz
526 : REAL(KIND=dp), DIMENSION(0:3) :: ac, e0, e1, f, r
527 : REAL(KIND=dp), DIMENSION(0:3, 0:3) :: z
528 :
529 1537894 : order_ = order
530 1537894 : calc = .FALSE.
531 :
532 1537894 : IF (order_ > 0) THEN
533 5240326 : calc(0:order_) = .TRUE.
534 : ELSE
535 0 : order_ = -1*order_
536 0 : calc(order_) = .TRUE.
537 : END IF
538 :
539 1537894 : rho = a + b
540 :
541 1537894 : CALL calc_fx(a, b, f(0:order_), order_)
542 1537894 : CALL calc_rs(rho, r(0))
543 1537894 : CALL calc_g(r(0), -1, ac(0:order_), order_)
544 1537894 : CALL calc_g(r(0), 0, e0(0:order_), order_)
545 1537894 : CALL calc_g(r(0), 1, e1(0:order_), order_)
546 1537894 : CALL calc_z(a, b, z(0:order_, 0:order_), order_)
547 :
548 : !! calculate first partial derivatives
549 1537894 : IF (order_ >= 1) THEN
550 1537894 : r(1) = (-1.0_dp/3.0_dp)*r(0)/rho
551 : tr = e0(1) &
552 : + fpp*ac(1)*f(0) &
553 : - fpp*ac(1)*f(0)*z(0, 0)**4 &
554 1537894 : + (e1(1) - e0(1))*f(0)*z(0, 0)**4
555 : tz = fpp*ac(0)*f(1) &
556 : - fpp*ac(0)*f(1)*z(0, 0)**4 &
557 : - fpp*ac(0)*f(0)*4.0_dp*z(0, 0)**3 &
558 : + (e1(0) - e0(0))*f(1)*z(0, 0)**4 &
559 1537894 : + (e1(0) - e0(0))*f(0)*4.0_dp*z(0, 0)**3
560 : END IF
561 :
562 : !! calculate second partial derivatives
563 1537894 : IF (order_ >= 2) THEN
564 626644 : r(2) = (-4.0_dp/3.0_dp)*r(1)/rho
565 : trr = e0(2) &
566 : + fpp*ac(2)*f(0) &
567 : - fpp*ac(2)*f(0)*z(0, 0)**4 &
568 626644 : + (e1(2) - e0(2))*f(0)*z(0, 0)**4
569 : trz = fpp*ac(1)*f(1) &
570 : - fpp*ac(1)*f(1)*z(0, 0)**4 &
571 : - fpp*ac(1)*f(0)*4.0_dp*z(0, 0)**3 &
572 : + (e1(1) - e0(1))*f(1)*z(0, 0)**4 &
573 626644 : + (e1(1) - e0(1))*f(0)*4.0_dp*z(0, 0)**3
574 : tzz = fpp*ac(0)*f(2) &
575 : - fpp*ac(0)*f(2)*z(0, 0)**4 &
576 : - fpp*ac(0)*f(1)*8.0_dp*z(0, 0)**3 &
577 : - fpp*ac(0)*f(0)*12.0_dp*z(0, 0)**2 &
578 : + (e1(0) - e0(0))*f(2)*z(0, 0)**4 &
579 : + (e1(0) - e0(0))*f(1)*8.0_dp*z(0, 0)**3 &
580 626644 : + (e1(0) - e0(0))*f(0)*12.0_dp*z(0, 0)**2
581 : END IF
582 :
583 : !! calculate third derivatives
584 1537894 : IF (order_ >= 3) THEN
585 :
586 0 : r(3) = (-7.0_dp/3.0_dp)*r(2)/rho
587 :
588 : trrr = e0(3) &
589 : + fpp*ac(3)*f(0) &
590 : - fpp*ac(3)*f(0)*z(0, 0)**4 &
591 0 : + (e1(3) - e0(3))*f(0)*z(0, 0)**4
592 :
593 : trrz = fpp*ac(2)*f(1) &
594 : - fpp*ac(2)*f(1)*z(0, 0)**4 &
595 : - fpp*ac(2)*f(0)*4.0_dp*z(0, 0)**3 &
596 : + (e1(2) - e0(2))*f(1)*z(0, 0)**4 &
597 0 : + (e1(2) - e0(2))*f(0)*4.0_dp*z(0, 0)**3
598 :
599 : trzz = fpp*ac(1)*f(2) &
600 : - fpp*ac(1)*f(2)*z(0, 0)**4 &
601 : - fpp*ac(1)*f(1)*8.0_dp*z(0, 0)**3 &
602 : - fpp*ac(1)*f(0)*12.0_dp*z(0, 0)**2 &
603 : + (e1(1) - e0(1))*f(2)*z(0, 0)**4 &
604 : + (e1(1) - e0(1))*f(1)*8.0_dp*z(0, 0)**3 &
605 0 : + (e1(1) - e0(1))*f(0)*12.0_dp*z(0, 0)**2
606 :
607 : tzzz = fpp*ac(0)*f(3) &
608 : - fpp*ac(0)*f(3)*z(0, 0)**4 &
609 : - fpp*ac(0)*f(2)*12.0_dp*z(0, 0)**3 &
610 : - fpp*ac(0)*f(1)*36.0_dp*z(0, 0)**2 &
611 : - fpp*ac(0)*f(0)*24.0_dp*z(0, 0) &
612 : + (e1(0) - e0(0))*f(3)*z(0, 0)**4 &
613 : + (e1(0) - e0(0))*f(2)*12.0_dp*z(0, 0)**3 &
614 : + (e1(0) - e0(0))*f(1)*36.0_dp*z(0, 0)**2 &
615 0 : + (e1(0) - e0(0))*f(0)*24.0_dp*z(0, 0)
616 : END IF
617 :
618 1537894 : m = 0
619 1537894 : IF (calc(0)) THEN
620 : ed(m) = e0(0) &
621 : + fpp*ac(0)*f(0)*(1.0_dp - z(0, 0)**4) &
622 1537894 : + (e1(0) - e0(0))*f(0)*z(0, 0)**4
623 1537894 : m = m + 1
624 : END IF
625 1537894 : IF (calc(1)) THEN
626 1537894 : ed(m) = tr*r(1) + tz*z(1, 0)
627 1537894 : ed(m + 1) = tr*r(1) + tz*z(0, 1)
628 1537894 : m = m + 2
629 : END IF
630 1537894 : IF (calc(2)) THEN
631 : ed(m) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(1, 0) &
632 626644 : + tr*r(2) + tzz*z(1, 0)**2 + tz*z(2, 0)
633 : ed(m + 1) = trr*r(1)**2 + trz*r(1)*(z(0, 1) + z(1, 0)) &
634 626644 : + tr*r(2) + tzz*z(1, 0)*z(0, 1) + tz*z(1, 1)
635 : ed(m + 2) = trr*r(1)**2 + 2.0_dp*trz*r(1)*z(0, 1) &
636 626644 : + tr*r(2) + tzz*z(0, 1)**2 + tz*z(0, 2)
637 626644 : m = m + 3
638 : END IF
639 1537894 : IF (calc(3)) THEN
640 : ed(m) = &
641 : trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(1, 0) &
642 : + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(1, 0) + tr*r(3) &
643 : + 3.0_dp*trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**3 &
644 : + 3.0_dp*trz*r(1)*z(2, 0) &
645 0 : + 3.0_dp*tzz*z(1, 0)*z(2, 0) + tz*z(3, 0)
646 : ed(m + 1) = &
647 : trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(1, 0) + z(0, 1)) &
648 : + 2.0_dp*trzz*r(1)*z(1, 0)*z(0, 1) &
649 : + 2.0_dp*trz*(r(2)*z(1, 0) + r(1)*z(1, 1)) &
650 : + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(0, 1) + tr*r(3) &
651 : + trzz*r(1)*z(1, 0)**2 + tzzz*z(1, 0)**2*z(0, 1) &
652 : + 2.0_dp*tzz*z(1, 0)*z(1, 1) &
653 0 : + trz*r(1)*z(2, 0) + tzz*z(2, 0)*z(0, 1) + tz*z(2, 1)
654 : ed(m + 2) = &
655 : trrr*r(1)**3 + trrz*r(1)**2*(2.0_dp*z(0, 1) + z(1, 0)) &
656 : + 2.0_dp*trzz*r(1)*z(0, 1)*z(1, 0) &
657 : + 2.0_dp*trz*(r(2)*z(0, 1) + r(1)*z(1, 1)) &
658 : + 3.0_dp*trr*r(2)*r(1) + trz*r(2)*z(1, 0) + tr*r(3) &
659 : + trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**2*z(1, 0) &
660 : + 2.0_dp*tzz*z(0, 1)*z(1, 1) &
661 0 : + trz*r(1)*z(0, 2) + tzz*z(0, 2)*z(1, 0) + tz*z(1, 2)
662 : ed(m + 3) = &
663 : trrr*r(1)**3 + 3.0_dp*trrz*r(1)**2*z(0, 1) &
664 : + 3.0_dp*trr*r(1)*r(2) + 3.0_dp*trz*r(2)*z(0, 1) + tr*r(3) &
665 : + 3.0_dp*trzz*r(1)*z(0, 1)**2 + tzzz*z(0, 1)**3 &
666 : + 3.0_dp*trz*r(1)*z(0, 2) &
667 0 : + 3.0_dp*tzz*z(0, 1)*z(0, 2) + tz*z(0, 3)
668 : END IF
669 :
670 1537894 : END SUBROUTINE pw_lsd_ed_loc
671 :
672 : END MODULE xc_perdew_wang
|
