Compact Formulae for Three-Center Nuclear Attraction Integrals Over Exponential Type Functions

TL;DR

This paper presents new analytical formulas for three-center nuclear attraction integrals over exponential type functions, significantly improving computational efficiency in molecular orbital calculations.

Contribution

It introduces compact analytical expressions and an efficient evaluation method for three-center nuclear attraction integrals, enhancing accuracy and speed in quantum chemistry computations.

Findings

Substantial efficiency gains over existing methods

Effective use of rational minimax approximants

Accurate double precision evaluations

Abstract

In any ab initio molecular orbital (MO) calculations, the major task involves the computation of the so-called molecular multi-center integrals. Multi-center integral calculations is a very challenging mathematical problem in nature. Quantum mechanics only determines which integrals we evaluate, but the techniques employed for their evaluations are entirely mathematical. The three-center nuclear attraction integrals occur in a very large number even for small molecules and are among of the most difficult molecular integrals to compute efficiently to a high pre-determined accuracy. In the present contribution, we report analytical expressions for the three-center nuclear attraction integrals over exponential type functions. We describe how to compute the formula to obtain an efficient evaluation in double precision arithmetic. This requires the rational minimax approximants that minimize the maximum error on the interval of evaluation. The numerical tests show a substantial gain in efficiency over the state-of-the-art.

Figures

Tables

Table 1. Table 1: Evaluation of ℐ ​ ( s ) ℐ 𝑠 {\cal I}(s)  ( 30 ) involved in ( 28 ). s = 0.05 𝑠 0.05 s=0.05 , ζ 1 = 1.5 subscript 𝜁 1 1.5 \zeta_{1}=1.5 , ζ 2 = 2.0 subscript 𝜁 2 2.0 \zeta_{2}=2.0 , R 1 = 2.0 subscript 𝑅 1 2.0 R_{1}=2.0 , R 2 = 5.5 subscript 𝑅 2 5.5 R_{2}=5.5 and v = 3.225 𝑣 3.225 v=3.225

$\nu -\frac{1}{2}$	$n_{\gamma }$	$n_x$	$\lambda$	$ℐ{\left(s\right)}^{†}$	$ℐ{\left(s\right)}^{\text{Matlab}}$
7	11	3	1	0.151 189 722 612 165 (-1)	0.151 189 722 613 836 (-1)
7	11	4	0	-0.770 700 245 226 897 (-1)	-0.770 700 245 226 550 (-1)
7	11	4	2	0.911 341 847 656 817 (-1)	0.911 341 847 657 070 (-1)
6	13	3	1	0.862 532 316 505 739 (-3)	0.862 532 316 513 426 (-3)
6	13	4	0	-0.412 000 772 378 986 (-2)	-0.412 000 772 377 796 (-2)
6	13	4	2	0.492 236 336 705 101 (-2)	0.492 236 336 704 074 (-2)
7	15	3	1	0.106 814 986 690 961 (-1)	0.106 814 986 691 068 (-1)
9	17	3	1	0.266 323 983 838 913 ( 1)	0.266 323 983 839 022 ( 1)
9	19	3	1	0.184 547 358 116 701 ( 1)	0.184 547 358 116 765 ( 1)
9	19	4	0	-0.762 928 846 920 085 ( 1)	-0.762 928 846 919 988 ( 1)
9	19	5	1	-0.347 485 191 071 318 ( 2)	-0.347 485 191 071 172 ( 2)
10	21	3	1	0.257 290 058 890 616 ( 2)	0.257 290 058 890 630 ( 2)
10	21	4	0	-0.101 363 984 175 823 ( 3)	-0.101 363 984 175 821 ( 3)
10	21	4	2	0.125 297 943 142 392 ( 3)	0.125 297 943 142 393 ( 3)
10	21	5	1	-0.440 284 382 122 921 ( 3)	-0.440 284 382 122 792 ( 3)
10	19	4	0	-0.163 020 781 236 490 ( 3)	-0.163 020 781 236 467 ( 3)
10	19	5	1	-0.746 484 060 054 242 ( 3)	-0.746 484 060 053 884 ( 3)
11	23	3	1	0.372 866 539 760 214 ( 3)	0.372 866 539 760 223 ( 3)
11	21	4	0	-0.235 942 964 388 561 ( 4)	-0.235 942 964 388 555 ( 4)
11	23	4	2	0.174 686 408 404 681 ( 4)	0.174 686 408 404 684 ( 4)
11	23	5	1	-0.579 875 446 713 486 ( 4)	-0.579 875 446 713 279 ( 4)
11	23	5	3	0.850 707 087 650 976 ( 4)	0.850 707 087 650 815 ( 4)

