ai_eri_debug.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Calculation of Coulomb integrals over Cartesian Gaussian-type functions
!>      (electron repulsion integrals, ERIs).
!> \par Literature
!>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
!> \par History
!>      none
!> \author J. Hutter (07.2009)
! **************************************************************************************************
MODULE ai_eri_debug
 
   USE gamma,                           ONLY: fgamma_ref
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_eri_debug'
 
   INTEGER, PARAMETER            :: lmax = 5
 
   REAL(dp)                      :: xa, xb, xc, xd
   REAL(dp), DIMENSION(3)        :: a, b, c, d
   REAL(dp), DIMENSION(3)        :: p, q, w
   REAL(dp)                      :: xsi, eta, rho, t
 
   REAL(dp), DIMENSION(0:4*lmax) :: fm, i0m
 
   PRIVATE
 
CONTAINS
 
! **************************************************************************************************
!> \brief   Calculation of the primitive two-center Coulomb integrals over
!>          Cartesian Gaussian-type functions.
!> \param ya ...
!> \param yb ...
!> \param yc ...
!> \param yd ...
!> \param rA ...
!> \param rB ...
!> \param rC ...
!> \param rD ...
!> \date    07.2009
!> \author  J. Hutter
!> \version 1.0
! **************************************************************************************************
   SUBROUTINE init_os(ya, yb, yc, yd, rA, rB, rC, rD)
      REAL(dp)                                           :: ya, yb, yc, yd
      REAL(dp), DIMENSION(3)                             :: ra, rb, rc, rd
 
      REAL(dp)                                           :: eab, ecd, kab, kcd
 
      xa = ya
      xb = yb
      xc = yc
      xd = yd
      a = ra
      b = rb
      c = rc
      d = rd
 
      xsi = xa + xb
      eta = xc + xd
 
      p = (xa*a + xb*b)/xsi
      q = (xc*c + xd*d)/eta
      w = (xsi*p + eta*q)/(xsi + eta)
 
      rho = xsi*eta/(xsi + eta)
 
      t = rho*sum((p - q)**2)
 
      fm = fgamma_ref(4*lmax, t)
 
      eab = -xa*xb/xsi*sum((a - b)**2)
      kab = sqrt(2._dp)*pi**1.25_dp/xsi*exp(eab)
 
      ecd = -xc*xd/eta*sum((c - d)**2)
      kcd = sqrt(2._dp)*pi**1.25_dp/eta*exp(ecd)
 
      i0m = kab*kcd/sqrt(xsi + eta)*fm
 
   END SUBROUTINE init_os
 
! **************************************************************************************************
 
