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