Table 2. Table 2: Evaluation of ℐ ​ ( s ) ℐ 𝑠 {\cal I}(s)  ( 30 ) involved in ( 28 ). s = 0.99 𝑠 0.99 s=0.99 , ζ 1 = 1.5 subscript 𝜁 1 1.5 \zeta_{1}=1.5 , ζ 2 = 2.0 subscript 𝜁 2 2.0 \zeta_{2}=2.0 , R 1 = 2.0 subscript 𝑅 1 2.0 R_{1}=2.0 , R 2 = 5.5 subscript 𝑅 2 5.5 R_{2}=5.5 and v = 1.9450 𝑣 1.9450 v=1.9450

$\nu -\frac{1}{2}$	$n_{\gamma }$	$n_x$	$\lambda$	$ℐ{\left(s\right)}^{†}$	$ℐ{\left(s\right)}^{\text{Matlab}}$
8	9	5	1	-0.136 578 999 110 210 (-3)	-0.136 578 513 775 021 (-3)
8	9	5	3	0.157 429 717 614 749 (-3)	0.157 428 199 599 963 (-3)
9	11	4	0	-0.133 767 585 018 979 (-3)	-0.133 767 051 011 091 (-3)
9	11	5	1	-0.238 345 682 495 741 (-2)	-0.238 345 138 495 788 (-2)
9	11	5	3	0.275 872 683 252 550 (-2)	0.275 872 094 789 520 (-2)
9	11	4	2	0.145 980 032 943 987 (-3)	0.145 979 955 050 279 (-3)
9	13	5	1	-0.189 782 159 585 373 (-2)	-0.189 782 160 850 882 (-2)
9	17	3	1	0.331 772 864 261 456 (-5)	0.331 772 945 380 560 (-5)
9	19	3	1	0.201 688 439 270 122 (-5)	0.201 688 456 114 013 (-5)
9	19	4	0	-0.301 438 469 081 214 (-4)	-0.301 438 463 736 758 (-4)
9	19	5	1	-0.476 213 479 931 115 (-3)	-0.476 213 471 529 263 (-3)
10	21	3	1	0.246 169 752 226 837 (-4)	0.246 169 761 624 060 (-4)
10	21	4	0	-0.358 276 851 579 279 (-3)	-0.358 276 848 288 597 (-3)
10	21	4	2	0.396 246 479 172 108 (-3)	0.396 246 476 283 224 (-3)
10	21	5	1	-0.551 474 883 342 359 (-2)	-0.551 474 878 125 191 (-2)
10	19	4	0	-0.651 321 577 489 068 (-3)	-0.651 321 566 692 786 (-3)
10	19	5	1	-0.103 223 435 438 069 (-1)	-0.103 223 433 320 636 (-1)
11	23	3	1	0.307 741 810 026 843 (-3)	0.307 741 813 898 588 (-3)
11	21	4	0	-0.843 566 525 213 554 (-2)	-0.843 566 520 507 011 (-2)
11	23	4	2	0.483 636 533 415 352 (-2)	0.483 636 533 188 303 (-2)
11	23	5	1	-0.654 156 086 743 768 (-1)	-0.654 156 083 528 505 (-1)
11	23	5	3	0.778 484 244 434 089 (-1)	0.778 484 240 745 456 (-1)

Equations

$$B_{n,l}^m(\zeta ,\text{vv}r)=\frac{(\zeta r{)}^l}{2^{n+l}(n+l)!}{\hat{k}}_{n-\frac{1}{2}}(\zeta r)Y_l^m({\theta }_{\text{vv}r},{\phi }_{\text{vv}r}),$$

$$Y_l^m(\theta ,\varphi )={\mathrm{i}}^{m+|m|}\sqrt{\frac{(2l+1)(l-|m|)!}{4\pi (l+|m|)!}}{P}_l^{|m|}(\cos \theta )e^{\mathrm{i}m\varphi },$$

$$P_l^m(x)=(1-x^2{)}^{m/2}\frac{{\mathrm{d}}^{l+m}}{\mathrm{d}{x}^{l+m}}\left(\frac{(x^2-1{)}^l}{2^{l}l!}\right).$$

$${\hat{k}}_{n+\frac{1}{2}}(z)=2^n(1/2{)}_n{e}^{-z}{}_1{F}_1(-n;-2n;2z),$$

