Double Exponential Transformation For Computing Three-Centre

arXiv:1901.00260v2 [math.NA] 3 Jan 2019

Nuclear Attraction Integrals

Jordan Lovrod and Hassan Safouhi

∗

Mathematical Division

Campus Saint-Jean, University of Alberta

8406 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 structure

calculations. Even for relatively small systems, millions of these laborious calculations need to be

executed. Highly eﬃcient and accurate methods for evaluati ng molecular integrals are therefore all

the more vital in order to perform the calculations necessary for large systems. When using a basis

set of B functions, an analytical expression for the three-centre nuclear attraction integrals can be

derived via the Fourier transform method. However, due to the presence of the highly oscillatory

semi-inﬁnite spherical Bessel integral, the analytical exp ression still remains problematic. By apply-

ing the S transformation, the spherical Bessel integral can be converted into a much more favorable

sine integral. In the present work, we then apply two types of double exponential transformations

to the resulting sine integral, which leads to a highly eﬃcient 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 the

numerical s ection illustrate the advantages of the method.

Keywords.

Numerical integration; Os cillatory integrals; Molecular multi-centre integrals; double exponential

transformation; Trapezoidal rule; Bessel functions.

∗

Corresponding author: hsafouhi@ualberta.ca.

The corresponding author acknowledges the ﬁnancial support for this research by the Natural Sciences and Engineering

Research Council of Canada (NSERC) - Grant RGPIN-2016-04317.

1 Introduction

The Schrö dinger equation, ﬁrst 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]. Solutions

to 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 calculations

that use the positions of the nuclei and the number of electrons to solve the Schrödinger equation,

still present signiﬁcant computational diﬃculties. Although the Schrödinger equation can be solved

analytically 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 solving

the Schrödinger equation extremely laborious, if not impossible.

Three-centre nuclear attraction integrals are the rate determining step of ab initio and density

functional theory (DFT) molecular structure calculations. Three-centre integrals, one of the most

common types of molecular multi-centre integrals, arise when each of the three nuclei occupy a

diﬀerent 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 Schrodinger

equation for the hydrogen atom or a hydrogen-like ion. This technique is often referred to as the

LCAO-MO method. The complexity of the three-centre nuclear attraction integrals, as well as

that of other multi-centre integrals, can be mitigated by strategically choosing the basis of atomic

orbitals with which to perform the calculations. It has been shown that a suitable atomic orbital

basis s hould satisfy two conditions for analytical solutions to the Schrödinger equation: exponential

