(git:47b47fb)
Loading...
Searching...
No Matches
xc_pade.F
Go to the documentation of this file.
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! **************************************************************************************************
19MODULE 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
72CONTAINS
73
74! **************************************************************************************************
75!> \brief ...
76!> \param cutoff ...
77!> \param debug ...
78! **************************************************************************************************
79 SUBROUTINE pade_init(cutoff, debug)
80
81 REAL(kind=dp), INTENT(IN) :: cutoff
82 LOGICAL, INTENT(IN), OPTIONAL :: debug
83
84 eps_rho = cutoff
85 CALL set_util(cutoff)
86
87 CALL cite_reference(goedecker1996)
88
89 IF (PRESENT(debug)) THEN
90 debug_flag = debug
91 ELSE
92 debug_flag = .false.
93 END IF
94
95 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 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 IF (PRESENT(reference)) THEN
113 reference = "S. Goedecker, M. Teter and J. Hutter," &
114 //" Phys. Rev. B 54, 1703 (1996)"
115 END IF
116 IF (PRESENT(shortform)) THEN
117 shortform = "S. Goedecker et al., PRB 54, 1703 (1996)"
118 END IF
119
120 IF (PRESENT(needs)) THEN
121 IF (.NOT. PRESENT(lsd)) THEN
122 cpabort("Arguments mismatch.")
123 END IF
124 IF (lsd) THEN
125 needs%rho_spin = .true.
126 ELSE
127 needs%rho = .true.
128 END IF
129 END IF
130
131 IF (PRESENT(max_deriv)) max_deriv = 3
132
133 END SUBROUTINE pade_info
134
135! **************************************************************************************************
136!> \brief ...
137!> \param deriv_set ...
138!> \param rho_set ...
139!> \param order ...
140! **************************************************************************************************
141 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rs
150 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
151 POINTER :: e_0, e_r, e_rr, e_rrr
152 TYPE(xc_derivative_type), POINTER :: deriv
153
154 calc = .false.
155 calc(0:order) = .true.
156
157 n = product(rho_set%local_bounds(2, :) - rho_set%local_bounds(1, :) + [1, 1, 1])
158 ALLOCATE (rs(n))
159
160 CALL calc_rs_pw(rho_set%rho, rs, n)
161 IF (calc(0) .AND. calc(1)) THEN
162 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
163 allocate_deriv=.true.)
164 CALL xc_derivative_get(deriv, deriv_data=e_0)
165 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
166 allocate_deriv=.true.)
167 CALL xc_derivative_get(deriv, deriv_data=e_r)
168 CALL pade_lda_01(n, rho_set%rho, rs, e_0, e_r)
169 ELSE IF (calc(0)) THEN
170 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
171 allocate_deriv=.true.)
172 CALL xc_derivative_get(deriv, deriv_data=e_0)
173 CALL pade_lda_0(n, rho_set%rho, rs, e_0)
174 ELSE IF (calc(1)) THEN
175 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
176 allocate_deriv=.true.)
177 CALL xc_derivative_get(deriv, deriv_data=e_r)
178 CALL pade_lda_1(n, rho_set%rho, rs, e_r)
179 END IF
180 IF (calc(2)) THEN
181 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
182 allocate_deriv=.true.)
183 CALL xc_derivative_get(deriv, deriv_data=e_rr)
184 CALL pade_lda_2(n, rho_set%rho, rs, e_rr)
185 END IF
186 IF (calc(3)) THEN
187 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
188 allocate_deriv=.true.)
189 CALL xc_derivative_get(deriv, deriv_data=e_rrr)
190 CALL pade_lda_3(n, rho_set%rho, rs, e_rrr)
191 END IF
192
193 DEALLOCATE (rs)
194
195 END SUBROUTINE pade_lda_pw_eval
196
197! **************************************************************************************************
198!> \brief ...
199!> \param deriv_set ...
200!> \param rho_set ...
201!> \param order ...
202! **************************************************************************************************
203 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 POINTER :: e_0, e_ra, e_rara, e_rarara, e_rararb, &
214 e_rarb, e_rarbrb, e_rb, e_rbrb, &
215 e_rbrbrb
216 REAL(kind=dp), DIMENSION(4) :: fx
217 TYPE(xc_derivative_type), POINTER :: deriv
218
219 calc = .false.
220 calc(0:order) = .true.
221
222 IF (calc(0)) THEN
223 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
224 allocate_deriv=.true.)
225 CALL xc_derivative_get(deriv, deriv_data=e_0)
226 END IF
227 IF (calc(1)) THEN
228 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
229 allocate_deriv=.true.)
230 CALL xc_derivative_get(deriv, deriv_data=e_ra)
231 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
232 allocate_deriv=.true.)
233 CALL xc_derivative_get(deriv, deriv_data=e_rb)
234 END IF
235 IF (calc(2)) THEN
236 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
237 allocate_deriv=.true.)
238 CALL xc_derivative_get(deriv, deriv_data=e_rara)
239 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
240 allocate_deriv=.true.)
241 CALL xc_derivative_get(deriv, deriv_data=e_rarb)
242 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
243 allocate_deriv=.true.)
244 CALL xc_derivative_get(deriv, deriv_data=e_rbrb)
245 END IF
246 IF (calc(3)) THEN
247 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
248 allocate_deriv=.true.)
249 CALL xc_derivative_get(deriv, deriv_data=e_rarara)
250 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
251 allocate_deriv=.true.)
252 CALL xc_derivative_get(deriv, deriv_data=e_rararb)
253 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
254 allocate_deriv=.true.)
255 CALL xc_derivative_get(deriv, deriv_data=e_rarbrb)
256 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
257 allocate_deriv=.true.)
258 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!$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 END SUBROUTINE pade_lsd_pw_eval
296
297! **************************************************************************************************
298!> \brief ...
299!> \param n ...
300!> \param rho ...
301!> \param rs ...
302!> \param pot ...
303! **************************************************************************************************
304 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!$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 END SUBROUTINE pade_lda_0
325
326! **************************************************************************************************
327!> \brief ...
328!> \param n ...
329!> \param rho ...
330!> \param rs ...
331!> \param pot ...
332! **************************************************************************************************
333 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!$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 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 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!$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 END SUBROUTINE pade_lda_01
401
402! **************************************************************************************************
403!> \brief ...
404!> \param n ...
405!> \param rho ...
406!> \param rs ...
407!> \param pot ...
408! **************************************************************************************************
409 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!$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 END SUBROUTINE pade_lda_2
444
445! **************************************************************************************************
446!> \brief ...
447!> \param n ...
448!> \param rho ...
449!> \param rs ...
450!> \param pot ...
451! **************************************************************************************************
452 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!$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 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 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 rhoab = rhoa + rhob
513
514 IF (rhoab > eps_rho) THEN
515
516 fa0 = a0 + fx(1)*da0
517 fa1 = a1 + fx(1)*da1
518 fa2 = a2 + fx(1)*da2
519 fa3 = a3 + fx(1)*da3
520 fb1 = b1 + fx(1)*db1
521 fb2 = b2 + fx(1)*db2
522 fb3 = b3 + fx(1)*db3
523 fb4 = b4 + fx(1)*db4
524
525 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
526 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
527
528 pot0 = pot0 - p/q*rhoab
529
530 END IF
531
532 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 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 rhoab = rhoa + rhob
553
554 IF (rhoab > eps_rho) THEN
555
556 fa0 = a0 + fx(1)*da0
557 fa1 = a1 + fx(1)*da1
558 fa2 = a2 + fx(1)*da2
559 fa3 = a3 + fx(1)*da3
560 fb1 = b1 + fx(1)*db1
561 fb2 = b2 + fx(1)*db2
562 fb3 = b3 + fx(1)*db3
563 fb4 = b4 + fx(1)*db4
564
565 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
566 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
567 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
568 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
569 4.0_dp*fb4*rs)*rs)*rs
570 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
571 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
572
573 dr = (dpv*q - p*dq)/(q*q)
574 dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
575 dc = f13*rs*dr - p/q
576
577 pota = pota + dc - dx*rhob
578 potb = potb + dc + dx*rhoa
579
580 END IF
581
582 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 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 rhoab = rhoa + rhob
604
605 IF (rhoab > eps_rho) THEN
606
607 fa0 = a0 + fx(1)*da0
608 fa1 = a1 + fx(1)*da1
609 fa2 = a2 + fx(1)*da2
610 fa3 = a3 + fx(1)*da3
611 fb1 = b1 + fx(1)*db1
612 fb2 = b2 + fx(1)*db2
613 fb3 = b3 + fx(1)*db3
614 fb4 = b4 + fx(1)*db4
615
616 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
617 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
618 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
619 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
620 4.0_dp*fb4*rs)*rs)*rs
621 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
622 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
623
624 dr = (dpv*q - p*dq)/(q*q)
625 dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
626 dc = f13*rs*dr - p/q
627
628 pot0 = pot0 - p/q*rhoab
629 pota = pota + dc - dx*rhob
630 potb = potb + dc + dx*rhoa
631
632 END IF
633
634 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 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 rhoab = rhoa + rhob
658
659 IF (rhoab > eps_rho) THEN
660
661 fa0 = a0 + fx(1)*da0
662 fa1 = a1 + fx(1)*da1
663 fa2 = a2 + fx(1)*da2
664 fa3 = a3 + fx(1)*da3
665 fb1 = b1 + fx(1)*db1
666 fb2 = b2 + fx(1)*db2
667 fb3 = b3 + fx(1)*db3
668 fb4 = b4 + fx(1)*db4
669
670 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
671 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
672
673 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
674 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
675 4.0_dp*fb4*rs)*rs)*rs
676
677 d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
678 d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
679
680 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
681 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
682
683 dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
684 dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
685 4.0_dp*db4*rs)*rs)*rs
686
687 dr = (dpv*q - p*dq)/(q*q)
688 drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
689 dx = (xp*q - p*xq)/(q*q)
690 dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
691 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 or = 1.0_dp/rhoab
694 yt = rhob*or
695 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 - 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 + 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 - 4.0_dp*dx*fx(3)*xt*xt*or
709
710 END IF
711
712 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 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 IF (.NOT. debug_flag) cpabort("Routine not tested")
736
737 rhoab = rhoa + rhob
738
739 IF (rhoab > eps_rho) THEN
740
741 fa0 = a0 + fx(1)*da0
742 fa1 = a1 + fx(1)*da1
743 fa2 = a2 + fx(1)*da2
744 fa3 = a3 + fx(1)*da3
745 fb1 = b1 + fx(1)*db1
746 fb2 = b2 + fx(1)*db2
747 fb3 = b3 + fx(1)*db3
748 fb4 = b4 + fx(1)*db4
749
750 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
751 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
752
753 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
754 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
755 4.0_dp*fb4*rs)*rs)*rs
756
757 d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
758 d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
759
760 d3p = 6.0_dp*fa3
761 d3q = 6.0_dp*fb3 + 24.0_dp*fb4*rs
762
763 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
764 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
765
766 dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
767 dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
768 4.0_dp*db4*rs)*rs)*rs
769
770 d2xp = 2.0_dp*da2 + 6.0_dp*da3*rs
771 d2xq = 2.0_dp*db2 + (6.0_dp*db3 + 12.0_dp*db4*rs)*rs
772
773 dr = (dpv*q - p*dq)/(q*q)
774 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 6.0_dp*p*dq*dq*dq + 6.0_dp*p*dq*d2q*q - p*d3q*q*q)/(q*q*q*q)
777 dx = (xp*q - p*xq)/(q*q)
778 dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
779 dxxx = 6.0_dp*xq*(q*xp*xq - p*xq*xq)/(q*q*q*q)
780 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 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 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 or = 1.0_dp/rhoab
788 yt = rhob*or
789 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 dxr*or*or*yt*rs*(-8.0_dp/3.0_dp*fx(2) + 4.0_dp*fx(3)*yt)
795 potaab = potaab + 0.0_dp
796 potabb = potabb + 0.0_dp
797 potbbb = potbbb + 0.0_dp
798
799 END IF
800
801 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 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 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!$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 END SUBROUTINE pade_fxc_eval
846
847END MODULE xc_pade
848
xc_functionals_utilities::calc_fx
Definition xc_functionals_utilities.F:37
xc_functionals_utilities::calc_rs
Definition xc_functionals_utilities.F:41
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::goedecker1996
integer, save, public goedecker1996
Definition bibliography.F:36
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
pw_types
Definition pw_types.F:24
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_functionals_utilities
Utility routines for the functional calculations.
Definition xc_functionals_utilities.F:15
xc_functionals_utilities::calc_rs_pw
subroutine, public calc_rs_pw(rho, rs, n)
...
Definition xc_functionals_utilities.F:117
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition xc_functionals_utilities.F:54
xc_pade
Calculate the LDA functional in the Pade approximation Literature: S. Goedecker, M....
Definition xc_pade.F:19
xc_pade::pade_fxc_eval
subroutine, public pade_fxc_eval(rho_a, rho_b, fxc_aa, fxc_ab, fxc_bb)
...
Definition xc_pade.F:812
xc_pade::pade_init
subroutine, public pade_init(cutoff, debug)
...
Definition xc_pade.F:80
xc_pade::pade_lsd_pw_eval
subroutine, public pade_lsd_pw_eval(deriv_set, rho_set, order)
...
Definition xc_pade.F:204
xc_pade::pade_info
subroutine, public pade_info(reference, shortform, lsd, needs, max_deriv)
...
Definition xc_pade.F:106
xc_pade::pade_lda_pw_eval
subroutine, public pade_lda_pw_eval(deriv_set, rho_set, order)
...
Definition xc_pade.F:142
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78