r12_fundamentals.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 Radial fundamentals for Gaussian two-electron integrals.
!>        Entry m+1 stores order m. Private kernels provide unnormalized coefficients;
!>        the dispatcher returns their common multiplicative normalization separately.
!>        Gaussian-product attenuation and geometry prefactors belong to the caller.
! **************************************************************************************************
MODULE r12_fundamentals
   USE gamma,                           ONLY: fgamma => fgamma_0
   USE input_constants,                 ONLY: &
        do_potential_tshpsc, do_potential_coulomb, do_potential_gaussian, do_potential_long, &
        do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
        do_potential_truncated
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
   USE t_c_g0,                          ONLY: t_c_g0_n
   USE t_sh_p_s_c,                      ONLY: trunc_cs_poly_n20
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: compute_r12_fundamentals
 
CONTAINS
 
! **************************************************************************************************
!> \brief Compute the unnormalized Coulomb Boys coefficients.
!> \param F coefficients of orders zero through m_max
!> \param m_max maximum order, nonnegative
!> \param T reduced squared distance between Gaussian product centers
! **************************************************************************************************
   SUBROUTINE compute_coulomb_fundamentals(F, m_max, T)
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: F(m_max + 1)
      REAL(dp), INTENT(IN)                               :: T
 
      CALL fgamma(m_max, t, f)
   END SUBROUTINE compute_coulomb_fundamentals
 
! **************************************************************************************************
!> \brief Compute radial coefficients for erf(omega*r)/r.
!> \param F coefficients of orders zero through m_max
!> \param m_max maximum order, nonnegative
!> \param rho reduced Gaussian exponent, positive
!> \param T reduced squared distance between Gaussian product centers
!> \param omega range-separation parameter in inverse length units
! **************************************************************************************************
   SUBROUTINE compute_long_range_fundamentals(F, m_max, rho, T, omega)
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: F(m_max + 1)
      REAL(dp), INTENT(IN)                               :: rho, T, omega
 
      INTEGER                                            :: m
      REAL(dp)                                           :: omega2, omega_corr2, scaled_T, weight
 
      omega2 = omega**2
      omega_corr2 = omega2/(omega2 + rho)
      scaled_t = t*omega_corr2
      CALL fgamma(m_max, scaled_t, f)
      weight = sqrt(omega_corr2)
      DO m = 1, m_max + 1
         f(m) = f(m)*weight
         weight = weight*omega_corr2
      END DO
   END SUBROUTINE compute_long_range_fundamentals
 
! **************************************************************************************************
!> \brief Compute radial coefficients for erfc(omega*r)/r.
!> \param F coefficients of orders zero through m_max
!> \param m_max maximum order, nonnegative
!> \param rho reduced Gaussian exponent, positive
!> \param T reduced squared distance between Gaussian product centers
!> \param omega range-separation parameter in inverse length units
!> \param work scratch space for the long-range coefficients
! **************************************************************************************************
   SUBROUTINE compute_short_range_fundamentals(F, m_max, rho, T, omega, work)
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: F(m_max + 1)
      REAL(dp), INTENT(IN)                               :: rho, T, omega
      REAL(dp), INTENT(INOUT)                            :: work(m_max + 1)
 
      CALL compute_coulomb_fundamentals(f, m_max, t)
      CALL compute_long_range_fundamentals(work, m_max, rho, t, omega)
      f = f - work
   END SUBROUTINE compute_short_range_fundamentals
 
! **************************************************************************************************
!> \brief Compute radial coefficients for Coulomb truncated at cutoff_radius.
!> \param F coefficients of orders zero through m_max
!> \param m_max maximum order, nonnegative
!> \param rho reduced Gaussian exponent, positive
!> \param T reduced squared distance between Gaussian product centers
!> \param cutoff_radius truncation radius
!> \note The truncated-Coulomb tables must be initialized before this call.
! **************************************************************************************************
   SUBROUTINE compute_truncated_fundamentals(F, m_max, rho, T, cutoff_radius)
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: F(m_max + 1)
      REAL(dp), INTENT(IN)                               :: rho, T, cutoff_radius
 
      LOGICAL                                            :: use_gamma
      REAL(dp)                                           :: R
 
      r = cutoff_radius*sqrt(rho)
      CALL t_c_g0_n(f, use_gamma, r, t, m_max)
      IF (use_gamma) CALL compute_coulomb_fundamentals(f, m_max, t)
   END SUBROUTINE compute_truncated_fundamentals
 