decay at inﬁnity [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 conditions

for analytical solutions to the Schrödinger equation. Generally, ETFs are a better suited basis than

Gaussian type functions (GTFs) [4, 5], which are AOs that behave like exp(−x

). That being said

however, GTFs are commonly us ed for chemical calculations [6 ]. Unfortunately, these f unctions

display neither exponential decay at inﬁnity nor a cusp at the origin. In order to compensate, a

relatively large basis set of GTFs is required. ETFs are more favorable than GTFs therefore, since

a smaller basis set can be used to obtain results that are just as accurate [7].

The most common ETFs are the Slater type f unctions (STFs) [8]. Despite the relative simplicity of

their analytical expression, however, the fact that multi-centre integrals over STFs can be extremely

problematic 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 by

Shavitt due to their remarkably simple Weierstrass transform [9]. B functions involve reduced Bessel

functions and, in fact, they can be expressed as linear combinations of STFs [10]. Compared to

other 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 G rotendorst et al in 1988 [18], is one of the most

successful approaches yet for the evaluation of multi-centre integrals. This particular method is

especially suitable for a basis of B functions. In particular, the Fourier transformation method

allows for the development of an analytical expression for three-centre nuclear attraction integrals.

The challenge, however, is that these analytical expressions involve semi-inﬁnite spherical Bessel

integrals, 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 a

spherical Bess el function rather than a simple trigonometric function. Moreover, as λ and v become

large, the zeros of the integrand become closer and the oscillations become stronger, thus further

complicating the numerical evaluation of the integral. It is possible to rewrite such integrals as

slowly convergent inﬁnite series of integrals of alternating sign. It may then seem appropriate to

apply series convergence acceleration techniq ues, 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 prevent

this approach from rendering accurate results.

As shown in previous work by Safouhi [21], through a reformalized integration by parts with respect

to xdx, the semi-inﬁnite spherical Bessel integral involved in the analytical expression of the three-

centre nuclear attraction integrals can be transformed into an integral involving the simple sine

function. This transformation, called the S transformation, results in a s ine integral that has

equidistant zeros, thus making it much more numerically favorable than the initial spherical Bessel

integral [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 the

asymptotic limit. While this method is highly accurate and eﬃcient, there are some ranges of

parameters 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 s implify the integrand to a sine function, a double exponential

transformation (or DE transformation) is applied. The DE transformations introduced by Ooura et

al in 1991 [23] and reﬁned by Ooura et al in 1999 [24], map the interval (0, ∞) onto (−∞, ∞) and are

such that the transformed integrand decays double exponentially at both inﬁnities. The trapezoidal

formula with an equal mesh size is then applied to the integral, and the resulting summation can

be truncated at s ome moderate positive and negative values. The ﬁnal result is a highly eﬃcient

quadrature formula in which relatively few function evaluations are required in order to obtain

highly accurate approximations. Both the original and the reﬁned DE transformations are applied

to the spherical Bessel integral involved in the three-centre molecular integral, and their eﬃciency

for this particular problem is compared in the discussion section. The numerical results illustrate

the high accuracy of these algorithms applied to three-center nuclear attraction integrals over B

functions with a miscellany of diﬀerent parameters.

2 General deﬁnitions and properties

The spherical Bessel function j

(x) of order λ ∈ N

is deﬁned by [25, 26]:

(x) = (−1)



xdx



(x) = (−1)



xdx





sin x



. (1)

The reduced Bessel function

(z) is deﬁned by [9, 27]:

(z) =

(z) = z

−z

j=0

(n + j)!

j! (n − j)!

(2 z)

, (2)

where n ∈ N, and where K

(z) is the modiﬁed Bessel function of the second kind [28].

The B functions are deﬁned by [10, 27]:

n,l

(ζ, ~r) =

(ζr)

n+l

(n + l)!

n−

(ζr)Y

(θ

, ϕ

), (3)

where n, l, m are the quantum numb ers and Y

(θ, ϕ) denotes the surface spherical harmonic, which

is deﬁned using the Condon-Shortley phase convention [29]:

(θ, ϕ) = i

m+|m|



(2l + 1)(l −|m|)!

4π(l +|m|)!



1/2

|m|

(cos θ)e

imϕ

, (4)

where P

(x) is the associated Legendre polynomial of the lth degree and mth order.

The Fourier transform

n,l

(ζ, ~p) of B

n,l

(ζ, ~r) (3) is given by [14]

n,l

2n+l−1

(−i|p|)

(ζ

+|p|

)

n+l+1

(θ

, ϕ

). (5)

The Gaunt coeﬃcients are deﬁned by [30–33]:





θ =0

2π

ϕ=0

(θ, ϕ)

∗

(θ, ϕ)Y

(θ, ϕ) sin(θ)dθdϕ. (6)

The Gaunt coeﬃcients linearize the product of two spherical harmonics:

(θ, ϕ)

∗

(θ, ϕ) =

l=l

min



, m

|l, m

− m



−m

(θ, ϕ), (7)

where the subscript l = l

min,2

implies that the summation index l runs in steps of 2 f rom l

min

+ l

, and where the constant l

min

is given by:

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:





R −

# —





∗

R −

# —

OC|



R −

# —



R, (9)

where A, B and C are three arbitrary points of the euclidian space E

, while O is the origin of the

ﬁxed coordinate sys tem. n

and n

stand for the principal quantum numb ers.

After performing a translation of vector

# —

OA, we can write the integral I

as follows:

(ζ

, ζ

) =

(ζ

, ~r)

∗



~r −



(ζ

, ~r −

)d~r, (10)

where ~r =

#—

R −

# —

OA,

# —

AC and

# —

AB.

3 Double exponential transformation for the computation of s emi-

inﬁnite integrals

For more details on the DE transformations and their applications, we refer the readers to the work

done by Ooura et al [23, 24]. The double exponential (DE) transformations present a particularly

eﬃcient method for computing semi-inﬁnite integrals of the form:

∞

(x) sin(vx)dx, (11)

where f

(x) is a slowly decaying function and v is a constant.

We will denote the semi-inﬁnite integral we wish to evaluate by:

I =

∞

(x)dx. (12)

Supp ose that f

is a function such that for some θ, λ ∈ R and for all n ∈ N, n ≥ n

, where n

some large natural number:

(λn + θ) = 0. (13)

In other words, after some large x value, f

(x) has inﬁnite equidistant zeros occurring every λ, with

a shift of θ with respect to the origin. Now provided that f

(x) is of the form of the integrand

in (11), it can easily be determined that λ = π v and θ = 0.

Now let φ(t) be a function such that:

φ(t) ∼ t as t → ∞ and φ(t) ∼ 0 as t → −∞. (14)

Let M be a large positive constant. Applying the variable transformation x = Mφ(t) to the integral

I, we get:

I =

−1

(∞)

−1

(0)

(Mφ(t))Mφ

′

(t)dt (15)

= M

∞

−∞

(Mφ(t))φ

′

(t)dt. (16)

Now applying the trap ezoidal rule with a mesh size of h and a shift of

with respect to the origin,

we obtain the approximation:

I ≈ Mh

∞

n=−∞

Mφ



nh +



′



nh +



. (17)

Let M and h be such that:

Mh = λ, (18)

and consider f



Mφ



nh +





for large n:

Mφ



nh +



∼ f



nh +



= f

(Mnh + θ)

= 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 the

summation in (17) at some negative N

−

∈ Z.

As a result, we obtain the approximation:

I ≈ M h

n=N

−

Mφ



nh +



′



nh +



, (20)

for some N

−

, N

∈ Z.

3.1 A DE transformation

One transformation that satisﬁes the conditions in (14) is given by [23]:

(t) =

1 − exp(−K sinh t)

, (21)

where K is s ome positive constant.

−1 −0.5 0 0.5 1 1.5 2

0.5

1.5

(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) ∼











0 as t → −∞

as t → 0

t as t → ∞.

(22)

Similarly, φ

′

(t) also encounters a singularity when t = 0. By the same method, it can be easily

found that:

′

(t) ∼











0 as t → −∞

as t → 0

1 as t → ∞.

(23)

3.2 A more robust DE transformation

A second transformation that satisﬁes the conditions in (14) is given by [24]:

(t) =

1 − exp(−2x − α(1 − e

−t

) − β(e

− 1))

(24)

where α, β are constants such that:

β = O(1), α = O((M log M )

−1/2

) and 0 ≤ α ≤ β ≤ 1. (25)

Some assignments that satisfy the constraints in (25 ) are given by [24]:

α =

1 + M log(1 + M)/4π

and β =

. (26)

−2 −1 0 1 2 3

1.5

2.5

α+β+2

(t)

Figure 2: The graph of φ

(t) given by (24), with the assignments for α and β in (26), where M = 4.5.

In a similar fashion to φ

(t), φ

(t) and φ

′

(t) both have singularities when t = 0. By Taylor

expansion, it can be determined that:

(t) ∼











0 as t → −∞

2 + α + β

as t → 0

t as t → ∞,

(27)

and that:

′

(t) ∼











0 as t → −∞

(α + β + 2)

+ (α − β)

2(α + β + 2)

as t → 0

1 as t → ∞.

(28)

4 The three-centre nuclear attraction integrals

Using the Fourier transform method, an analytical expressi on f or three-centre nuclear attraction

integrals over B functions (10) can be found [17, 18]:

8(4π)

(−1)

(2l

+ 1)!!(2l

+ 1)!!(n

+ l

+ n

+ l

+ 1)!ζ

−1

+ l

)!(n

+ l

′

=−l

′

(−1)

′



′

− l

′

− m

′



(2l

′

+ 1)!![2(l

−l

′

) + 1]!!

′

=−l

′

(−1)

′



′

− l

′

− m

′



(2l

′

+ 1)!![2(l

−l

′

) + 1]!!

′

l=l

′

min



′

|lm

′

− m

′



′

−m

′

(θ

, ϕ

)

−l

′

−l

′

λ=l

′′

min

(−i)



− l

′

− m

′

− l

′

− m

′

|λµ



∆l

j=0



∆l



(−1)

−j+1

+ n

+ l

− j + 1)!

s=0

−l

′

(1 − s)

−l

′

(θ

, ϕ

)

+∞

x=0



γ(s, x)





γ(s, x)



(vx)dx

ds, (29)

where the constant l

min

is deﬁned in (8) and:

γ(s, x) =

(1 − s)ζ

+ sζ

+ s(1 − s)x

~v = (1 − s)

−

, v = k~vk and R

= k

= l

− l

′

+ l

− l

′

= 2(n

+ l

+ n

+ l

) − (l

′

+ l

′

) − l + 1

ν = n

+ n

+ l

− l − j +

µ = (m

− m

′

) − (m

− m

′

)

∆l =



′

+ l

′

− l)/2



