(git:374b731)
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-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! **************************************************************************************************
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)) &
122 cpabort("Arguments mismatch.")
123 IF (lsd) THEN
124 needs%rho_spin = .true.
125 ELSE
126 needs%rho = .true.
127 END IF
128 END IF
129
130 IF (PRESENT(max_deriv)) max_deriv = 3
131
132 END SUBROUTINE pade_info
133
134! **************************************************************************************************
135!> \brief ...
136!> \param deriv_set ...
137!> \param rho_set ...
138!> \param order ...
139! **************************************************************************************************
140 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rs
149 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
150 POINTER :: e_0, e_r, e_rr, e_rrr
151 TYPE(xc_derivative_type), POINTER :: deriv
152
153 calc = .false.
154 IF (order >= 0) calc(0:order) = .true.
155 IF (order < 0) calc(-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 IF (order >= 0) calc(0:order) = .true.
221 IF (order < 0) calc(-order) = .true.
222
223 IF (calc(0)) THEN
224 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
225 allocate_deriv=.true.)
226 CALL xc_derivative_get(deriv, deriv_data=e_0)
227 END IF
228 IF (calc(1)) THEN
229 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
230 allocate_deriv=.true.)
231 CALL xc_derivative_get(deriv, deriv_data=e_ra)
232 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
233 allocate_deriv=.true.)
234 CALL xc_derivative_get(deriv, deriv_data=e_rb)
235 END IF
236 IF (calc(2)) THEN
237 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
238 allocate_deriv=.true.)
239 CALL xc_derivative_get(deriv, deriv_data=e_rara)
240 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
241 allocate_deriv=.true.)
242 CALL xc_derivative_get(deriv, deriv_data=e_rarb)
243 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
244 allocate_deriv=.true.)
245 CALL xc_derivative_get(deriv, deriv_data=e_rbrb)
246 END IF
247 IF (calc(3)) THEN
248 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
249 allocate_deriv=.true.)
250 CALL xc_derivative_get(deriv, deriv_data=e_rarara)
251 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhob], &
252 allocate_deriv=.true.)
253 CALL xc_derivative_get(deriv, deriv_data=e_rararb)
254 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob, deriv_rhob], &
255 allocate_deriv=.true.)
256 CALL xc_derivative_get(deriv, deriv_data=e_rarbrb)
257 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
258 allocate_deriv=.true.)
259 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!$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 END SUBROUTINE pade_lsd_pw_eval
297
298! **************************************************************************************************
299!> \brief ...
300!> \param n ...
301!> \param rho ...
302!> \param rs ...
303!> \param pot ...
304! **************************************************************************************************
305 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!$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 END SUBROUTINE pade_lda_0
326
327! **************************************************************************************************
328!> \brief ...
329!> \param n ...
330!> \param rho ...
331!> \param rs ...
332!> \param pot ...
333! **************************************************************************************************
334 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!$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 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 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!$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 END SUBROUTINE pade_lda_01
402
403! **************************************************************************************************
404!> \brief ...
405!> \param n ...
406!> \param rho ...
407!> \param rs ...
408!> \param pot ...
409! **************************************************************************************************
410 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!$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 END SUBROUTINE pade_lda_2
445
446! **************************************************************************************************
447!> \brief ...
448!> \param n ...
449!> \param rho ...
450!> \param rs ...
451!> \param pot ...
452! **************************************************************************************************
453 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!$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 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 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 rhoab = rhoa + rhob
514
515 IF (rhoab > eps_rho) THEN
516
517 fa0 = a0 + fx(1)*da0
518 fa1 = a1 + fx(1)*da1
519 fa2 = a2 + fx(1)*da2
520 fa3 = a3 + fx(1)*da3
521 fb1 = b1 + fx(1)*db1
522 fb2 = b2 + fx(1)*db2
523 fb3 = b3 + fx(1)*db3
524 fb4 = b4 + fx(1)*db4
525
526 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
527 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
528
529 pot0 = pot0 - p/q*rhoab
530
531 END IF
532
533 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 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 rhoab = rhoa + rhob
554
555 IF (rhoab > eps_rho) THEN
556
557 fa0 = a0 + fx(1)*da0
558 fa1 = a1 + fx(1)*da1
559 fa2 = a2 + fx(1)*da2
560 fa3 = a3 + fx(1)*da3
561 fb1 = b1 + fx(1)*db1
562 fb2 = b2 + fx(1)*db2
563 fb3 = b3 + fx(1)*db3
564 fb4 = b4 + fx(1)*db4
565
566 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
567 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
568 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
569 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
570 4.0_dp*fb4*rs)*rs)*rs
571 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
572 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
573
574 dr = (dpv*q - p*dq)/(q*q)
575 dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
576 dc = f13*rs*dr - p/q
577
578 pota = pota + dc - dx*rhob
579 potb = potb + dc + dx*rhoa
580
581 END IF
582
583 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 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 rhoab = rhoa + rhob
605
606 IF (rhoab > eps_rho) THEN
607
608 fa0 = a0 + fx(1)*da0
609 fa1 = a1 + fx(1)*da1
610 fa2 = a2 + fx(1)*da2
611 fa3 = a3 + fx(1)*da3
612 fb1 = b1 + fx(1)*db1
613 fb2 = b2 + fx(1)*db2
614 fb3 = b3 + fx(1)*db3
615 fb4 = b4 + fx(1)*db4
616
617 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
618 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
619 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
620 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
621 4.0_dp*fb4*rs)*rs)*rs
622 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
623 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
624
625 dr = (dpv*q - p*dq)/(q*q)
626 dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx(2)/rhoab
627 dc = f13*rs*dr - p/q
628
629 pot0 = pot0 - p/q*rhoab
630 pota = pota + dc - dx*rhob
631 potb = potb + dc + dx*rhoa
632
633 END IF
634
635 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 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 rhoab = rhoa + rhob
659
660 IF (rhoab > eps_rho) THEN
661
662 fa0 = a0 + fx(1)*da0
663 fa1 = a1 + fx(1)*da1
664 fa2 = a2 + fx(1)*da2
665 fa3 = a3 + fx(1)*da3
666 fb1 = b1 + fx(1)*db1
667 fb2 = b2 + fx(1)*db2
668 fb3 = b3 + fx(1)*db3
669 fb4 = b4 + fx(1)*db4
670
671 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
672 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
673
674 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
675 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
676 4.0_dp*fb4*rs)*rs)*rs
677
678 d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
679 d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
680
681 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
682 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
683
684 dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
685 dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
686 4.0_dp*db4*rs)*rs)*rs
687
688 dr = (dpv*q - p*dq)/(q*q)
689 drr = (d2p*q*q - p*q*d2q - 2.0_dp*dpv*q*dq + 2.0_dp*p*dq*dq)/(q*q*q)
690 dx = (xp*q - p*xq)/(q*q)
691 dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
692 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 or = 1.0_dp/rhoab
695 yt = rhob*or
696 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 - 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 + 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 - 4.0_dp*dx*fx(3)*xt*xt*or
710
711 END IF
712
713 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 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 IF (.NOT. debug_flag) cpabort("Routine not tested")
737
738 rhoab = rhoa + rhob
739
740 IF (rhoab > eps_rho) THEN
741
742 fa0 = a0 + fx(1)*da0
743 fa1 = a1 + fx(1)*da1
744 fa2 = a2 + fx(1)*da2
745 fa3 = a3 + fx(1)*da3
746 fb1 = b1 + fx(1)*db1
747 fb2 = b2 + fx(1)*db2
748 fb3 = b3 + fx(1)*db3
749 fb4 = b4 + fx(1)*db4
750
751 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
752 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
753
754 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
755 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
756 4.0_dp*fb4*rs)*rs)*rs
757
758 d2p = 2.0_dp*fa2 + 6.0_dp*fa3*rs
759 d2q = 2.0_dp*fb2 + (6.0_dp*fb3 + 12.0_dp*fb4*rs)*rs
760
761 d3p = 6.0_dp*fa3
762 d3q = 6.0_dp*fb3 + 24.0_dp*fb4*rs
763
764 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
765 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
766
767 dxp = da1 + (2.0_dp*da2 + 3.0_dp*da3*rs)*rs
768 dxq = db1 + (2.0_dp*db2 + (3.0_dp*db3 + &
769 4.0_dp*db4*rs)*rs)*rs
770
771 d2xp = 2.0_dp*da2 + 6.0_dp*da3*rs
772 d2xq = 2.0_dp*db2 + (6.0_dp*db3 + 12.0_dp*db4*rs)*rs
773
774 dr = (dpv*q - p*dq)/(q*q)
775 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 6.0_dp*p*dq*dq*dq + 6.0_dp*p*dq*d2q*q - p*d3q*q*q)/(q*q*q*q)
778 dx = (xp*q - p*xq)/(q*q)
779 dxx = 2.0_dp*xq*(p*xq - xp*q)/(q*q*q)
780 dxxx = 6.0_dp*xq*(q*xp*xq - p*xq*xq)/(q*q*q*q)
781 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 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 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 or = 1.0_dp/rhoab
789 yt = rhob*or
790 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 dxr*or*or*yt*rs*(-8.0_dp/3.0_dp*fx(2) + 4.0_dp*fx(3)*yt)
796 potaab = potaab + 0.0_dp
797 potabb = potabb + 0.0_dp
798 potbbb = potbbb + 0.0_dp
799
800 END IF
801
802 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 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 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!$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 END SUBROUTINE pade_fxc_eval
847
848END MODULE xc_pade
849
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:28
bibliography::goedecker1996
integer, save, public goedecker1996
Definition bibliography.F:43
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:66
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:813
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:141
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:49
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