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