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