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        7744 :    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        7744 :       exponent = 0.0_dp
      65             : 
      66        7744 :       IF (radius < 1.0E-6_dp) RETURN
      67        7744 :       IF (threshold < 1.0E-12_dp) RETURN
      68             : 
      69        7744 :       exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2
      70             : 
      71        7744 :    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   147266682 :    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   147266682 :       IF (l .LT. 0) THEN
     111           0 :          CPABORT("The angular momentum quantum number is negative")
     112             :       END IF
     113             : 
     114   147266682 :       IF (alpha .EQ. 0.0_dp) THEN
     115           0 :          CPABORT("The Gaussian function exponent is zero")
     116             :       ELSE
     117   147266682 :          a = ABS(alpha)
     118             :       END IF
     119             : 
     120   147266682 :       IF (threshold .NE. 0.0_dp) THEN
     121   147266682 :          t = ABS(threshold)
     122             :       ELSE
     123           0 :          CPABORT("The requested threshold is zero")
     124             :       END IF
     125             : 
     126   147266682 :       radius = 0.0_dp
     127   147266682 :       IF (PRESENT(rlow)) radius = rlow
     128   147266682 :       IF (prefactor .EQ. 0.0_dp) RETURN
     129             : 
     130             :       ! MAX: facilitate early exit
     131   147017774 :       r = MAX(SQRT(0.5_dp*REAL(l, dp)/a), radius)
     132             : 
     133   147017774 :       d = ABS(prefactor); g(1) = d
     134   147017774 :       IF (l .NE. 0) THEN
     135    97423362 :          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   147017774 :       IF (g(1) .LT. t) RETURN ! early exit
     141             : 
     142    77167164 :       radius = r*next + step
     143   158105486 :       DO i = 1, itermax
     144   158105486 :          g(1) = d*EXP(REAL(-a*radius*radius, KIND=prec_exp))*radius**l
     145   158105486 :          IF (g(1) .LT. t) EXIT
     146   158105486 :          r = radius; radius = r*next + step
     147             :       END DO
     148             : 
     149             :       ! consider absolute and relative accuracy (interval width)
     150    77167164 :       IF (PRESENT(epsabs)) THEN
     151    68894772 :          eps = epsabs
     152     8272392 :       ELSE IF (.NOT. PRESENT(epsrel)) THEN
     153             :          eps = mineps
     154             :       ELSE
     155             :          eps = HUGE(eps)
     156             :       END IF
     157    68894772 :       IF (PRESENT(epsrel)) eps = MIN(eps, epsrel*r)
     158             : 
     159    77167164 :       dr = 0.0_dp
     160  1357031996 :       DO i = i + 1, itermax
     161  1357031996 :          rd = radius - r
     162             :          ! check if finished or no further progress
     163  1357031996 :          IF ((rd .LT. eps) .OR. (rd .EQ. dr)) RETURN
     164  3839594496 :          s = r + rd*contract ! interval contraction
     165  3839594496 :          g = d*EXP(REAL(-a*s*s, KIND=prec_exp))*s**l
     166  2255899043 :          DO j = 1, SIZE(contract)
     167  2255899043 :             IF (g(j) .LT. t) THEN
     168   912958894 :                radius = s(j) ! interval [r, sj)
     169   912958894 :                EXIT
     170             :             ELSE
     171   976034211 :                r = s(j) ! interval [sj, radius)
     172             :             END IF
     173             :          END DO
     174  1279864832 :          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    36909918 :    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    36909918 :       epsin_local = 1.0E-2_dp
     227    36909918 :       IF (PRESENT(epsabs)) epsin_local = epsabs
     228             : 
     229    36909918 :       IF (PRESENT(pab)) THEN
     230       56352 :          prefactor_local = cutoff
     231      112862 :          DO lxa = 0, la_max
     232      169372 :          DO lxb = 0, lb_max
     233      169767 :             DO lya = 0, la_max - lxa
     234      170004 :             DO lyb = 0, lb_max - lxb
     235      170557 :                DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya
     236      170873 :                DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
     237      228252 :                   la = (/lxa, lya, lza/)
     238      228252 :                   lb = (/lxb, lyb, lzb/)
     239       57063 :                   ico = coset(lxa, lya, lza)
     240       57063 :                   jco = coset(lxb, lyb, lzb)
     241      114126 :                   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       56352 :          prefactor_local = prefactor*prefactor_local
     249             :       ELSE
     250    36853566 :          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   147639672 :       rad_a = SQRT(SUM((ra - rp)**2))
     259   147639672 :       rad_b = SQRT(SUM((rb - rp)**2))
     260             : 
     261   113018464 :       polycoef(0:la_max + lb_max) = 0.0_dp
     262    94048686 :       DO lxa = 0, la_max
     263   183436302 :       DO lxb = 0, lb_max
     264   336295678 :          coef(0:la_max + lb_max) = 0.0_dp
     265    89387616 :          bini = 1.0_dp
     266    89387616 :          s1 = 1.0_dp
     267   219084457 :          DO i = 0, lxa
     268             :             binj = 1.0_dp
     269             :             s2 = 1.0_dp
     270   318229186 :             DO j = 0, lxb
     271   188532345 :                coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
     272   188532345 :                binj = (binj*(lxb - j))/(j + 1)
     273   318229186 :                s2 = s2*(rad_b)
     274             :             END DO
     275   129696841 :             bini = (bini*(lxa - i))/(i + 1)
     276   219084457 :             s1 = s1*(rad_a)
     277             :          END DO
     278   314674223 :          DO i = 0, lxa + lxb
     279   257535455 :             polycoef(i) = MAX(polycoef(i), coef(i))
     280             :          END DO
     281             :       END DO
     282             :       END DO
     283             : 
     284   113018464 :       polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
     285    36909918 :       radius = 0.0_dp
     286   113018464 :       DO i = 0, la_max + lb_max
     287   113018464 :          radius = MAX(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
     288             :       END DO
     289             : 
     290    36909918 :    END FUNCTION exp_radius_very_extended
     291             : 
     292             : ! **************************************************************************************************
     293             : !> \brief ...
     294             : !> \param alpha ...
     295             : !> \param l ...
     296             : !> \return ...
     297             : ! **************************************************************************************************
     298     2206788 :    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     2206788 :       IF ((l/2)*2 == l) THEN
     310             :          !even l:
     311             :          gaussint_sph = ROOTPI*0.5_dp**(l/2 + 2)*dfac(l + 1) &
     312     2206788 :                         /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     2206788 :    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     4988184 :    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     4988184 :       t1 = 0._dp
     343     4988184 :       t2 = 0._dp
     344     4988184 :       t3 = 0._dp
     345     4988184 :       t4 = 0._dp
     346     4988184 :       imod = MODULO(m, 4)
     347     2177544 :       SELECT CASE (imod)
     348             :       CASE (0)
     349    53617512 :          DO i2 = 1, n
     350    53617512 :             DO i1 = 1, m, 4
     351   890969952 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     352   890969952 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     353   890969952 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     354   890969952 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     355             :             END DO
     356             :          END DO
     357             :       CASE (1)
     358     2358628 :          mminus3 = m - 3
     359    63702520 :          DO i2 = 1, n
     360    61343892 :             DO i1 = 1, mminus3, 4
     361   454680860 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     362   454680860 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     363   454680860 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     364   454680876 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     365             :             END DO
     366    63702520 :             t1 = t1 + A(m, i2)*B(m, i2)
     367             :          END DO
     368             :       CASE (2)
     369       47204 :          mminus3 = m - 3
     370     1559020 :          DO i2 = 1, n
     371     1511816 :             DO i1 = 1, mminus3, 4
     372    33343552 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     373    33343552 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     374    33343552 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     375    33343552 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     376             :             END DO
     377     1511816 :             t1 = t1 + A(m - 1, i2)*B(m - 1, i2)
     378     1559020 :             t2 = t2 + A(m, i2)*B(m, i2)
     379             :          END DO
     380             :       CASE (3)
     381      404808 :          mminus3 = m - 3
     382    11250144 :          DO i2 = 1, n
     383     6261960 :             DO i1 = 1, mminus3, 4
     384    55247556 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     385    55247556 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     386    55247556 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     387    55261620 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     388             :             END DO
     389     6261960 :             t1 = t1 + A(m - 2, i2)*B(m - 2, i2)
     390     6261960 :             t2 = t2 + A(m - 1, i2)*B(m - 1, i2)
     391     6666768 :             t3 = t3 + A(m, i2)*B(m, i2)
     392             :          END DO
     393             :       END SELECT
     394     4988184 :       trace_r_AxB = t1 + t2 + t3 + t4
     395             : 
     396     4988184 :    END FUNCTION trace_r_AxB
     397             : 
     398             : END MODULE ao_util
     399             :