$$(x{)}_n=\frac{\Gamma (x+n)}{\Gamma (x)}.$$

$${\hat{k}}_{\nu }(z)=\sqrt{\frac{2}{\pi }}{z}^{\nu }{K}_{\nu }(z),$$

$${\hat{k}}_{\nu +1}(z)$$

$$\frac{K_{\nu }(xz)}{z^{\nu }}=\frac{1}{2}{\int }_0^{\infty }\frac{e^{-\frac{x}{2}(t+\frac{z^2}{t})}}{t^{\nu +1}}\mathrm{d}t$$

$${\left(\frac{\mathrm{d}}{z\mathrm{d}z}\right)}^m\left(z^{\pm \nu }{K}_{\nu }(z)\right)$$

$${\chi }_{n,l}^m(\zeta ,\text{vv}r)=\mathcal{N}(\zeta ,n)r^{n-1}{e}^{-\zeta r}{Y}_l^m({\theta }_{\text{vv}r},{\phi }_{\text{vv}r}),$$

$$\mathcal{N}(\zeta ,n)=\sqrt{\frac{(2\zeta {)}^{2n+1}}{(2n)!}}.$$

$${\chi }_{n,l}^m(\zeta ,\text{vv}r)=\frac{\mathcal{N}(\zeta ,n)}{{\zeta }^{n-1}}\sum _{p=\tilde{p}}^{n-l}\frac{(-1{)}^{n-l-p}{2}^{2p+2l-n}(l+p)!}{(2p-n+l)!(n-l-p)!}{B}_{p,l}^m(\zeta ,\text{vv}r),$$

$$\tilde{p}=\left\{\begin{array}[]{llll}\displaystyle\frac{n-l}{2}&\textrm{if}&\quad n-l&\textrm{is even}\\ \displaystyle\frac{n-l+1}{2}&\textrm{if}&\quad n-l&\textrm{is odd}.\end{array}\right.$$

$$\bar{f}(\vec{k})=(2\pi {)}^{-3/2}{\int }_{\vec{r}}{e}^{-i\vec{k}\cdot \vec{r}}f(\vec{r})d\vec{r}\text{and}f(\vec{r})=(2\pi {)}^{-3/2}{\int }_{\vec{k}}{e}^{\text{leavevmode}i\vec{k}\cdot \vec{r}}\bar{f}(\vec{k})d\vec{k}.$$

$$e^{\pm \mathrm{i}\vec{p}\cdot \vec{r}}=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}4\pi (\pm \mathrm{i}{)}^l{j}_l(|\vec{p}||\vec{r}|)Y_l^m({\theta }_{\vec{r}},{\varphi }_{\vec{r}})[Y_l^m({\theta }_{\vec{p}},{\varphi }_{\vec{p}}){]}^{\ast },$$

$$j_n(x)=(-x{)}^n{\left(\frac{\mathrm{d}}{x\mathrm{d}x}\right)}^n\left(\frac{\sin (x)}{x}\right).$$

$$j_n(x)=\sqrt{\frac{\pi }{2x}}{J}_{n+1/2}(x).$$

$$j_{n+1}(x)=\frac{2n+1}{x}{j}_n(x)-j_{n-1}(x)\text{and}{J}_{\eta +1}(x)=\frac{2\eta }{x}{J}_{\eta }(x)-J_{\eta -1}(x).$$

$$J_{\eta }(z)$$

$${\left(\frac{\mathrm{d}}{z\mathrm{d}z}\right)}^m\left(z^{\eta }{J}_{\eta }(z)\right)$$

$${\bar{B}}_{n,l}^m(\zeta ,\vec{p})=\sqrt{\frac{2}{\pi }}{\zeta }^{2n+l-1}\frac{(-\mathrm{i}p{)}^l}{({\zeta }^2+p^2{)}^{n+l+1}}{Y}_l^m({\theta }_{\vec{p}},{\varphi }_{\vec{p}}).$$

$$\frac{1}{|\vec{r}-{\vec{R}}_1|}=\frac{1}{2{\pi }^2}{\int }_{\vec{p}}\frac{e^{-\mathrm{i}\vec{p}\cdot (\vec{r}-{\vec{R}}_1)}}{p^2}{\mathrm{d}}^3\vec{p}.$$

$$⟨l_1{m}_1|l_2{m}_2|l_3{m}_3⟩={\int }_0^{2\pi }{\int }_0^{\pi }[Y_{l_1}^{m_1}(\theta ,\varphi ){]}^{\ast }{Y}_{l_2}^{m_2}(\theta ,\varphi )Y_{l_3}^{m_3}(\theta ,\varphi )\sin \theta \mathrm{d}\theta \mathrm{d}\varphi .$$

