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 All kind of helpful little routines
10 : !> \par History
11 : !> - Cleaned (22.04.2021, MK)
12 : !> \author CJM & JGH
13 : ! **************************************************************************************************
14 : MODULE ao_util
15 :
16 : USE kinds, ONLY: dp,&
17 : sp
18 : USE mathconstants, ONLY: dfac,&
19 : fac,&
20 : rootpi
21 : USE orbital_pointers, ONLY: coset
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ao_util'
29 :
30 : ! Public subroutines
31 :
32 : PUBLIC :: exp_radius, &
33 : exp_radius_very_extended, &
34 : gauss_exponent, &
35 : gaussint_sph, &
36 : trace_r_AxB
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief The exponent of a primitive Gaussian function for a given radius
42 : !> and threshold is calculated.
43 : !> \param l ...
44 : !> \param radius ...
45 : !> \param threshold ...
46 : !> \param prefactor ...
47 : !> \return ...
48 : !> \date 07.03.1999
49 : !> \par Variables
50 : !> - exponent : Exponent of the primitive Gaussian function.
51 : !> - l : Angular momentum quantum number l.
52 : !> - prefactor: Prefactor of the Gaussian function (e.g. a contraction
53 : !> coefficient).
54 : !> - radius : Calculated radius of the Gaussian function.
55 : !> - threshold: Threshold for radius.
56 : !> \author MK
57 : !> \version 1.0
58 : ! **************************************************************************************************
59 8348 : FUNCTION gauss_exponent(l, radius, threshold, prefactor) RESULT(exponent)
60 : INTEGER, INTENT(IN) :: l
61 : REAL(KIND=dp), INTENT(IN) :: radius, threshold, prefactor
62 : REAL(KIND=dp) :: exponent
63 :
64 8348 : exponent = 0.0_dp
65 :
66 8348 : IF (radius < 1.0E-6_dp) RETURN
67 8348 : IF (threshold < 1.0E-12_dp) RETURN
68 :
69 8348 : exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2
70 :
71 8348 : END FUNCTION gauss_exponent
72 :
73 : ! **************************************************************************************************
74 : !> \brief The radius of a primitive Gaussian function for a given threshold
75 : !> is calculated.
76 : !> g(r) = prefactor*r**l*exp(-alpha*r**2) - threshold = 0
77 : !> \param l Angular momentum quantum number l.
78 : !> \param alpha Exponent of the primitive Gaussian function.
79 : !> \param threshold Threshold for function g(r).
80 : !> \param prefactor Prefactor of the Gaussian function (e.g. a contraction
81 : !> coefficient).
82 : !> \param epsabs Absolute tolerance (radius)
83 : !> \param epsrel Relative tolerance (radius)
84 : !> \param rlow Optional lower bound radius, e.g., when searching maximum radius
85 : !> \return Calculated radius of the Gaussian function
86 : !> \par Variables
87 : !> - g : function g(r)
88 : !> - prec_exp: single/double precision exponential function
89 : !> - itermax : Maximum number of iterations
90 : !> - contract: Golden Section Search (GSS): [0.38, 0.62]
91 : !> Equidistant sampling: [0.2, 0.4, 0.6, 0.8]
92 : !> Bisection (original approach): [0.5]
93 : !> Array size may trigger vectorization
94 : ! **************************************************************************************************
95 147211270 : FUNCTION exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow) RESULT(radius)
96 : INTEGER, INTENT(IN) :: l
97 : REAL(KIND=dp), INTENT(IN) :: alpha, threshold, prefactor
98 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: epsabs, epsrel, rlow
99 : REAL(KIND=dp) :: radius
100 :
101 : INTEGER, PARAMETER :: itermax = 5000, prec_exp = sp
102 : REAL(KIND=dp), PARAMETER :: contract(*) = [0.38, 0.62], &
103 : mineps = 1.0E-12, next = 2.0, &
104 : step = 1.0
105 :
106 : INTEGER :: i, j
107 : REAL(KIND=dp) :: a, d, dr, eps, r, rd, t
108 : REAL(KIND=dp), DIMENSION(SIZE(contract)) :: g, s
109 :
110 147211270 : IF (l < 0) THEN
111 0 : CPABORT("The angular momentum quantum number is negative")
112 : END IF
113 :
114 147211270 : IF (alpha == 0.0_dp) THEN
115 0 : CPABORT("The Gaussian function exponent is zero")
116 : ELSE
117 147211270 : a = ABS(alpha)
118 : END IF
119 :
120 147211270 : IF (threshold /= 0.0_dp) THEN
121 147211270 : t = ABS(threshold)
122 : ELSE
123 0 : CPABORT("The requested threshold is zero")
124 : END IF
125 :
126 147211270 : radius = 0.0_dp
127 147211270 : IF (PRESENT(rlow)) radius = rlow
128 147211270 : IF (prefactor == 0.0_dp) RETURN
129 :
130 : ! MAX: facilitate early exit
131 146949582 : r = MAX(SQRT(0.5_dp*REAL(l, dp)/a), radius)
132 :
133 146949582 : d = ABS(prefactor); g(1) = d
134 146949582 : IF (l /= 0) THEN
135 97009619 : g(1) = g(1)*EXP(REAL(-a*r*r, KIND=prec_exp))*r**l
136 : END IF
137 : ! original approach may return radius=0
138 : ! with g(r) != g(radius)
139 : !radius = r
140 146949582 : IF (g(1) < t) RETURN ! early exit
141 :
142 77799520 : radius = r*next + step
143 159163635 : DO i = 1, itermax
144 159163635 : g(1) = d*EXP(REAL(-a*radius*radius, KIND=prec_exp))*radius**l
145 159163635 : IF (g(1) < t) EXIT
146 159163635 : r = radius; radius = r*next + step
147 : END DO
148 :
149 : ! consider absolute and relative accuracy (interval width)
150 77799520 : IF (PRESENT(epsabs)) THEN
151 69781890 : eps = epsabs
152 8017630 : ELSE IF (.NOT. PRESENT(epsrel)) THEN
153 : eps = mineps
154 : ELSE
155 : eps = HUGE(eps)
156 : END IF
157 69781890 : IF (PRESENT(epsrel)) eps = MIN(eps, epsrel*r)
158 :
159 77799520 : dr = 0.0_dp
160 1373645598 : DO i = i + 1, itermax
161 1373645598 : rd = radius - r
162 : ! check if finished or no further progress
163 1373645598 : IF ((rd < eps) .OR. (rd == dr)) RETURN
164 3887538234 : s = r + rd*contract ! interval contraction
165 3887538234 : g = d*EXP(REAL(-a*s*s, KIND=prec_exp))*s**l
166 2296708027 : DO j = 1, SIZE(contract)
167 2296708027 : IF (g(j) < t) THEN
168 919571189 : radius = s(j) ! interval [r, sj)
169 919571189 : EXIT
170 : ELSE
171 1000861949 : r = s(j) ! interval [sj, radius)
172 : END IF
173 : END DO
174 1295846078 : dr = rd
175 : END DO
176 0 : IF (i >= itermax) THEN
177 0 : CPABORT("Maximum number of iterations reached")
178 : END IF
179 :
180 : END FUNCTION exp_radius
181 :
182 : ! **************************************************************************************************
183 : !> \brief computes the radius of the Gaussian outside of which it is smaller
184 : !> than eps
185 : !> \param la_min ...
186 : !> \param la_max ...
187 : !> \param lb_min ...
188 : !> \param lb_max ...
189 : !> \param pab ...
190 : !> \param o1 ...
191 : !> \param o2 ...
192 : !> \param ra ...
193 : !> \param rb ...
194 : !> \param rp ...
195 : !> \param zetp ...
196 : !> \param eps ...
197 : !> \param prefactor ...
198 : !> \param cutoff ...
199 : !> \param epsabs ...
200 : !> \return ...
201 : !> \par History
202 : !> 03.2007 new version that assumes that the Gaussians origante from spherical
203 : !> Gaussians
204 : !> \note
205 : !> can optionally screen by the maximum element of the pab block
206 : ! **************************************************************************************************
207 37550351 : FUNCTION exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, &
208 : zetp, eps, prefactor, cutoff, epsabs) RESULT(radius)
209 :
210 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max
211 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: pab
212 : INTEGER, OPTIONAL :: o1, o2
213 : REAL(KIND=dp), INTENT(IN) :: ra(3), rb(3), rp(3), zetp, eps, &
214 : prefactor, cutoff
215 : REAL(KIND=dp), OPTIONAL :: epsabs
216 : REAL(KIND=dp) :: radius
217 :
218 : INTEGER :: i, ico, j, jco, la(3), lb(3), lxa, lxb, &
219 : lya, lyb, lza, lzb
220 : REAL(KIND=dp) :: bini, binj, coef(0:20), epsin_local, &
221 : polycoef(0:60), prefactor_local, &
222 : rad_a, rad_b, s1, s2
223 :
224 : ! get the local prefactor, we'll now use the largest density matrix element of the block to screen
225 :
226 37550351 : epsin_local = 1.0E-2_dp
227 37550351 : IF (PRESENT(epsabs)) epsin_local = epsabs
228 :
229 37550351 : IF (PRESENT(pab)) THEN
230 56897 : prefactor_local = cutoff
231 113952 : DO lxa = 0, la_max
232 171007 : DO lxb = 0, lb_max
233 171402 : DO lya = 0, la_max - lxa
234 171639 : DO lyb = 0, lb_max - lxb
235 172192 : DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya
236 172508 : DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
237 230432 : la = [lxa, lya, lza]
238 230432 : lb = [lxb, lyb, lzb]
239 57608 : ico = coset(lxa, lya, lza)
240 57608 : jco = coset(lxb, lyb, lzb)
241 115216 : prefactor_local = MAX(ABS(pab(o1 + ico, o2 + jco)), prefactor_local)
242 : END DO
243 : END DO
244 : END DO
245 : END DO
246 : END DO
247 : END DO
248 56897 : prefactor_local = prefactor*prefactor_local
249 : ELSE
250 37493454 : prefactor_local = prefactor*MAX(1.0_dp, cutoff)
251 : END IF
252 :
253 : !
254 : ! assumes that we can compute the radius for the case where
255 : ! the Gaussians a and b are both on the z - axis, but at the same
256 : ! distance as the original a and b
257 : !
258 150201404 : rad_a = SQRT(SUM((ra - rp)**2))
259 150201404 : rad_b = SQRT(SUM((rb - rp)**2))
260 :
261 114831115 : polycoef(0:la_max + lb_max) = 0.0_dp
262 95609215 : DO lxa = 0, la_max
263 186267666 : DO lxb = 0, lb_max
264 340390602 : coef(0:la_max + lb_max) = 0.0_dp
265 90658451 : bini = 1.0_dp
266 90658451 : s1 = 1.0_dp
267 222034054 : DO i = 0, lxa
268 : binj = 1.0_dp
269 : s2 = 1.0_dp
270 322056698 : DO j = 0, lxb
271 190681095 : coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
272 190681095 : binj = (binj*(lxb - j))/(j + 1)
273 322056698 : s2 = s2*(rad_b)
274 : END DO
275 131375603 : bini = (bini*(lxa - i))/(i + 1)
276 222034054 : s1 = s1*(rad_a)
277 : END DO
278 318912616 : DO i = 0, lxa + lxb
279 260853752 : polycoef(i) = MAX(polycoef(i), coef(i))
280 : END DO
281 : END DO
282 : END DO
283 :
284 114831115 : polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
285 37550351 : radius = 0.0_dp
286 114831115 : DO i = 0, la_max + lb_max
287 114831115 : radius = MAX(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
288 : END DO
289 :
290 37550351 : END FUNCTION exp_radius_very_extended
291 :
292 : ! **************************************************************************************************
293 : !> \brief ...
294 : !> \param alpha ...
295 : !> \param l ...
296 : !> \return ...
297 : ! **************************************************************************************************
298 2589870 : FUNCTION gaussint_sph(alpha, l)
299 :
300 : ! calculates the radial integral over a spherical Gaussian
301 : ! of the form
302 : ! r**(2+l) * exp(-alpha * r**2)
303 : !
304 :
305 : REAL(dp), INTENT(IN) :: alpha
306 : INTEGER, INTENT(IN) :: l
307 : REAL(dp) :: gaussint_sph
308 :
309 2589870 : IF ((l/2)*2 == l) THEN
310 : !even l:
311 : gaussint_sph = ROOTPI*0.5_dp**(l/2 + 2)*dfac(l + 1) &
312 2589870 : /SQRT(alpha)**(l + 3)
313 : ELSE
314 : !odd l:
315 0 : gaussint_sph = 0.5_dp*fac((l + 1)/2)/SQRT(alpha)**(l + 3)
316 : END IF
317 :
318 2589870 : END FUNCTION gaussint_sph
319 :
320 : ! **************************************************************************************************
321 : !> \brief ...
322 : !> \param A ...
323 : !> \param lda ...
324 : !> \param B ...
325 : !> \param ldb ...
326 : !> \param m ...
327 : !> \param n ...
328 : !> \return ...
329 : ! **************************************************************************************************
330 5106107 : PURE FUNCTION trace_r_AxB(A, lda, B, ldb, m, n)
331 :
332 : INTEGER, INTENT(in) :: lda
333 : REAL(dp), INTENT(in) :: A(lda, *)
334 : INTEGER, INTENT(in) :: ldb
335 : REAL(dp), INTENT(in) :: B(ldb, *)
336 : INTEGER, INTENT(in) :: m, n
337 : REAL(dp) :: trace_r_AxB
338 :
339 : INTEGER :: i1, i2, imod, mminus3
340 : REAL(dp) :: t1, t2, t3, t4
341 :
342 5106107 : t1 = 0._dp
343 5106107 : t2 = 0._dp
344 5106107 : t3 = 0._dp
345 5106107 : t4 = 0._dp
346 5106107 : imod = MODULO(m, 4)
347 2224760 : SELECT CASE (imod)
348 : CASE (0)
349 55479224 : DO i2 = 1, n
350 55479224 : DO i1 = 1, m, 4
351 925393152 : t1 = t1 + A(i1, i2)*B(i1, i2)
352 925393152 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
353 925393152 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
354 925393152 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
355 : END DO
356 : END DO
357 : CASE (1)
358 2392010 : mminus3 = m - 3
359 65113604 : DO i2 = 1, n
360 62721594 : DO i1 = 1, mminus3, 4
361 469858508 : t1 = t1 + A(i1, i2)*B(i1, i2)
362 469858508 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
363 469858508 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
364 469858524 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
365 : END DO
366 65113604 : t1 = t1 + A(m, i2)*B(m, i2)
367 : END DO
368 : CASE (2)
369 57276 : mminus3 = m - 3
370 2024404 : DO i2 = 1, n
371 1967128 : DO i1 = 1, mminus3, 4
372 77744144 : t1 = t1 + A(i1, i2)*B(i1, i2)
373 77744144 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
374 77744144 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
375 77744144 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
376 : END DO
377 1967128 : t1 = t1 + A(m - 1, i2)*B(m - 1, i2)
378 2024404 : t2 = t2 + A(m, i2)*B(m, i2)
379 : END DO
380 : CASE (3)
381 432061 : mminus3 = m - 3
382 11596966 : DO i2 = 1, n
383 6490859 : DO i1 = 1, mminus3, 4
384 55793131 : t1 = t1 + A(i1, i2)*B(i1, i2)
385 55793131 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
386 55793131 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
387 55818571 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
388 : END DO
389 6490859 : t1 = t1 + A(m - 2, i2)*B(m - 2, i2)
390 6490859 : t2 = t2 + A(m - 1, i2)*B(m - 1, i2)
391 6922920 : t3 = t3 + A(m, i2)*B(m, i2)
392 : END DO
393 : END SELECT
394 5106107 : trace_r_AxB = t1 + t2 + t3 + t4
395 :
396 5106107 : END FUNCTION trace_r_AxB
397 :
398 : END MODULE ao_util
399 :
|
