(git:ccc2433)
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 ! **************************************************************************************************
19 MODULE xc_pade
20  USE bibliography, ONLY: goedecker1996,&
21  cite_reference
22  USE kinds, ONLY: dp
23  USE pw_types, ONLY: pw_r3d_rs_type
24  USE xc_derivative_desc, ONLY: deriv_rho,&
25  deriv_rhoa,&
26  deriv_rhob
27  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
28  xc_dset_get_derivative
29  USE xc_derivative_types, ONLY: xc_derivative_get,&
30  xc_derivative_type
31  USE xc_functionals_utilities, ONLY: calc_fx,&
32  calc_rs,&
33  calc_rs_pw,&
34  set_util
35  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
36  USE xc_rho_set_types, ONLY: xc_rho_set_type
37 #include "../base/base_uses.f90"
38 
39  IMPLICIT NONE
40 
41  PRIVATE
42 
43  REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
44  f23 = 2.0_dp*f13, &
45  f43 = 4.0_dp*f13
46 
47  REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429e+0_dp, &
48  a1 = 0.2217058676663745e+1_dp, &
49  a2 = 0.7405551735357053e+0_dp, &
50  a3 = 0.1968227878617998e-1_dp, &
51  b1 = 1.0000000000000000e+0_dp, &
52  b2 = 0.4504130959426697e+1_dp, &
53  b3 = 0.1110667363742916e+1_dp, &
54  b4 = 0.2359291751427506e-1_dp
55 
56  REAL(KIND=dp), PARAMETER :: da0 = 0.119086804055547e+0_dp, &
57  da1 = 0.6157402568883345e+0_dp, &
58  da2 = 0.1574201515892867e+0_dp, &
59  da3 = 0.3532336663397157e-2_dp, &
60  db1 = 0.0000000000000000e+0_dp, &
61  db2 = 0.2673612973836267e+0_dp, &
62  db3 = 0.2052004607777787e+0_dp, &
63  db4 = 0.4200005045691381e-2_dp
64 
65  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pade'
66 
67  PUBLIC :: pade_lda_pw_eval, pade_lsd_pw_eval, pade_info, pade_init, pade_fxc_eval
68 
69  REAL(KIND=dp) :: eps_rho
70  LOGICAL :: debug_flag
71 
72 CONTAINS
73 
74 ! **************************************************************************************************
75 !> \brief ...
76 !> \param cutoff ...
77 !> \param debug ...
78 ! **************************************************************************************************
79  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 
848 END MODULE xc_pade
849 
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