22 #include "../base/base_uses.f90"
28 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ao_util'
60 INTEGER,
INTENT(IN) :: l
61 REAL(kind=
dp),
INTENT(IN) :: radius, threshold, prefactor
62 REAL(kind=
dp) :: exponent
66 IF (radius < 1.0e-6_dp)
RETURN
67 IF (threshold < 1.0e-12_dp)
RETURN
69 exponent = log(abs(prefactor)*radius**l/threshold)/radius**2
95 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
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, &
107 REAL(kind=
dp) :: a, d, dr, eps, r, rd, t
108 REAL(kind=
dp),
DIMENSION(SIZE(contract)) :: g, s
111 cpabort(
"The angular momentum quantum number is negative")
114 IF (alpha .EQ. 0.0_dp)
THEN
115 cpabort(
"The Gaussian function exponent is zero")
120 IF (threshold .NE. 0.0_dp)
THEN
123 cpabort(
"The requested threshold is zero")
127 IF (
PRESENT(rlow)) radius = rlow
128 IF (prefactor .EQ. 0.0_dp)
RETURN
131 r = max(sqrt(0.5_dp*real(l,
dp)/a), radius)
133 d = abs(prefactor); g(1) = d
135 g(1) = g(1)*exp(real(-a*r*r, kind=prec_exp))*r**l
140 IF (g(1) .LT. t)
RETURN
142 radius = r*next + step
144 g(1) = d*exp(real(-a*radius*radius, kind=prec_exp))*radius**l
145 IF (g(1) .LT. t)
EXIT
146 r = radius; radius = r*next + step
150 IF (
PRESENT(epsabs))
THEN
152 ELSE IF (.NOT.
PRESENT(epsrel))
THEN
157 IF (
PRESENT(epsrel)) eps = min(eps, epsrel*r)
160 DO i = i + 1, itermax
163 IF ((rd .LT. eps) .OR. (rd .EQ. dr))
RETURN
165 g = d*exp(real(-a*s*s, kind=prec_exp))*s**l
166 DO j = 1,
SIZE(contract)
167 IF (g(j) .LT. t)
THEN
176 IF (i .GE. itermax)
THEN
177 cpabort(
"Maximum number of iterations reached")
207 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)
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, &
215 REAL(kind=
dp),
OPTIONAL :: epsabs
216 REAL(kind=
dp) :: radius
218 INTEGER :: i, ico, j, jco, la(3), lb(3), lxa, lxb, &
220 REAL(kind=
dp) :: bini, binj, coef(0:20), epsin_local, &
221 polycoef(0:60), prefactor_local, &
226 epsin_local = 1.0e-2_dp
227 IF (
PRESENT(epsabs)) epsin_local = epsabs
229 IF (
PRESENT(pab))
THEN
230 prefactor_local = cutoff
233 DO lya = 0, la_max - lxa
234 DO lyb = 0, lb_max - lxb
235 DO lza = max(la_min - lxa - lya, 0), la_max - lxa - lya
236 DO lzb = max(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
237 la = (/lxa, lya, lza/)
238 lb = (/lxb, lyb, lzb/)
239 ico =
coset(lxa, lya, lza)
240 jco =
coset(lxb, lyb, lzb)
241 prefactor_local = max(abs(pab(o1 + ico, o2 + jco)), prefactor_local)
248 prefactor_local = prefactor*prefactor_local
250 prefactor_local = prefactor*max(1.0_dp, cutoff)
258 rad_a = sqrt(sum((ra - rp)**2))
259 rad_b = sqrt(sum((rb - rp)**2))
261 polycoef(0:la_max + lb_max) = 0.0_dp
264 coef(0:la_max + lb_max) = 0.0_dp
271 coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
272 binj = (binj*(lxb - j))/(j + 1)
275 bini = (bini*(lxa - i))/(i + 1)
279 polycoef(i) = max(polycoef(i), coef(i))
284 polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
286 DO i = 0, la_max + lb_max
287 radius = max(radius,
exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
305 REAL(
dp),
INTENT(IN) :: alpha
306 INTEGER,
INTENT(IN) :: l
309 IF ((l/2)*2 == l)
THEN
312 /sqrt(alpha)**(l + 3)
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
339 INTEGER :: i1, i2, imod, mminus3
340 REAL(
dp) :: t1, t2, t3, t4
351 t1 = t1 + a(i1, i2)*b(i1, i2)
352 t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
353 t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
354 t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
360 DO i1 = 1, mminus3, 4
361 t1 = t1 + a(i1, i2)*b(i1, i2)
362 t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
363 t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
364 t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
366 t1 = t1 + a(m, i2)*b(m, i2)
371 DO i1 = 1, mminus3, 4
372 t1 = t1 + a(i1, i2)*b(i1, i2)
373 t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
374 t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
375 t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
377 t1 = t1 + a(m - 1, i2)*b(m - 1, i2)
378 t2 = t2 + a(m, i2)*b(m, i2)
383 DO i1 = 1, mminus3, 4
384 t1 = t1 + a(i1, i2)*b(i1, i2)
385 t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
386 t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
387 t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
389 t1 = t1 + a(m - 2, i2)*b(m - 2, i2)
390 t2 = t2 + a(m - 1, i2)*b(m - 1, i2)
391 t3 = t3 + a(m, i2)*b(m, i2)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
All kind of helpful little routines.
real(dp) function, public gaussint_sph(alpha, l)
...
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
pure real(dp) function, public trace_r_axb(A, lda, B, ldb, m, n)
...
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
real(kind=dp) function, public gauss_exponent(l, radius, threshold, prefactor)
The exponent of a primitive Gaussian function for a given radius and threshold is calculated.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public sp
Definition of mathematical constants and functions.
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:, :, :), allocatable, public coset