The semi-inﬁnite integral in (29), which will from now on be referred to as I(s), is deﬁned by:

I(s) =

+∞



γ(s, x)





γ(s, x)



(vx)dx. (30)

The challenge of evaluating (29) is mainly caused by the presence of the s emi-inﬁnite integral I(s)

(30), whose integrand is highly oscillatory due to the presence of the spherical Bessel function. This

is especially the case when λ and v are large. Notice also that when the value of s is close to 0 or 1,

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:



γ(s, x)





γ(s, x)



becomes constant, and the integrand can be reduced to x

(vx). As a result, the oscillations of

(vx) cannot be restrained, and the evaluation of I(s) becomes most laborious. In Figure 3, the

highly oscillatory nature of the integrand can be observed.

0 5 10 15 20 25 30

−4

−2

Integrand of I(s) (10

−2

)

Figure 3: The integrand of I(s) where s = 0.01, ν =

/2, n

= 7, n

= 2, λ = 2, R

= 9, ζ

= 2, R

= 3.5, and

= 1 (see row 9 of Table 1).

4.1 Application of the 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 twi ce continuously diﬀerentiable functions deﬁned

on the i nterval [0, ∞[. Now if g(x) is such that



xdx



λ−1

g(x)) for l = 0, 1, 2, . . . , λ are deﬁnied

and:

lim

x→0

l−λ+1



xdx





λ−1

g(x)



λ−l−1

(x) = 0

lim

x→∞

l−λ+1



xdx





λ−1

g(x)



λ−l−1

(x) = 0, (31)

holds true for all l = 0, 1, 2, . . . , λ − 1, then:

∞

f(x)dx =

∞



xdx



λ−1

g(x))

