Computing Four-Center Two-Electron Coulomb Integrals Using Exponential Transformations and Trapezoidal Rule

The numerical evaluations of the four-center two-electron Coulomb integrals are among the most time-consuming computations involved in molecular electronic structure calculations. In the present paper we extend the double exponential (DE) transform method, previously developed for the numerical evaluation of threecenter one-electron molecular integrals [J. Lovrod, H. Safouhi, Molecular Physics (2019) DOI:10.1030/0026867.2019.1619854], to four-center two-electron integrals. The fast convergence properties analyzed in the numerical section illustrate the advantages of the new approach.


Introduction
The double exponential (DE) transformation method for the numerical evaluation of three-center oneelectron molecular integrals [1] is extended in the present paper for the evaluation of the notoriously difficult four-center molecular integrals.
The analytical expression of a four-center molecular integral, which can be derived via the Fourier transform method, contains a semi-infinite integral J(s, t), the integrand of which is slowly decaying and involves the spherical Bessel function j λ (vx). In particular, when λ and/or v are large, the oscillations of the integrand can become huge, thus complicating the accurate approximation of the integral by any of the existing numerical programming platforms. By applying the S transformation [2] followed by the DE transformation [1], we obtain a bi-infinite integral, which can be approximated by a trapezoidal rule. Relatively few collocation points were required in order to obtain accurate approximations.

General definitions and properties
The spherical Bessel function j λ (x) of order λ ∈ N 0 , [3,4], is given by e-mail: jlovrod@ualberta.ca e-mail: hsafouhi@ualberta.ca The reduced Bessel functionk n+ 1 2 (z) with half-integral order [5,6] is given bŷ where n ∈ N, and where K n+ 1 2 (z) is the modified Bessel function of the second kind [7]. The B functions [6,8] are defined by where n, l, m are the quantum numbers, Y m l (θ, ϕ) denotes the surface spherical harmonic, which is defined using the Condon-Shortley phase convention [9] and P m l (x) is the associated Legendre polynomial of the lth degree and mth order.
The Fourier transformB m n,l (ζ, #p ) of B m n,l (ζ, #r ) (3) is given by [10]: The four center two-electron Coulomb integrals are given by: where n, l and m are the quantum numbers, A, B, C and D are arbitrary points in the Euclidean space, and O stands for the origin of the fixed coordinate system. In the case where A = C, we obtain the expression of three-center exchange integrals. Two-center exchange integrals correspond to the case where A = C and B = D.
By performing the translations of vectors # -OA and # -OD, we can write the four center two-electron Coulomb integrals as follows: where

Four-center two-electron Coulomb integrals
The analytical expression of the four-center two-electron Coulomb integrals in terms of B functions can be derived via the Fourier transform method [11,12]. The problematic semi-infinite oscillatory integral involved in the obtained analytical expression is given by: where s, t ∈ [0, 1], n γ 12 , n γ 34 , n x , λ are positive integers, ν 12 and ν 34 are half integers and where: The evaluation of J(s, t) is extremely difficult due to the spherical Bessel function j λ (vx), particularly when λ and v are large. The rapid oscillations of the integrand of J(s, t) (7) can be observed in Figure 1.

Application of the S transformation
The S transformation [2], converts the integrals involving spherical Bessel functions into sine integrals by first repeatedly differentiating with respect to x followed by dividing by x. Using this transformation, J(s, t) Eq. (7) becomes Now by letting the integral in (9) can be rewritten as: Since g(x) in Eq. (10) is a slowly decaying analytic function on (0, ∞) and v is constant, the semi-infinite integral after the S transformation (11) is a suitable candidate for a DE transformation.

DE transformation
The DE transformations (see Ooura and Mori [13,14] for details) provide an efficient way to approximate integrals of the form: where g 0 (x) is a slowly decaying analytic function on (0, +∞), and v is a constant. Let φ(τ) be such that φ(τ) ∼ τ as τ → ∞ and φ(τ) ∼ 0 as τ → −∞, and let M be some large positive constant. The variable transformation x = Mφ(τ) results in: The trapezoidal formula with a constant mesh size h: In our application, this summation is not symmetric. For future reference, we define: as well as Two possible DE transformation functions are [13,14]: where K = 6, α = β/ 1 + M log(1 + M)/4π, and β = 1/4.

Application to the spherical Bessel integral
Since g(x) (10) is slowly decaying and analytic on (0, ∞), the integral in (11) is of the same form as the integral in (12). It follows that J(s, t) (11) can be approximated using the DE transformation. Applying the transformation x = Mφ(τ) to (11) results in: where g is defined in (10).

Numerical discussion
The approximation (18) was implemented in Python with the symbolic computation package SymPy, and calculations were completed using IEEE 754 double precision. To ensure that our integrator produced correct approximations, we compared our results with the output of a MATLAB built-in numerical integration function that uses global adaptive quadrature. In our implementation, we truncated the summation at N + > 0 when S φ N + /I φ N + < 10 −15 and at N − < 0 when S φ N − /I φ N − < 10 −15 . As can be seen in Figure 1 and Figure 2, the S transformation converts the spherical Bessel integrals into a simpler sine integral, but the integrand has slow convergence ( Figure 2). The convergence is greatly improved by applying a DE variable transformation x = Mφ(τ) (Figure 3). Because the transformation induces rapid convergence to zero, sinc quadrature (trapezoidal rule) renders highly efficient approximation, and very few collocation points are required to obtain accurate results when using either transformation φ 1 (τ) or φ 2 (τ) (17). Usually, less collocation points are needed when using the transformation φ 2 (τ) to approximate J(s, t) to a given degree of accuracy (Figure 4).
• n φ 1 and n φ 2 represent the total number of collocation points used for the summation in (18).
• ε φ 1 and ε φ 2 are the approximate relative errors as in (20) using transformation φ 1 and φ 2 , respectively. accuracy, while still small enough to keep the calculation times low. The need for a higher value M for the transformation φ 1 (τ) in comparison with the lower value M for the transformation φ 2 (τ) is a feature similar to the results reported in [1] for the three-center molecular integrals. We calculate the relative approximate error of the approximation by: where S φ N + and I φ N + are defined in (15)  n φ 1 (τ ) φ 2 (τ ) Figure 4. The number of collocation points n needed to complete the approximation J(s, t) depending on the natural logarithm of the relative error ε with the variable assignments from the first row of Table 1.

Conclusion
The computation of the multicenter two-electron Coulomb integrals involves the highest degree of difficulty in the molecular electronic structure calculations. The present paper proposed an efficient solution for the computation of the four-center two-electron Coulomb integrals expressed in terms of B functions. To this aim, the starting expression of such an integral is successively transformed by a procedure involving the following steps: (i) simplification using translations; (ii) use of the Fourier transform method yielding an expression in terms of a highly oscillatory semi-infinite integral, the integrand of which involves a spherical Bessel function; (iii) use of the S transform of Safouhi [2] which simplifies the problem by transforming the spherical Bessel integral into a semi-infinite sine integral with slowly decaying oscillatory integrand; (iv) use of the DE transformation of Ooura and Mori [13,14] which results in a bi-infinite integral with fast decaying integrand at both infinities; (v) Sinc quadrature approximating the integral by a bi-infinite sum which can be truncated on both sides. As demonstrated in the numerical section, relatively few collocation points are needed to approximate such integrals to high accuracy.