$$[Y_{l_1}^{m_1}(\theta ,\varphi ){]}^{\ast }{Y}_{l_2}^{m_2}(\theta ,\varphi )=\sum _{l=l_{\min },2}^{l_1+l_2}⟨l_2{m}_2|l_1{m}_1|l{m}_2-m_1⟩Y_l^{m_2-m_1}(\theta ,\varphi ),$$

$$l_{\rm min}=\left\{\begin{array}[]{lccr}\max(|l_{1}-l_{2}|,|m_{2}-m_{1}|)&{\rm if}&l_{1}+l_{2}+\max(|l_{1}-l_{2}|,|m_{2}-m_{1}|)&{\rm is\leavevmode\nobreak\ even}\\ \max(|l_{1}-l_{2}|,|m_{2}-m_{1}|)+1&{\rm if}&l_{1}+l_{2}+\max(|l_{1}-l_{2}|,|m_{2}-m_{1}|)&{\rm is\leavevmode\nobreak\ odd.}\\ \end{array}\right.$$

$${\mathcal{I}}_{n_1,l_1,m_1}^{n_2,l_2,m_2}={\int }_{\vec{R}}{\left[B_{n_1,l_1}^{m_1}({\zeta }_1,\vec{R}-\overrightarrow{OA})\right]}^{\ast }\frac{1}{|\vec{R}-\overrightarrow{OC}|}{B}_{n_2,l_2}^{m_2}({\zeta }_2,\vec{R}-\overrightarrow{OB}){\mathrm{d}}^3\vec{R},$$

$${\mathcal{I}}_{n_1,l_1,m_1}^{n_2,l_2,m_2}=\frac{1}{2{\pi }^2}{\int }_{\vec{x}}\frac{e^{\mathrm{i}\vec{x}\cdot {\vec{R}}_1}}{x^2}{⟨B_{n_1,l_1}^{m_1}({\zeta }_1,\vec{r})\left|e^{-\mathrm{i}\vec{x}\cdot \vec{r}}\right|B_{n_2,l_2}^{m_2}({\zeta }_2,\vec{r}-{\vec{R}}_2)⟩}_{\vec{r}}{\mathrm{d}}^3\vec{x},$$

$${⟨B_{n_1,l_1}^{m_1}({\zeta }_1,\vec{r})\left|e^{-\mathrm{i}\vec{x}\cdot \vec{r}}\right|B_{n_2,l_2}^{m_2}({\zeta }_2,\vec{r}-{\vec{R}}_2)⟩}_{\vec{r}}={\int }_{\vec{r}}{\left[B_{n_1,l_1}^{m_1}({\zeta }_1,\vec{r})\right]}^{\ast }{e}^{-\mathrm{i}\vec{x}\cdot \vec{r}}{B}_{n_2,l_2}^{m_2}({\zeta }_2,\vec{r}-{\vec{R}}_2){\mathrm{d}}^3\vec{r},$$

$${\mathcal{I}}_{n_1,l_1,m_1}^{n_2,l_2,m_2}$$

$$[{\gamma }_{12}(s,x){]}^2$$

$$n_{\gamma }$$

Full text

Compact Formulae for Three-Center Nuclear Attraction Integrals Over Exponential Type Functions

Richard M. Slevinsky*§* and Hassan Safouhi*†*111The corresponding author (HS) acknowledges the financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) - Grant RGPIN-2016-04317

*§**Department of Mathematics, University of Manitoba

186 Dysart Rd, Winnipeg (MB), Canada R3T 2M8

Richard.Slevinsky@umanitoba.ca

†Campus Saint-Jean, University of Alberta

8406, 91 Street, Edmonton (AB), Canada T6C 4G9

hsafouhi@ualberta.ca*

Abstract.   In any ab initio molecular orbital (MO) calculations, the major task involves the computation of the so-called molecular multi-center integrals. Multi-center integral calculations is a very challenging mathematical problem in nature. Quantum mechanics only determines which integrals we evaluate, but the techniques employed for their evaluations are entirely mathematical. The three-center nuclear attraction integrals occur in a very large number even for small molecules and are among of the most difficult molecular integrals to compute efficiently.

In the present contribution, we report analytical expressions for the three-center nuclear attraction integrals over exponential type functions. We describe how to compute the formula to obtain an efficient evaluation in double precision arithmetic. This requires the rational minimax approximants that minimize the maximum error on the interval of evaluation.

