Calculation of STOs electron repulsion integrals by ellipsoidal expansion and large-order approximations

Abstract

For general two-electron two-centre integrals over Slater-type orbitals (STOs), the use of the Neumann expansion for the Coulomb interaction potential yields infinite series in terms of few basic functions. In many important cases the number of terms necessary to achieve convergence by a straightforward summation is large and one is forced to calculate the basic integrals of high order. We present a systematic approach to calculation of the higher-order terms in the Neumann series by large-order expansions of the basic integrals. The final expressions are shown to be transparent and straightforward to implement, and all auxiliary quantities can be calculated analytically. Moreover, numerical stability and computational efficiency are also discussed. Results of the present work can be used to speed up calculations of the STOs integral files, but also to study convergence of the Neumann expansion and develop appropriate convergence accelerators.

Appendices

Appendix 1: Auxiliary integrals

Virtually all final working formulae obtained in the present paper are given in terms of the basic integrals $E_n(z)$ and $a_n(z)$. They are defined through the integral representations

$$\begin{aligned} E_n(z)&= \int _1^\infty dt\, \frac{e^{-z t}}{t^n}, \end{aligned}$$

(61)

$$\begin{aligned} a_n(z)&= \int _0^1 dt\, t^n e^{-z t}, \end{aligned}$$

(62)

where n is an arbitrary integer in the former and a nonnegative integer in the latter. Calculation of $a_n$ is most easily carried out with help of the Miller algorithm [49] as discussed by Harris [22]. The integral $E_n$ is usually called the generalised exponential integral. Computation of $E_n$ differs depending on the sign of n. For a negative integer n the following recursion is completely stable in the upward direction

$$\begin{aligned} E_n(z) = -\frac{n}{z} E_{n-1}(z)+\frac{e^{-z}}{z}. \end{aligned}$$

(63)

For positive n and $z<1$ one uses the series expansion

$$\begin{aligned} E_n(z) = \frac{(-z)^{n-1}}{(n-1)!}\left[ \varPsi (n)-\log z\right] -\sum _{\begin{array}{c} k=0 \\ k\ne n-1 \end{array}}^\infty \frac{(-z)^k}{k!\,(1-n+k)}, \end{aligned}$$

(64)

where $\varPsi (n)$ is the digamma function at integer argument. The above infinite summations converge to the machine precision in, at most, few tens of terms. Finally, for positive n and $z>1$ the continued fraction (CF) formula can be applied

$$\begin{aligned} E_n(z) = e^{-z} \left( \frac{1}{z+}\frac{p}{1+}\frac{1}{z+}\frac{p+1}{1+}\frac{2}{z+}\cdots \right) . \end{aligned}$$

(65)

To evaluate the CF one can use the Lentz algorithm [43]. The only inconvenience is that consecutive numerators and denominators in the Lentz scheme grow very quickly with the number of terms retained in Eq. (65). Therefore, it is necessary to rescale them from time to time by a small number to avoid numerical overflows. Let us also recall the leading terms of the large-order asymptotic expansions for $E_n(z)$ and $a_n(z)$ which read

$$\begin{aligned} E_n(z)&= \frac{e^{-z}}{n+z} + {\mathcal {O}}\left( \frac{1}{n^2}\right) , \end{aligned}$$

(66)

$$\begin{aligned} a_n(z)&= \frac{e^{-z}}{n-z} + {\mathcal {O}}\left( \frac{1}{n^2}\right) , \end{aligned}$$

(67)

where $n>z$.

Appendix 2: Large-order asymptotic formulae for $i_\mu (\alpha )$ and $k_\mu (\alpha )$

According to the work of Sidi et al. [47, 48] the modified spherical Bessel functions posses the following large-order expansions

$$\begin{aligned} i_\mu (z)&= \frac{\sqrt{\pi }}{2}\frac{(z/2)^\mu }{\varGamma (\mu +\frac{3}{2})}\sum _{m=0}^\infty \frac{b_m(z)}{(\mu +1/2)^m}, \end{aligned}$$

(68)

$$\begin{aligned} k_\mu (z)&= \frac{\sqrt{\pi }}{4}\frac{\varGamma (\mu +\frac{1}{2})}{(z/2)^{\mu +1}}\sum _{m=0}^\infty \frac{b_m(z)}{(\mu +1/2)^m}(-1)^m, \end{aligned}$$

(69)

as $\mu \rightarrow \infty$ at a fixed z, where

$$\begin{aligned} b_m(z)=\sum _{k=1}^m (-1)^{m-k} \frac{S_{mk}}{k!} (z/2)^{2k}, \end{aligned}$$

(70)

for $m>0$ and $b_0(z)=1$. The quantities $S_{mk}$ in the above expression are the Stirling numbers of the second kind [43] defined recursively as

$$\begin{aligned} S_{m0}&=\delta _{m0},\;\;\; S_{m1}=1,\;\;\; S_{mm}=1, \end{aligned}$$

(71)

$$\begin{aligned} S_{mk}&=S_{m-1,k-1}+k S_{m-1,k}, \end{aligned}$$

(72)

and we additionally adapt the convention $S_{mk}=0$ for $m<k$ or $m<0$.