ai_overlap_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 Two-center overlap 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_overlap_debug
 
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap_debug'
 
   INTEGER, PARAMETER            :: lmax = 5
 
   REAL(dp)                      :: xa, xb
   REAL(dp), DIMENSION(3)        :: a, b
   REAL(dp), DIMENSION(3)        :: p
   REAL(dp)                      :: xsi, zeta, ss
 
   PRIVATE
   PUBLIC :: init_os_overlap2, os_overlap2
 
CONTAINS
 
! **************************************************************************************************
!> \brief   Calculation of overlap integrals over
!>          Cartesian Gaussian-type functions.
!> \param ya ...
!> \param yb ...
!> \param rA ...
!> \param rB ...
! **************************************************************************************************
   SUBROUTINE init_os_overlap2(ya, yb, rA, rB)
      REAL(dp)                                           :: ya, yb
      REAL(dp), DIMENSION(3)                             :: ra, rb
 
      xa = ya
      xb = yb
      a = ra
      b = rb
 
      xsi = xa + xb
      zeta = xa*xb/xsi
 
      p = (xa*a + xb*b)/xsi
 
      ss = (pi/xsi)**(3._dp/2._dp)*exp(-zeta*sum((a - b)**2))
 
   END SUBROUTINE init_os_overlap2
 
! **************************************************************************************************
 
! **************************************************************************************************
!> \brief ...
!> \param an ...
!> \param bn ...
!> \return ...
! **************************************************************************************************
   RECURSIVE FUNCTION os_overlap2(an, bn) RESULT(IAB)
      INTEGER, DIMENSION(3)                              :: an, bn
      REAL(dp)                                           :: iab
 
      INTEGER, DIMENSION(3), PARAMETER                   :: i1 = (/1, 0, 0/), i2 = (/0, 1, 0/), &
                                                            i3 = (/0, 0, 1/)
 
      iab = 0._dp
      IF (any(an < 0)) RETURN
      IF (any(bn < 0)) RETURN
 
      IF (sum(an + bn) == 0) THEN
         iab = ss
         RETURN
      END IF
 
      IF (bn(1) > 0) THEN
         iab = os_overlap2(an + i1, bn - i1) + (a(1) - b(1))*os_overlap2(an, bn - i1)
      ELSEIF (bn(2) > 0) THEN
         iab = os_overlap2(an + i2, bn - i2) + (a(2) - b(2))*os_overlap2(an, bn - i2)
      ELSEIF (bn(3) > 0) THEN
         iab = os_overlap2(an + i3, bn - i3) + (a(3) - b(3))*os_overlap2(an, bn - i3)
      ELSE
         IF (an(1) > 0) THEN
            iab = (p(1) - a(1))*os_overlap2(an - i1, bn) + &
                  0.5_dp*(an(1) - 1)/xsi*os_overlap2(an - i1 - i1, bn)
         ELSEIF (an(2) > 0) THEN
            iab = (p(2) - a(2))*os_overlap2(an - i2, bn) + &
                  0.5_dp*(an(2) - 1)/xsi*os_overlap2(an - i2 - i2, bn)
         ELSEIF (an(3) > 0) THEN
            iab = (p(3) - a(3))*os_overlap2(an - i3, bn) + &
                  0.5_dp*(an(3) - 1)/xsi*os_overlap2(an - i3 - i3, bn)
         ELSE
            cpabort("I(0000)")
         END IF
      END IF
 
   END FUNCTION os_overlap2
 
! **************************************************************************************************
 
END MODULE ai_overlap_debug