Line data Source code
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 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 7704 : 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 7704 : exponent = 0.0_dp
65 :
66 7704 : IF (radius < 1.0E-6_dp) RETURN
67 7704 : IF (threshold < 1.0E-12_dp) RETURN
68 :
69 7704 : exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2
70 :
71 7704 : 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 140345794 : 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 140345794 : IF (l .LT. 0) THEN
111 0 : CPABORT("The angular momentum quantum number is negative")
112 : END IF
113 :
114 140345794 : IF (alpha .EQ. 0.0_dp) THEN
115 0 : CPABORT("The Gaussian function exponent is zero")
116 : ELSE
117 140345794 : a = ABS(alpha)
118 : END IF
119 :
120 140345794 : IF (threshold .NE. 0.0_dp) THEN
121 140345794 : t = ABS(threshold)
122 : ELSE
123 0 : CPABORT("The requested threshold is zero")
124 : END IF
125 :
126 140345794 : radius = 0.0_dp
127 140345794 : IF (PRESENT(rlow)) radius = rlow
128 140345794 : IF (prefactor .EQ. 0.0_dp) RETURN
129 :
130 : ! MAX: facilitate early exit
131 140038947 : r = MAX(SQRT(0.5_dp*REAL(l, dp)/a), radius)
132 :
133 140038947 : d = ABS(prefactor); g(1) = d
134 140038947 : IF (l .NE. 0) THEN
135 91680162 : 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 140038947 : IF (g(1) .LT. t) RETURN ! early exit
141 :
142 75599422 : radius = r*next + step
143 156291713 : DO i = 1, itermax
144 156291713 : g(1) = d*EXP(REAL(-a*radius*radius, KIND=prec_exp))*radius**l
145 156291713 : IF (g(1) .LT. t) EXIT
146 156291713 : r = radius; radius = r*next + step
147 : END DO
148 :
149 : ! consider absolute and relative accuracy (interval width)
150 75599422 : IF (PRESENT(epsabs)) THEN
151 68211457 : eps = epsabs
152 7387965 : ELSE IF (.NOT. PRESENT(epsrel)) THEN
153 : eps = mineps
154 : ELSE
155 0 : eps = HUGE(eps)
156 : END IF
157 75599422 : IF (PRESENT(epsrel)) eps = MIN(eps, epsrel*r)
158 :
159 75599422 : dr = 0.0_dp
160 1325939211 : DO i = i + 1, itermax
161 1325939211 : rd = radius - r
162 : ! check if finished or no further progress
163 1325939211 : IF ((rd .LT. eps) .OR. (rd .EQ. dr)) RETURN
164 3751019367 : s = r + rd*contract ! interval contraction
165 3751019367 : g = d*EXP(REAL(-a*s*s, KIND=prec_exp))*s**l
166 2232247307 : DO j = 1, SIZE(contract)
167 2232247307 : IF (g(j) .LT. t) THEN
168 881271290 : radius = s(j) ! interval [r, sj)
169 881271290 : EXIT
170 : ELSE
171 981907518 : r = s(j) ! interval [sj, radius)
172 : END IF
173 : END DO
174 1250339789 : dr = rd
175 : END DO
176 0 : IF (i .GE. 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 36680248 : 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 36680248 : epsin_local = 1.0E-2_dp
227 36680248 : IF (PRESENT(epsabs)) epsin_local = epsabs
228 :
229 36680248 : IF (PRESENT(pab)) THEN
230 54629 : prefactor_local = cutoff
231 109416 : DO lxa = 0, la_max
232 164203 : DO lxb = 0, lb_max
233 164598 : DO lya = 0, la_max - lxa
234 164835 : DO lyb = 0, lb_max - lxb
235 165388 : DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya
236 165704 : DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
237 221360 : la = (/lxa, lya, lza/)
238 221360 : lb = (/lxb, lyb, lzb/)
239 55340 : ico = coset(lxa, lya, lza)
240 55340 : jco = coset(lxb, lyb, lzb)
241 110680 : 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 54629 : prefactor_local = prefactor*prefactor_local
249 : ELSE
250 36625619 : 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 146720992 : rad_a = SQRT(SUM((ra - rp)**2))
259 146720992 : rad_b = SQRT(SUM((rb - rp)**2))
260 :
261 111834332 : polycoef(0:la_max + lb_max) = 0.0_dp
262 93204368 : DO lxa = 0, la_max
263 181128286 : DO lxb = 0, lb_max
264 328593272 : coef(0:la_max + lb_max) = 0.0_dp
265 87923918 : bini = 1.0_dp
266 87923918 : s1 = 1.0_dp
267 214931611 : DO i = 0, lxa
268 : binj = 1.0_dp
269 : s2 = 1.0_dp
270 310691423 : DO j = 0, lxb
271 183683730 : coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
272 183683730 : binj = (binj*(lxb - j))/(j + 1)
273 310691423 : s2 = s2*(rad_b)
274 : END DO
275 127007693 : bini = (bini*(lxa - i))/(i + 1)
276 214931611 : s1 = s1*(rad_a)
277 : END DO
278 308744674 : DO i = 0, lxa + lxb
279 252220554 : polycoef(i) = MAX(polycoef(i), coef(i))
280 : END DO
281 : END DO
282 : END DO
283 :
284 111834332 : polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
285 36680248 : radius = 0.0_dp
286 111834332 : DO i = 0, la_max + lb_max
287 111834332 : radius = MAX(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
288 : END DO
289 :
290 36680248 : END FUNCTION exp_radius_very_extended
291 :
292 : ! **************************************************************************************************
293 : !> \brief ...
294 : !> \param alpha ...
295 : !> \param l ...
296 : !> \return ...
297 : ! **************************************************************************************************
298 2202716 : 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 2202716 : IF ((l/2)*2 == l) THEN
310 : !even l:
311 : gaussint_sph = ROOTPI*0.5_dp**(l/2 + 2)*dfac(l + 1) &
312 2202716 : /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 2202716 : 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 4967716 : 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 4967716 : t1 = 0._dp
343 4967716 : t2 = 0._dp
344 4967716 : t3 = 0._dp
345 4967716 : t4 = 0._dp
346 4967716 : imod = MODULO(m, 4)
347 2147744 : SELECT CASE (imod)
348 : CASE (0)
349 52807104 : DO i2 = 1, n
350 52807104 : DO i1 = 1, m, 4
351 887379296 : t1 = t1 + A(i1, i2)*B(i1, i2)
352 887379296 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
353 887379296 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
354 887379296 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
355 : END DO
356 : END DO
357 : CASE (1)
358 2371400 : mminus3 = m - 3
359 64307696 : DO i2 = 1, n
360 61936296 : DO i1 = 1, mminus3, 4
361 461328440 : t1 = t1 + A(i1, i2)*B(i1, i2)
362 461328440 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
363 461328440 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
364 461328456 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
365 : END DO
366 64307696 : t1 = t1 + A(m, i2)*B(m, i2)
367 : END DO
368 : CASE (2)
369 49852 : mminus3 = m - 3
370 1891636 : DO i2 = 1, n
371 1841784 : DO i1 = 1, mminus3, 4
372 67764000 : t1 = t1 + A(i1, i2)*B(i1, i2)
373 67764000 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
374 67764000 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
375 67764000 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
376 : END DO
377 1841784 : t1 = t1 + A(m - 1, i2)*B(m - 1, i2)
378 1891636 : t2 = t2 + A(m, i2)*B(m, i2)
379 : END DO
380 : CASE (3)
381 398720 : mminus3 = m - 3
382 11005364 : DO i2 = 1, n
383 6037648 : DO i1 = 1, mminus3, 4
384 52876732 : t1 = t1 + A(i1, i2)*B(i1, i2)
385 52876732 : t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
386 52876732 : t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
387 52890268 : t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
388 : END DO
389 6037648 : t1 = t1 + A(m - 2, i2)*B(m - 2, i2)
390 6037648 : t2 = t2 + A(m - 1, i2)*B(m - 1, i2)
391 6436368 : t3 = t3 + A(m, i2)*B(m, i2)
392 : END DO
393 : END SELECT
394 4967716 : trace_r_AxB = t1 + t2 + t3 + t4
395 :
396 4967716 : END FUNCTION trace_r_AxB
397 :
398 : END MODULE ao_util
399 :
|
