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 LDA functional in the Pade approximation
10 : !> Literature: S. Goedecker, M. Teter and J. Hutter,
11 : !> Phys. Rev. B 54, 1703 (1996)
12 : !> \note
13 : !> Order of derivatives is: LDA 0; 1; 2; 3;
14 : !> LSD 0; a b; aa ab bb; aaa aab abb bbb;
15 : !> \par History
16 : !> JGH (26.02.2003) : OpenMP enabled
17 : !> \author JGH (15.02.2002)
18 : ! **************************************************************************************************
19 : MODULE xc_pade
20 : USE bibliography, ONLY: Goedecker1996,&
21 : cite_reference
22 : USE kinds, ONLY: dp
23 : USE pw_types, ONLY: pw_r3d_rs_type
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_rs_pw,&
34 : set_util
35 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
36 : USE xc_rho_set_types, ONLY: xc_rho_set_type
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
44 : f23 = 2.0_dp*f13, &
45 : f43 = 4.0_dp*f13
46 :
47 : REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, &
48 : a1 = 0.2217058676663745E+1_dp, &
49 : a2 = 0.7405551735357053E+0_dp, &
50 : a3 = 0.1968227878617998E-1_dp, &
51 : b1 = 1.0000000000000000E+0_dp, &
52 : b2 = 0.4504130959426697E+1_dp, &
53 : b3 = 0.1110667363742916E+1_dp, &
54 : b4 = 0.2359291751427506E-1_dp
55 :
56 : REAL(KIND=dp), PARAMETER :: da0 = 0.119086804055547E+0_dp, &
57 : da1 = 0.6157402568883345E+0_dp, &
58 : da2 = 0.1574201515892867E+0_dp, &
59 : da3 = 0.3532336663397157E-2_dp, &
60 : db1 = 0.0000000000000000E+0_dp, &
61 : db2 = 0.2673612973836267E+0_dp, &
62 : db3 = 0.2052004607777787E+0_dp, &
63 : db4 = 0.4200005045691381E-2_dp
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pade'
66 :
67 : PUBLIC :: pade_lda_pw_eval, pade_lsd_pw_eval, pade_info, pade_init, pade_fxc_eval
68 :
69 : REAL(KIND=dp) :: eps_rho
70 : LOGICAL :: debug_flag
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param cutoff ...
77 : !> \param debug ...
78 : ! **************************************************************************************************
79 93488 : SUBROUTINE pade_init(cutoff, debug)
80 :
81 : REAL(KIND=dp), INTENT(IN) :: cutoff
82 : LOGICAL, INTENT(IN), OPTIONAL :: debug
83 :
84 93488 : eps_rho = cutoff
85 93488 : CALL set_util(cutoff)
86 :
87 93488 : CALL cite_reference(Goedecker1996)
88 :
89 93488 : IF (PRESENT(debug)) THEN
90 0 : debug_flag = debug
91 : ELSE
92 93488 : debug_flag = .FALSE.
93 : END IF
94 :
95 93488 : END SUBROUTINE pade_init
96 :
97 : ! **************************************************************************************************
98 : !> \brief ...
99 : !> \param reference ...
100 : !> \param shortform ...
101 : !> \param lsd ...
102 : !> \param needs ...
103 : !> \param max_deriv ...
104 : ! **************************************************************************************************
105 88283 : SUBROUTINE pade_info(reference, shortform, lsd, needs, max_deriv)
106 :
107 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
108 : LOGICAL, INTENT(IN), OPTIONAL :: lsd
109 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
110 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
111 :
112 88283 : IF (PRESENT(reference)) THEN
113 : reference = "S. Goedecker, M. Teter and J. Hutter," &
114 573 : //" Phys. Rev. B 54, 1703 (1996)"
115 : END IF
116 88283 : IF (PRESENT(shortform)) THEN
117 573 : shortform = "S. Goedecker et al., PRB 54, 1703 (1996)"
118 : END IF
119 :
120 88283 : IF (PRESENT(needs)) THEN
121 87710 : IF (.NOT. PRESENT(lsd)) THEN
122 0 : CPABORT("Arguments mismatch.")
123 : END IF
124 87710 : IF (lsd) THEN
125 14129 : needs%rho_spin = .TRUE.
126 : ELSE
127 73581 : needs%rho = .TRUE.
128 : END IF
129 : END IF
130 :
131 88283 : IF (PRESENT(max_deriv)) max_deriv = 3
132 :
133 88283 : END SUBROUTINE pade_info
134 :
135 : ! **************************************************************************************************
136 : !> \brief ...
137 : !> \param deriv_set ...
138 : !> \param rho_set ...
139 : !> \param order ...
140 : ! **************************************************************************************************
141 79337 : SUBROUTINE pade_lda_pw_eval(deriv_set, rho_set, order)
142 :
143 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
144 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
145 : INTEGER, INTENT(IN), OPTIONAL :: order
146 :
147 : INTEGER :: n
148 : LOGICAL :: calc(0:4)
149 79337 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rs
150 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
151 79337 : POINTER :: e_0, e_r, e_rr, e_rrr
152 : TYPE(xc_derivative_type), POINTER :: deriv
153 :
154 79337 : calc = .FALSE.
155 243195 : IF (order >= 0) calc(0:order) = .TRUE.
156 79337 : IF (order < 0) calc(-order) = .TRUE.
157 :
158 317348 : n = PRODUCT(rho_set%local_bounds(2, :) - rho_set%local_bounds(1, :) + [1, 1, 1])
159 238011 : ALLOCATE (rs(n))
160 :
161 79337 : CALL calc_rs_pw(rho_set%rho, rs, n)
162 79337 : IF (calc(0) .AND. calc(1)) THEN
163 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
164 76311 : allocate_deriv=.TRUE.)
165 76311 : CALL xc_derivative_get(deriv, deriv_data=e_0)
166 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
167 76311 : allocate_deriv=.TRUE.)
168 76311 : CALL xc_derivative_get(deriv, deriv_data=e_r)
169 76311 : CALL pade_lda_01(n, rho_set%rho, rs, e_0, e_r)
170 3026 : ELSE IF (calc(0)) THEN
171 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
172 3026 : allocate_deriv=.TRUE.)
173 3026 : CALL xc_derivative_get(deriv, deriv_data=e_0)
174 3026 : CALL pade_lda_0(n, rho_set%rho, rs, e_0)
175 0 : ELSE IF (calc(1)) THEN
176 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
177 0 : allocate_deriv=.TRUE.)
178 0 : CALL xc_derivative_get(deriv, deriv_data=e_r)
179 0 : CALL pade_lda_1(n, rho_set%rho, rs, e_r)
180 : END IF
181 79337 : IF (calc(2)) THEN
182 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
183 8210 : allocate_deriv=.TRUE.)
184 8210 : CALL xc_derivative_get(deriv, deriv_data=e_rr)
185 8210 : CALL pade_lda_2(n, rho_set%rho, rs, e_rr)
186 : END IF
187 79337 : IF (calc(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_rrr)
191 0 : CALL pade_lda_3(n, rho_set%rho, rs, e_rrr)
192 : END IF
193 :
194 79337 : DEALLOCATE (rs)
195 :
196 79337 : END SUBROUTINE pade_lda_pw_eval
197 :
198 : ! **************************************************************************************************
199 : !> \brief ...
200 : !> \param deriv_set ...
201 : !> \param rho_set ...
202 : !> \param order ...
203 : ! **************************************************************************************************
204 14149 : SUBROUTINE pade_lsd_pw_eval(deriv_set, rho_set, order)
205 :
206 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
207 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
208 : INTEGER, INTENT(IN), OPTIONAL :: order
209 :
210 : INTEGER :: i, j, k
211 : LOGICAL :: calc(0:4)
212 : REAL(KIND=dp) :: rhoa, rhob, rs
213 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
214 14149 : POINTER :: e_0, e_ra, e_rara, e_rarara, e_rararb, &
215 14149 : e_rarb, e_rarbrb, e_rb, e_rbrb, &
216 14149 : e_rbrbrb
217 : REAL(KIND=dp), DIMENSION(4) :: fx
218 : TYPE(xc_derivative_type), POINTER :: deriv
219 :
220 14149 : calc = .FALSE.
221 42021 : IF (order >= 0) calc(0:order) = .TRUE.
222 14149 : IF (order < 0) calc(-order) = .TRUE.
223 :
224 14149 : IF (calc(0)) THEN
225 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
226 14149 : allocate_deriv=.TRUE.)
227 14149 : CALL xc_derivative_get(deriv, deriv_data=e_0)
228 : END IF
229 14149 : IF (calc(1)) THEN
230 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
231 12729 : allocate_deriv=.TRUE.)
232 12729 : CALL xc_derivative_get(deriv, deriv_data=e_ra)
233 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
234 12729 : allocate_deriv=.TRUE.)
235 12729 : CALL xc_derivative_get(deriv, deriv_data=e_rb)
236 : END IF
237 14149 : IF (calc(2)) THEN
238 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
239 994 : allocate_deriv=.TRUE.)
240 994 : CALL xc_derivative_get(deriv, deriv_data=e_rara)
241 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
242 994 : allocate_deriv=.TRUE.)
243 994 : CALL xc_derivative_get(deriv, deriv_data=e_rarb)
244 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
245 994 : allocate_deriv=.TRUE.)
246 994 : CALL xc_derivative_get(deriv, deriv_data=e_rbrb)
247 : END IF
248 14149 : IF (calc(3)) THEN
249 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
250 0 : allocate_deriv=.TRUE.)
251 0 : CALL xc_derivative_get(deriv, deriv_data=e_rarara)
252 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
253 0 : allocate_deriv=.TRUE.)
254 0 : CALL xc_derivative_get(deriv, deriv_data=e_rararb)
255 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
256 0 : allocate_deriv=.TRUE.)
257 0 : CALL xc_derivative_get(deriv, deriv_data=e_rarbrb)
258 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
259 0 : allocate_deriv=.TRUE.)
260 0 : CALL xc_derivative_get(deriv, deriv_data=e_rbrbrb)
261 : END IF
262 :
263 : !$OMP PARALLEL DO PRIVATE(i,j,k,fx,rhoa,rhob,rs) DEFAULT(NONE)&
264 14149 : !$OMP SHARED(rho_set,order,e_0,e_ra,e_rb,calc,e_rara,e_rarb,e_rbrb,e_rarara,e_rararb,e_rarbrb,e_rbrbrb)
265 : DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
266 : DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
267 : DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
268 :
269 : rhoa = rho_set%rhoa(i, j, k)
270 : rhob = rho_set%rhob(i, j, k)
271 : fx(1) = rhoa + rhob
272 :
273 : CALL calc_rs(fx(1), rs)
274 : CALL calc_fx(rhoa, rhob, fx, ABS(order))
275 :
276 : IF (calc(0) .AND. calc(1)) THEN
277 : CALL pade_lsd_01(rhoa, rhob, rs, fx, &
278 : e_0(i, j, k), e_ra(i, j, k), e_rb(i, j, k))
279 : ELSE IF (calc(0)) THEN
280 : CALL pade_lsd_0(rhoa, rhob, rs, fx, e_0(i, j, k))
281 : ELSE IF (calc(1)) THEN
282 : CALL pade_lsd_1(rhoa, rhob, rs, fx, &
283 : e_ra(i, j, k), e_rb(i, j, k))
284 : END IF
285 : IF (calc(2)) THEN
286 : CALL pade_lsd_2(rhoa, rhob, rs, fx, &
287 : e_rara(i, j, k), e_rarb(i, j, k), e_rbrb(i, j, k))
288 : END IF
289 : IF (calc(3)) THEN
290 : CALL pade_lsd_3(rhoa, rhob, rs, fx, &
291 : e_rarara(i, j, k), e_rararb(i, j, k), e_rarbrb(i, j, k), e_rbrbrb(i, j, k))
292 : END IF
293 : END DO
294 : END DO
295 : END DO
296 :
297 14149 : END SUBROUTINE pade_lsd_pw_eval
298 :
299 : ! **************************************************************************************************
300 : !> \brief ...
301 : !> \param n ...
302 : !> \param rho ...
303 : !> \param rs ...
304 : !> \param pot ...
305 : ! **************************************************************************************************
306 3026 : SUBROUTINE pade_lda_0(n, rho, rs, pot)
307 :
308 : INTEGER, INTENT(IN) :: n
309 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs
310 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
311 :
312 : INTEGER :: ip
313 : REAL(KIND=dp) :: epade, p, q
314 :
315 : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade) DEFAULT(NONE)&
316 3026 : !$OMP SHARED(n,rho,eps_rho,pot,rs)
317 : DO ip = 1, n
318 : IF (rho(ip) > eps_rho) THEN
319 : p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
320 : q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
321 : epade = -p/q
322 : pot(ip) = pot(ip) + epade*rho(ip)
323 : END IF
324 : END DO
325 :
326 3026 : END SUBROUTINE pade_lda_0
327 :
328 : ! **************************************************************************************************
329 : !> \brief ...
330 : !> \param n ...
331 : !> \param rho ...
332 : !> \param rs ...
333 : !> \param pot ...
334 : ! **************************************************************************************************
335 0 : SUBROUTINE pade_lda_1(n, rho, rs, pot)
336 :
337 : INTEGER, INTENT(IN) :: n
338 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs
339 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
340 :
341 : INTEGER :: ip
342 : REAL(KIND=dp) :: depade, dpv, dq, epade, p, q
343 :
344 : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade,dpv,dq,depade) DEFAULT(NONE)&
345 0 : !$OMP SHARED(n,rho,eps_rho,rs,pot)
346 :
347 : DO ip = 1, n
348 : IF (rho(ip) > eps_rho) THEN
349 :
350 : p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
351 : q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
352 : epade = -p/q
353 :
354 : dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
355 : dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
356 : depade = f13*rs(ip)*(dpv*q - p*dq)/(q*q)
357 :
358 : pot(ip) = pot(ip) + epade + depade
359 :
360 : END IF
361 : END DO
362 :
363 0 : END SUBROUTINE pade_lda_1
364 :
365 : ! **************************************************************************************************
366 : !> \brief ...
367 : !> \param n ...
368 : !> \param rho ...
369 : !> \param rs ...
370 : !> \param pot0 ...
371 : !> \param pot1 ...
372 : ! **************************************************************************************************
373 76311 : SUBROUTINE pade_lda_01(n, rho, rs, pot0, pot1)
374 :
375 : INTEGER, INTENT(IN) :: n
376 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs
377 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot0, pot1
378 :
379 : INTEGER :: ip
380 : REAL(KIND=dp) :: depade, dpv, dq, epade, p, q
381 :
382 : !$OMP PARALLEL DO PRIVATE(ip,p,q,epade,dpv,dq,depade) DEFAULT(NONE)&
383 76311 : !$OMP SHARED(n,rho,eps_rho,pot0,pot1)
384 :
385 : DO ip = 1, n
386 : IF (rho(ip) > eps_rho) THEN
387 :
388 : p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
389 : q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
390 : epade = -p/q
391 :
392 : dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
393 : dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
394 : depade = f13*rs(ip)*(dpv*q - p*dq)/(q*q)
395 :
396 : pot0(ip) = pot0(ip) + epade*rho(ip)
397 : pot1(ip) = pot1(ip) + epade + depade
398 :
399 : END IF
400 : END DO
401 :
402 76311 : END SUBROUTINE pade_lda_01
403 :
404 : ! **************************************************************************************************
405 : !> \brief ...
406 : !> \param n ...
407 : !> \param rho ...
408 : !> \param rs ...
409 : !> \param pot ...
410 : ! **************************************************************************************************
411 8210 : SUBROUTINE pade_lda_2(n, rho, rs, pot)
412 :
413 : INTEGER, INTENT(IN) :: n
414 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs
415 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
416 :
417 : INTEGER :: ip
418 : REAL(KIND=dp) :: d2p, d2q, dpv, dq, p, q, rsr, t1, t2, t3
419 :
420 : !$OMP PARALLEL DO PRIVATE(ip,p,q,dpv,dq,d2p,d2q,rsr,t1,t2,t3) DEFAULT(NONE)&
421 8210 : !$OMP SHARED(n,rho,eps_rho,rs)
422 :
423 : DO ip = 1, n
424 : IF (rho(ip) > eps_rho) THEN
425 :
426 : p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
427 : q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
428 :
429 : dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
430 : dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
431 :
432 : d2p = 2.0_dp*a2 + 6.0_dp*a3*rs(ip)
433 : d2q = 2.0_dp*b2 + (6.0_dp*b3 + 12.0_dp*b4*rs(ip))*rs(ip)
434 :
435 : rsr = rs(ip)/rho(ip)
436 : t1 = (p*dq - dpv*q)/(q*q)
437 : t2 = (d2p*q - p*d2q)/(q*q)
438 : t3 = (p*dq*dq - dpv*q*dq)/(q*q*q)
439 :
440 : pot(ip) = pot(ip) - f13*(f23*t1 + f13*t2*rs(ip) + f23*t3*rs(ip))*rsr
441 :
442 : END IF
443 : END DO
444 :
445 8210 : END SUBROUTINE pade_lda_2
446 :
447 : ! **************************************************************************************************
448 : !> \brief ...
449 : !> \param n ...
450 : !> \param rho ...
451 : !> \param rs ...
452 : !> \param pot ...
453 : ! **************************************************************************************************
454 0 : SUBROUTINE pade_lda_3(n, rho, rs, pot)
455 :
456 : INTEGER, INTENT(IN) :: n
457 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs
458 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: pot
459 :
460 : INTEGER :: ip
461 : REAL(KIND=dp) :: ab1, ab2, ab3, d2p, d2q, d3p, d3q, dpv, &
462 : dq, p, q, rsr1, rsr2, rsr3
463 :
464 : !$OMP PARALLEL DO PRIVATE(ip,p,q,dpv,dq,d2p,d2q,d3p,d3q,ab1,ab2,ab3,rsr1,rsr2,rsr3) DEFAULT(NONE)&
465 0 : !$OMP SHARED(n,rho,eps_rho,rs,pot)
466 :
467 : DO ip = 1, n
468 : IF (rho(ip) > eps_rho) THEN
469 :
470 : p = a0 + (a1 + (a2 + a3*rs(ip))*rs(ip))*rs(ip)
471 : q = (b1 + (b2 + (b3 + b4*rs(ip))*rs(ip))*rs(ip))*rs(ip)
472 :
473 : dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs(ip))*rs(ip)
474 : dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs(ip))*rs(ip))*rs(ip)
475 :
476 : d2p = 2.0_dp*a2 + 6.0_dp*a3*rs(ip)
477 : d2q = 2.0_dp*b2 + (6.0_dp*b3 + 12.0_dp*b4*rs(ip))*rs(ip)
478 :
479 : d3p = 6.0_dp*a3
480 : d3q = 6.0_dp*b3 + 24.0_dp*b4*rs(ip)
481 :
482 : ab1 = (dpv*q - p*dq)/(q*q)
483 : ab2 = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
484 : ab3 = (d3p*q*q - p*q*d3q - 3.0_dp*dpv*q*d2q + 3.0_dp*p*dq*d2q)/(q*q*q)
485 : ab3 = ab3 - 3.0_dp*ab2*dq/q
486 : rsr1 = rs(ip)/(rho(ip)*rho(ip))
487 : rsr2 = f13*f13*rs(ip)*rsr1
488 : rsr3 = f13*rs(ip)*rsr2
489 : rsr1 = -f23*f23*f23*rsr1
490 : pot(ip) = pot(ip) + rsr1*ab1 + rsr2*ab2 + rsr3*ab3
491 :
492 : END IF
493 : END DO
494 :
495 0 : END SUBROUTINE pade_lda_3
496 :
497 : ! **************************************************************************************************
498 : !> \brief ...
499 : !> \param rhoa ...
500 : !> \param rhob ...
501 : !> \param rs ...
502 : !> \param fx ...
503 : !> \param pot0 ...
504 : ! **************************************************************************************************
505 50766496 : SUBROUTINE pade_lsd_0(rhoa, rhob, rs, fx, pot0)
506 :
507 : REAL(KIND=dp), INTENT(IN) :: rhoa, rhob, rs
508 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fx
509 : REAL(KIND=dp), INTENT(INOUT) :: pot0
510 :
511 : REAL(KIND=dp) :: fa0, fa1, fa2, fa3, fb1, fb2, fb3, fb4, &
512 : p, q, rhoab
513 :
514 50766496 : rhoab = rhoa + rhob
515 :
516 50766496 : IF (rhoab > eps_rho) THEN
517 :
518 48174092 : fa0 = a0 + fx(1)*da0
519 48174092 : fa1 = a1 + fx(1)*da1
520 48174092 : fa2 = a2 + fx(1)*da2
521 48174092 : fa3 = a3 + fx(1)*da3
522 48174092 : fb1 = b1 + fx(1)*db1
523 48174092 : fb2 = b2 + fx(1)*db2
524 48174092 : fb3 = b3 + fx(1)*db3
525 48174092 : fb4 = b4 + fx(1)*db4
526 :
527 48174092 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
528 48174092 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
529 :
530 48174092 : pot0 = pot0 - p/q*rhoab
531 :
532 : END IF
533 :
534 50766496 : END SUBROUTINE pade_lsd_0
535 :
536 : ! **************************************************************************************************
537 : !> \brief ...
538 : !> \param rhoa ...
539 : !> \param rhob ...
540 : !> \param rs ...
541 : !> \param fx ...
542 : !> \param pota ...
543 : !> \param potb ...
544 : ! **************************************************************************************************
545 0 : SUBROUTINE pade_lsd_1(rhoa, rhob, rs, fx, pota, potb)
546 :
547 : REAL(KIND=dp), INTENT(IN) :: rhoa, rhob, rs
548 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fx
549 : REAL(KIND=dp), INTENT(INOUT) :: pota, potb
550 :
551 : REAL(KIND=dp) :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
552 : fb1, fb2, fb3, fb4, p, q, rhoab, xp, xq
553 :
554 0 : rhoab = rhoa + rhob
555 :
556 0 : IF (rhoab > eps_rho) THEN
557 :
558 0 : fa0 = a0 + fx(1)*da0
559 0 : fa1 = a1 + fx(1)*da1
560 0 : fa2 = a2 + fx(1)*da2
561 0 : fa3 = a3 + fx(1)*da3
562 0 : fb1 = b1 + fx(1)*db1
563 0 : fb2 = b2 + fx(1)*db2
564 0 : fb3 = b3 + fx(1)*db3
565 0 : fb4 = b4 + fx(1)*db4
566 :
567 0 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
568 0 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
569 0 : dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
570 : dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
571 0 : 4.0_dp*fb4*rs)*rs)*rs
572 0 : xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
573 0 : xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
574 :
575 0 : dr = (dpv*q - p*dq)/(q*q)
576 0 : dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
577 0 : dc = f13*rs*dr - p/q
578 :
579 0 : pota = pota + dc - dx*rhob
580 0 : potb = potb + dc + dx*rhoa
581 :
582 : END IF
583 :
584 0 : END SUBROUTINE pade_lsd_1
585 :
586 : ! **************************************************************************************************
587 : !> \brief ...
588 : !> \param rhoa ...
589 : !> \param rhob ...
590 : !> \param rs ...
591 : !> \param fx ...
592 : !> \param pot0 ...
593 : !> \param pota ...
594 : !> \param potb ...
595 : ! **************************************************************************************************
596 472504503 : SUBROUTINE pade_lsd_01(rhoa, rhob, rs, fx, pot0, pota, potb)
597 :
598 : REAL(KIND=dp), INTENT(IN) :: rhoa, rhob, rs
599 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fx
600 : REAL(KIND=dp), INTENT(INOUT) :: pot0, pota, potb
601 :
602 : REAL(KIND=dp) :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
603 : fb1, fb2, fb3, fb4, p, q, rhoab, xp, xq
604 :
605 472504503 : rhoab = rhoa + rhob
606 :
607 472504503 : IF (rhoab > eps_rho) THEN
608 :
609 437842889 : fa0 = a0 + fx(1)*da0
610 437842889 : fa1 = a1 + fx(1)*da1
611 437842889 : fa2 = a2 + fx(1)*da2
612 437842889 : fa3 = a3 + fx(1)*da3
613 437842889 : fb1 = b1 + fx(1)*db1
614 437842889 : fb2 = b2 + fx(1)*db2
615 437842889 : fb3 = b3 + fx(1)*db3
616 437842889 : fb4 = b4 + fx(1)*db4
617 :
618 437842889 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
619 437842889 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
620 437842889 : dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
621 : dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
622 437842889 : 4.0_dp*fb4*rs)*rs)*rs
623 437842889 : xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
624 437842889 : xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
625 :
626 437842889 : dr = (dpv*q - p*dq)/(q*q)
627 437842889 : dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
628 437842889 : dc = f13*rs*dr - p/q
629 :
630 437842889 : pot0 = pot0 - p/q*rhoab
631 437842889 : pota = pota + dc - dx*rhob
632 437842889 : potb = potb + dc + dx*rhoa
633 :
634 : END IF
635 :
636 472504503 : END SUBROUTINE pade_lsd_01
637 :
638 : ! **************************************************************************************************
639 : !> \brief ...
640 : !> \param rhoa ...
641 : !> \param rhob ...
642 : !> \param rs ...
643 : !> \param fx ...
644 : !> \param potaa ...
645 : !> \param potab ...
646 : !> \param potbb ...
647 : ! **************************************************************************************************
648 30697597 : SUBROUTINE pade_lsd_2(rhoa, rhob, rs, fx, potaa, potab, potbb)
649 :
650 : REAL(KIND=dp), INTENT(IN) :: rhoa, rhob, rs
651 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fx
652 : REAL(KIND=dp), INTENT(INOUT) :: potaa, potab, potbb
653 :
654 : REAL(KIND=dp) :: d2p, d2q, dpv, dq, dr, drr, dx, dxp, &
655 : dxq, dxr, dxx, fa0, fa1, fa2, fa3, &
656 : fb1, fb2, fb3, fb4, or, p, q, rhoab, &
657 : xp, xq, xt, yt
658 :
659 30697597 : rhoab = rhoa + rhob
660 :
661 30697597 : IF (rhoab > eps_rho) THEN
662 :
663 30552303 : fa0 = a0 + fx(1)*da0
664 30552303 : fa1 = a1 + fx(1)*da1
665 30552303 : fa2 = a2 + fx(1)*da2
666 30552303 : fa3 = a3 + fx(1)*da3
667 30552303 : fb1 = b1 + fx(1)*db1
668 30552303 : fb2 = b2 + fx(1)*db2
669 30552303 : fb3 = b3 + fx(1)*db3
670 30552303 : fb4 = b4 + fx(1)*db4
671 :
672 30552303 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
673 30552303 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
674 :
675 30552303 : dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
676 : dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
677 30552303 : 4.0_dp*fb4*rs)*rs)*rs
678 :
679 30552303 : d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
680 30552303 : d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
681 :
682 30552303 : xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
683 30552303 : xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
684 :
685 30552303 : dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
686 : dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
687 30552303 : 4.0_dp*db4*rs)*rs)*rs
688 :
689 30552303 : dr = (dpv*q - p*dq)/(q*q)
690 30552303 : drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
691 30552303 : dx = (xp*q - p*xq)/(q*q)
692 30552303 : dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
693 30552303 : dxr = (dxp*q*q + dpv*xq*q - xp*dq*q - p*dxq*q - 2.0_dp*dpv*q*xq + 2.0_dp*p*dq*xq)/(q*q*q)
694 :
695 30552303 : or = 1.0_dp/rhoab
696 30552303 : yt = rhob*or
697 30552303 : xt = rhoa*or
698 :
699 : potaa = potaa + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
700 : + f43*rs*fx(2)*dxr*yt*or &
701 : - 4.0_dp*fx(2)*fx(2)*dxx*yt*yt*or &
702 30552303 : - 4.0_dp*dx*fx(3)*yt*yt*or
703 : potab = potab + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
704 : + f23*rs*fx(2)*dxr*(yt - xt)*or &
705 : + 4.0_dp*fx(2)*fx(2)*dxx*xt*yt*or &
706 30552303 : + 4.0_dp*dx*fx(3)*xt*yt*or
707 : potbb = potbb + f23*f13*dr*rs*or - f13*f13*drr*rs*rs*or &
708 : - f43*rs*fx(2)*dxr*xt*or &
709 : - 4.0_dp*fx(2)*fx(2)*dxx*xt*xt*or &
710 30552303 : - 4.0_dp*dx*fx(3)*xt*xt*or
711 :
712 : END IF
713 :
714 30697597 : END SUBROUTINE pade_lsd_2
715 :
716 : ! **************************************************************************************************
717 : !> \brief ...
718 : !> \param rhoa ...
719 : !> \param rhob ...
720 : !> \param rs ...
721 : !> \param fx ...
722 : !> \param potaaa ...
723 : !> \param potaab ...
724 : !> \param potabb ...
725 : !> \param potbbb ...
726 : ! **************************************************************************************************
727 0 : SUBROUTINE pade_lsd_3(rhoa, rhob, rs, fx, potaaa, potaab, potabb, potbbb)
728 :
729 : REAL(KIND=dp), INTENT(IN) :: rhoa, rhob, rs
730 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fx
731 : REAL(KIND=dp), INTENT(INOUT) :: potaaa, potaab, potabb, potbbb
732 :
733 : REAL(KIND=dp) :: d2p, d2q, d2xp, d2xq, d3p, d3q, dpv, dq, dr, drr, drrr, dx, dxp, dxq, dxr, &
734 : dxrr, dxx, dxxr, dxxx, fa0, fa1, fa2, fa3, fb1, fb2, fb3, fb4, or, p, q, rhoab, xp, xq, &
735 : xt, yt
736 :
737 0 : IF (.NOT. debug_flag) CPABORT("Routine not tested")
738 :
739 0 : rhoab = rhoa + rhob
740 :
741 0 : IF (rhoab > eps_rho) THEN
742 :
743 0 : fa0 = a0 + fx(1)*da0
744 0 : fa1 = a1 + fx(1)*da1
745 0 : fa2 = a2 + fx(1)*da2
746 0 : fa3 = a3 + fx(1)*da3
747 0 : fb1 = b1 + fx(1)*db1
748 0 : fb2 = b2 + fx(1)*db2
749 0 : fb3 = b3 + fx(1)*db3
750 0 : fb4 = b4 + fx(1)*db4
751 :
752 0 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
753 0 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
754 :
755 0 : dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
756 : dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
757 0 : 4.0_dp*fb4*rs)*rs)*rs
758 :
759 0 : d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
760 0 : d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
761 :
762 0 : d3p = 6.0_dp*fa3
763 0 : d3q = 6.0_dp*fb3 + 24.0_dp*fb4*rs
764 :
765 0 : xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
766 0 : xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
767 :
768 0 : dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
769 : dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
770 0 : 4.0_dp*db4*rs)*rs)*rs
771 :
772 0 : d2xp = 2.0_dp*da2 + 6.0_dp*da3*rs
773 0 : d2xq = 2.0_dp*db2 + (6.0_dp*db3 + 12.0_dp*db4*rs)*rs
774 :
775 0 : dr = (dpv*q - p*dq)/(q*q)
776 0 : drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
777 : drrr = (d3p*q*q*q - 3.0_dp*d2p*dq*q*q + 6.0_dp*dpv*dq*dq*q - 3.0_dp*dpv*d2q*q*q - &
778 0 : 6.0_dp*p*dq*dq*dq + 6.0_dp*p*dq*d2q*q - p*d3q*q*q)/(q*q*q*q)
779 0 : dx = (xp*q - p*xq)/(q*q)
780 0 : dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
781 0 : dxxx = 6.0_dp*xq*(q*xp*xq - p*xq*xq)/(q*q*q*q)
782 0 : dxr = (dxp*q*q + dpv*xq*q - xp*dq*q - p*dxq*q - 2.0_dp*dpv*q*xq + 2.0_dp*p*dq*xq)/(q*q*q)
783 : dxxr = 2.0_dp*(2.0_dp*dxq*q*p*xq - dxq*q*q*xp + xq*xq*q*dpv - xq*q*q*dxp + &
784 0 : 2.0_dp*xq*q*xp*dq - 3.0_dp*xq*xq*dq*p)/(q*q*q*q)
785 : dxrr = (q*q*q*d2xp - 2.0_dp*q*q*dxp*dq - q*q*xp*d2q - q*q*d2p*xq - &
786 : 2.0_dp*q*q*dpv*dxq - q*q*p*d2xq + 4.0_dp*dq*q*dpv*xq + 4.0_dp*dq*q*p*dxq + &
787 0 : 2.0_dp*dq*dq*q*xp - 6.0_dp*dq*dq*p*xq + 2.0_dp*d2q*q*p*xq)/(q*q*q*q)
788 :
789 0 : or = 1.0_dp/rhoab
790 0 : yt = rhob*or
791 0 : xt = rhoa*or
792 :
793 : potaaa = potaaa + 8.0_dp/27.0_dp*dr*rs*or*or + &
794 : 1.0_dp/9.0_dp*drr*rs*rs*or*or + &
795 : 1.0_dp/27.0_dp*drrr*rs**3*or*or + &
796 0 : dxr*or*or*yt*rs*(-8.0_dp/3.0_dp*fx(2) + 4.0_dp*fx(3)*yt)
797 0 : potaab = potaab + 0.0_dp
798 0 : potabb = potabb + 0.0_dp
799 0 : potbbb = potbbb + 0.0_dp
800 :
801 : END IF
802 :
803 0 : END SUBROUTINE pade_lsd_3
804 :
805 : ! **************************************************************************************************
806 : !> \brief ...
807 : !> \param rho_a ...
808 : !> \param rho_b ...
809 : !> \param fxc_aa ...
810 : !> \param fxc_ab ...
811 : !> \param fxc_bb ...
812 : ! **************************************************************************************************
813 2 : SUBROUTINE pade_fxc_eval(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb)
814 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
815 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: fxc_aa, fxc_ab, fxc_bb
816 :
817 : INTEGER :: i, j, k
818 : INTEGER, DIMENSION(2, 3) :: bo
819 : REAL(KIND=dp) :: eaa, eab, ebb, rhoa, rhob, rs
820 : REAL(KIND=dp), DIMENSION(4) :: fx
821 :
822 20 : bo(1:2, 1:3) = rho_a%pw_grid%bounds_local(1:2, 1:3)
823 : !$OMP PARALLEL DO PRIVATE(i,j,k,fx,rhoa,rhob,rs,eaa,eab,ebb) DEFAULT(NONE)&
824 2 : !$OMP SHARED(bo,rho_a,rho_b,fxc_aa,fxc_ab,fxc_bb)
825 : DO k = bo(1, 3), bo(2, 3)
826 : DO j = bo(1, 2), bo(2, 2)
827 : DO i = bo(1, 1), bo(2, 1)
828 :
829 : rhoa = rho_a%array(i, j, k)
830 : rhob = rho_b%array(i, j, k)
831 : fx(1) = rhoa + rhob
832 :
833 : CALL calc_rs(fx(1), rs)
834 : CALL calc_fx(rhoa, rhob, fx, 2)
835 :
836 : eaa = 0.0_dp; eab = 0.0_dp; ebb = 0.0_dp
837 : CALL pade_lsd_2(rhoa, rhob, rs, fx, eaa, eab, ebb)
838 :
839 : fxc_aa%array(i, j, k) = eaa
840 : fxc_ab%array(i, j, k) = eab
841 : fxc_bb%array(i, j, k) = ebb
842 :
843 : END DO
844 : END DO
845 : END DO
846 :
847 2 : END SUBROUTINE pade_fxc_eval
848 :
849 : END MODULE xc_pade
850 :
|
