generic_shg_integrals_init.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation
!>        of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of
!>        the contraction matrices
!> \par History
!>      created [08.2016]
!> \author Dorothea Golze
! **************************************************************************************************
MODULE generic_shg_integrals_init
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fac,&
                                              ifac,&
                                              pi
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init'
 
   PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, contraction_matrix_shg_rx2m
 
! **************************************************************************************************
 
CONTAINS
 
! **************************************************************************************************
!> \brief contraction matrix for SHG integrals
!> \param basis ...
!> \param scon_shg contraction matrix
! **************************************************************************************************
   SUBROUTINE contraction_matrix_shg(basis, scon_shg)
 
      TYPE(gto_basis_set_type), POINTER                  :: basis
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: scon_shg
 
      INTEGER                                            :: ipgf, iset, ishell, l, maxpgf, maxshell, &
                                                            nset
      INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
      REAL(kind=dp)                                      :: aif, gcc
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: norm
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: zet
 
      nset = basis%nset
      npgf => basis%npgf
      nshell => basis%nshell
      zet => basis%zet
 
      maxpgf = SIZE(basis%gcc, 1)
      maxshell = SIZE(basis%gcc, 2)
      ALLOCATE (norm(basis%nset, maxshell))
      ALLOCATE (scon_shg(maxpgf, maxshell, nset))
      scon_shg = 0.0_dp
 
      CALL basis_norm_shg(basis, norm)
 
      DO iset = 1, nset
         DO ishell = 1, nshell(iset)
            l = basis%l(ishell, iset)
            DO ipgf = 1, npgf(iset)
               aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
               gcc = basis%gcc(ipgf, ishell, iset)
               scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
            END DO
         END DO
      END DO
 
      DEALLOCATE (norm)
 
   END SUBROUTINE contraction_matrix_shg
 
!***************************************************************************************************
!> \brief normalization solid harmonic Gaussians (SHG)
!> \param basis ...
!> \param norm ...
! **************************************************************************************************
   SUBROUTINE basis_norm_shg(basis, norm)
 
      TYPE(gto_basis_set_type), POINTER                  :: basis
      REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: norm
 
      INTEGER                                            :: ipgf, iset, ishell, jpgf, l
      REAL(kind=dp)                                      :: aai, aaj, cci, ccj, expa, ppl
 
      norm = 0._dp
 
      DO iset = 1, basis%nset
         DO ishell = 1, basis%nshell(iset)
            l = basis%l(ishell, iset)
            expa = 0.5_dp*real(2*l + 3, dp)
            ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1)
            ppl = ppl/(2._dp**real(2*l + 1, dp))
            ppl = ppl/real(2*l + 1, dp)
            DO ipgf = 1, basis%npgf(iset)
               cci = basis%gcc(ipgf, ishell, iset)
               aai = basis%zet(ipgf, iset)
               DO jpgf = 1, basis%npgf(iset)
                  ccj = basis%gcc(jpgf, ishell, iset)
                  aaj = basis%zet(jpgf, iset)
                  norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa
               END DO
            END DO
            norm(iset, ishell) = 1.0_dp/sqrt(norm(iset, ishell))
         END DO
      END DO
 
   END SUBROUTINE basis_norm_shg
 
