Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Radial fundamentals for Gaussian two-electron integrals.
10 : !> Entry m+1 stores order m. Private kernels provide unnormalized coefficients;
11 : !> the dispatcher returns their common multiplicative normalization separately.
12 : !> Gaussian-product attenuation and geometry prefactors belong to the caller.
13 : ! **************************************************************************************************
14 : MODULE r12_fundamentals
15 : USE gamma, ONLY: fgamma => fgamma_0
16 : USE input_constants, ONLY: &
17 : do_potential_TShPSC, do_potential_coulomb, do_potential_gaussian, do_potential_long, &
18 : do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
19 : do_potential_truncated
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: pi
22 : USE t_c_g0, ONLY: t_c_g0_n
23 : USE t_sh_p_s_c, ONLY: trunc_CS_poly_n20
24 : #include "./base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 : PRIVATE
28 :
29 : PUBLIC :: compute_r12_fundamentals
30 :
31 : CONTAINS
32 :
33 : ! **************************************************************************************************
34 : !> \brief Compute the unnormalized Coulomb Boys coefficients.
35 : !> \param F coefficients of orders zero through m_max
36 : !> \param m_max maximum order, nonnegative
37 : !> \param T reduced squared distance between Gaussian product centers
38 : ! **************************************************************************************************
39 180168881 : SUBROUTINE compute_coulomb_fundamentals(F, m_max, T)
40 : INTEGER, INTENT(IN) :: m_max
41 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
42 : REAL(dp), INTENT(IN) :: T
43 :
44 180168881 : CALL fgamma(m_max, T, F)
45 180168881 : END SUBROUTINE compute_coulomb_fundamentals
46 :
47 : ! **************************************************************************************************
48 : !> \brief Compute radial coefficients for erf(omega*r)/r.
49 : !> \param F coefficients of orders zero through m_max
50 : !> \param m_max maximum order, nonnegative
51 : !> \param rho reduced Gaussian exponent, positive
52 : !> \param T reduced squared distance between Gaussian product centers
53 : !> \param omega range-separation parameter in inverse length units
54 : ! **************************************************************************************************
55 31147533 : SUBROUTINE compute_long_range_fundamentals(F, m_max, rho, T, omega)
56 : INTEGER, INTENT(IN) :: m_max
57 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
58 : REAL(dp), INTENT(IN) :: rho, T, omega
59 :
60 : INTEGER :: m
61 : REAL(dp) :: omega2, omega_corr2, scaled_T, weight
62 :
63 31147533 : omega2 = omega**2
64 31147533 : omega_corr2 = omega2/(omega2 + rho)
65 31147533 : scaled_T = T*omega_corr2
66 31147533 : CALL fgamma(m_max, scaled_T, F)
67 31147533 : weight = SQRT(omega_corr2)
68 132199456 : DO m = 1, m_max + 1
69 101051923 : F(m) = F(m)*weight
70 132199456 : weight = weight*omega_corr2
71 : END DO
72 31147533 : END SUBROUTINE compute_long_range_fundamentals
73 :
74 : ! **************************************************************************************************
75 : !> \brief Compute radial coefficients for erfc(omega*r)/r.
76 : !> \param F coefficients of orders zero through m_max
77 : !> \param m_max maximum order, nonnegative
78 : !> \param rho reduced Gaussian exponent, positive
79 : !> \param T reduced squared distance between Gaussian product centers
80 : !> \param omega range-separation parameter in inverse length units
81 : !> \param work scratch space for the long-range coefficients
82 : ! **************************************************************************************************
83 16285894 : SUBROUTINE compute_short_range_fundamentals(F, m_max, rho, T, omega, work)
84 : INTEGER, INTENT(IN) :: m_max
85 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
86 : REAL(dp), INTENT(IN) :: rho, T, omega
87 : REAL(dp), INTENT(INOUT) :: work(m_max + 1)
88 :
89 16285894 : CALL compute_coulomb_fundamentals(F, m_max, T)
90 16285894 : CALL compute_long_range_fundamentals(work, m_max, rho, T, omega)
91 76850163 : F = F - work
92 16285894 : END SUBROUTINE compute_short_range_fundamentals
93 :
94 : ! **************************************************************************************************
95 : !> \brief Compute radial coefficients for Coulomb truncated at cutoff_radius.
96 : !> \param F coefficients of orders zero through m_max
97 : !> \param m_max maximum order, nonnegative
98 : !> \param rho reduced Gaussian exponent, positive
99 : !> \param T reduced squared distance between Gaussian product centers
100 : !> \param cutoff_radius truncation radius
101 : !> \note The truncated-Coulomb tables must be initialized before this call.
102 : ! **************************************************************************************************
103 263527871 : SUBROUTINE compute_truncated_fundamentals(F, m_max, rho, T, cutoff_radius)
104 : INTEGER, INTENT(IN) :: m_max
105 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
106 : REAL(dp), INTENT(IN) :: rho, T, cutoff_radius
107 :
108 : LOGICAL :: use_gamma
109 : REAL(dp) :: R
110 :
111 263527871 : R = cutoff_radius*SQRT(rho)
112 263527871 : CALL t_c_g0_n(F, use_gamma, R, T, m_max)
113 263527871 : IF (use_gamma) CALL compute_coulomb_fundamentals(F, m_max, T)
114 263527871 : END SUBROUTINE compute_truncated_fundamentals
115 :
116 : ! **************************************************************************************************
117 : !> \brief Compute radial coefficients for exp(-alpha*r*r).
118 : !> \param F coefficients of orders zero through m_max
119 : !> \param m_max maximum order, nonnegative
120 : !> \param rho reduced Gaussian exponent, positive
121 : !> \param T reduced squared distance between Gaussian product centers
122 : !> \param alpha Gaussian operator exponent in inverse squared length units
123 : ! **************************************************************************************************
124 301982 : SUBROUTINE compute_gaussian_fundamentals(F, m_max, rho, T, alpha)
125 : INTEGER, INTENT(IN) :: m_max
126 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
127 : REAL(dp), INTENT(IN) :: rho, T, alpha
128 :
129 : INTEGER :: m
130 : REAL(dp) :: scaled_T, weight
131 :
132 301982 : scaled_T = -alpha*T/(rho + alpha)
133 301982 : weight = 1.0_dp
134 1104172 : DO m = 1, m_max + 1
135 802190 : F(m) = EXP(scaled_T)*weight
136 1104172 : weight = weight*alpha/(rho + alpha)
137 : END DO
138 301982 : END SUBROUTINE compute_gaussian_fundamentals
139 :
140 : ! **************************************************************************************************
141 : !> \brief Compose radial kernels before applying their common normalization.
142 : !> \param F coefficients of orders zero through m_max, excluding normalization;
143 : !> MIX_LG includes its distinct component normalizations in F
144 : !> \param normalization common multiplicative normalization; one for MIX_LG
145 : !> \param m_max maximum order, nonnegative
146 : !> \param potential_type operator selector from input_constants
147 : !> \param rho reduced Gaussian exponent, positive
148 : !> \param T reduced squared distance between Gaussian product centers
149 : !> \param work reusable scratch space for composing radial kernels
150 : !> \param omega range-separation parameter; Gaussian operators use alpha=omega**2
151 : !> \param cutoff_radius cutoff parameter for truncated Coulomb and TShPSC
152 : !> \param scale_coulomb Coulomb coefficient in mixed operators
153 : !> \param scale_longrange long-range coefficient in mixed operators
154 : !> \param scale_gaussian Gaussian coefficient in mixed operators
155 : !> \note Spatial screening and table initialization belong to the caller.
156 : ! **************************************************************************************************
157 430754328 : SUBROUTINE compute_r12_fundamentals(F, normalization, m_max, potential_type, rho, T, work, &
158 : omega, cutoff_radius, scale_coulomb, scale_longrange, scale_gaussian)
159 : REAL(dp), INTENT(OUT) :: normalization
160 : INTEGER, INTENT(IN) :: m_max
161 : REAL(dp), INTENT(OUT) :: F(m_max + 1)
162 : INTEGER, INTENT(IN) :: potential_type
163 : REAL(dp), INTENT(IN) :: rho, T
164 : REAL(dp), INTENT(INOUT) :: work(m_max + 1)
165 : REAL(dp), INTENT(IN) :: omega, cutoff_radius, scale_coulomb, &
166 : scale_longrange, scale_gaussian
167 :
168 : REAL(dp) :: gaussian_normalization, R
169 :
170 430754328 : normalization = 2.0_dp*pi/rho
171 574041742 : SELECT CASE (potential_type)
172 : CASE (do_potential_coulomb)
173 143287414 : CALL compute_coulomb_fundamentals(F, m_max, T)
174 : CASE (do_potential_truncated)
175 255158713 : CALL compute_truncated_fundamentals(F, m_max, rho, T, cutoff_radius)
176 : CASE (do_potential_TShPSC)
177 1160668 : R = cutoff_radius*SQRT(rho)
178 6054616 : F = 0.0_dp
179 1160668 : CALL trunc_CS_poly_n20(F, R, T, m_max)
180 : CASE (do_potential_short)
181 16285894 : CALL compute_short_range_fundamentals(F, m_max, rho, T, omega, work)
182 : CASE (do_potential_long)
183 2077323 : CALL compute_long_range_fundamentals(F, m_max, rho, T, omega)
184 : CASE (do_potential_gaussian)
185 0 : CALL compute_gaussian_fundamentals(F, m_max, rho, T, omega**2)
186 0 : normalization = (pi/(rho + omega**2))**1.5_dp
187 : CASE (do_potential_mix_cl)
188 4113176 : CALL compute_coulomb_fundamentals(F, m_max, T)
189 4113176 : CALL compute_long_range_fundamentals(work, m_max, rho, T, omega)
190 14752877 : F = F*scale_coulomb + work*scale_longrange
191 : CASE (do_potential_mix_cl_trunc)
192 8369158 : CALL compute_truncated_fundamentals(F, m_max, rho, T, cutoff_radius)
193 8369158 : CALL compute_coulomb_fundamentals(work, m_max, T)
194 32066925 : F = F*(scale_coulomb + scale_longrange) - work*scale_longrange
195 8369158 : CALL compute_long_range_fundamentals(work, m_max, rho, T, omega)
196 32066925 : F = F + work*scale_longrange
197 : CASE (do_potential_mix_lg)
198 301982 : CALL compute_long_range_fundamentals(F, m_max, rho, T, omega)
199 301982 : CALL compute_gaussian_fundamentals(work, m_max, rho, T, omega**2)
200 301982 : gaussian_normalization = (pi/(rho + omega**2))**1.5_dp
201 1104172 : F = work*gaussian_normalization*scale_gaussian + F*normalization*scale_longrange
202 301982 : normalization = 1.0_dp
203 : CASE DEFAULT
204 430754328 : CPABORT("Unsupported operator for radial r12 fundamentals")
205 : END SELECT
206 430754328 : END SUBROUTINE compute_r12_fundamentals
207 :
208 : END MODULE r12_fundamentals
|
