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