! **************************************************************************************************
!> \brief Compute radial coefficients for exp(-alpha*r*r).
!> \param F coefficients of orders zero through m_max
!> \param m_max maximum order, nonnegative
!> \param rho reduced Gaussian exponent, positive
!> \param T reduced squared distance between Gaussian product centers
!> \param alpha Gaussian operator exponent in inverse squared length units
! **************************************************************************************************
   SUBROUTINE compute_gaussian_fundamentals(F, m_max, rho, T, alpha)
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: F(m_max + 1)
      REAL(dp), INTENT(IN)                               :: rho, T, alpha
 
      INTEGER                                            :: m
      REAL(dp)                                           :: scaled_T, weight
 
      scaled_t = -alpha*t/(rho + alpha)
      weight = 1.0_dp
      DO m = 1, m_max + 1
         f(m) = exp(scaled_t)*weight
         weight = weight*alpha/(rho + alpha)
      END DO
   END SUBROUTINE compute_gaussian_fundamentals
 
! **************************************************************************************************
!> \brief Compose radial kernels before applying their common normalization.
!> \param F coefficients of orders zero through m_max, excluding normalization;
!>        MIX_LG includes its distinct component normalizations in F
!> \param normalization common multiplicative normalization; one for MIX_LG
!> \param m_max maximum order, nonnegative
!> \param potential_type operator selector from input_constants
!> \param rho reduced Gaussian exponent, positive
!> \param T reduced squared distance between Gaussian product centers
!> \param work reusable scratch space for composing radial kernels
!> \param omega range-separation parameter; Gaussian operators use alpha=omega**2
!> \param cutoff_radius cutoff parameter for truncated Coulomb and TShPSC
!> \param scale_coulomb Coulomb coefficient in mixed operators
!> \param scale_longrange long-range coefficient in mixed operators
!> \param scale_gaussian Gaussian coefficient in mixed operators
!> \note Spatial screening and table initialization belong to the caller.
! **************************************************************************************************
   SUBROUTINE compute_r12_fundamentals(F, normalization, m_max, potential_type, rho, T, work, &
                                       omega, cutoff_radius, scale_coulomb, scale_longrange, scale_gaussian)
      REAL(dp), INTENT(OUT)                              :: normalization
      INTEGER, INTENT(IN)                                :: m_max
      REAL(dp), INTENT(OUT)                              :: f(m_max + 1)
      INTEGER, INTENT(IN)                                :: potential_type
      REAL(dp), INTENT(IN)                               :: rho, t
      REAL(dp), INTENT(INOUT)                            :: work(m_max + 1)
      REAL(dp), INTENT(IN)                               :: omega, cutoff_radius, scale_coulomb, &
                                                            scale_longrange, scale_gaussian
 
      REAL(dp)                                           :: gaussian_normalization, r
 
      normalization = 2.0_dp*pi/rho
      SELECT CASE (potential_type)
      CASE (do_potential_coulomb)
         CALL compute_coulomb_fundamentals(f, m_max, t)
      CASE (do_potential_truncated)
         CALL compute_truncated_fundamentals(f, m_max, rho, t, cutoff_radius)
      CASE (do_potential_tshpsc)
         r = cutoff_radius*sqrt(rho)
         f = 0.0_dp
         CALL trunc_cs_poly_n20(f, r, t, m_max)
      CASE (do_potential_short)
         CALL compute_short_range_fundamentals(f, m_max, rho, t, omega, work)
      CASE (do_potential_long)
         CALL compute_long_range_fundamentals(f, m_max, rho, t, omega)
      CASE (do_potential_gaussian)
         CALL compute_gaussian_fundamentals(f, m_max, rho, t, omega**2)
         normalization = (pi/(rho + omega**2))**1.5_dp
      CASE (do_potential_mix_cl)
         CALL compute_coulomb_fundamentals(f, m_max, t)
         CALL compute_long_range_fundamentals(work, m_max, rho, t, omega)
         f = f*scale_coulomb + work*scale_longrange
      CASE (do_potential_mix_cl_trunc)
         CALL compute_truncated_fundamentals(f, m_max, rho, t, cutoff_radius)
         CALL compute_coulomb_fundamentals(work, m_max, t)
         f = f*(scale_coulomb + scale_longrange) - work*scale_longrange
         CALL compute_long_range_fundamentals(work, m_max, rho, t, omega)
         f = f + work*scale_longrange
      CASE (do_potential_mix_lg)
         CALL compute_long_range_fundamentals(f, m_max, rho, t, omega)
         CALL compute_gaussian_fundamentals(work, m_max, rho, t, omega**2)
         gaussian_normalization = (pi/(rho + omega**2))**1.5_dp
         f = work*gaussian_normalization*scale_gaussian + f*normalization*scale_longrange
         normalization = 1.0_dp
      CASE DEFAULT
         cpabort("Unsupported operator for radial r12 fundamentals")
      END SELECT
   END SUBROUTINE compute_r12_fundamentals
 
END MODULE r12_fundamentals