sin(x)dx.

In short, this transformation, called the S transformation, can simplify spherical Bessel integrals

into integrals involving the sine function. As such, the integral I(s) can be rewritten as:

I(s) =

λ+1

+∞







xdx



+λ−1



γ(s, x)





γ(s, x)







sin(vx)dx. (32)

2 4 6 8 10 12 14 16

−0.2

−0.1

0.1

0.2

Integrand of I(s) after S transformation

Figure 4: The integrand of I(s) in (32) after applying the S transformation, where s = 0.01, ν =

/2, n

= 7, n

= 2,

λ = 2, R

= 9, ζ

= 2, R

= 3.5, and ζ

= 1 (see row 9 of Table 1). The integrand prior to the transformation can

be 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) =



xdx



+λ−1



γ(s, x)





γ(s, x)



. (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 variable

transformation, such as φ

(t) (21) or φ

(t) (24).

Applying the DE transformation φ(t) to the equation f or I(s) in (32), we obtain the approximation:

I(s) =

λ+1

+∞

x=0

f(x) sin(vx)dx

≈

λ+1

n=N

−

Mφ



nh +



′



nh +



= I

(s). (34)

5 Numerical results and discussion

Table 1 lists values for I(s) obtained using the the approximation (34) with the transformations

(t) (21) and φ

(t) (24). In this ﬁrst table, we restrict the variable assignments to values that

can be handled by a MATLAB built-in numerical integration function that uses global adaptive

quadrature set to an accuracy of 15 correct digits. In Table 2, we restrict the variable assignm ents

to values that cannot be handled by the MATLAB integration function. For more problematic

variable assignments, particularly when λ and/or v are large, the MATLAB automatic integrator

fails to provide results to the same accuracy eﬃciently. We theref ore opted to use the same algorithm

in Python which, with the help of the symbolic computation package SymPy, was able to complete

the approximations signiﬁcantly faster. The relative errors corresponding to the approximations

using each transformation are listed in Table 2 as ε

and ε

, respectively.

As shown by Ooura et al [23, 24], the relative error for the approximation I

(s) using a DE

transformation can be written:

≈ exp



−



, (35)

where c is a constant. Ooura et al argue that we can assume that the relative error is approximately

given by [24]:

≈ exp



−A



= exp



−AM

v π



, (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.

5.1 Eﬀect of collocation points

For the variable assignments used in Figure 6 for the transformation φ

(t) (21), in our present

integrator, an optimal value for M was found to be M ≈ 54.25338. Over the course of our calcula-

tions for Figure 6, therefore, we assigned said value to M in order to preserve the behaviour of the

approximation that was conducted by the integrator.

In our integrator, we ﬁrst took the summation in (34) over positive n values, until the addition

of higher terms stopped signiﬁcantly contributing to the overall approximation. Similarly, we then

took the summation over negative n values, again truncating the summation once the contributions

of each term became less signiﬁcant than the predetermined error tolerance. Using this method, for

the 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 the

portion of the integral below zero is often higher than the number of collocation points needed to

approximate the portion of the integral ab ove zero. This is due to the asymmetry of the transformed

integrand which is represented by the summation in (34), that is to say of the integrand of I(s) =

∞

−∞



Mφ(t)



· Mφ

′

(t)dt. This asymmetric behavior can be observed in Figure 5.

−6 −4 −2 0 2 4 6

−0.5

0.5

Integrand of I(s) =

∞

−∞



Mφ(t)



· Mφ

′

(t)dt

Figure 5: The integrand of I(s) =

∞

−∞



Mφ(t)



· Mφ

′

(t)dt where s = 0.01, ν =

/2, n

= 5, n

= 0, λ = 0,

= 6.31, ζ

= 1, R

= 2.0, and ζ

= 1 (see row 2 of Table 1).

Now in order see a pattern in our data regarding the eﬀect of the collocation points on the accuracy

of the approximation, for our ﬁrst point on the graph, a moderate number of median collocation

points was taken in order to ensure s uﬃcient prox imity to the exact result. More precisely, for the

ﬁrst point on the graph in Figure 6, we took the summation in (34) over 70 collocation points, with

a lower bound of -62 and an upper bound of 7. From then on, the upper and lower bounds of the

summation were each extended by one. We proceeded to extend the bounds of the summation in

this manner – increasing the number of collocation points by two at a time until we reached 138

collocation points, which was the necessary amount of collocation points originally determined by

our integrator (see Table 1).

For the variable assignments used φ

(t) (24), our present integrator found M ≈ 21.70135 to be an

optimal value for M; we therefore as signed 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 was

done for transformation φ

(t), for our ﬁrst point of data, we took the summation over 45 collocation

points, with a lower bound of -37 and an upper bound of 7. From then on, the upper and lower

bounds of the summation were each extended by one until we reached a total of 97 collocation

points, 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 of

collocation points increase. Whether our integrator us es transformation φ

(t) (21) or transformation

(t) (24), the behaviour of the absolute error for the approximation of I(s) based on the number

of collocation points is strikingly similar. In fact, after a certain number of collocation points, both

graphs 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 ass ignments f or the evaluation of I(s); with the variable ass ignments

we 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 signiﬁcant in determining the eﬃciency

of our integrator. The consequences of M on our integrator’s eﬃciency will be addressed in the

following subsection.

70 80 90 100 110 120

n for the transformation φ

(t)

Absolute error (10

−3

)

50 60 70 80 90

n for the transformation φ

(t)

Figure 6: Absolute error dependence of I

(s) ( 34) on the total number of collocat ion points using transformations

(t) (21) and φ

(t) (24), and t he variable assignments from the second row of Table 1. The absolute error was found

by comparing the results of the approximation with the result obtained by a MATLAB b uilt-in numerical integration

function that uses global adaptive quadrature.

5.2 Eﬀect of M

In Figure 7, it can easily be seen that the absolute error converges to zero much more quick ly using

(t) (24). In other words, for a given set of variable assignments, a smaller mesh size is required

using transformation φ

(t) (21) to obtain results that are as accurate as the results using φ

(t) (24).

The challenge is to ﬁnd 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 excessively

large, as quite a number of collocation points in the summation (34) would be needed to compensate

for the unnecessarily small mesh size. Yet if M is too small, as shown in Figure 7, the error of the

approximation will not approach zero. This is what is meant by ﬁnding an “optimal value” for M.

Our present integrator ﬁnds an optimal value for M by ﬁrst letting M be a relatively small positive

constant, say M

, based on the relative error tolerance input by the user and according to the

following formula [24]:

−π

log



√



, (37)

and ε

denotes the relative error tolerance. For the transformation φ

(t) (21), we use A = 2, and for

the transformation φ

(t) (24), we use A = 5. For details regarding the selection of each subsequent

, we direct the reader to the algorithm in [24].

The integrator then calculates I

(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 result

due to the fact that the value assigned to M is too small to expect the error to converge to zero.

Rather, I

(s) serves as a value to which a second, more accurate approximation can be compared

for 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

. The

ensuing approximation, I

(s), often times meets the user’s relative error tolerance requirements,

at which p oint the integrator outputs the value of I

(s) as its ﬁnal result. Otherwise, the process

must be repeated with another value assigned to M, until a result of suﬃcient accuracy is found.

This could take up to a maximum of four diﬀerent assignments to M, each of which is larger than

the previous.

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

Absolute error (10

−5

)

(t)

Figure 7: Absolute error dependence of I

(s) ( 34) on the total number of collocat ion points using transformations

(t) (21) and φ

(t) (24), and t he variable assignments from the second row of Table 1. The absolute error was found

by comparing the results of the approximation with the result obtained by a MATLAB b uilt-in numerical integration

function that uses global adaptive quadrature.

Clearly, with this type of integrator, it is ideal to ﬁnd an op timal value for M on the s econd attempt,

and highly undesirable to not ﬁnd an optimal value until the fourth attempt, because the approxi-

mation (34) would need to be calculated four times over, meaning that numerous time-consuming

computations would have to be executed. It is for this reason that M also holds signiﬁcant weight

while comparing the eﬃciency of our algorithm using transformations φ

(t) (21) and φ

(t) (24). For

instance, by solely comparing the collocation points in Table 1, it may appear as though φ

(t) (24) is

without fail the more eﬃcient transformation to use. That being said, since the value n

, which we

are using to denote the total number of values that are assigned to M before successfully completing

the approximation for I(s), also gives valuable insight into the eﬃciency of each transformation. A

high n

value signiﬁes that our integrator had to complete the approximation (34) several times

before ﬁnding an optimal value for M and successfully approximating the integral, thus multiplying

the 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 ) is

truncated with 81 collocation points using transformation φ

(t) (21), and only 72 collocation points

using transformation φ

(t) (24), the values for n

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 our

present integrator is more eﬃcient in its approximation using transformation φ

(t) (21). Generally

speaking, however, the two-parameter transformation, φ

(t) (24) lends itself to a more eﬃcient

approximation of the semi-inﬁnite spherical Bessel integrals.

6 Conclusion

The presence of the semi-inﬁnite spherical Bessel integral in the analytical expression for three-centre

nuclear attraction integrals causes the rapid and accurate numerical evaluation of the integrals

to be extremely challenging. The S transformation, introduced by Safouhi, is able to transform

the spherical Bessel integrals into considerably simpler sine integrals. The zeros of the integrand

become eq uidistant, which signiﬁcantly increases our ability to apply various extrapolation and

transformation methods. Nevertheless, due to its slow convergence, the semi-inﬁnite sine integrals

still pose substantial computational diﬃculties.

In the present work, we presented an eﬃcient way of treating the arduous sine integrals that result

from the S transformation. Utilizing the slowly decaying, highly oscillatory nature of the integrands,

we applied two distinct double exponential transformations. These transformations each rendered

highly eﬃcient quadrature formulae that spanned from negative to positive inﬁnity. That being

said, we were able to truncate the summations at relatively small positive and negative integers.

It was f ound that remarkably few collocation points were required to obtain suﬃciently accurate

results. This was especially the case for the second transformation. In our numerical section, it can

be seen that by using the second transformation, the semi-inﬁnite spherical Bessel integral can be

approximated with a high predetermined accuracy in under 100 collocation points. The numerical

tests served to further verify the precision and demonstrate the signiﬁcance of the method.

7 Numerical tables

Table 1: Computation of I (s).

s ν n

λ R

I(s)

Matlab

max

0.99

/2 1 0 0 24.00 1.5 2.0 1 .113 874 170 637 205( 0) 142 45 2 .1(-16) 93 34 2 .1(-16)

0.01

/2 5 0 0 6.31 1.0 2.0 1 .638 243 453 884 443( 0) 138 41 2 .1(-16) 97 33 2 .1(-16)

0.99

/2 5 0 0 4.50 2.0 1.5 1 .701 581 269 512 308( 0) 139 42 2 .1(-16) 97 33 2 .1(-16)

0.99

/2 5 1 0 3.00 1.5 3.5 2 .242 778 918 544 382(-3) 78 41 2 .1(-16) 78 37 4 .1(-16)

0.99

/2 9 1 1 6.00 2.0 3.5 1 .183 138 910 224 195( 1) 139 42 2 .1(-16) 97 33 2 .1(-16)

0.01

/2 9 2 1 8.50 2.0 3.5 2 .248 336 723 989 982(-3) 81 44 2 .1(-16) 72 35 4 .5(-18)

0.01

/2 3 2 1 3.00 2.0 5.0 1 .285 091 100 421 789(-2) 85 41 3 .1(-16) 81 37 3 .1(-16)

0.99

/2 5 2 2 4.00 2.5 5.5 1 .112 567 767 257 153( 0) 139 42 2 .1(-16) 92 33 2 .1(-16)

0.01

/2 7 2 2 9.00 2.0 3.5 1 .183 269 571 025 263(-2) 141 44 2 .1(-16) 97 33 2 .1(-16)

0.01

/2 9 3 2 4.00 2.5 5.5 1 .167 566 737 865 368(-1) 91 37 3 .1(-16) 92 36 3 .1(-16)

• I(s)

Matlab

are obtained by m eans of a MATLAB built-in numerical integration function that uses global adaptive quadrature.

• ε

stand for the relative errors obtained by using the approximation (34) with the t ransformation φ

(t) (21).

• ε

stand for the relative errors obtained by using the approximation (34) with the t ransformation φ

(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

and n

represent the total number of values that are assigned to M before successfully completing th e 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 compu tation package, SymPy

• The error tolerance of ε = 10

−15

was used in our calculation.

• Calculations were completed using IEEE 754 double precision

• Numbers in parentheses represent powers of 10.

Table 2: Computation of I (s).

s ν n

λ R

I(s)

max

I(s)

max

0.99

/2 9 2 1 35 2.5 3.5 0.5 .737 829 982 455 112( 0) 82 45 2 .1(-16) .737 829 982 455 148( 0) 84 41 4 .1(-16)

0.01

/2 3 2 2 45 1.5 2.0 1.0 .103 851 188 766 556(-2) 142 45 2 .1(-16) .103 851 188 766 626(-2) 93 34 2 .1(-16)

0.99

/2 5 3 3 50 1.5 2.0 1.0 .317 648 310 603 089(-1) 141 44 2 .1(-16) .317 648 310 603 304(-1) 93 34 2 .1(-16)

0.01

/2 6 4 3 55 1.5 2.0 1.0 .989 261 101 060 361(-3) 83 46 2 .1(-16) .989 261 101 060 357(-3) 86 42 4 .5(-20)

0.99

/2 6 4 4 55 1.5 2.0 1.0 .366 551 897 499 993(-1) 141 44 2 .1(-16) .366 551 897 500 241(-1) 93 34 2 .5(-16)

0.01

/2 9 5 4 55 1.5 2.0 1.0 .698 018 626 122 216(-3) 83 46 2 .1(-16) .698 018 626 122 213(-3) 82 40 4 .2(-19)

0.99

/2 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

/2 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

/2 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

/2 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

/2 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 m eans of a MATLAB built-in numerical integration function that uses global adaptive quadrature.

• ε

stand for the relative errors obtained by using the approximation (34) with the t ransformation φ

(t) (21).

• ε

stand for the relative errors obtained by using the approximation (34) with the t ransformation φ

(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

and n

represent the total number of values that are assigned to M before successfully completing th e 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 compu tation package, SymPy

• The error tolerance of ε = 10

−15

was used in our calculation.

• Calculations were completed using IEEE 754 double precision

• Numbers in parentheses represent powers of 10.

References

[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. Graﬃ

(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 states

of any molecular system. Proc. R. Soc. Lond. Series A, M ath. & Phys. Sciences., 200:542–554,

1950.

[5] S.F. Boys. Electronic wave functions. II. a calculation for the ground state of the Beryllium

atom. 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-electron

wave 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 ﬁrst 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. Shav itt. The Gaussian f unction 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 integrals

and 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 Bes sel functions

of 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 functions

and 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-electron

repulsion integrals involv ing Slater functions. J. Chem. Phys., 40:3083–3086, 1964.

[17] H.P. Trivedi and E.O. Steinborn. Fourier transform of a two-center product of exponential-type

orbitals. 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-electron

multicenter 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

) 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 s ine, spherical Bessel and reduced Bessel functions f or improving

convergence of semi-inﬁnite very oscillatory integrals: The evaluation of three-center nuclear

attraction integrals over B functions. J. Phys. A: Math. Gen., 34:2801–2818, 2001.

[22] L. Berlu and H. Safouhi. An extremely eﬃcient and rapid algorithm for a numerical evaluation

of 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 the

inﬁnite 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, Fifth

edition, 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 ﬁelds represented by spherical-harmonics ex-

pansions for molecular calculations. III. Translations of reduced Bessel functions, Slater-typ e

s-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 coeﬃcients of real spheri-

cal harmonics and their relation to Gaunt coeﬃcients. 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. Eﬃcient evaluation of vector translation coeﬃcients in multiparticle light-s cattering

theories. J. Comput. Phys., 139:137–165, 1998.