Keywords:

Three-center nuclear attraction integrals. Molecular integrals. Exponential type functions. Bessel functions. Molecular structure calculations.

1 Introduction

Molecular electronic structure theory is a highly interdisciplinary research topic whose progress depends crucially on the progress in numerical analysis, mainly numerical integration, and only a very few mathematicians were interested in this challenging mathematical problem. In this regards, Sack has played a prominent role with his work in the 1960s on addition theorems [1, 2, 3, 4]. Hellmann’s book [5] has a section dealing with the evaluation of these extremely difficult integrals. The pursuit of accurate and efficient algorithms for the numerical evaluation of molecular integrals for the purpose of electronic structure calculations has led to a substantial body of work. While molecular integrals with Gaussian type functions (GTFs) [6, 7] as a basis have been for a long time the most easily calculated, their theoretical deficiency both in the long and short ranges [8, 9] has led to a renewed interest in an exponential type function (ETF) basis. ETFs are better suited than GTFs to represent electron wave functions near the nucleus and at long range, provided that multi-center integrals using such functions could be computed efficiently. Examples of such a basis include the Slater type functions (STFs) [10], the Laguerre type functions, and the $B$ functions [11, 12, 13]. We note that many researchers hope that the next generation of ab initio programs will be based on the usage of ETFs. Indeed much effort is being made to develop efficient molecular algorithms for integrals over conventional ETFs (STFs or $B$ functions) [14, 15, 16, 17, 18, 19] and references therein.

The $B$ function basis has been at the forefront of recent developments [20, 21, 22, 23, 24, 25], as their simple Fourier transform leads to compact analytical expressions for the integrals. Unfortunately, through the expansions, some of the resulting semi-infinite integrals have integrands composed of the Bessel functions. As the most complicated part of the analytical expressions, the semi-infinite integrals have become the bottleneck of the calculations in the $B$ function approach. Since the identification of this computational problem, a large body of work is devoted to dealing with this bottleneck. The most common element to most algorithms [26, 27, 28, 29, 30, 31] involves an intense integration-then-extrapolation approach, which is usually characterized by subdividing the integral between the oscillatory Bessel function’s zeros, by integrating by a quadrature, and using the resulting sequence to estimate the remainder of the integral. Of course, many improvements to the general extrapolation procedure are documented and indeed some of the most recent of these improvements are considered state-of-the-art for multi-center integrals. While these methods and their refinements are generally highly accurate and efficient, there are some ranges of parameters where either failure is inevitable or the computation becomes extremely heavy.

In the comparative study [32] of the most popular extrapolation methods and sequence transformations for computing semi-infinite integrals, the authors conclude that having asymptotic series representations for integrals and applying sequence transformations to accelerate their convergence or to sum their divergence leads to the most efficient algorithms for computing the integrals. However, when such asymptotic series representations do not exist, refinements to either the numerical steepest descent method [33, 34, 35, 36, 37, 38] or the extrapolation methods [39, 40, 41, 42, 43] must be made to obtain a desirable outcome. This conclusion does not particularly challenge any preconceived notions. However, it does emphasize that the ultimate goal in computing semi-infinite integrals is to find analytical expressions. Several examples of this approach have been documented in the literature on molecular integrals in the $B$ function basis [22, 44]. The expressions that are obtained have greatly simplified the calculation of one- and two-center integrals. The pursuit of such analytical expressions stopped in the 1990s at the three-center integrals because of the increased complexity of the integrands.

In this work, we report a significant breakthrough for the three-center nuclear attraction integrals. This breakthrough takes the form of an analytical expression for the semi-infinite integrals. In section 2, we introduce the general definitions and properties required for the discussion. In section 3, we prove our main result. In section 4, we hold a numerical discussion on its performance and compare it to existing methods based on extrapolation. Our comparison leads us to believe that the analytical expression is useful in a numerical setting and performs extremely well. Indeed, the decrease in calculation time is substantial compared with the existing methods.

Numerical experiments conducted for the integrals in (36), show that the use of the analytic expressions is $10^{3}$ times more efficient than the state-of-the-art [S7].

