Double Exponential Transformation For Computing Three-Centre Nuclear Attraction Integrals
aa r X i v : . [ m a t h . NA ] J a n Double Exponential Transformation For Computing Three-CentreNuclear Attraction Integrals
Jordan Lovrod and Hassan Safouhi ∗ Mathematical DivisionCampus Saint-Jean, University of Alberta8406 91 Street, Edmonton (AB) T6C 4G9, Canada
Abstract .Three-centre nuclear attraction integrals, which arise in density functional and ab initio calcula-tions, are one of the most time-consuming computations involved in molecular electronic structurecalculations. Even for relatively small systems, millions of these laborious calculations need to beexecuted. Highly efficient and accurate methods for evaluating molecular integrals are therefore allthe more vital in order to perform the calculations necessary for large systems. When using a basisset of B functions, an analytical expression for the three-centre nuclear attraction integrals can bederived via the Fourier transform method. However, due to the presence of the highly oscillatorysemi-infinite spherical Bessel integral, the analytical expression still remains problematic. By apply-ing the S transformation, the spherical Bessel integral can be converted into a much more favorablesine integral. In the present work, we then apply two types of double exponential transformationsto the resulting sine integral, which leads to a highly efficient and accurate quadrature formulae.This method facilitates the approximation of the molecular integrals to a high predetermined accu-racy, while still keeping the calculation times low. The fast convergence properties analyzed in thenumerical section illustrate the advantages of the method. Keywords .Numerical integration; Oscillatory integrals; Molecular multi-centre integrals; double exponentialtransformation; Trapezoidal rule; Bessel functions. ∗ Corresponding author: [email protected] corresponding author acknowledges the financial support for this research by the Natural Sciences and EngineeringResearch Council of Canada (NSERC) - Grant RGPIN-2016-04317. Introduction
The Schrödinger equation, first published in Schrödinger’s 1926 paper, was a ground-breaking devel-opment that provided new insight into the behavior of the electrons in a given system [1]. Solutionsto the Schrödinger equation yield valuable details regarding the electron placement, total energy,and other properties of the system. Nevertheless, ab initio calculations, that is to say calculationsthat use the positions of the nuclei and the number of electrons to solve the Schrödinger equation,still present significant computational difficulties. Although the Schrödinger equation can be solvedanalytically for hydrogen and other atomic nuclei with only one electron, the inter-electron interac-tion involved in systems comprised of more complex atoms makes the process of analytically solvingthe Schrödinger equation extremely laborious, if not impossible.Three-centre nuclear attraction integrals are the rate determining step of ab initio and densityfunctional theory (DFT) molecular structure calculations. Three-centre integrals, one of the mostcommon types of molecular multi-centre integrals, arise when each of the three nuclei occupy adifferent position in space.A common approach for calculating molecular orbitals (MOs) is to build each MO as a linear combi-nation of atomic orbitals (LCAO), where atomic orbitals (AOs) are the solutions to the Schrodingerequation for the hydrogen atom or a hydrogen-like ion. This technique is often referred to as theLCAO-MO method. The complexity of the three-centre nuclear attraction integrals, as well asthat of other multi-centre integrals, can be mitigated by strategically choosing the basis of atomicorbitals with which to perform the calculations. It has been shown that a suitable atomic orbitalbasis should satisfy two conditions for analytical solutions to the Schrödinger equation: exponentialdecay at infinity [2] and a cusp at the origin [3].One choice for the basis functions are the exponential type functions (ETFs) which behave like exp( − x ) . ETFs are appropriate wave functions in that they satisfy the aforementioned conditionsfor analytical solutions to the Schrödinger equation. Generally, ETFs are a better suited basis thanGaussian type functions (GTFs) [4, 5], which are AOs that behave like exp( − x ) . That being saidhowever, GTFs are commonly used for chemical calculations [6]. Unfortunately, these functionsdisplay neither exponential decay at infinity nor a cusp at the origin. In order to compensate, arelatively large basis set of GTFs is required. ETFs are more favorable than GTFs therefore, sincea smaller basis set can be used to obtain results that are just as accurate [7].The most common ETFs are the Slater type functions (STFs) [8]. Despite the relative simplicity oftheir analytical expression, however, the fact that multi-centre integrals over STFs can be extremelyproblematic for polyatomic molecules has averted the use of STFs as a basis set. B functions are another class of ETFs, and their use for this particular purpose was proposed byShavitt due to their remarkably simple Weierstrass transform [9]. B functions involve reduced Besselfunctions and, in fact, they can be expressed as linear combinations of STFs [10]. Compared toother ETFs, B functions are much better adapted to the evaluation of multi-centre integrals [10–13].This is in part due to their exceptionally straightforward Fourier transform [14, 15].The Fourier transformation method introduced by Bonham et al in 1964 [16], generalized by Trivedi et al in 1983 [17], and further generalized by Grotendorst et al in 1988 [18], is one of the mostsuccessful approaches yet for the evaluation of multi-centre integrals. This particular method isespecially suitable for a basis of B functions. In particular, the Fourier transformation methodallows for the development of an analytical expression for three-centre nuclear attraction integrals.The challenge, however, is that these analytical expressions involve semi-infinite spherical Bessel2ntegrals, which are highly oscillatory due to the presence of the spherical Bessel function j λ ( vx ) .Approximating oscillatory integrals can be problematic, especially when the oscillatory part is aspherical Bessel function rather than a simple trigonometric function. Moreover, as λ and v becomelarge, the zeros of the integrand become closer and the oscillations become stronger, thus furthercomplicating the numerical evaluation of the integral. It is possible to rewrite such integrals asslowly convergent infinite series of integrals of alternating sign. It may then seem appropriate toapply series convergence acceleration techniques, such as Wynn’s Epsilon Method [19] or Levin’s u transform [20]. In the case where λ and v are large, however, the excessive calculation times preventthis approach from rendering accurate results.As shown in previous work by Safouhi [21], through a reformalized integration by parts with respectto x d x , the semi-infinite spherical Bessel integral involved in the analytical expression of the three-centre nuclear attraction integrals can be transformed into an integral involving the simple sinefunction. This transformation, called the S transformation, results in a sine integral that hasequidistant zeros, thus making it much more numerically favorable than the initial spherical Besselintegral [21, 22].The S transformation supplies remarkable theoretical and computational power in computing spher-ical Bessel integrals. It applies considerable pressure on the envelope of oscillatory integrals in theasymptotic limit. While this method is highly accurate and efficient, there are some ranges ofparameters where either failure is inevitable or the computation becomes extremely heavy.In this work, after utilizing the Fourier transformation method to obtain an analytical expression,followed by the S transformation to simplify the integrand to a sine function, a double exponentialtransformation (or DE transformation) is applied. The DE transformations introduced by Ooura etal in 1991 [23] and refined by Ooura et al in 1999 [24], map the interval (0 , ∞ ) onto ( −∞ , ∞ ) and aresuch that the transformed integrand decays double exponentially at both infinities. The trapezoidalformula with an equal mesh size is then applied to the integral, and the resulting summation canbe truncated at some moderate positive and negative values. The final result is a highly efficientquadrature formula in which relatively few function evaluations are required in order to obtainhighly accurate approximations. Both the original and the refined DE transformations are appliedto the spherical Bessel integral involved in the three-centre molecular integral, and their efficiencyfor this particular problem is compared in the discussion section. The numerical results illustratethe high accuracy of these algorithms applied to three-center nuclear attraction integrals over B functions with a miscellany of different parameters. The spherical Bessel function j λ ( x ) of order λ ∈ N is defined by [25, 26]: j λ ( x ) = ( − λ x λ (cid:18) d x d x (cid:19) λ j ( x ) = ( − λ x λ (cid:18) d x d x (cid:19) λ (cid:18) sin xx (cid:19) . (1)The reduced Bessel function ˆ k n + ( z ) is defined by [9, 27]: ˆ k n + ( z ) = r π z n + K n + ( z ) = z n e − z n X j =0 ( n + j )! j ! ( n − j )! 1(2 z ) j , (2)where n ∈ N , and where K n + ( z ) is the modified Bessel function of the second kind [28].3he B functions are defined by [10, 27]: B mn,l ( ζ, ~r ) = ( ζr ) l n + l ( n + l )! ˆ k n − ( ζr ) Y ml ( θ ~r , ϕ ~r ) , (3)where n, l, m are the quantum numbers and Y ml ( θ, ϕ ) denotes the surface spherical harmonic, whichis defined using the Condon-Shortley phase convention [29]: Y ml ( θ, ϕ ) = i m + | m | (cid:20) (2 l + 1)( l −| m | )!4 π ( l + | m | )! (cid:21) / P | m | l (cos θ ) e i mϕ , (4)where P ml ( x ) is the associated Legendre polynomial of the l th degree and m th order.The Fourier transform ¯ B mn,l ( ζ, ~p ) of B mn,l ( ζ, ~r ) (3) is given by [14] ¯ B mn,l = r π ζ n + l − ( − i | p | ) l ( ζ + | p | ) n + l +1 Y ml ( θ ~p , ϕ ~p ) . (5)The Gaunt coefficients are defined by [30–33]: (cid:10) l m | l m | l m (cid:11) = Z πθ =0 Z πϕ =0 h Y m l ( θ, ϕ ) i ∗ Y m l ( θ, ϕ ) Y m l ( θ, ϕ ) sin( θ )d θ d ϕ. (6)The Gaunt coefficients linearize the product of two spherical harmonics: h Y m l ( θ, ϕ ) i ∗ Y m l ( θ, ϕ ) = l + l X l = l min , (cid:10) l , m | l , m | l, m − m (cid:11) Y m − m l ( θ, ϕ ) , (7)where the subscript l = l min , implies that the summation index l runs in steps of 2 from l min to l + l , and where the constant l min is given by: l min = ( max( | l − l | , | m − m | ) , l + l + max( | l − l | , | m − m | ) even, max( | l − l | , | m − m | ) + 1 , l + l + max( | l − l | , | m − m | ) odd. (8)The three-centre nuclear attraction integral over B functions is given by: I n ,l ,m n ,l ,m = Z (cid:20) B m n ,l (cid:16) ζ , ~R − OA (cid:17)(cid:21) ∗ | ~R − OC | B m n ,l (cid:16) ζ , ~R − OB (cid:17) d ~R, (9)where A , B and C are three arbitrary points of the euclidian space E , while O is the origin of thefixed coordinate system. n and n stand for the principal quantum numbers.After performing a translation of vector OA , we can write the integral I n ,l ,m n ,l ,m as follows: I n ,l ,m n ,l ,m ( ζ , ζ , ~R , ~R ) = Z h B m n ,l ( ζ , ~r ) i ∗ (cid:12)(cid:12)(cid:12) ~r − ~R (cid:12)(cid:12)(cid:12) B m n ,l ( ζ , ~r − ~R )d ~r, (10)where ~r = R − OA , ~R = AC and ~R = AB . 4 Double exponential transformation for the computation of semi-infinite integrals
For more details on the DE transformations and their applications, we refer the readers to the workdone by Ooura et al [23, 24]. The double exponential (DE) transformations present a particularlyefficient method for computing semi-infinite integrals of the form: Z ∞ f ( x ) sin( vx )d x, (11)where f ( x ) is a slowly decaying function and v is a constant.We will denote the semi-infinite integral we wish to evaluate by: I = Z ∞ f ( x )d x. (12)Suppose that f is a function such that for some θ, λ ∈ R and for all n ∈ N , n ≥ n , where n issome large natural number: f ( λn + θ ) = 0 . (13)In other words, after some large x value, f ( x ) has infinite equidistant zeros occurring every λ , witha shift of θ with respect to the origin. Now provided that f ( x ) is of the form of the integrandin (11), it can easily be determined that λ = π v and θ = 0 .Now let φ ( t ) be a function such that: φ ( t ) ∼ t as t → ∞ and φ ( t ) ∼ as t → −∞ . (14)Let M be a large positive constant. Applying the variable transformation x = M φ ( t ) to the integral I , we get: I = Z φ − ( ∞ ) φ − (0) f ( M φ ( t )) M φ ′ ( t )d t (15) = M Z ∞−∞ f ( M φ ( t )) φ ′ ( t )d t. (16)Now applying the trapezoidal rule with a mesh size of h and a shift of θM with respect to the origin,we obtain the approximation: I ≈ M h ∞ X n = −∞ f M φ (cid:18) nh + θM (cid:19)! φ ′ (cid:18) nh + θM (cid:19) . (17)Let M and h be such that: M h = λ, (18)and consider f (cid:18) M φ (cid:16) nh + θM (cid:17)(cid:19) for large n : f M φ (cid:18) nh + θM (cid:19)! ∼ f M (cid:18) nh + θM (cid:19)! f ( M nh + θ )= f ( λn + θ )= 0 . (19)Given the equation (19), we can truncate the summation in (17) at some positive N + ∈ Z . Fur-thermore, since by (14), φ ( t ) tends to zero as t → −∞ , it follows that we can also truncate thesummation in (17) at some negative N − ∈ Z .As a result, we obtain the approximation: I ≈ M h N + X n = N − f M φ (cid:18) nh + θM (cid:19)! φ ′ (cid:18) nh + θM (cid:19) , (20)for some N − , N + ∈ Z . One transformation that satisfies the conditions in (14) is given by [23]: φ ( t ) = t − exp( − K sinh t ) , (21)where K is some positive constant. − − . . . . . K t φ ( t ) Figure 1:
The graph of φ ( t ) given by (21), where K = 6 . Notice that φ ( t ) encounters a singularity when t = 0 . By Taylor expansion, the behaviour of φ ( t ) about zero can be determined. More precisely: φ ( t ) ∼ as t → −∞ K as t → t as t → ∞ . (22)Similarly, φ ′ ( t ) also encounters a singularity when t = 0 . By the same method, it can be easily6ound that: φ ′ ( t ) ∼ as t → −∞ as t → as t → ∞ . (23) A second transformation that satisfies the conditions in (14) is given by [24]: φ ( t ) = t − exp( − x − α (1 − e − t ) − β ( e t − (24)where α, β are constants such that: β = O (1) , α = O (( M log M ) − / ) and ≤ α ≤ β ≤ . (25)Some assignments that satisfy the constraints in (25) are given by [24]: α = β p M log(1 + M ) / π and β = 14 . (26) − − . . α + β +2 t φ ( t ) Figure 2:
The graph of φ ( t ) given by (24), with the assignments for α and β in (26), where M = 4 . . In a similar fashion to φ ( t ) , φ ( t ) and φ ′ ( t ) both have singularities when t = 0 . By Taylorexpansion, it can be determined that: φ ( t ) ∼ as t → −∞
12 + α + β as t → t as t → ∞ , (27)and that: φ ′ ( t ) ∼ as t → −∞ ( α + β + 2) + ( α − β )2( α + β + 2) as t → as t → ∞ . (28)7 The three-centre nuclear attraction integrals
Using the Fourier transform method, an analytical expression for three-centre nuclear attractionintegrals over B functions (10) can be found [17, 18]: I n ,l ,m n ,l ,m = 8(4 π ) ( − l + l (2 l + 1)!!(2 l + 1)!!( n + l + n + l + 1)! ζ n + l − ζ n + l − ( n + l )!( n + l )! × l X l ′ =0 l ′ X m ′ = − l ′ i l + l ′ ( − l ′ (cid:10) l m | l ′ m ′ | l − l ′ m − m ′ (cid:11) (2 l ′ + 1)!![2( l − l ′ ) + 1]!! × l X l ′ =0 l ′ X m ′ = − l ′ i l + l ′ ( − l ′ (cid:10) l m | l ′ m ′ | l − l ′ m − m ′ (cid:11) (2 l ′ + 1)!![2( l − l ′ ) + 1]!! × l ′ + l ′ X l = l ′ min , (cid:10) l ′ m ′ | l ′ m ′ | lm ′ − m ′ (cid:11) R l Y m ′ − m ′ l ( θ ~R , ϕ ~R ) × l − l ′ + l − l ′ X λ = l ′′ min , ( − i ) λ (cid:10) l − l ′ m − m ′ | l − l ′ m − m ′ | λµ (cid:11) × ∆ l X j =0 (cid:18) ∆ lj (cid:19) ( − j n + n + l + l − j +1 ( n + n + l + l − j + 1)! × Z s =0 s n + l + l − l ′ (1 − s ) n + l + l − l ′ Y µλ ( θ ~v , ϕ ~v ) × "Z + ∞ x =0 x n x ˆ k ν (cid:2) R γ ( s, x ) (cid:3)(cid:2) γ ( s, x ) (cid:3) n γ j λ ( vx )d x d s, (29)where the constant l min is defined in (8) and: γ ( s, x ) = q (1 − s ) ζ + sζ + s (1 − s ) x ~v = (1 − s ) ~R − ~R , v = k ~v k and R = k ~R k n x = l − l ′ + l − l ′ n γ = 2( n + l + n + l ) − ( l ′ + l ′ ) − l + 1 ν = n + n + l + l − l − j + 12 µ = ( m − m ′ ) − ( m − m ′ )∆ l = (cid:2) ( l ′ + l ′ − l ) / (cid:3) . The semi-infinite integral in (29), which will from now on be referred to as I ( s ) , is defined by: I ( s ) = Z + ∞ x n x ˆ k ν (cid:2) R γ ( s, x ) (cid:3)(cid:2) γ ( s, x ) (cid:3) n γ j λ ( vx )d x. (30)The challenge of evaluating (29) is mainly caused by the presence of the semi-infinite integral I ( s ) (30), whose integrand is highly oscillatory due to the presence of the spherical Bessel function. This8s especially the case when λ and v are large. Notice also that when the value of s is close to or ,the oscillation of the integrand become sharp. Indeed, if we make the substitution s = 0 or s = 1 in the integrand, the exponentially decaying part of the integrand: ˆ k ν (cid:2) R γ ( s, x ) (cid:3)(cid:2) γ ( s, x ) (cid:3) n γ , becomes constant, and the integrand can be reduced to x n x j λ ( vx ) . As a result, the oscillations of j λ ( vx ) cannot be restrained, and the evaluation of I ( s ) becomes most laborious. In Figure 3, thehighly oscillatory nature of the integrand can be observed. − − x I n t e g r a nd o f I ( s ) ( − ) Figure 3:
The integrand of I ( s ) where s = 0 . , ν = / , n γ = 7 , n x = 2 , λ = 2 , R = 9 , ζ = 2 , R = 3 . , and ζ = 1 (see row 9 of Table 1). S transformation We shall now reiterate a theorem which is stated and proven by Safouhi [21].
Theorem 1. [21] Suppose that f ( x ) is a function integrable on [0, ∞ [ and of the form: f ( x ) = g ( x ) j λ ( x ) , where g ( x ) ∈ C ([0, ∞ [) , which is the set of all twice continuously differentiable functions definedon the interval [0, ∞ [ . Now if g ( x ) is such that (cid:16) dxdx (cid:17) l ( x λ − g ( x )) for l = 0 , , , . . . , λ are definiedand: lim x → x l − λ +1 (cid:18) dxdx (cid:19) (cid:16) x λ − g ( x ) (cid:17) j λ − l − ( x ) = 0lim x →∞ x l − λ +1 (cid:18) dxdx (cid:19) (cid:16) x λ − g ( x ) (cid:17) j λ − l − ( x ) = 0 , (31) holds true for all l = 0 , , , . . . , λ − , then: Z ∞ f ( x )d x = Z ∞ "(cid:18) dxdx (cid:19) λ ( x λ − g ( x )) sin( x )d x.
9n short, this transformation, called the S transformation, can simplify spherical Bessel integralsinto integrals involving the sine function. As such, the integral I ( s ) can be rewritten as: I ( s ) = 1 v λ +1 Z + ∞ (cid:18) dxdx (cid:19) λ x n x + λ − ˆ k ν (cid:2) R γ ( s, x ) (cid:3)(cid:2) γ ( s, x ) (cid:3) n γ ! sin( vx ) dx. (32) − . − . . . x I n t e g r a nd o f I ( s ) a f t e r S tr a n s f o r m a t i o n Figure 4:
The integrand of I ( s ) in (32) after applying the S transformation, where s = 0 . , ν = / , n γ = 7 , n x = 2 , λ = 2 , R = 9 , ζ = 2 , R = 3 . , and ζ = 1 (see row 9 of Table 1). The integrand prior to the transformation canbe observed in Figure 3. Now, let the part of the integrand not involving the sine function be referred to as f ( x ) , i.e. let: f ( x ) = (cid:18) dxdx (cid:19) λ x n x + λ − ˆ k ν (cid:2) R γ ( s, x ) (cid:3)(cid:2) γ ( s, x ) (cid:3) n γ ! . (33)Notice that f ( x ) is a decaying function and is therefore of the same form as the function f ( x ) ,which was introduced in (11). The integrand of I ( s ) is hence a suitable candidate for a DE variabletransformation, such as φ ( t ) (21) or φ ( t ) (24).Applying the DE transformation φ ( t ) to the equation for I ( s ) in (32), we obtain the approximation: I ( s ) = 1 v λ +1 Z + ∞ x =0 f ( x ) sin( vx ) dx ≈ M hv λ +1 N + X n = N − f M φ (cid:18) nh + θM (cid:19)! φ ′ (cid:18) nh + θM (cid:19) = I M ( s ) . (34) Table 1 lists values for I ( s ) obtained using the the approximation (34) with the transformations φ ( t ) (21) and φ ( t ) (24). In this first table, we restrict the variable assignments to values thatcan be handled by a MATLAB built-in numerical integration function that uses global adaptivequadrature set to an accuracy of correct digits. In Table 2, we restrict the variable assignments10o values that cannot be handled by the MATLAB integration function. For more problematicvariable assignments, particularly when λ and/or v are large, the MATLAB automatic integratorfails to provide results to the same accuracy efficiently. We therefore opted to use the same algorithmin Python which, with the help of the symbolic computation package SymPy, was able to completethe approximations significantly faster. The relative errors corresponding to the approximationsusing each transformation are listed in Table 2 as ε φ M and ε φ M , respectively.As shown by Ooura et al [23, 24], the relative error for the approximation I M ( s ) using a DEtransformation can be written: ε M ≈ exp (cid:18) − ch (cid:19) , (35)where c is a constant. Ooura et al argue that we can assume that the relative error is approximatelygiven by [24]: ε M ≈ exp (cid:18) − Ah (cid:19) = exp (cid:18) − AMv π (cid:19) , (36)where A is a constant. For the transformation φ ( t ) (21), we use A = 2 , and for the transformation φ ( t ) (24), we use A = 5 . For the variable assignments used in Figure 6 for the transformation φ ( t ) (21), in our presentintegrator, an optimal value for M was found to be M ≈ . . Over the course of our calcula-tions for Figure 6, therefore, we assigned said value to M in order to preserve the behaviour of theapproximation that was conducted by the integrator.In our integrator, we first took the summation in (34) over positive n values, until the additionof higher terms stopped significantly contributing to the overall approximation. Similarly, we thentook the summation over negative n values, again truncating the summation once the contributionsof each term became less significant than the predetermined error tolerance. Using this method, forthe variable assignments in Figure 6 and while using transformation φ ( t ) , the summation in (34)was truncated at -96 and at 41.In particular, it is worth noting that the number of collocation points needed to approximate theportion of the integral below zero is often higher than the number of collocation points needed toapproximate the portion of the integral above zero. This is due to the asymmetry of the transformedintegrand which is represented by the summation in (34), that is to say of the integrand of I ( s ) = R ∞−∞ f (cid:0) M φ ( t ) (cid:1) · M φ ′ ( t )d t . This asymmetric behavior can be observed in Figure 5.11 − − − . . x I n t e g r a nd o f I ( s ) = R ∞ − ∞ f (cid:0) M φ ( t ) (cid:1) · M φ ′ ( t ) d t Figure 5:
The integrand of I ( s ) = R ∞−∞ f (cid:0) Mφ ( t ) (cid:1) · Mφ ′ ( t )d t where s = 0 . , ν = / , n γ = 5 , n x = 0 , λ = 0 , R = 6 . , ζ = 1 , R = 2 . , and ζ = 1 (see row 2 of Table 1). Now in order see a pattern in our data regarding the effect of the collocation points on the accuracyof the approximation, for our first point on the graph, a moderate number of median collocationpoints was taken in order to ensure sufficient proximity to the exact result. More precisely, for thefirst point on the graph in Figure 6, we took the summation in (34) over 70 collocation points, witha lower bound of -62 and an upper bound of 7. From then on, the upper and lower bounds of thesummation were each extended by one. We proceeded to extend the bounds of the summation inthis manner – increasing the number of collocation points by two at a time until we reached 138collocation points, which was the necessary amount of collocation points originally determined byour integrator (see Table 1).For the variable assignments used φ ( t ) (24), our present integrator found M ≈ . to be anoptimal value for M ; we therefore assigned said value to M to collect our data points. Furthermore,while using transformation φ ( t ) (24) and the variable assignments in Figure 6 to approximate I ( s ) ,our present integrator truncated the summation in (34) at -63 and at 33. In a similar way as wasdone for transformation φ ( t ) , for our first point of data, we took the summation over 45 collocationpoints, with a lower bound of -37 and an upper bound of 7. From then on, the upper and lowerbounds of the summation were each extended by one until we reached a total of 97 collocationpoints, which was the number of collocation points evaluated by our integrator.Notice that the graphs in Figure 6 both rapidly decrease and approach zero as the number ofcollocation points increase. Whether our integrator uses transformation φ ( t ) (21) or transformation φ ( t ) (24), the behaviour of the absolute error for the approximation of I ( s ) based on the numberof collocation points is strikingly similar. In fact, after a certain number of collocation points, bothgraphs follow a pattern of double exponential decay.With the particular variable assignments used in Figure 6, it is clear that using transformation φ ( t ) (24), less collocation points are needed for the absolute error to converge to zero. This is represen-tative of all possible variable assignments for the evaluation of I ( s ) ; with the variable assignmentswe selected, there are always more collocation points evaluated when using transformation φ ( t ) (21) (see Table 1 ).That being said however, the parameter M is perhaps just as significant in determining the efficiencyof our integrator. The consequences of M on our integrator’s efficiency will be addressed in thefollowing subsection. 12 n for the transformation φ ( t ) A b s o l u t ee rr o r( − )
50 60 70 80 900246810 n for the transformation φ ( t ) Figure 6:
Absolute error dependence of I M ( s ) (34) on the total number of collocation points using transformations φ ( t ) (21) and φ ( t ) (24), and the variable assignments from the second row of Table 1. The absolute error was foundby comparing the results of the approximation with the result obtained by a MATLAB built-in numerical integrationfunction that uses global adaptive quadrature. M In Figure 7, it can easily be seen that the absolute error converges to zero much more quickly using φ ( t ) (24). In other words, for a given set of variable assignments, a smaller mesh size is requiredusing transformation φ ( t ) (21) to obtain results that are as accurate as the results using φ ( t ) (24).The challenge is to find a value for M that is large enough such that the error converges to zero,yet small enough so as to not make h , the mesh size, prohibitively small. Clearly, given (18), as M increases in size, h decreases. For this reason, it would be counterproductive to let M be excessivelylarge, as quite a number of collocation points in the summation (34) would be needed to compensatefor the unnecessarily small mesh size. Yet if M is too small, as shown in Figure 7, the error of theapproximation will not approach zero. This is what is meant by finding an “optimal value” for M .Our present integrator finds an optimal value for M by first letting M be a relatively small positiveconstant, say M , based on the relative error tolerance input by the user and according to thefollowing formula [24]: M = − πA log (cid:0) √ ε (cid:1) , (37)and ε denotes the relative error tolerance. For the transformation φ ( t ) (21), we use A = 2 , and forthe transformation φ ( t ) (24), we use A = 5 . For details regarding the selection of each subsequent M i , we direct the reader to the algorithm in [24].The integrator then calculates I M ( s ) using equation (34), with the value of M assigned to M .The result is not a viable approximation in the sense that it is quite far away from the desired resultdue to the fact that the value assigned to M is too small to expect the error to converge to zero.Rather, I M ( s ) serves as a value to which a second, more accurate approximation can be comparedfor the purpose of determining a relative error.For the second evaluation of (34), we assign a larger and more suitable value to M , say M . Theensuing approximation, I M ( s ) , often times meets the user’s relative error tolerance requirements,at which point the integrator outputs the value of I M ( s ) as its final result. Otherwise, the processmust be repeated with another value assigned to M , until a result of sufficient accuracy is found.This could take up to a maximum of four different assignments to M , each of which is larger thanthe previous. 13 M A b s o l u t ee rr o r( − ) φ ( t ) φ ( t ) Figure 7:
Absolute error dependence of I M ( s ) (34) on the total number of collocation points using transformations φ ( t ) (21) and φ ( t ) (24), and the variable assignments from the second row of Table 1. The absolute error was foundby comparing the results of the approximation with the result obtained by a MATLAB built-in numerical integrationfunction that uses global adaptive quadrature. Clearly, with this type of integrator, it is ideal to find an optimal value for M on the second attempt,and highly undesirable to not find an optimal value until the fourth attempt, because the approxi-mation (34) would need to be calculated four times over, meaning that numerous time-consumingcomputations would have to be executed. It is for this reason that M also holds significant weightwhile comparing the efficiency of our algorithm using transformations φ ( t ) (21) and φ ( t ) (24). Forinstance, by solely comparing the collocation points in Table 1, it may appear as though φ ( t ) (24) iswithout fail the more efficient transformation to use. That being said, since the value n M , which weare using to denote the total number of values that are assigned to M before successfully completingthe approximation for I ( s ) , also gives valuable insight into the efficiency of each transformation. Ahigh n M value signifies that our integrator had to complete the approximation (34) several timesbefore finding an optimal value for M and successfully approximating the integral, thus multiplyingthe total number of calculations that needed to be computed by the integrator.Observing the sixth row in Table 1, for example, it can be seen that although the summation in (34) istruncated with 81 collocation points using transformation φ ( t ) (21), and only 72 collocation pointsusing transformation φ ( t ) (24), the values for n M are found to be 2 and 4, using transformation φ ( t ) (21) and φ ( t ) (24), respectively. In this particular case, therefore, it can be argued that ourpresent integrator is more efficient in its approximation using transformation φ ( t ) (21). Generallyspeaking, however, the two-parameter transformation, φ ( t ) (24) lends itself to a more efficientapproximation of the semi-infinite spherical Bessel integrals. The presence of the semi-infinite spherical Bessel integral in the analytical expression for three-centrenuclear attraction integrals causes the rapid and accurate numerical evaluation of the integralsto be extremely challenging. The S transformation, introduced by Safouhi, is able to transformthe spherical Bessel integrals into considerably simpler sine integrals. The zeros of the integrandbecome equidistant, which significantly increases our ability to apply various extrapolation andtransformation methods. Nevertheless, due to its slow convergence, the semi-infinite sine integralsstill pose substantial computational difficulties. 14n the present work, we presented an efficient way of treating the arduous sine integrals that resultfrom the S transformation. Utilizing the slowly decaying, highly oscillatory nature of the integrands,we applied two distinct double exponential transformations. These transformations each renderedhighly efficient quadrature formulae that spanned from negative to positive infinity. That beingsaid, we were able to truncate the summations at relatively small positive and negative integers.It was found that remarkably few collocation points were required to obtain sufficiently accurateresults. This was especially the case for the second transformation. In our numerical section, it canbe seen that by using the second transformation, the semi-infinite spherical Bessel integral can beapproximated with a high predetermined accuracy in under 100 collocation points. The numericaltests served to further verify the precision and demonstrate the significance of the method. able 1: Computation of I ( s ) . s ν n γ n x λ R ζ R ζ I ( s ) Matlab n φ max φ n φ M ε φ M n φ max φ n φ M ε φ M / / / / / / / / / / • I ( s ) Matlab are obtained by means of a MATLAB built-in numerical integration function that uses global adaptive quadrature.• ε φ M stand for the relative errors obtained by using the approximation (34) with the transformation φ ( t ) (21).• ε φ M stand for the relative errors obtained by using the approximation (34) with the transformation φ ( t ) (24).• n φ and n φ represent the number of collocation points needed to approximate the integral using (34) with transformations φ ( t ) (21) and φ ( t ) (24),respectively.• n φ M and n φ M represent the total number of values that are assigned to M before successfully completing the approximation (34) with transformations φ ( t ) (21) and φ ( t ) (24), respectively.• max φ and max φ represent the upper limit of the summation in (34) with transformations φ ( t ) (21) and φ ( t ) (24), respectively.• The approximation (34) was implemented in Python with symbolic computation package, SymPy• The error tolerance of ε = 10 − was used in our calculation.• Calculations were completed using IEEE 754 double precision• Numbers in parentheses represent powers of 10. able 2: Computation of I ( s ) . s ν n γ n x λ R ζ R ζ I ( s ) φ n φ max φ n φ M ε φ M I ( s ) φ n φ max φ n φ M ε φ M / / / / / / /
21 6 5 55 1.5 2.0 1.5 .564 022 921 688 577(-2) 81 44 2 .1(-16) .564 022 921 688 556(-2) 84 41 4 .1(-16)0.01 /
29 6 6 60 1.5 2.0 1.0 .414 131 181 249 046( 0) 142 45 2 .1(-16) .414 131 181 249 327( 0) 93 34 2 .1(-16)0.01 /
25 7 6 55 2.0 3.0 2.0 .284 952 750 728 949(-2) 82 45 2 .1(-16) .284 952 750 728 952(-2) 82 40 4 .1(-16)0.01 /
23 7 7 55 2.5 2.0 1.0 .107 097 615 409 941(-1) 143 45 2 .1(-16) .107 097 615 410 016(-1) 93 34 2 .1(-16)0.01 /
33 7 7 65 2.0 2.0 1.0 .167 421 970 712 946(-1) 143 45 2 .1(-16) .167 421 970 713 064(-1) 93 34 2 .1(-16) • I ( s ) Matlab are obtained by means of a MATLAB built-in numerical integration function that uses global adaptive quadrature.• ε φ M stand for the relative errors obtained by using the approximation (34) with the transformation φ ( t ) (21).• ε φ M stand for the relative errors obtained by using the approximation (34) with the transformation φ ( t ) (24).• n φ and n φ represent the number of collocation points needed to approximate the integral using (34) with transformations φ ( t ) (21) and φ ( t ) (24),respectively.• n φ M and n φ M represent the total number of values that are assigned to M before successfully completing the approximation (34) with transformations φ ( t ) (21) and φ ( t ) (24), respectively.• max φ and max φ represent the upper limit of the summation in (34) with transformations φ ( t ) (21) and φ ( t ) (24), respectively.• The approximation (34) was implemented in Python with symbolic computation package, SymPy• The error tolerance of ε = 10 − was used in our calculation.• Calculations were completed using IEEE 754 double precision• Numbers in parentheses represent powers of 10. eferences [1] E. Schrödinger. An undulatory theory of the mechanics of atoms and molecules. Phys. Rev. ,28:1049–1070, 1926.[2] S. Agmon.
Bounds on exponential decay of eigenfunctions of Schrödinger operators . in S. Graffi(editor), Schrödinger operators. Springer-Verlag, Berlin, 1985.[3] T. Kato. On the eigenfunctions of many-particle systems in quantum mechanics.
Commun.Pure Appl. Math. , 10:151–177, 1957.[4] S.F. Boys. Electronic wave functions. I. a general method of calculation for the stationary statesof any molecular system.
Proc. R. Soc. Lond. Series A, Math. & Phys. Sciences. , 200:542–554,1950.[5] S.F. Boys. Electronic wave functions. II. a calculation for the ground state of the Berylliumatom.
Proc. R. Soc. Lond. Series A, Math. & Phys. Sciences. , 201:125–137, 1950.[6] S.F. Boys. The integral formulae for the variational solution of the molecular many-electronwave equation in terms of gaussian functions with direct electronic correlation.
Proc. R. Soc.Lond. Series A, Math. & Phys. Sciences. , 258:402–411, 1960.[7] R.M. Slevinsky, T. Temga, M. Mouattamid, and H. Safouhi. One- and two-center ETF-integrals of first order in relativistic calculation of NMR parameters.
J. Phys. A: Math. Theor. ,43:225202, 2010.[8] J.C. Slater. Analytic atomic wave functions.
Phys. Rev. , 42:33–43, 1932.[9] I. Shavitt.
The Gaussian function in calculation of statistical mechanics and quantum me-chanics, Methods in Computational Physics. 2. Quantum Mechanics . edited by B. Alder, S.Fernbach, M. Rotenberg, Academic Press, New York, 1963.[10] E. Filter and E.O. Steinborn. Extremely compact formulas for molecular one-electron integralsand Coulomb integrals over Slater-type orbitals.
Phys. Rev. A. , 18:1–11, 1978.[11] E. Filter and E.O. Steinborn. The three-dimensional convolution of reduced Bessel functionsof physical interest.
J. Math. Phys. , 19:79–84, 1978.[12] E.J. Weniger and E.O. Steinborn. Numerical properties of the convolution theorems of B functions. Phys. Rev. A. , 28:2026–2041, 1983.[13] E.J. Weniger. The spherical tensor gradient operator.
Collect. Czech. Chem. Commun. ,70:1125–1271, 2005.[14] E.J. Weniger and E.O. Steinborn. The Fourier transforms of some exponential-type functionsand their relevance to multicenter problems.
J. Chem. Phys. , 78:6121–6132, 1983.[15] A.W. Niukkanen. Fourier transforms of atomic orbitals. I. reduction to fourdimensional har-monics and quadratic transformations.
Int. J. Quantum Chem. , 25:941–955, 1984.[16] R.A. Bonham, J.L. Peacher, and H.L. Cox. On the calculation of multicenter two-electronrepulsion integrals involving Slater functions.
J. Chem. Phys. , 40:3083–3086, 1964.1817] H.P. Trivedi and E.O. Steinborn. Fourier transform of a two-center product of exponential-typeorbitals. application to one- and two-electron multicenter integrals.
Phys. Rev. A. , 27:670–679,1983.[18] J. Grotendorst and E.O. Steinborn. Numerical evaluation of molecular one- and two-electronmulticenter integrals with exponential-type orbitals via the Fourier-transform method.
Phys.Rev. A. , 38:3857–3876, 1988.[19] P. Wynn. On a device for computing the e m ( S n ) transformation. Math. Tables Aids Comput. ,10:91–96, 1956.[20] D. Levin. Developement of non-linear transformations for improving convergence of sequences.
Int. J. Comput. Math. , B3:371–388, 1973.[21] H. Safouhi. The properties of sine, spherical Bessel and reduced Bessel functions for improvingconvergence of semi-infinite very oscillatory integrals: The evaluation of three-center nuclearattraction integrals over B functions. J. Phys. A: Math. Gen. , 34:2801–2818, 2001.[22] L. Berlu and H. Safouhi. An extremely efficient and rapid algorithm for a numerical evaluationof three-center nuclear attraction integrals over Slater type functions.
J. Phys. A: Math. Gen. ,36:11791–11805, 2003.[23] T. Ooura and M. Mori. The double exponential formula for oscillatory functions over theinfinite interval.
J. Comput. Appl. Math. , 38:353–360, 1991.[24] T. Ooura and M. Mori. A robust double exponential formula for fourier-type integrals.
J.Comput. Appl. Math. , 112:229–241, 1999.[25] G.B. Arfken and H.J. Weber.
Mathematical methods for Physicists . Academic Press, Fifthedition, 1995.[26] M. Abramowitz and I.A. Stegun.
Handbook of Mathematical Functions . Dover, New York,1965.[27] E.O. Steinborn and E. Filter. Translations of fields represented by spherical-harmonics ex-pansions for molecular calculations. III. Translations of reduced Bessel functions, Slater-types-orbitals, and other functions.
Theor. Chim. Acta. , 38:273–281, 1975.[28] G.N. Watson.
A Treatise on the Theory of Bessel Functions . Cambridge University Press,Second Edition, Cambridge, England, 1944.[29] E.U. Condon and G.H. Shortley.
The theory of atomic spectra . Cambridge University Press,Cambridge, England, 1951.[30] J.A. Gaunt. The triplets of helium.
Phil. Trans. Roy. Soc. , A. 228:151–196, 1929.[31] H.H.H. Homeier and E.O. Steinborn. Some properties of the coupling coefficients of real spheri-cal harmonics and their relation to Gaunt coefficients.
J. Mol. Struct. (THEOCEM) , 368:31–37,1996.[32] E.J. Weniger and E.O. Steinborn. Programs for the coupling of spherical harmonics.
Comput.Phys.Commun. , 25:149–157, 1982.[33] Yu-Lin Xu. Efficient evaluation of vector translation coefficients in multiparticle light-scatteringtheories.