(git:1f285aa)
xc_functionals_utilities.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 Utility routines for the functional calculations
10 !> \par History
11 !> JGH (20.02.2001) : Added setup routine
12 !> JGH (26.02.2003) : OpenMP enabled
13 !> \author JGH (15.02.2002)
14 ! **************************************************************************************************
15 MODULE xc_functionals_utilities
16 
17  USE kinds, ONLY: dp
18  USE mathconstants, ONLY: pi
19 #include "../base/base_uses.f90"
20 
21  IMPLICIT NONE
22 
23  PRIVATE
24 
25  REAL(KIND=dp), PARAMETER :: rsfac = 0.6203504908994000166680065_dp ! (4*pi/3)^(-1/3)
26  REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
27  f23 = 2.0_dp*f13, &
28  f43 = 4.0_dp*f13, &
29  f53 = 5.0_dp*f13
30  REAL(KIND=dp), PARAMETER :: fxfac = 1.923661050931536319759455_dp ! 1/(2^(4/3) - 2)
31 
32  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_functionals_utilities'
33 
34  PUBLIC :: set_util, calc_rs, calc_fx, &
35  calc_wave_vector, calc_z, calc_rs_pw, calc_srs_pw
36 
37  INTERFACE calc_fx
38  MODULE PROCEDURE calc_fx_array, calc_fx_single
39  END INTERFACE
40 
41  INTERFACE calc_rs
42  MODULE PROCEDURE calc_rs_array, calc_rs_single
43  END INTERFACE
44 
45  REAL(KIND=dp), SAVE :: eps_rho
46 
47 CONTAINS
48 
49 ! **************************************************************************************************
50 !> \brief ...
51 !> \param cutoff ...
52 ! **************************************************************************************************
53  SUBROUTINE set_util(cutoff)
54 
55  REAL(kind=dp) :: cutoff
56 
57  eps_rho = cutoff
58 
59  END SUBROUTINE set_util
60 
61 ! **************************************************************************************************
62 !> \brief ...
63 !> \param rho ...
64 !> \param rs ...
65 ! **************************************************************************************************
66  ELEMENTAL SUBROUTINE calc_rs_single(rho, rs)
67 
68 ! rs parameter : f*rho**(-1/3)
69 
70  REAL(kind=dp), INTENT(IN) :: rho
71  REAL(kind=dp), INTENT(OUT) :: rs
72 
73  IF (rho < eps_rho) THEN
74  rs = 0.0_dp
75  ELSE
76  rs = rsfac*rho**(-f13)
77  END IF
78 
79  END SUBROUTINE calc_rs_single
80 
81 ! **************************************************************************************************
82 !> \brief ...
83 !> \param rho ...
84 !> \param rs ...
85 ! **************************************************************************************************
86  SUBROUTINE calc_rs_array(rho, rs)
87 
88 ! rs parameter : f*rho**(-1/3)
89 
90  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rho
91  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: rs
92 
93  INTEGER :: k
94 
95  IF (SIZE(rs) < SIZE(rho)) THEN
96  cpabort("Size of array rs too small.")
97  END IF
98 
99 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(rs,eps_rho,rho)
100  DO k = 1, SIZE(rs)
101  IF (rho(k) < eps_rho) THEN
102  rs(k) = 0.0_dp
103  ELSE
104  rs(k) = rsfac*rho(k)**(-f13)
105  END IF
106  END DO
107 
108  END SUBROUTINE calc_rs_array
109 
110 ! **************************************************************************************************
111 !> \brief ...
112 !> \param rho ...
113 !> \param rs ...
114 !> \param n ...
115 ! **************************************************************************************************
116  SUBROUTINE calc_rs_pw(rho, rs, n)
117 
118 ! rs parameter : f*rho**(-1/3)
119 
120  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho
121  REAL(kind=dp), DIMENSION(*), INTENT(OUT) :: rs
122  INTEGER, INTENT(IN) :: n
123 
124  INTEGER :: k
125 
126 !$OMP PARALLEL DO PRIVATE(k) SHARED(n,rs,rho,eps_rho) DEFAULT(NONE)
127  DO k = 1, n
128  IF (rho(k) < eps_rho) THEN
129  rs(k) = 0.0_dp
130  ELSE
131  rs(k) = rsfac*rho(k)**(-f13)
132  END IF
133  END DO
134 
135  END SUBROUTINE calc_rs_pw
136 
137 ! **************************************************************************************************
138 !> \brief ...
139 !> \param rho ...
140 !> \param x ...
141 !> \param n ...
142 ! **************************************************************************************************
143  SUBROUTINE calc_srs_pw(rho, x, n)
144 
145 ! rs parameter : f*rho**(-1/3)
146 ! x = sqrt(rs)
147 
148  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho
149  REAL(kind=dp), DIMENSION(*), INTENT(OUT) :: x
150  INTEGER, INTENT(in) :: n
151 
152  INTEGER :: ip
153 
154  CALL calc_rs_pw(rho, x, n)
155 
156 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) SHARED(x,n)
157  DO ip = 1, n
158  x(ip) = sqrt(x(ip))
159  END DO
160 
161  END SUBROUTINE calc_srs_pw
162 
163 ! **************************************************************************************************
164 !> \brief ...
165 !> \param tag ...
166 !> \param rho ...
167 !> \param grho ...
168 !> \param s ...
169 ! **************************************************************************************************
170  SUBROUTINE calc_wave_vector(tag, rho, grho, s)
171 
172 ! wave vector s = |nabla rho| / (2(3pi^2)^1/3 * rho^4/3)
173 
174  CHARACTER(len=*), INTENT(IN) :: tag
175  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho
176  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: s
177 
178  INTEGER :: ip, n
179  REAL(kind=dp) :: fac
180 
181 ! TAGS: U: total density, spin wave vector
182 ! R: spin density, total density wave vector
183 
184  fac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
185  IF (tag(1:1) == "u" .OR. tag(1:1) == "U") fac = fac*(2.0_dp)**f13
186  IF (tag(1:1) == "r" .OR. tag(1:1) == "R") fac = fac*(2.0_dp)**f13
187 
188  n = SIZE(s) !FM it was sizederiv_rho
189  !FM IF ( n > SIZE(s) ) &
190  !FM CPABORT("Incompatible array sizes.")
191  !FM IF ( n > SIZE(grho) ) &
192  !FM CPABORT("Incompatible array sizes.")
193 
194 !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) SHARED(rho,eps_rho,s,fac,grho,n)
195  DO ip = 1, n
196  IF (rho(ip) < eps_rho) THEN
197  s(ip) = 0.0_dp
198  ELSE
199  s(ip) = fac*grho(ip)*rho(ip)**(-f43)
200  END IF
201  END DO
202 
203  END SUBROUTINE calc_wave_vector
204 
205 ! **************************************************************************************************
206 !> \brief ...
207 !> \param n ...
208 !> \param rhoa ...
209 !> \param rhob ...
210 !> \param fx ...
211 !> \param m ...
212 ! **************************************************************************************************
213  SUBROUTINE calc_fx_array(n, rhoa, rhob, fx, m)
214 
215 ! spin interpolation function and derivatives
216 !
217 ! f(x) = ( (1+x)^(4/3) + (1-x)^(4/3) - 2 ) / (2^(4/3)-2)
218 ! df(x) = (4/3)( (1+x)^(1/3) - (1-x)^(1/3) ) / (2^(4/3)-2)
219 ! d2f(x) = (4/9)( (1+x)^(-2/3) + (1-x)^(-2/3) ) / (2^(4/3)-2)
220 ! d3f(x) = (-8/27)( (1+x)^(-5/3) - (1-x)^(-5/3) ) / (2^(4/3)-2)
221 !
222 
223  INTEGER, INTENT(IN) :: n
224  REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob
225  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: fx
226  INTEGER, INTENT(IN) :: m
227 
228  INTEGER :: ip
229  REAL(kind=dp) :: rhoab, x
230 
231 ! number of points
232 ! order of derivatives
233 ! *** Parameters ***
234 
235  IF (m > 3) cpabort("Order too high.")
236 !! IF (.NOT.ASSOCIATED(fx)) THEN
237 !! ALLOCATE(fx(n,m+1), STAT=ierr)
238 !! IF (ierr /= 0) CALL stop_memory(routineP, "fx", n*m)
239 !! ELSE
240  IF (SIZE(fx, 1) < n) cpabort("SIZE(fx,1) too small")
241  IF (SIZE(fx, 2) < m) cpabort("SIZE(fx,2) too small")
242 !! END IF
243 
244 !$OMP PARALLEL DO PRIVATE(ip,x,rhoab) DEFAULT(NONE) SHARED(fx,m,eps_rho,n)
245  DO ip = 1, n
246  rhoab = rhoa(ip) + rhob(ip)
247  IF (rhoab < eps_rho) THEN
248  fx(ip, 1:m) = 0.0_dp
249  ELSE
250  x = (rhoa(ip) - rhob(ip))/rhoab
251  IF (x < -1.0_dp) THEN
252  IF (m >= 0) fx(ip, 1) = 1.0_dp
253  IF (m >= 1) fx(ip, 2) = -f43*fxfac*2.0_dp**f13
254  IF (m >= 2) fx(ip, 3) = f13*f43*fxfac/2.0_dp**f23
255  IF (m >= 3) fx(ip, 4) = f23*f13*f43*fxfac/2.0_dp**f53
256  ELSE IF (x > 1.0_dp) THEN
257  IF (m >= 0) fx(ip, 1) = 1.0_dp
258  IF (m >= 1) fx(ip, 2) = f43*fxfac*2.0_dp**f13
259  IF (m >= 2) fx(ip, 3) = f13*f43*fxfac/2.0_dp**f23
260  IF (m >= 3) fx(ip, 4) = -f23*f13*f43*fxfac/2.0_dp**f53
261  ELSE
262  IF (m >= 0) &
263  fx(ip, 1) = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
264  IF (m >= 1) &
265  fx(ip, 2) = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
266  IF (m >= 2) &
267  fx(ip, 3) = ((1.0_dp + x)**(-f23) + (1.0_dp - x)**(-f23))* &
268  fxfac*f43*f13
269  IF (m >= 3) &
270  fx(ip, 4) = ((1.0_dp + x)**(-f53) - (1.0_dp - x)**(-f53))* &
271  fxfac*f43*f13*(-f23)
272  END IF
273  END IF
274  END DO
275 
276  END SUBROUTINE calc_fx_array
277 
278 ! **************************************************************************************************
279 !> \brief ...
280 !> \param rhoa ...
281 !> \param rhob ...
282 !> \param fx ...
283 !> \param m ...
284 ! **************************************************************************************************
285  PURE SUBROUTINE calc_fx_single(rhoa, rhob, fx, m)
286 
287 ! spin interpolation function and derivatives
288 !
289 ! f(x) = ( (1+x)^(4/3) + (1-x)^(4/3) - 2 ) / (2^(4/3)-2)
290 ! df(x) = (4/3)( (1+x)^(1/3) - (1-x)^(1/3) ) / (2^(4/3)-2)
291 ! d2f(x) = (4/9)( (1+x)^(-2/3) + (1-x)^(-2/3) ) / (2^(4/3)-2)
292 ! d3f(x) = (-8/27)( (1+x)^(-5/3) - (1-x)^(-5/3) ) / (2^(4/3)-2)
293 !
294 
295  REAL(kind=dp), INTENT(IN) :: rhoa, rhob
296  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: fx
297  INTEGER, INTENT(IN) :: m
298 
299  REAL(kind=dp) :: rhoab, x
300 
301  rhoab = rhoa + rhob
302  IF (rhoab < eps_rho) THEN
303  fx(1:m) = 0.0_dp
304  ELSE
305  x = (rhoa - rhob)/rhoab
306  IF (x < -1.0_dp) THEN
307  IF (m >= 0) fx(1) = 1.0_dp
308  IF (m >= 1) fx(2) = -f43*fxfac*2.0_dp**f13
309  IF (m >= 2) fx(3) = f13*f43*fxfac/2.0_dp**f23
310  IF (m >= 3) fx(4) = f23*f13*f43*fxfac/2.0_dp**f53
311  ELSE IF (x > 1.0_dp) THEN
312  IF (m >= 0) fx(1) = 1.0_dp
313  IF (m >= 1) fx(2) = f43*fxfac*2.0_dp**f13
314  IF (m >= 2) fx(3) = f13*f43*fxfac/2.0_dp**f23
315  IF (m >= 3) fx(4) = -f23*f13*f43*fxfac/2.0_dp**f53
316  ELSE
317  IF (m >= 0) &
318  fx(1) = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
319  IF (m >= 1) &
320  fx(2) = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
321  IF (m >= 2) &
322  fx(3) = ((1.0_dp + x)**(-f23) + (1.0_dp - x)**(-f23))* &
323  fxfac*f43*f13
324  IF (m >= 3) &
325  fx(4) = ((1.0_dp + x)**(-f53) - (1.0_dp - x)**(-f53))* &
326  fxfac*f43*f13*(-f23)
327  END IF
328  END IF
329 
330  END SUBROUTINE calc_fx_single
331 
332 ! **************************************************************************************************
333 !> \brief ...
334 !> \param a ...
335 !> \param b ...
336 !> \param z ...
337 !> \param order ...
338 ! **************************************************************************************************
339  PURE SUBROUTINE calc_z(a, b, z, order)
340 
341  REAL(kind=dp), INTENT(IN) :: a, b
342  REAL(kind=dp), DIMENSION(0:, 0:), INTENT(OUT) :: z
343  INTEGER, INTENT(IN) :: order
344 
345  REAL(kind=dp) :: c, d
346 
347  c = a + b
348 
349  z(0, 0) = (a - b)/c
350  IF (order >= 1) THEN
351  d = c*c
352  z(1, 0) = 2.0_dp*b/d
353  z(0, 1) = -2.0_dp*a/d
354  END IF
355  IF (order >= 2) THEN
356  d = d*c
357  z(2, 0) = -4.0_dp*b/d
358  z(1, 1) = 2.0_dp*(a - b)/d
359  z(0, 2) = 4.0_dp*a/d
360  END IF
361  IF (order >= 3) THEN
362  d = d*c
363  z(3, 0) = 12.0_dp*b/d
364  z(2, 1) = -4.0_dp*(a - 2.0_dp*b)/d
365  z(1, 2) = -4.0_dp*(2.0_dp*a - b)/d
366  z(0, 3) = -12.0_dp*a/d
367  END IF
368 
369  END SUBROUTINE calc_z
370 
371 END MODULE xc_functionals_utilities
372 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
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::calc_wave_vector
subroutine, public calc_wave_vector(tag, rho, grho, s)
...
Definition: xc_functionals_utilities.F:171
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition: xc_functionals_utilities.F:54
xc_functionals_utilities::calc_z
pure subroutine, public calc_z(a, b, z, order)
...
Definition: xc_functionals_utilities.F:340
xc_functionals_utilities::calc_srs_pw
subroutine, public calc_srs_pw(rho, x, n)
...
Definition: xc_functionals_utilities.F:144