ai_overlap3_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 Three-center integrals over Cartesian Gaussian-type functions
!> \par Literature
!>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
!> \par History
!>      none
!> \author Dorothea Golze
! **************************************************************************************************
MODULE ai_overlap3_debug
 
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3_debug'
 
   INTEGER, PARAMETER            :: lmax = 5
 
   REAL(dp)                      :: xa, xb, xc
   REAL(dp), DIMENSION(3)        :: a, b, c
   REAL(dp), DIMENSION(3)        :: p, g
   REAL(dp)                      :: xsi, zeta, sss
 
   PRIVATE
   PUBLIC :: init_os_overlap3, os_overlap3
 
CONTAINS
 
! **************************************************************************************************
!> \brief   Calculation of three-center integrals over
!>          Cartesian Gaussian-type functions.
!> \param ya ...
!> \param yb ...
!> \param yc ...
!> \param rA ...
!> \param rB ...
!> \param rC ...
! **************************************************************************************************
   SUBROUTINE init_os_overlap3(ya, yb, yc, rA, rB, rC)
      REAL(dp)                                           :: ya, yb, yc
      REAL(dp), DIMENSION(3)                             :: ra, rb, rc
 
      REAL(dp)                                           :: fpc, ss
 
      xa = ya
      xb = yb
      xc = yc
      a = ra
      b = rb
      c = rc
 
      xsi = xa + xb
      zeta = xa*xb/xsi
 
      p = (xa*a + xb*b)/xsi
      g = (xsi*p + xc*c)/(xsi + xc)
 
      ss = (pi/xsi)**(3._dp/2._dp)*exp(-zeta*sum((a - b)**2))
 
      fpc = exp(-xsi*xc/(xsi + xc)*sum((p - c)**2))
      sss = (xsi/(xsi + xc))**(3._dp/2._dp)*ss*fpc
 
   END SUBROUTINE init_os_overlap3
 
! **************************************************************************************************
 
! **************************************************************************************************
!> \brief ...
!> \param an ...
!> \param cn ...
!> \param bn ...
!> \return ...
! **************************************************************************************************
   RECURSIVE FUNCTION os_overlap3(an, cn, bn) RESULT(IACB)
      INTEGER, DIMENSION(3)                              :: an, cn, bn
      REAL(dp)                                           :: iacb
 
      INTEGER, DIMENSION(3), PARAMETER                   :: i1 = (/1, 0, 0/), i2 = (/0, 1, 0/), &
                                                            i3 = (/0, 0, 1/)
 
      iacb = 0._dp
      IF (any(an < 0)) RETURN
      IF (any(bn < 0)) RETURN
      IF (any(cn < 0)) RETURN
 
      IF (sum(an + cn + bn) == 0) THEN
         iacb = sss
         RETURN
      END IF
 
      IF (bn(1) > 0) THEN
         iacb = os_overlap3(an, cn + i1, bn - i1) + (c(1) - b(1))*os_overlap3(an, cn, bn - i1)
      ELSEIF (bn(2) > 0) THEN
         iacb = os_overlap3(an, cn + i2, bn - i2) + (c(2) - b(2))*os_overlap3(an, cn, bn - i2)
      ELSEIF (bn(3) > 0) THEN
         iacb = os_overlap3(an, cn + i3, bn - i3) + (c(3) - b(3))*os_overlap3(an, cn, bn - i3)
      ELSE
         IF (cn(1) > 0) THEN
            iacb = os_overlap3(an + i1, cn - i1, bn) + (a(1) - c(1))*os_overlap3(an, cn - i1, bn)
         ELSEIF (cn(2) > 0) THEN
            iacb = os_overlap3(an + i2, cn - i2, bn) + (a(2) - c(2))*os_overlap3(an, cn - i2, bn)
         ELSEIF (cn(3) > 0) THEN
            iacb = os_overlap3(an + i3, cn - i3, bn) + (a(3) - c(3))*os_overlap3(an, cn - i3, bn)
         ELSE
            IF (an(1) > 0) THEN
               iacb = (g(1) - a(1))*os_overlap3(an - i1, cn, bn) + &
                      0.5_dp*(an(1) - 1)/(xsi + xc)*os_overlap3(an - i1 - i1, cn, bn)
            ELSEIF (an(2) > 0) THEN
               iacb = (g(2) - a(2))*os_overlap3(an - i2, cn, bn) + &
                      0.5_dp*(an(2) - 1)/(xsi + xc)*os_overlap3(an - i2 - i2, cn, bn)
            ELSEIF (an(3) > 0) THEN
               iacb = (g(3) - a(3))*os_overlap3(an - i3, cn, bn) + &
                      0.5_dp*(an(3) - 1)/(xsi + xc)*os_overlap3(an - i3 - i3, cn, bn)
            ELSE
               cpabort("I(0000)")
            END IF
         END IF
      END IF
 
   END FUNCTION os_overlap3
 
! **************************************************************************************************
 
END MODULE ai_overlap3_debug