ao_util.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief All kind of helpful little routines
!> \par History
!>      - Cleaned (22.04.2021, MK)
!> \author CJM & JGH
! **************************************************************************************************
MODULE ao_util
 
   USE kinds,                           ONLY: dp,&
                                              sp
   USE mathconstants,                   ONLY: dfac,&
                                              fac,&
                                              rootpi
   USE orbital_pointers,                ONLY: coset
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ao_util'
 
   ! Public subroutines
 
   PUBLIC :: exp_radius, &
             exp_radius_very_extended, &
             gauss_exponent, &
             gaussint_sph, &
             trace_r_axb
 
CONTAINS
 
! **************************************************************************************************
!> \brief  The exponent of a primitive Gaussian function for a given radius
!>         and threshold is calculated.
!> \param l ...
!> \param radius ...
!> \param threshold ...
!> \param prefactor ...
!> \return ...
!> \date   07.03.1999
!> \par Variables
!>      - exponent : Exponent of the primitive Gaussian function.
!>      - l        : Angular momentum quantum number l.
!>      - prefactor: Prefactor of the Gaussian function (e.g. a contraction
!>                   coefficient).
!>      - radius   : Calculated radius of the Gaussian function.
!>      - threshold: Threshold for radius.
!> \author MK
!> \version 1.0
! **************************************************************************************************
   FUNCTION gauss_exponent(l, radius, threshold, prefactor) RESULT(exponent)
      INTEGER, INTENT(IN)                                :: l
      REAL(kind=dp), INTENT(IN)                          :: radius, threshold, prefactor
      REAL(kind=dp)                                      :: exponent
 
      exponent = 0.0_dp
 
      IF (radius < 1.0e-6_dp) RETURN
      IF (threshold < 1.0e-12_dp) RETURN
 
      exponent = log(abs(prefactor)*radius**l/threshold)/radius**2
 
   END FUNCTION gauss_exponent
 
! **************************************************************************************************
!> \brief  The radius of a primitive Gaussian function for a given threshold
!>         is calculated.
!>               g(r) = prefactor*r**l*exp(-alpha*r**2) - threshold = 0
!> \param l Angular momentum quantum number l.
!> \param alpha Exponent of the primitive Gaussian function.
!> \param threshold Threshold for function g(r).
!> \param prefactor Prefactor of the Gaussian function (e.g. a contraction
!>                     coefficient).
!> \param epsabs Absolute tolerance (radius)
!> \param epsrel Relative tolerance (radius)
!> \param rlow Optional lower bound radius, e.g., when searching maximum radius
!> \return Calculated radius of the Gaussian function
!> \par Variables
!>        - g       : function g(r)
!>        - prec_exp: single/double precision exponential function
!>        - itermax : Maximum number of iterations
!>        - contract: Golden Section Search (GSS): [0.38, 0.62]
!>                    Equidistant sampling: [0.2, 0.4, 0.6, 0.8]
!>                    Bisection (original approach): [0.5]
!>                    Array size may trigger vectorization
! **************************************************************************************************
   FUNCTION exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow) RESULT(radius)
      INTEGER, INTENT(IN)                                :: l
      REAL(kind=dp), INTENT(IN)                          :: alpha, threshold, prefactor
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: epsabs, epsrel, rlow
      REAL(kind=dp)                                      :: radius
 
      INTEGER, PARAMETER                                 :: itermax = 5000, prec_exp = sp
      REAL(kind=dp), PARAMETER :: contract(*) = [0.38_dp, 0.62_dp], mineps = 1.0e-12_dp, &
         next = 2.0_dp, step = 1.0_dp
 
      INTEGER                                            :: i, j
      REAL(kind=dp)                                      :: a, d, dr, eps, r, rd, t
      REAL(kind=dp), DIMENSION(SIZE(contract))           :: g, s
 
      IF (l < 0) THEN
         cpabort("The angular momentum quantum number is negative")
      END IF
 
      IF (alpha == 0.0_dp) THEN
         cpabort("The Gaussian function exponent is zero")
      ELSE
         a = abs(alpha)
      END IF
 
      IF (threshold /= 0.0_dp) THEN
         t = abs(threshold)
      ELSE
         cpabort("The requested threshold is zero")
      END IF
 
      radius = 0.0_dp
      IF (PRESENT(rlow)) radius = rlow
      IF (prefactor == 0.0_dp) RETURN
 
      ! MAX: facilitate early exit
      r = max(sqrt(0.5_dp*real(l, dp)/a), radius)
 
      d = abs(prefactor); g(1) = d
      IF (l /= 0) THEN
         g(1) = g(1)*exp(real(-a*r*r, kind=prec_exp))*r**l
      END IF
      ! original approach may return radius=0
      ! with g(r) != g(radius)
      !radius = r
      IF (g(1) < t) RETURN ! early exit
 
      radius = r*next + step
      DO i = 1, itermax
         g(1) = d*exp(real(-a*radius*radius, kind=prec_exp))*radius**l
         IF (g(1) < t) EXIT
         r = radius; radius = r*next + step
      END DO
 
      ! consider absolute and relative accuracy (interval width)
      IF (PRESENT(epsabs)) THEN
         eps = epsabs
      ELSE IF (.NOT. PRESENT(epsrel)) THEN
         eps = mineps
      ELSE
         eps = huge(eps)
      END IF
      IF (PRESENT(epsrel)) eps = min(eps, epsrel*r)
 
      dr = 0.0_dp
      DO i = i + 1, itermax
         rd = radius - r
         ! check if finished or no further progress
         IF ((rd < eps) .OR. (rd == dr)) RETURN
         s = r + rd*contract ! interval contraction
         g = d*exp(real(-a*s*s, kind=prec_exp))*s**l
         DO j = 1, SIZE(contract)
            IF (g(j) < t) THEN
               radius = s(j) ! interval [r, sj)
               EXIT
            ELSE
               r = s(j) ! interval [sj, radius)
            END IF
         END DO
         dr = rd
      END DO
      IF (i >= itermax) THEN
         cpabort("Maximum number of iterations reached")
      END IF
 
   END FUNCTION exp_radius
 
