Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 86272 : SUBROUTINE set_util(cutoff)
54 :
55 : REAL(KIND=dp) :: cutoff
56 :
57 86272 : eps_rho = cutoff
58 :
59 86272 : END SUBROUTINE set_util
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param rho ...
64 : !> \param rs ...
65 : ! **************************************************************************************************
66 522427096 : 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 522427096 : IF (rho < eps_rho) THEN
74 32226076 : rs = 0.0_dp
75 : ELSE
76 490201020 : rs = rsfac*rho**(-f13)
77 : END IF
78 :
79 522427096 : END SUBROUTINE calc_rs_single
80 :
81 : ! **************************************************************************************************
82 : !> \brief ...
83 : !> \param rho ...
84 : !> \param rs ...
85 : ! **************************************************************************************************
86 0 : 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 0 : IF (SIZE(rs) < SIZE(rho)) THEN
96 0 : CPABORT("Size of array rs too small.")
97 : END IF
98 :
99 0 : !$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 0 : END SUBROUTINE calc_rs_array
109 :
110 : ! **************************************************************************************************
111 : !> \brief ...
112 : !> \param rho ...
113 : !> \param rs ...
114 : !> \param n ...
115 : ! **************************************************************************************************
116 68320 : 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 68320 : !$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 68320 : END SUBROUTINE calc_rs_pw
136 :
137 : ! **************************************************************************************************
138 : !> \brief ...
139 : !> \param rho ...
140 : !> \param x ...
141 : !> \param n ...
142 : ! **************************************************************************************************
143 827 : 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 827 : CALL calc_rs_pw(rho, x, n)
155 :
156 827 : !$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 827 : END SUBROUTINE calc_srs_pw
162 :
163 : ! **************************************************************************************************
164 : !> \brief ...
165 : !> \param tag ...
166 : !> \param rho ...
167 : !> \param grho ...
168 : !> \param s ...
169 : ! **************************************************************************************************
170 1656 : 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 1656 : fac = 1.0_dp/(2.0_dp*(3.0_dp*pi*pi)**f13)
185 1656 : IF (tag(1:1) == "u" .OR. tag(1:1) == "U") fac = fac*(2.0_dp)**f13
186 1656 : IF (tag(1:1) == "r" .OR. tag(1:1) == "R") fac = fac*(2.0_dp)**f13
187 :
188 1656 : 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 1656 : !$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 1656 : 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 0 : 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 0 : 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 0 : IF (SIZE(fx, 1) < n) CPABORT("SIZE(fx,1) too small")
241 0 : IF (SIZE(fx, 2) < m) CPABORT("SIZE(fx,2) too small")
242 : !! END IF
243 :
244 0 : !$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 0 : END SUBROUTINE calc_fx_array
277 :
278 : ! **************************************************************************************************
279 : !> \brief ...
280 : !> \param rhoa ...
281 : !> \param rhob ...
282 : !> \param fx ...
283 : !> \param m ...
284 : ! **************************************************************************************************
285 511560443 : 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 511560443 : rhoab = rhoa + rhob
302 511560443 : IF (rhoab < eps_rho) THEN
303 61950487 : fx(1:m) = 0.0_dp
304 : ELSE
305 479334367 : x = (rhoa - rhob)/rhoab
306 479334367 : IF (x < -1.0_dp) THEN
307 84983 : IF (m >= 0) fx(1) = 1.0_dp
308 84983 : IF (m >= 1) fx(2) = -f43*fxfac*2.0_dp**f13
309 84983 : IF (m >= 2) fx(3) = f13*f43*fxfac/2.0_dp**f23
310 84983 : IF (m >= 3) fx(4) = f23*f13*f43*fxfac/2.0_dp**f53
311 479249384 : ELSE IF (x > 1.0_dp) THEN
312 203936 : IF (m >= 0) fx(1) = 1.0_dp
313 203936 : IF (m >= 1) fx(2) = f43*fxfac*2.0_dp**f13
314 203936 : IF (m >= 2) fx(3) = f13*f43*fxfac/2.0_dp**f23
315 203936 : IF (m >= 3) fx(4) = -f23*f13*f43*fxfac/2.0_dp**f53
316 : ELSE
317 479045448 : IF (m >= 0) &
318 479045448 : fx(1) = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
319 479045448 : IF (m >= 1) &
320 430076310 : fx(2) = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
321 479045448 : IF (m >= 2) &
322 : fx(3) = ((1.0_dp + x)**(-f23) + (1.0_dp - x)**(-f23))* &
323 30381237 : fxfac*f43*f13
324 479045448 : IF (m >= 3) &
325 : fx(4) = ((1.0_dp + x)**(-f53) - (1.0_dp - x)**(-f53))* &
326 0 : fxfac*f43*f13*(-f23)
327 : END IF
328 : END IF
329 :
330 511560443 : END SUBROUTINE calc_fx_single
331 :
332 : ! **************************************************************************************************
333 : !> \brief ...
334 : !> \param a ...
335 : !> \param b ...
336 : !> \param z ...
337 : !> \param order ...
338 : ! **************************************************************************************************
339 1887206 : 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 1887206 : c = a + b
348 :
349 1887206 : z(0, 0) = (a - b)/c
350 1887206 : IF (order >= 1) THEN
351 1712550 : d = c*c
352 1712550 : z(1, 0) = 2.0_dp*b/d
353 1712550 : z(0, 1) = -2.0_dp*a/d
354 : END IF
355 1887206 : IF (order >= 2) THEN
356 626644 : d = d*c
357 626644 : z(2, 0) = -4.0_dp*b/d
358 626644 : z(1, 1) = 2.0_dp*(a - b)/d
359 626644 : z(0, 2) = 4.0_dp*a/d
360 : END IF
361 1887206 : IF (order >= 3) THEN
362 0 : d = d*c
363 0 : z(3, 0) = 12.0_dp*b/d
364 0 : z(2, 1) = -4.0_dp*(a - 2.0_dp*b)/d
365 0 : z(1, 2) = -4.0_dp*(2.0_dp*a - b)/d
366 0 : z(0, 3) = -12.0_dp*a/d
367 : END IF
368 :
369 1887206 : END SUBROUTINE calc_z
370 :
371 : END MODULE xc_functionals_utilities
372 :
|
