(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
15MODULE 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
47CONTAINS
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
371END MODULE xc_functionals_utilities
372
xc_functionals_utilities::calc_fx
Definition xc_functionals_utilities.F:37
xc_functionals_utilities::calc_rs
Definition xc_functionals_utilities.F:41
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
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
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