! **************************************************************************************************
!> \brief ...
!> \param an ...
!> \param bn ...
!> \param cn ...
!> \param dn ...
!> \param mi ...
!> \return ...
! **************************************************************************************************
   RECURSIVE FUNCTION os(an, bn, cn, dn, mi) RESULT(IABCD)
      INTEGER, DIMENSION(3)                              :: an, bn, cn, dn
      INTEGER, OPTIONAL                                  :: mi
      REAL(dp)                                           :: iabcd
 
      INTEGER, DIMENSION(3), PARAMETER                   :: i1 = (/1, 0, 0/), i2 = (/0, 1, 0/), &
                                                            i3 = (/0, 0, 1/)
 
      INTEGER                                            :: m
 
      m = 0
      IF (PRESENT(mi)) m = mi
 
      iabcd = 0._dp
      IF (any(an < 0)) RETURN
      IF (any(bn < 0)) RETURN
      IF (any(cn < 0)) RETURN
      IF (any(dn < 0)) RETURN
 
      IF (sum(an + bn + cn + dn) == 0) THEN
         iabcd = i0m(m)
         RETURN
      END IF
 
      IF (dn(1) > 0) THEN
         iabcd = os(an, bn, cn + i1, dn - i1) - (d(1) - c(1))*os(an, bn, cn, dn - i1)
      ELSEIF (dn(2) > 0) THEN
         iabcd = os(an, bn, cn + i2, dn - i2) - (d(2) - c(2))*os(an, bn, cn, dn - i2)
      ELSEIF (dn(3) > 0) THEN
         iabcd = os(an, bn, cn + i3, dn - i3) - (d(3) - c(3))*os(an, bn, cn, dn - i3)
      ELSE
         IF (bn(1) > 0) THEN
            iabcd = os(an + i1, bn - i1, cn, dn) - (b(1) - a(1))*os(an, bn - i1, cn, dn)
         ELSEIF (bn(2) > 0) THEN
            iabcd = os(an + i2, bn - i2, cn, dn) - (b(2) - a(2))*os(an, bn - i2, cn, dn)
         ELSEIF (bn(3) > 0) THEN
            iabcd = os(an + i3, bn - i3, cn, dn) - (b(3) - a(3))*os(an, bn - i3, cn, dn)
         ELSE
            IF (cn(1) > 0) THEN
               iabcd = ((q(1) - c(1)) + xsi/eta*(p(1) - a(1)))*os(an, bn, cn - i1, dn) + &
                       0.5_dp*an(1)/eta*os(an - i1, bn, cn - i1, dn) + &
                       0.5_dp*(cn(1) - 1)/eta*os(an, bn, cn - i1 - i1, dn) - &
                       xsi/eta*os(an + i1, bn, cn - i1, dn)
            ELSEIF (cn(2) > 0) THEN
               iabcd = ((q(2) - c(2)) + xsi/eta*(p(2) - a(2)))*os(an, bn, cn - i2, dn) + &
                       0.5_dp*an(2)/eta*os(an - i2, bn, cn - i2, dn) + &
                       0.5_dp*(cn(2) - 1)/eta*os(an, bn, cn - i2 - i2, dn) - &
                       xsi/eta*os(an + i2, bn, cn - i2, dn)
            ELSEIF (cn(3) > 0) THEN
               iabcd = ((q(3) - c(3)) + xsi/eta*(p(3) - a(3)))*os(an, bn, cn - i3, dn) + &
                       0.5_dp*an(3)/eta*os(an - i3, bn, cn - i3, dn) + &
                       0.5_dp*(cn(3) - 1)/eta*os(an, bn, cn - i3 - i3, dn) - &
                       xsi/eta*os(an + i3, bn, cn - i3, dn)
            ELSE
               IF (an(1) > 0) THEN
                  iabcd = (p(1) - a(1))*os(an - i1, bn, cn, dn, m) + &
                          (w(1) - p(1))*os(an - i1, bn, cn, dn, m + 1) + &
                          0.5_dp*(an(1) - 1)/xsi*os(an - i1 - i1, bn, cn, dn, m) - &
                          0.5_dp*(an(1) - 1)/xsi*rho/xsi*os(an - i1 - i1, bn, cn, dn, m + 1)
               ELSEIF (an(2) > 0) THEN
                  iabcd = (p(2) - a(2))*os(an - i2, bn, cn, dn, m) + &
                          (w(2) - p(2))*os(an - i2, bn, cn, dn, m + 1) + &
                          0.5_dp*(an(2) - 1)/xsi*os(an - i2 - i2, bn, cn, dn, m) - &
                          0.5_dp*(an(2) - 1)/xsi*rho/xsi*os(an - i2 - i2, bn, cn, dn, m + 1)
               ELSEIF (an(3) > 0) THEN
                  iabcd = (p(3) - a(3))*os(an - i3, bn, cn, dn, m) + &
                          (w(3) - p(3))*os(an - i3, bn, cn, dn, m + 1) + &
                          0.5_dp*(an(3) - 1)/xsi*os(an - i3 - i3, bn, cn, dn, m) - &
                          0.5_dp*(an(3) - 1)/xsi*rho/xsi*os(an - i3 - i3, bn, cn, dn, m + 1)
               ELSE
                  cpabort("I(0000)")
               END IF
            END IF
         END IF
      END IF
 
   END FUNCTION os
 
! **************************************************************************************************
 
END MODULE ai_eri_debug