! **************************************************************************************************
!> \brief computes the radius of the Gaussian outside of which it is smaller
!>      than eps
!> \param la_min ...
!> \param la_max ...
!> \param lb_min ...
!> \param lb_max ...
!> \param pab ...
!> \param o1 ...
!> \param o2 ...
!> \param ra ...
!> \param rb ...
!> \param rp ...
!> \param zetp ...
!> \param eps ...
!> \param prefactor ...
!> \param cutoff ...
!> \param epsabs ...
!> \return ...
!> \par History
!>      03.2007 new version that assumes that the Gaussians origante from spherical
!>              Gaussians
!> \note
!>      can optionally screen by the maximum element of the pab block
! **************************************************************************************************
   FUNCTION exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, &
                                     zetp, eps, prefactor, cutoff, epsabs) RESULT(radius)
 
      INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max
      REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: pab
      INTEGER, OPTIONAL                                  :: o1, o2
      REAL(kind=dp), INTENT(IN)                          :: ra(3), rb(3), rp(3), zetp, eps, &
                                                            prefactor, cutoff
      REAL(kind=dp), OPTIONAL                            :: epsabs
      REAL(kind=dp)                                      :: radius
 
      INTEGER                                            :: i, ico, j, jco, la(3), lb(3), lxa, lxb, &
                                                            lya, lyb, lza, lzb
      REAL(kind=dp)                                      :: bini, binj, coef(0:20), epsin_local, &
                                                            polycoef(0:60), prefactor_local, &
                                                            rad_a, rad_b, s1, s2
 