! **************************************************************************************************
!> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis
!>        at the same atom
!> \param orb_basis orbital basis
!> \param ri_basis   ...
!> \param orb_index index for orbital basis
!> \param ri_index index for ri basis
!> \param scon_mix mixed contraction matrix
! **************************************************************************************************
   SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix)
 
      TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
      INTEGER, DIMENSION(:, :, :), POINTER               :: orb_index, ri_index
      REAL(kind=dp), DIMENSION(:, :, :, :), POINTER      :: scon_mix
 
      INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, &
         lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, &
         nset_orb, nset_ri
      INTEGER, DIMENSION(:), POINTER                     :: npgf_orb, npgf_ri, nshell_orb, nshell_ri
      REAL(kind=dp)                                      :: cjf, const, const1, const2, gcc_orb, &
                                                            gcc_ri, prefac, scon1, scon2, zet
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: shg_fac
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: norm_orb, norm_ri, zet_orb, zet_ri
 
      nset_orb = orb_basis%nset
      npgf_orb => orb_basis%npgf
      nshell_orb => orb_basis%nshell
      zet_orb => orb_basis%zet
      nset_ri = ri_basis%nset
      npgf_ri => ri_basis%npgf
      nshell_ri => ri_basis%nshell
      zet_ri => ri_basis%zet
 
      maxpgf_orb = SIZE(orb_basis%gcc, 1)
      maxshell_orb = SIZE(orb_basis%gcc, 2)
      ALLOCATE (norm_orb(nset_orb, maxshell_orb))
      maxpgf_ri = SIZE(ri_basis%gcc, 1)
      maxshell_ri = SIZE(ri_basis%gcc, 2)
      ALLOCATE (norm_ri(nset_ri, maxshell_ri))
 
      CALL basis_norm_shg(orb_basis, norm_orb)
      CALL basis_norm_shg(ri_basis, norm_ri)
 
      ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb))
      ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri))
 
      !** index orbital basis set
      nf_orb = 0
      DO iset = 1, nset_orb
         DO ishell = 1, nshell_orb(iset)
            DO ipgf = 1, npgf_orb(iset)
               nf_orb = nf_orb + 1
               orb_index(ipgf, ishell, iset) = nf_orb
            END DO
         END DO
      END DO
 
      !** index ri basis set
      nf_ri = 0
      DO iset = 1, nset_ri
         DO ishell = 1, nshell_ri(iset)
            DO ipgf = 1, npgf_ri(iset)
               nf_ri = nf_ri + 1
               ri_index(ipgf, ishell, iset) = nf_ri
            END DO
         END DO
      END DO
 
      lmax_orb = maxval(orb_basis%lmax)
      lmax_ri = maxval(ri_basis%lmax)
      nl_max = int((lmax_orb + lmax_ri)/2) + 1
      ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max))
      scon_mix = 0.0_dp
 
      ALLOCATE (shg_fac(0:nl_max - 1))
      shg_fac(0) = 1.0_dp
 
      DO iset = 1, nset_orb
         DO ishell = 1, nshell_orb(iset)
            l1 = orb_basis%l(ishell, iset)
            const1 = sqrt(1.0_dp/real(2*l1 + 1, dp))
            DO jset = 1, nset_ri
               DO jshell = 1, nshell_ri(jset)
                  l2 = ri_basis%l(jshell, jset)
                  const2 = sqrt(1.0_dp/real(2*l2 + 1, dp))
                  nl = int((l1 + l2)/2)
                  IF (l1 == 0 .OR. l2 == 0) nl = 0
                  DO il = 0, nl
                     l = l1 + l2 - 2*il
                     const = const1*const2*2.0_dp*sqrt(pi*real(2*l + 1, dp))
                     DO iil = 1, il
                        shg_fac(iil) = fac(l + iil - 1)*ifac(l)*real(l, dp) &
                                       *fac(il)/fac(il - iil)/fac(iil)
                     END DO
                     DO ipgf = 1, npgf_orb(iset)
                        forb = orb_index(ipgf, ishell, iset)
                        gcc_orb = orb_basis%gcc(ipgf, ishell, iset)
                        scon1 = norm_orb(iset, ishell)*gcc_orb
                        DO jpgf = 1, npgf_ri(jset)
                           fri = ri_index(jpgf, jshell, jset)
                           gcc_ri = ri_basis%gcc(jpgf, jshell, jset)
                           scon2 = norm_ri(jset, jshell)*gcc_ri
                           zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset)
                           cjf = 1.0_dp/((2._dp*zet)**l)
                           prefac = const*cjf*scon1*scon2
                           DO iil = 0, il
                              scon_mix(iil + 1, fri, forb, il + 1) = prefac*shg_fac(iil)/zet**iil
                           END DO
                        END DO
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
 
      DEALLOCATE (norm_orb, norm_ri, shg_fac)
 
   END SUBROUTINE contraction_matrix_shg_mix
 
! **************************************************************************************************
!> \brief ...
!> \param basis ...
!> \param m ...
!> \param scon_shg ...
!> \param scon_rx2m ...
! **************************************************************************************************
   SUBROUTINE contraction_matrix_shg_rx2m(basis, m, scon_shg, scon_rx2m)
 
      TYPE(gto_basis_set_type), POINTER                  :: basis
      INTEGER, INTENT(IN)                                :: m
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_shg
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: scon_rx2m
 
      INTEGER                                            :: ipgf, iset, ishell, j, l, maxpgf, &
                                                            maxshell, nset
      INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: shg_fac
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: zet
 
      npgf => basis%npgf
      nshell => basis%nshell
      zet => basis%zet
      nset = basis%nset
 
      maxpgf = SIZE(basis%gcc, 1)
      maxshell = SIZE(basis%gcc, 2)
      ALLOCATE (scon_rx2m(maxpgf, m + 1, maxshell, nset))
      scon_rx2m = 0.0_dp
      ALLOCATE (shg_fac(0:m))
      shg_fac(0) = 1.0_dp
 
      DO iset = 1, nset
         DO ishell = 1, nshell(iset)
            l = basis%l(ishell, iset)
            DO j = 1, m
               shg_fac(j) = fac(l + j - 1)*ifac(l)*real(l, dp) &
                            *fac(m)/fac(m - j)/fac(j)
            END DO
            DO ipgf = 1, npgf(iset)
               DO j = 0, m
                  scon_rx2m(ipgf, j + 1, ishell, iset) = scon_shg(ipgf, ishell, iset) &
                                                         *shg_fac(j)/zet(ipgf, iset)**j
               END DO
            END DO
         END DO
      END DO
 
      DEALLOCATE (shg_fac)
 
   END SUBROUTINE contraction_matrix_shg_rx2m
 
END MODULE generic_shg_integrals_init