Since the nuclear magnetic resonance (NMR) integrals, such as first order integrals in the relativistic calculations of the nuclear shielding tensor [45, 46], are also expressed in terms of semi-infinite integrals similar to three-center nuclear attraction integrals, the proposed analytic development will also have a major impact in computing NMR properties of molecular systems. In addition, on the algorithmic side, existing density functional theory based programs, such as ADF (Amsterdam Density Functional) will benefit from the results of this research, since even in such an approach, integrals similar to three-center nuclear attraction are needed. The proposed methods could also be implemented in the existing frameworks, such as Slater type orbital package (STOP) and SMILE (Slater molecular integrals for large electronic systems) where STFs are expanded in terms of GTFs.

2 General definitions and properties

The functions $B_{n,l}^{m}(\zeta,\vv{r})$ are defined by [12, 13]:

[TABLE]

where $n,l,$ and $m$ are the quantum numbers such that $n=1,2,\ldots,$ $l=0,1,\ldots,n-1,$ and $m=-l,-l+1,\ldots l-1,l$ and where $Y_{l}^{m}(\theta,\phi)$ denotes the surface spherical harmonic, and is defined explicitly using the Condon-Shortley phase convention as [47]:

[TABLE]

where $P_{l}^{m}(x)$ is the associated Legendre polynomial of $l^{\rm th}$ degree and $m^{\rm th}$ order:

[TABLE]

The radial part of the $B$ functions is given by the reduced Bessel functions $\hat{k}_{n-\frac{1}{2}}(\zeta r)$. For half-integral orders, the reduced Bessel functions can be represented by an exponential multiplied by a terminating confluent hypergeometric function [12, 13]:

[TABLE]

where $(x)_{n}$ stands for the Pochhammer symbol, which may be defined in terms of the Gamma function $\Gamma(z)$ as [48]:

[TABLE]

Reduced Bessel functions are related to modified Bessel functions of the second kind by [49]:

[TABLE]

Reduced and modified Bessel functions satisfy the following three-term recurrence relations [48]:

[TABLE]

The integral representation of the modified Bessel function is given by [48]:

[TABLE]

The modified Bessel functions of the second kind satisfy the following property for some $m\in\mathbb{N}_{0}$ [48]:

[TABLE]

The normalized STFs are defined by [10]:

[TABLE]

where $N(\zeta,n)$ stands for the normalization factor and it is given by:

[TABLE]

STFs can be expressed as finite linear combinations of $B$ functions [13]:

[TABLE]

where:

[TABLE]

A given function $f(\vec{r})$ and its Fourier transform $\bar{f}(\vec{k})$ are connected by the symmetric relationships [50]:

[TABLE]

Consequently, the Rayleigh expansion of the plane wavefunctions is useful [51]:

[TABLE]

where the spherical Bessel function $j_{n}(x)$ of order $n\in\mathbb{N}_{0}$ is defined by [48]:

[TABLE]

Spherical Bessel functions are related to Bessel functions of the first kind $J_{n+1/2}(x)$ by [48]:

[TABLE]

Spherical Bessel functions and Bessel functions of the first kind satisfy the three-term recurrence relations [48]:

[TABLE]

The Bessel functions satisfy the following properties [48]:

[TABLE]

Due to the Rayleigh expansion, the $B$ functions have a relatively simple Fourier transform [22]:

[TABLE]

The Fourier integral representation of the Coulomb operator is given by [52]:

[TABLE]

The Gaunt coefficients are defined as [53, 54, 55, 56]:

[TABLE]

These coefficients linearize the product of two spherical harmonics:

[TABLE]

where the subscript $l=l_{\rm min},2$ in the summation implies that the summation index $l$ runs in steps of $2$ from $l_{\rm min}$ to $l_{1}+l_{2}$, and the constant $l_{\rm min}$ is given by:

[TABLE]

The three-center nuclear attraction integrals over $B$ functions are given by:

[TABLE]

where $A$, $B$, and $C$ are three arbitrary points of $\mathbb{R}^{3}$, and $O$ is the origin. By performing a translation of $\overrightarrow{OA}$ and by substituting the integral representation of the Coulomb operator in the above equation, the integrals can be re-written as:

[TABLE]

where:

[TABLE]

and where $\vec{r}=\vec{R}-\overrightarrow{OA}$, $\vec{R}_{1}=\overrightarrow{OC}$, and $\vec{R}_{2}=\overrightarrow{AB}$.

Using all of the above properties, the three-center nuclear attraction integrals over $B$ functions have the representation [57, 58]:

[TABLE]

where :

[TABLE]

and where $v$ and $R_{2}$ stand for the modulus of $\vec{v}$ and $\vec{R}_{2}$ respectively.

The bottleneck in the numerical evaluation of this expression is the semi-infinite integral:

[TABLE]

which has varying degrees of oscillation and attenuation depending upon the values of the parameters.