! get the local prefactor, we'll now use the largest density matrix element of the block to screen
 
      epsin_local = 1.0e-2_dp
      IF (PRESENT(epsabs)) epsin_local = epsabs
 
      IF (PRESENT(pab)) THEN
         prefactor_local = cutoff
         DO lxa = 0, la_max
         DO lxb = 0, lb_max
            DO lya = 0, la_max - lxa
            DO lyb = 0, lb_max - lxb
               DO lza = max(la_min - lxa - lya, 0), la_max - lxa - lya
               DO lzb = max(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
                  la = [lxa, lya, lza]
                  lb = [lxb, lyb, lzb]
                  ico = coset(lxa, lya, lza)
                  jco = coset(lxb, lyb, lzb)
                  prefactor_local = max(abs(pab(o1 + ico, o2 + jco)), prefactor_local)
               END DO
               END DO
            END DO
            END DO
         END DO
         END DO
         prefactor_local = prefactor*prefactor_local
      ELSE
         prefactor_local = prefactor*max(1.0_dp, cutoff)
      END IF
 
      !
      ! assumes that we can compute the radius for the case where
      ! the Gaussians a and b are both on the z - axis, but at the same
      ! distance as the original a and b
      !
      rad_a = sqrt(sum((ra - rp)**2))
      rad_b = sqrt(sum((rb - rp)**2))
 
      polycoef(0:la_max + lb_max) = 0.0_dp
      DO lxa = 0, la_max
      DO lxb = 0, lb_max
         coef(0:la_max + lb_max) = 0.0_dp
         bini = 1.0_dp
         s1 = 1.0_dp
         DO i = 0, lxa
            binj = 1.0_dp
            s2 = 1.0_dp
            DO j = 0, lxb
               coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
               binj = (binj*(lxb - j))/(j + 1)
               s2 = s2*(rad_b)
            END DO
            bini = (bini*(lxa - i))/(i + 1)
            s1 = s1*(rad_a)
         END DO
         DO i = 0, lxa + lxb
            polycoef(i) = max(polycoef(i), coef(i))
         END DO
      END DO
      END DO
 
      polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
      radius = 0.0_dp
      DO i = 0, la_max + lb_max
         radius = max(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
      END DO
 
   END FUNCTION exp_radius_very_extended
 
! **************************************************************************************************
!> \brief ...
!> \param alpha ...
!> \param l ...
!> \return ...
! **************************************************************************************************
   FUNCTION gaussint_sph(alpha, l)
 
      !  calculates the radial integral over a spherical Gaussian
      !  of the form
      !     r**(2+l) * exp(-alpha * r**2)
      !
 
      REAL(dp), INTENT(IN)                               :: alpha
      INTEGER, INTENT(IN)                                :: l
      REAL(dp)                                           :: gaussint_sph
 
      IF ((l/2)*2 == l) THEN
         !even l:
         gaussint_sph = rootpi*0.5_dp**(l/2 + 2)*dfac(l + 1) &
                        /sqrt(alpha)**(l + 3)
      ELSE
         !odd l:
         gaussint_sph = 0.5_dp*fac((l + 1)/2)/sqrt(alpha)**(l + 3)
      END IF
 
   END FUNCTION gaussint_sph
 
! **************************************************************************************************
!> \brief ...
!> \param A ...
!> \param lda ...
!> \param B ...
!> \param ldb ...
!> \param m ...
!> \param n ...
!> \return ...
! **************************************************************************************************
   PURE FUNCTION trace_r_axb(A, lda, B, ldb, m, n)
 
      INTEGER, INTENT(in)                                :: lda
      REAL(dp), INTENT(in)                               :: a(lda, *)
      INTEGER, INTENT(in)                                :: ldb
      REAL(dp), INTENT(in)                               :: b(ldb, *)
      INTEGER, INTENT(in)                                :: m, n
      REAL(dp)                                           :: trace_r_axb
 
      INTEGER                                            :: i1, i2, imod, mminus3
      REAL(dp)                                           :: t1, t2, t3, t4
 
      t1 = 0._dp
      t2 = 0._dp
      t3 = 0._dp
      t4 = 0._dp
      imod = modulo(m, 4)
      SELECT CASE (imod)
      CASE (0)
         DO i2 = 1, n
            DO i1 = 1, m, 4
               t1 = t1 + a(i1, i2)*b(i1, i2)
               t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
               t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
               t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
            END DO
         END DO
      CASE (1)
         mminus3 = m - 3
         DO i2 = 1, n
            DO i1 = 1, mminus3, 4
               t1 = t1 + a(i1, i2)*b(i1, i2)
               t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
               t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
               t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
            END DO
            t1 = t1 + a(m, i2)*b(m, i2)
         END DO
      CASE (2)
         mminus3 = m - 3
         DO i2 = 1, n
            DO i1 = 1, mminus3, 4
               t1 = t1 + a(i1, i2)*b(i1, i2)
               t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
               t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
               t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
            END DO
            t1 = t1 + a(m - 1, i2)*b(m - 1, i2)
            t2 = t2 + a(m, i2)*b(m, i2)
         END DO
      CASE (3)
         mminus3 = m - 3
         DO i2 = 1, n
            DO i1 = 1, mminus3, 4
               t1 = t1 + a(i1, i2)*b(i1, i2)
               t2 = t2 + a(i1 + 1, i2)*b(i1 + 1, i2)
               t3 = t3 + a(i1 + 2, i2)*b(i1 + 2, i2)
               t4 = t4 + a(i1 + 3, i2)*b(i1 + 3, i2)
            END DO
            t1 = t1 + a(m - 2, i2)*b(m - 2, i2)
            t2 = t2 + a(m - 1, i2)*b(m - 1, i2)
            t3 = t3 + a(m, i2)*b(m, i2)
         END DO
      END SELECT
      trace_r_axb = t1 + t2 + t3 + t4
 
   END FUNCTION trace_r_axb
 
END MODULE ao_util