Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 12044 : 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 12044 : exponent = 0.0_dp
65 :
66 12044 : IF (radius < 1.0E-6_dp) RETURN
67 12044 : IF (threshold < 1.0E-12_dp) RETURN
68 :
69 12044 : exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2
70 :
71 12044 : 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 171089130 : 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_dp, 0.62_dp], mineps = 1.0E-12_dp, &
103 : next = 2.0_dp, step = 1.0_dp
104 :
105 : INTEGER :: i, j
106 : REAL(KIND=dp) :: a, d, dr, eps, r, rd, t
107 : REAL(KIND=dp), DIMENSION(SIZE(contract)) :: g, s
108 :
109 171089130 : IF (l < 0) THEN
110 0 : CPABORT("The angular momentum quantum number is negative")
111 : END IF
112 :
113 171089130 : IF (alpha == 0.0_dp) THEN
114 0 : CPABORT("The Gaussian function exponent is zero")
115 : ELSE
116 171089130 : a = ABS(alpha)
117 : END IF
118 :
119 171089130 : IF (threshold /= 0.0_dp) THEN
120 171089130 : t = ABS(threshold)
121 : ELSE
122 0 : CPABORT("The requested threshold is zero")
123 : END IF
124 :
125 171089130 : radius = 0.0_dp
126 171089130 : IF (PRESENT(rlow)) radius = rlow
127 171089130 : IF (prefactor == 0.0_dp) RETURN
128 :
129 : ! MAX: facilitate early exit
130 170779004 : r = MAX(SQRT(0.5_dp*REAL(l, dp)/a), radius)
131 :
132 170779004 : d = ABS(prefactor); g(1) = d
133 170779004 : IF (l /= 0) THEN
134 112100607 : g(1) = g(1)*EXP(REAL(-a*r*r, KIND=prec_exp))*r**l
135 : END IF
136 : ! original approach may return radius=0
137 : ! with g(r) != g(radius)
138 : !radius = r
139 170779004 : IF (g(1) < t) RETURN ! early exit
140 :
141 91269618 : radius = r*next + step
142 181938911 : DO i = 1, itermax
143 181938911 : g(1) = d*EXP(REAL(-a*radius*radius, KIND=prec_exp))*radius**l
144 181938911 : IF (g(1) < t) EXIT
145 181938911 : r = radius; radius = r*next + step
146 : END DO
147 :
148 : ! consider absolute and relative accuracy (interval width)
149 91269618 : IF (PRESENT(epsabs)) THEN
150 81242442 : eps = epsabs
151 10027176 : ELSE IF (.NOT. PRESENT(epsrel)) THEN
152 : eps = mineps
153 : ELSE
154 : eps = HUGE(eps)
155 : END IF
156 81242442 : IF (PRESENT(epsrel)) eps = MIN(eps, epsrel*r)
157 :
158 91269618 : dr = 0.0_dp
159 1643579665 : DO i = i + 1, itermax
160 1643579665 : rd = radius - r
161 : ! check if finished or no further progress
162 1643579665 : IF ((rd < eps) .OR. (rd == dr)) RETURN
163 4656930141 : s = r + rd*contract ! interval contraction
164 4656930141 : g = d*EXP(REAL(-a*s*s, KIND=prec_exp))*s**l
165 2742676988 : DO j = 1, SIZE(contract)
166 2742676988 : IF (g(j) < t) THEN
167 1104902836 : radius = s(j) ! interval [r, sj)
168 1104902836 : EXIT
169 : ELSE
170 1190366941 : r = s(j) ! interval [sj, radius)
171 : END IF
172 : END DO
173 1552310047 : dr = rd
174 : END DO
175 0 : IF (i >= itermax) THEN
176 0 : CPABORT("Maximum number of iterations reached")
177 : END IF
178 :
179 : END FUNCTION exp_radius
180 :
181 : ! **************************************************************************************************
182 : !> \brief computes the radius of the Gaussian outside of which it is smaller
183 : !> than eps
184 : !> \param la_min ...
185 : !> \param la_max ...
186 : !> \param lb_min ...
187 : !> \param lb_max ...
188 : !> \param pab ...
189 : !> \param o1 ...
190 : !> \param o2 ...
191 : !> \param ra ...
192 : !> \param rb ...
193 : !> \param rp ...
194 : !> \param zetp ...
195 : !> \param eps ...
196 : !> \param prefactor ...
197 : !> \param cutoff ...
198 : !> \param epsabs ...
199 : !> \return ...
200 : !> \par History
201 : !> 03.2007 new version that assumes that the Gaussians origante from spherical
202 : !> Gaussians
203 : !> \note
204 : !> can optionally screen by the maximum element of the pab block
205 : ! **************************************************************************************************
206 44611604 : FUNCTION exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, &
207 : zetp, eps, prefactor, cutoff, epsabs) RESULT(radius)
208 :
209 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max
210 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: pab
211 : INTEGER, OPTIONAL :: o1, o2
212 : REAL(KIND=dp), INTENT(IN) :: ra(3), rb(3), rp(3), zetp, eps, &
213 : prefactor, cutoff
214 : REAL(KIND=dp), OPTIONAL :: epsabs
215 : REAL(KIND=dp) :: radius
216 :
217 : INTEGER :: i, ico, j, jco, la(3), lb(3), lxa, lxb, &
218 : lya, lyb, lza, lzb
219 : REAL(KIND=dp) :: bini, binj, coef(0:20), epsin_local, &
220 : polycoef(0:60), prefactor_local, &
221 : rad_a, rad_b, s1, s2
222 :
223 : ! get the local prefactor, we'll now use the largest density matrix element of the block to screen
224 :
225 44611604 : epsin_local = 1.0E-2_dp
226 44611604 : IF (PRESENT(epsabs)) epsin_local = epsabs
227 :
228 44611604 : IF (PRESENT(pab)) THEN
229 65389 : prefactor_local = cutoff
230 130984 : DO lxa = 0, la_max
231 196579 : DO lxb = 0, lb_max
232 197142 : DO lya = 0, la_max - lxa
233 197499 : DO lyb = 0, lb_max - lxb
234 198412 : DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya
235 198968 : DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
236 266032 : la = [lxa, lya, lza]
237 266032 : lb = [lxb, lyb, lzb]
238 66508 : ico = coset(lxa, lya, lza)
239 66508 : jco = coset(lxb, lyb, lzb)
240 133016 : prefactor_local = MAX(ABS(pab(o1 + ico, o2 + jco)), prefactor_local)
241 : END DO
242 : END DO
243 : END DO
244 : END DO
245 : END DO
246 : END DO
247 65389 : prefactor_local = prefactor*prefactor_local
248 : ELSE
249 44546215 : prefactor_local = prefactor*MAX(1.0_dp, cutoff)
250 : END IF
251 :
252 : !
253 : ! assumes that we can compute the radius for the case where
254 : ! the Gaussians a and b are both on the z - axis, but at the same
255 : ! distance as the original a and b
256 : !
257 178446416 : rad_a = SQRT(SUM((ra - rp)**2))
258 178446416 : rad_b = SQRT(SUM((rb - rp)**2))
259 :
260 137798047 : polycoef(0:la_max + lb_max) = 0.0_dp
261 114160251 : DO lxa = 0, la_max
262 224410597 : DO lxb = 0, lb_max
263 416051116 : coef(0:la_max + lb_max) = 0.0_dp
264 110250346 : bini = 1.0_dp
265 110250346 : s1 = 1.0_dp
266 270347880 : DO i = 0, lxa
267 : binj = 1.0_dp
268 : s2 = 1.0_dp
269 393200129 : DO j = 0, lxb
270 233102595 : coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
271 233102595 : binj = (binj*(lxb - j))/(j + 1)
272 393200129 : s2 = s2*(rad_b)
273 : END DO
274 160097534 : bini = (bini*(lxa - i))/(i + 1)
275 270347880 : s1 = s1*(rad_a)
276 : END DO
277 387824551 : DO i = 0, lxa + lxb
278 318275904 : polycoef(i) = MAX(polycoef(i), coef(i))
279 : END DO
280 : END DO
281 : END DO
282 :
283 137798047 : polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
284 44611604 : radius = 0.0_dp
285 137798047 : DO i = 0, la_max + lb_max
286 137798047 : radius = MAX(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
287 : END DO
288 :
289 44611604 : END FUNCTION exp_radius_very_extended
290 :
291 : ! **************************************************************************************************
292 : !> \brief ...
293 : !> \param alpha ...
294 : !> \param l ...
295 : !> \return ...
296 : ! **************************************************************************************************
297 5522084 : FUNCTION gaussint_sph(alpha, l)
298 :
299 : ! calculates the radial integral over a spherical Gaussian
300 : ! of the form
301 : ! r**(2+l) * exp(-alpha * r**2)
302 : !
303 :
304 : REAL(dp), INTENT(IN) :: alpha
305 : INTEGER, INTENT(IN) :: l
306 : REAL(dp) :: gaussint_sph
307 :
308 5522084 : IF ((l/2)*2 == l) THEN
309 : !even l:
310 : gaussint_sph = ROOTPI*0.5_dp**(l/2 + 2)*dfac(l + 1) &
311 5522084 : /SQRT(alpha)**(l + 3)
312 : ELSE
313 : !odd l:
314 0 : gaussint_sph = 0.5_dp*fac((l + 1)/2)/SQRT(alpha)**(l + 3)
315 : END IF
316 :
317 5522084 : END FUNCTION gaussint_sph
318 :
319 : ! **************************************************************************************************
320 : !> \brief ...
321 : !> \param A ...
322 : !> \param lda ...
323 : !> \param B ...
324 : !> \param ldb ...
325 : !> \param m ...
326 : !> \param n ...
327 : !> \return ...
328 : ! **************************************************************************************************
329 11143803 : PURE FUNCTION trace_r_AxB(A, lda, B, ldb, m, n)
330 :
331 : INTEGER, INTENT(in) :: lda
332 : REAL(dp), INTENT(in) :: A(lda, *)
333 : INTEGER, INTENT(in) :: ldb
334 : REAL(dp), INTENT(in) :: B(ldb, *)
335 : INTEGER, INTENT(in) :: m, n
336 : REAL(dp) :: trace_r_AxB
337 :
338 : INTEGER :: i1, i2, imod, mminus3
339 : REAL(dp) :: t1, t2, t3, t4
340 :
341 11143803 : t1 = 0._dp
342 11143803 : t2 = 0._dp
343 11143803 : t3 = 0._dp
344 11143803 : t4 = 0._dp
345 11143803 : imod = MODULO(m, 4)
346 3028944 : SELECT CASE (imod)
347 : CASE (0)
348 91763168 : DO i2 = 1, n
349 91763168 : DO i1 = 1, m, 4
350 1954746800 : t1 = t1 + A(i1, i2)*B(i1, i2)
351 1954746800 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
352 1954746800 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
353 1954746800 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
354 : END DO
355 : END DO
356 : CASE (1)
357 7461046 : mminus3 = m - 3
358 201519148 : DO i2 = 1, n
359 194058102 : DO i1 = 1, mminus3, 4
360 1460848688 : t1 = t1 + A(i1, i2)*B(i1, i2)
361 1460848688 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
362 1460848688 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
363 1460848704 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
364 : END DO
365 201519148 : t1 = t1 + A(m, i2)*B(m, i2)
366 : END DO
367 : CASE (2)
368 83152 : mminus3 = m - 3
369 2065584 : DO i2 = 1, n
370 1982432 : DO i1 = 1, mminus3, 4
371 59361184 : t1 = t1 + A(i1, i2)*B(i1, i2)
372 59361184 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
373 59361184 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
374 59361184 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
375 : END DO
376 1982432 : t1 = t1 + A(m - 1, i2)*B(m - 1, i2)
377 2065584 : t2 = t2 + A(m, i2)*B(m, i2)
378 : END DO
379 : CASE (3)
380 570661 : mminus3 = m - 3
381 17656878 : DO i2 = 1, n
382 6513075 : DO i1 = 1, mminus3, 4
383 28146887 : t1 = t1 + A(i1, i2)*B(i1, i2)
384 28146887 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
385 28146887 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
386 28172327 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
387 : END DO
388 6513075 : t1 = t1 + A(m - 2, i2)*B(m - 2, i2)
389 6513075 : t2 = t2 + A(m - 1, i2)*B(m - 1, i2)
390 7083736 : t3 = t3 + A(m, i2)*B(m, i2)
391 : END DO
392 : END SELECT
393 11143803 : trace_r_AxB = t1 + t2 + t3 + t4
394 :
395 11143803 : END FUNCTION trace_r_AxB
396 :
397 : END MODULE ao_util
398 :
|