The semi-infinite integrals ${\cal I}(s)$ (30) involves spherical Bessel functions $j_{\lambda}(v\,x)$ and this is why their accurate and rapid numerical evaluation emerge as the critical calculation in the analytical expressions, in particular for values of $s$ close to [math] or $1$. If we make the substitution $s=0$ or $1$ in ${\cal I}_{s}$, $\gamma_{12}(s,x)$ becomes a constant and the integrand will be reduced to $x^{n_{x}}\,j_{\lambda}(v\,x)$, because the exponential decreasing part $\hat{k}_{\nu}$ of the integrand becomes a constant and consequently the rapid oscillations of $j_{\lambda}(v\,x)$ cannot be damped and suppressed as it can be seen from Figure (1). In addition, when the value of $v$ is large, the zeros of the integrands become closer and the oscillations become strong and the accurate numerical evaluation of the semi-infinite integrals becomes extremely challenging.

Originally, Gauss-Laguerre quadrature is used [57, 58], almost completely ignoring the possible effects of the oscillations.

The semi-infinite integrals can be transformed into an infinite series:

[TABLE]

where $j_{\lambda,v}^{0}\,=\,0$ and $j_{\lambda,v}^{n}=j_{\lambda+\frac{1}{2}}^{n}/v$ ($v\not=0$) for $n=1,2,\ldots$ are the successive positive zeros of $j_{\lambda}(vx)$, where $j_{\lambda+\frac{1}{2}}^{n}$ for $n=1,2,\ldots$ are the successive positive zeros of $j_{\lambda}(x)$.

Unfortunately, these infinite series are slowly convergent. The use of Levin’s $u$ transform [59] or the epsilon algorithm of Wynn [60], which are considered among the most powerful convergence accelerators, could improve convergence of the infinite series. Unfortunately, the calculation times are still prohibitive for a high pre-determined accuracy. Then, the extrapolation methods are implemented [26, 27, 28, 29, 30, 31], which use numerical integration of successive intervals to construct approximations the semi-infinite integral. While these methods are generally highly accurate, there are some ranges of parameters where the computation becomes extremely heavy.

3 The analytical development

We begin by considering integrals of the form [49]:

[TABLE]

where $\nu\in\mathbb{R}$, ${\rm Re}\,\eta>-1$, $\alpha>0$, $\beta>0$, and $|\arg z|<\frac{\pi}{2}$.

In these integrals, cylindrical Bessel functions are used instead of spherical Bessel functions. These integrals have an analytical expression in Watson’s Treatise on the Theory of Bessel functions [49, §13.47] that can be proved easily.

We begin by inserting the integral representation of the modified Bessel function (8) into the expression (32) and we reverse the order of integration:

[TABLE]

By using the convergent series for the Bessel function, the inner integrals are all readily obtainable [49, §6.22]:

[TABLE]

Then, the outer integral is evaluated by recognizing it as a modified Bessel function, and we obtain:

[TABLE]

3.1 The semi-infinite integral ${\cal I}(s)$

Let the integral ${\cal I}_{\lambda,r}^{\nu,\mu}(\alpha,\beta,z)$ be given by:

[TABLE]

The simpler case of ${\cal I}_{\lambda,r}^{\nu,\mu}(\alpha,\beta,z)$ is given by:

[TABLE]

We can also write:

[TABLE]

This simpler case is the essential starting point because identities can be used to increase the $r$ and $\mu$ parameters to obtain an expression for ${\cal I}_{\lambda,r}^{\nu,\mu}(\alpha,\beta,z)$. If $r,\mu\geq 0$, then by applying the following identities, which follow directly from equations (9) and (19):

[TABLE]

and:

[TABLE]

to the integral ${\cal I}_{\lambda+r,0}^{\nu-\mu,0}(\alpha,\beta,z)$ (39), we obtain:

[TABLE]

Upon changing the Bessel functions to their spherical and reduced counterparts in (36), we obtain:

[TABLE]

By using the analytic expression obtained in equation (36), we obtain:

[TABLE]

and:

[TABLE]

Then, upon inserting the expression for ${\cal I}_{\lambda+r,0}^{\nu-\mu,0}(\alpha,\beta,z)$ into (42), given the reflection formula $K_{\tau}(x)=K_{-\tau}(x)$, we obtain:

[TABLE]

where:

[TABLE]

Evidently, all that remains to do is to expand the derivations in (46), and the expression will be obtained. We start with the derivations with respect to $\beta$.

Using the identities:

[TABLE]

for some $\sigma\in\mathbb{N}_{0}$, the product rule yields:

[TABLE]

We continue with the derivations with respect to $\alpha$. Using the identities:

[TABLE]

for some $m\in\mathbb{N}_{0}$, the product rule yields:

[TABLE]

Then, combining these results and inserting them in (46), we obtain:

[TABLE]

Finally, replacing $\tau$ by its original value given by (47), we obtain:

[TABLE]

Now, the semi-infinite integrals ${\cal I}(s)$ (30) are formulated using the spherical and reduced Bessel functions, and they can be expressed as follows:

[TABLE]

where $\nu$, $n{\gamma}$, $n_{x}$, $v$, $\lambda$ and $s$ are given by (29). The integers $\mu$ and $r$, and the parameter $z$ are given by:

[TABLE]

Using the expressions of $\nu$, $n_{\gamma}$ and $n_{x}$ given by (29), the summation indices $\lambda$, $l$, $l_{1}^{\prime}$, $l_{2}^{\prime}$ and $j$ involed in (28) and defined by equation (24), one can easily verify that the integer $\mu\,\geq\,0$.

However, the integer $r$ can be negative, but this only happens when $\lambda=n_{x}$ leading to $r=-1$. In this case, the analytical development is no longer valid and a different method should be used to compute the corresponding semi-infinite integral ${\cal I}(s)$.

Using equation (54) and if $\lambda\not=n_{x}$, we obtain:

[TABLE]

For the case $r=-1$ corresponding to $n_{x}=\lambda$ in (30), one could use a most recent approach introduced in [61] and based on the $S$ transformation and the double exponential transformation. Other methods based on extrapolation methods and nonlinear transformations could also be used [26, 31].

4 Numerical Discussion

Since the order of the Bessel function is always an integer in (3.1), the accuracy and efficiency of this formula are directly related to that of the algorithm for the calculation of a sequence of Bessel functions $\{K_{\ell}(\cdot)\}_{\ell\in\mathbb{N}_{0}}$.

Given the recurrence relation (9) is stable in the upward direction, and given the reflection formula $K_{\nu}(z)=K_{-\nu}(z)$ an efficient algorithm would provide the seed values $K_{0}(z)$ and $K_{1}(z)$ and compute the remaining values by recurrence. In [62, §6.5], a very fast algorithm for these seed values is derived based on the use of rational minimax approximants. These approximants are rational functions with unknowns that are designed to minimize the maximum error on an interval. The algorithm achieves double precision accuracy by expressing $K_{0}(z)$ and $K_{1}(z)$ for $z\leq 1$ as:

[TABLE]

and for $z>1$ as:

[TABLE]

where $p_{n}^{(m)}(z)$ stands for a polynomial of degree $n$ in $z$ and where $m=0,1$ serves to distinguish between the polynomial for $K_{0}(z)$ and that for $K_{1}(z)$. The coefficients of these approximants are given in [63].

Tables 1 and 2 contain values obtained for the semi-infinite integrals ${\cal I}(s)$ (30). The calculations for values $s=0.05$ which is close to [math], and $s=0.99$ which is close to $1$.

In Tables 1 and 2, we listed values of the semi-infinite integrals ${\cal I}(s)$ obtained by a general MATLAB built-in numerical integration function that uses global adaptive quadrature, set to an accuracy of $15$ correct digits. However, the MATLAB built-in function is not always able to complete the integral approximation, particularly for large values of $\lambda$ and/or $v$. In contrast, the analytical expression is able to complete the computation to machine accuracy regardless of the values of the parameters.

To further illustrate the high efficiency of the analytic expressions, we computed the semi-infinite integrals ${\cal I}(s)$ using two existing methods proven to be accurate and efficient, namely the combination of $S$ (See [26]) and $\bar{D}$ transformations [64], and the transformation $\bar{D}$ with the $W$ algorithm of Sidi [65].

As can be seen from the numerical tables, the analytic expressions, which achieve machine accuracy, lead to a substantial gain in the calculation time. The gain in efficiency is of the order of $800$.

5 Conclusion

In this work, we report the analytical expression (3.1) for the semi-infinite integral bottleneck occurring in the three-center nuclear attraction integrals over $B$ functions. We describe how to compute the formula to obtain an efficient evaluation in double precision arithmetic. This requires the rational minimax approximants that minimize the maximum error on the interval of evaluation. The numerical tests show a substantial gain in efficiency over existing methods based on nonlinear transformations and extrapolation methods.

6 Numerical tables