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        7688 :    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        7688 :       exponent = 0.0_dp
      65             : 
      66        7688 :       IF (radius < 1.0E-6_dp) RETURN
      67        7688 :       IF (threshold < 1.0E-12_dp) RETURN
      68             : 
      69        7688 :       exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2
      70             : 
      71        7688 :    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   141684779 :    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   141684779 :       IF (l .LT. 0) THEN
     111           0 :          CPABORT("The angular momentum quantum number is negative")
     112             :       END IF
     113             : 
     114   141684779 :       IF (alpha .EQ. 0.0_dp) THEN
     115           0 :          CPABORT("The Gaussian function exponent is zero")
     116             :       ELSE
     117   141684779 :          a = ABS(alpha)
     118             :       END IF
     119             : 
     120   141684779 :       IF (threshold .NE. 0.0_dp) THEN
     121   141684779 :          t = ABS(threshold)
     122             :       ELSE
     123           0 :          CPABORT("The requested threshold is zero")
     124             :       END IF
     125             : 
     126   141684779 :       radius = 0.0_dp
     127   141684779 :       IF (PRESENT(rlow)) radius = rlow
     128   141684779 :       IF (prefactor .EQ. 0.0_dp) RETURN
     129             : 
     130             :       ! MAX: facilitate early exit
     131   141376484 :       r = MAX(SQRT(0.5_dp*REAL(l, dp)/a), radius)
     132             : 
     133   141376484 :       d = ABS(prefactor); g(1) = d
     134   141376484 :       IF (l .NE. 0) THEN
     135    92658797 :          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   141376484 :       IF (g(1) .LT. t) RETURN ! early exit
     141             : 
     142    76214871 :       radius = r*next + step
     143   157381055 :       DO i = 1, itermax
     144   157381055 :          g(1) = d*EXP(REAL(-a*radius*radius, KIND=prec_exp))*radius**l
     145   157381055 :          IF (g(1) .LT. t) EXIT
     146   157381055 :          r = radius; radius = r*next + step
     147             :       END DO
     148             : 
     149             :       ! consider absolute and relative accuracy (interval width)
     150    76214871 :       IF (PRESENT(epsabs)) THEN
     151    68705790 :          eps = epsabs
     152     7509081 :       ELSE IF (.NOT. PRESENT(epsrel)) THEN
     153             :          eps = mineps
     154             :       ELSE
     155             :          eps = HUGE(eps)
     156             :       END IF
     157    68705790 :       IF (PRESENT(epsrel)) eps = MIN(eps, epsrel*r)
     158             : 
     159    76214871 :       dr = 0.0_dp
     160  1334936780 :       DO i = i + 1, itermax
     161  1334936780 :          rd = radius - r
     162             :          ! check if finished or no further progress
     163  1334936780 :          IF ((rd .LT. eps) .OR. (rd .EQ. dr)) RETURN
     164  3776165727 :          s = r + rd*contract ! interval contraction
     165  3776165727 :          g = d*EXP(REAL(-a*s*s, KIND=prec_exp))*s**l
     166  2245000223 :          DO j = 1, SIZE(contract)
     167  2245000223 :             IF (g(j) .LT. t) THEN
     168   888023328 :                radius = s(j) ! interval [r, sj)
     169   888023328 :                EXIT
     170             :             ELSE
     171   986278314 :                r = s(j) ! interval [sj, radius)
     172             :             END IF
     173             :          END DO
     174  1258721909 :          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    36914190 :    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    36914190 :       epsin_local = 1.0E-2_dp
     227    36914190 :       IF (PRESENT(epsabs)) epsin_local = epsabs
     228             : 
     229    36914190 :       IF (PRESENT(pab)) THEN
     230       55346 :          prefactor_local = cutoff
     231      110850 :          DO lxa = 0, la_max
     232      166354 :          DO lxb = 0, lb_max
     233      166749 :             DO lya = 0, la_max - lxa
     234      166986 :             DO lyb = 0, lb_max - lxb
     235      167539 :                DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya
     236      167855 :                DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb
     237      224228 :                   la = (/lxa, lya, lza/)
     238      224228 :                   lb = (/lxb, lyb, lzb/)
     239       56057 :                   ico = coset(lxa, lya, lza)
     240       56057 :                   jco = coset(lxb, lyb, lzb)
     241      112114 :                   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       55346 :          prefactor_local = prefactor*prefactor_local
     249             :       ELSE
     250    36858844 :          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   147656760 :       rad_a = SQRT(SUM((ra - rp)**2))
     259   147656760 :       rad_b = SQRT(SUM((rb - rp)**2))
     260             : 
     261   112663025 :       polycoef(0:la_max + lb_max) = 0.0_dp
     262    93875863 :       DO lxa = 0, la_max
     263   182614980 :       DO lxb = 0, lb_max
     264   332430616 :          coef(0:la_max + lb_max) = 0.0_dp
     265    88739117 :          bini = 1.0_dp
     266    88739117 :          s1 = 1.0_dp
     267   217146951 :          DO i = 0, lxa
     268             :             binj = 1.0_dp
     269             :             s2 = 1.0_dp
     270   314451000 :             DO j = 0, lxb
     271   186043166 :                coef(lxa + lxb - i - j) = coef(lxa + lxb - i - j) + bini*binj*s1*s2
     272   186043166 :                binj = (binj*(lxb - j))/(j + 1)
     273   314451000 :                s2 = s2*(rad_b)
     274             :             END DO
     275   128407834 :             bini = (bini*(lxa - i))/(i + 1)
     276   217146951 :             s1 = s1*(rad_a)
     277             :          END DO
     278   311916098 :          DO i = 0, lxa + lxb
     279   254954425 :             polycoef(i) = MAX(polycoef(i), coef(i))
     280             :          END DO
     281             :       END DO
     282             :       END DO
     283             : 
     284   112663025 :       polycoef(0:la_max + lb_max) = polycoef(0:la_max + lb_max)*prefactor_local
     285    36914190 :       radius = 0.0_dp
     286   112663025 :       DO i = 0, la_max + lb_max
     287   112663025 :          radius = MAX(radius, exp_radius(i, zetp, eps, polycoef(i), epsin_local, rlow=radius))
     288             :       END DO
     289             : 
     290    36914190 :    END FUNCTION exp_radius_very_extended
     291             : 
     292             : ! **************************************************************************************************
     293             : !> \brief ...
     294             : !> \param alpha ...
     295             : !> \param l ...
     296             : !> \return ...
     297             : ! **************************************************************************************************
     298     2187848 :    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     2187848 :       IF ((l/2)*2 == l) THEN
     310             :          !even l:
     311             :          gaussint_sph = ROOTPI*0.5_dp**(l/2 + 2)*dfac(l + 1) &
     312     2187848 :                         /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     2187848 :    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     4965884 :    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     4965884 :       t1 = 0._dp
     343     4965884 :       t2 = 0._dp
     344     4965884 :       t3 = 0._dp
     345     4965884 :       t4 = 0._dp
     346     4965884 :       imod = MODULO(m, 4)
     347     2146344 :       SELECT CASE (imod)
     348             :       CASE (0)
     349    52783144 :          DO i2 = 1, n
     350    52783144 :             DO i1 = 1, m, 4
     351   886502656 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     352   886502656 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     353   886502656 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     354   886502656 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     355             :             END DO
     356             :          END DO
     357             :       CASE (1)
     358     2371160 :          mminus3 = m - 3
     359    64292816 :          DO i2 = 1, n
     360    61921656 :             DO i1 = 1, mminus3, 4
     361   461108840 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     362   461108840 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     363   461108840 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     364   461108856 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     365             :             END DO
     366    64292816 :             t1 = t1 + A(m, i2)*B(m, i2)
     367             :          END DO
     368             :       CASE (2)
     369       49116 :          mminus3 = m - 3
     370     1576884 :          DO i2 = 1, n
     371     1527768 :             DO i1 = 1, mminus3, 4
     372    33391824 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     373    33391824 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     374    33391824 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     375    33391824 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     376             :             END DO
     377     1527768 :             t1 = t1 + A(m - 1, i2)*B(m - 1, i2)
     378     1576884 :             t2 = t2 + A(m, i2)*B(m, i2)
     379             :          END DO
     380             :       CASE (3)
     381      399264 :          mminus3 = m - 3
     382    11007340 :          DO i2 = 1, n
     383     6041456 :             DO i1 = 1, mminus3, 4
     384    52880540 :                t1 = t1 + A(i1, i2)*B(i1, i2)
     385    52880540 :                t2 = t2 + A(i1 + 1, i2)*B(i1 + 1, i2)
     386    52880540 :                t3 = t3 + A(i1 + 2, i2)*B(i1 + 2, i2)
     387    52894076 :                t4 = t4 + A(i1 + 3, i2)*B(i1 + 3, i2)
     388             :             END DO
     389     6041456 :             t1 = t1 + A(m - 2, i2)*B(m - 2, i2)
     390     6041456 :             t2 = t2 + A(m - 1, i2)*B(m - 1, i2)
     391     6440720 :             t3 = t3 + A(m, i2)*B(m, i2)
     392             :          END DO
     393             :       END SELECT
     394     4965884 :       trace_r_AxB = t1 + t2 + t3 + t4
     395             : 
     396     4965884 :    END FUNCTION trace_r_AxB
     397             : 
     398             : END MODULE ao_util
     399             :