A fast algorithm for computing the Boys function

Beylkin, Gregory; Sharma, Sandeep

doi:10.1063/5.0062444

We present a new fast algorithm for computing the Boys function using a nonlinear approximation of the integrand via exponentials. The resulting algorithms evaluate the Boys function with real and complex valued arguments and are competitive with previously developed algorithms for the same purpose.

I. INTRODUCTION

The Boys function¹

F (n, z) = \int_{0}^{1} e^{- z t^{2}} t^{2 n} d t = \frac{1}{2} \int_{0}^{1} e^{- z s} s^{n - 1 / 2} d s

(1)

appears in problems of computing Gaussian integrals, and over the years, there were many algorithms proposed for its evaluation; see, e.g., Refs. 2–10. The Boys function is related to a number of special functions, for example, the error function and the incomplete Gamma function, and (for a pure imaginary argument) to the Fresnel integrals.

It is common (see, e.g., Refs. 5 and 7) to use recursion to compute the Boys function for different n. The recursion is obtained via integration by parts,

\begin{matrix} F (n, z) & = & - \frac{1}{2 z} \int_{0}^{1} \frac{d}{d s} (e^{- z s}) s^{n - 1 / 2} d s \\ = & \frac{n - 1 / 2}{z} F (n - 1, z) - \frac{1}{2 z} e^{- z}, \end{matrix}

(2)

and can be run starting with n = 1 so that we need to have the value $F (0, z)$ or starting from a large n = n_max and going to n = 1,

F (n - 1, z) = \frac{z}{n - 1 / 2} F (n, z) + \frac{1}{2 (n - 1 / 2)} e^{- z},

(3)

so that we need to have the value $F (n_{\max}, z)$ ⁠. In order to avoid a loss of accuracy, the choice of which recursion to use depends on the size z and n_max. Iterating recursion (2), the dominant term expressing $F (n, z)$ via $F (0, z)$ is $\prod_{j = 1}^{n} (j - 1 / 2) / z^{n}$ ⁠. We set

z^{*} = {(\prod_{j = 1}^{n} (j - 1 / 2))}^{1 / n}

and choose (2) when $|z| \geq z^{*}$ and (3) if $|z| < z^{*}$ ⁠. For example, if n_max = 18, then z* ≈ 6.75. We note that other choices of the parameter z* are possible.

At each step, recursions (2) and (3) require only three multiplications and one addition (since the coefficients can be computed in advance and stored), so it is hard to obtain a more efficient alternative if one needs to compute these functions for a range of n, 1 ≤ n ≤ n_max. In order to initialize these recursions, we need fast algorithms for computing $F (0, z)$ and $F (n_{\max}, z)$ ⁠. Computing $F (0, z)$ for real z is straightforward since

F (0, z) = \int_{0}^{1} e^{- z t^{2}} d t = \frac{\sqrt{π} Erf (\sqrt{z})}{2 \sqrt{z}} .

(4)

For a real argument, an optimized implementation of the error function Erf is available within programming languages. For a complex argument, we present an algorithm for computing $F (0, z)$ using a nonlinear approximation of the integrand following the approach in Ref. 11. We obtain a rational approximation of $F (0, z)$ with an additional exponential factor.

We note that, as a function of complex argument, the Boys function $F (0, z)$ can be highly oscillatory. In particular, if z is purely imaginary, then the Boys function is related to the Fresnel integrals,

S (y) = \int_{0}^{y} \sin (\frac{π}{2} t^{2}) d t, C (y) = \int_{0}^{y} \cos (\frac{π}{2} t^{2}) d t

so that

C (y) - i S (y) = \int_{0}^{y} e^{- i \frac{π}{2} t^{2}} d t = y \int_{0}^{1} e^{- i \frac{π}{2} y^{2} s^{2}} d s = y F (0, i \frac{π}{2} y^{2}) .

(5)

For computing $F (n_{\max}, z)$ ⁠, instead of tabulating this function as it is performed for a real argument in, e.g., Refs. 2, 5, 7, and 9, we use a nonlinear approximation of the integrand in (1) (see Ref. 11), leading to an approximation of the Boys function valid for the complex argument $R e (z) \geq 0$ with tight error estimates. For $R e (z) < 0$ ⁠, we compute $e^{z} F (n, z)$ instead of $F (n, z)$ ⁠. Based on these approximations, we develop two algorithms, for real and complex valued arguments. We refer to Refs. 3, 4, and 8 for previously developed algorithms for the Boys function with the complex argument. The complex argument appears in a number of problems, for example, in calculations with mixed Gaussian/plane wave bases in molecules and scattering problems,^12–16 in the context of complex scaling calculations of excited states,¹⁷ and in using gauge invariant basis functions for calculating magnetic properties.¹⁸

II. APPROXIMATION OF $F (0, z)$ FOR COMPLEX VALUED ARGUMENT

We have

F (0, z) = \int_{0}^{1} e^{- z t^{2}} d t = \frac{1}{2} \int_{0}^{1} e^{- z s} s^{- 1 / 2} d s

(6)

and use the integral

s^{- 1 / 2} = \frac{2}{\sqrt{π}} \int_{0}^{\infty} e^{- s t^{2}} d t

(7)

to obtain

F (0, z) = \frac{1}{\sqrt{π}} \int_{0}^{\infty} \frac{1 - e^{- (t^{2} + z)}}{t^{2} + z} d t = \frac{1}{\sqrt{π}} \int_{0}^{\infty} q (t^{2} + z) d t,

(8)

where

q (ξ) = (1 - e^{- ξ}) / ξ, ξ \in C,

is an analytic function. An algorithm for computing $F (0, z)$ is essentially a quadrature for the integral in (8). Note that if, instead, we were to use a quadrature to compute $F (0, z)$ via integrals in (6), then, for each z, we would need to evaluate as many exponentials as the number of quadrature terms. Importantly, when using (8), we need to evaluate e^−z only once and then use the result as a factor.

A. The case $R e (z) \geq 0$ ⁠.

We split integral (8) into three terms

\begin{align} F (0, z) & = \frac{1}{\sqrt{π}} \int_{0}^{\infty} \frac{1}{t^{2} + z} d t - \frac{e^{- z}}{\sqrt{π}} \int_{0}^{t_{\max}} \frac{e^{- t^{2}}}{t^{2} + z} d t \\ - \frac{e^{- z}}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{- t^{2}}}{t^{2} + z} d t \end{align}

(9)

and observe that the last term in (9) (without the factor e^−z) is estimated as

\begin{matrix} |\frac{1}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{- t^{2}}}{t^{2} + z} d t| & \leq & \frac{1}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{- t^{2}}}{|t^{2} + z|} d t \leq \frac{1}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{- t^{2}}}{t^{2}} d t \\ = & \frac{1}{\sqrt{π}} (\frac{e^{- t_{\max}^{2}}}{t_{\max}} - \sqrt{π} Erfc (t_{\max})) \\ = & ϵ_{t_{\max}} . \end{matrix}

(10)

Select t_max = e^7/4 to obtain $ε_{t_{\max}} \approx 5.9 \cdot 1 0^{- 18}$ ⁠. For the first term in (9), we have

\frac{1}{\sqrt{π}} \int_{0}^{\infty} \frac{1}{t^{2} + z} d t = \frac{1}{2} \sqrt{\frac{π}{z}} .

For $|z| \geq r_{0} = 0.35$ ⁠, we use quadrature (see the Appendix for details) to approximate the second term in (9) as

|\frac{1}{\sqrt{π}} \int_{0}^{t_{\max}} \frac{e^{- t^{2}}}{t^{2} + z} d t - \sum_{m = 1}^{M} \frac{w_{m} e^{- η_{m}}}{η_{m} + z}| \leq ε,

(11)

where M = 22 and nodes and weights are given in Table I. We note that it is possible to use the standard Gauss–Legendre quadrature on the interval $[0, t_{\max}]$ ⁠, but the number of terms, M, will be larger. As a result, we obtain approximation

|F (0, z) - (\frac{1}{2} \sqrt{\frac{π}{z}} - \frac{1}{2 \sqrt{π}} e^{- z} \sum_{m = 1}^{22} \frac{w_{m} e^{- η_{m}}}{η_{m} + z})| \leq 2 ε + ε_{t_{\max}}, |z| \geq r_{0} .

(12)

TABLE I.

The poles and weights in (12).

m	η_m	$w_{m} \cdot e^{- η_{m}}$	m	η_m	$w_{m} \cdot e^{- η_{m}}$
1	0.147 787 826 379 695 65 × 10⁻⁰²	0.866 431 027 201 416 54 × 10⁻⁰¹	12	0.125 395 022 879 192 93 × 10⁺⁰¹	0.574 448 042 214 302 23 × 10⁻⁰¹
2	0.133 172 764 137 258 17 × 10⁻⁰¹	0.857 720 608 434 394 68 × 10⁻⁰¹	13	0.172 446 342 335 733 95 × 10⁺⁰¹	0.420 819 943 469 454 42 × 10⁻⁰¹
3	0.370 635 914 520 525 41 × 10⁻⁰¹	0.839 350 436 829 178 76 × 10⁻⁰¹	14	0.237 152 482 627 818 63 × 10⁺⁰¹	0.258 385 394 482 232 72 × 10⁻⁰¹
4	0.727 525 124 228 827 62 × 10⁻⁰¹	0.809 661 970 413 229 21 × 10⁻⁰¹	15	0.326 137 969 960 783 55 × 10⁺⁰¹	0.124 450 241 572 555 60 × 10⁻⁰¹
5	0.120 236 941 228 785 68 × 10⁺⁰⁰	0.769 089 548 492 978 56 × 10⁻⁰¹	16	0.448 513 016 905 959 11 × 10⁺⁰¹	0.429 254 159 259 983 68 × 10⁻⁰²
6	0.179 574 293 958 937 73 × 10⁺⁰⁰	0.731 552 078 711 821 68 × 10⁻⁰¹	17	0.616 806 213 512 248 38 × 10⁺⁰¹	0.935 434 298 773 596 86 × 10⁻⁰³
7	0.253 534 046 984 087 27 × 10⁺⁰⁰	0.726 950 035 163 157 20 × 10⁻⁰¹	18	0.848 247 187 231 786 981 × 10⁺⁰¹	0.108 408 854 665 025 05 × 10⁻⁰³
8	0.350 388 652 780 721 95 × 10⁺⁰⁰	0.752 842 556 089 304 00 × 10⁻⁰¹	19	0.116 653 054 862 967 931 × 10⁺⁰²	0.527 186 796 676 167 36 × 10⁻⁰⁵
9	0.482 109 575 931 276 68 × 10⁺⁰⁰	0.770 943 953 645 196 33 × 10⁻⁰¹	20	0.160 424 171 322 883 281 × 10⁺⁰²	0.776 597 403 975 041 90 × 10⁻⁰⁷
10	0.663 028 993 158 374 16 × 10⁺⁰⁰	0.754 250 625 677 530 40 × 10⁻⁰¹	21	0.220 619 295 181 470 89 × 10⁺⁰²	0.221 381 724 226 800 93 × 10⁻⁰⁹
11	0.911 814 736 856 590 87 × 10⁺⁰⁰	0.689 686 192 650 315 33 × 10⁻⁰¹	22	0.303 401 120 947 083 07 × 10⁺⁰²	0.659 416 176 003 770 69 × 10⁻¹³

m	η_m	$w_{m} \cdot e^{- η_{m}}$	m	η_m	$w_{m} \cdot e^{- η_{m}}$
1	0.147 787 826 379 695 65 × 10⁻⁰²	0.866 431 027 201 416 54 × 10⁻⁰¹	12	0.125 395 022 879 192 93 × 10⁺⁰¹	0.574 448 042 214 302 23 × 10⁻⁰¹
2	0.133 172 764 137 258 17 × 10⁻⁰¹	0.857 720 608 434 394 68 × 10⁻⁰¹	13	0.172 446 342 335 733 95 × 10⁺⁰¹	0.420 819 943 469 454 42 × 10⁻⁰¹
3	0.370 635 914 520 525 41 × 10⁻⁰¹	0.839 350 436 829 178 76 × 10⁻⁰¹	14	0.237 152 482 627 818 63 × 10⁺⁰¹	0.258 385 394 482 232 72 × 10⁻⁰¹
4	0.727 525 124 228 827 62 × 10⁻⁰¹	0.809 661 970 413 229 21 × 10⁻⁰¹	15	0.326 137 969 960 783 55 × 10⁺⁰¹	0.124 450 241 572 555 60 × 10⁻⁰¹
5	0.120 236 941 228 785 68 × 10⁺⁰⁰	0.769 089 548 492 978 56 × 10⁻⁰¹	16	0.448 513 016 905 959 11 × 10⁺⁰¹	0.429 254 159 259 983 68 × 10⁻⁰²
6	0.179 574 293 958 937 73 × 10⁺⁰⁰	0.731 552 078 711 821 68 × 10⁻⁰¹	17	0.616 806 213 512 248 38 × 10⁺⁰¹	0.935 434 298 773 596 86 × 10⁻⁰³
7	0.253 534 046 984 087 27 × 10⁺⁰⁰	0.726 950 035 163 157 20 × 10⁻⁰¹	18	0.848 247 187 231 786 981 × 10⁺⁰¹	0.108 408 854 665 025 05 × 10⁻⁰³
8	0.350 388 652 780 721 95 × 10⁺⁰⁰	0.752 842 556 089 304 00 × 10⁻⁰¹	19	0.116 653 054 862 967 931 × 10⁺⁰²	0.527 186 796 676 167 36 × 10⁻⁰⁵
9	0.482 109 575 931 276 68 × 10⁺⁰⁰	0.770 943 953 645 196 33 × 10⁻⁰¹	20	0.160 424 171 322 883 281 × 10⁺⁰²	0.776 597 403 975 041 90 × 10⁻⁰⁷
10	0.663 028 993 158 374 16 × 10⁺⁰⁰	0.754 250 625 677 530 40 × 10⁻⁰¹	21	0.220 619 295 181 470 89 × 10⁺⁰²	0.221 381 724 226 800 93 × 10⁻⁰⁹
11	0.911 814 736 856 590 87 × 10⁺⁰⁰	0.689 686 192 650 315 33 × 10⁻⁰¹	22	0.303 401 120 947 083 07 × 10⁺⁰²	0.659 416 176 003 770 69 × 10⁻¹³

B. The case $R e (z) < 0$

In this case, we compute $e^{z} F (0, z)$ rather than $F (0, z)$ ⁠. Since the denominator in (8) can be zero, we cannot separate terms in $q (t^{2} + z)$ as in (9). Instead, we split integral (8) into two terms and obtain

e^{z} F (0, z) = \frac{e^{z}}{\sqrt{π}} \int_{0}^{t_{\max}} \frac{1 - e^{- (t^{2} + z)}}{t^{2} + z} d t + \frac{1}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{z} - e^{- t^{2}}}{t^{2} + z} d t .

(13)

The first term in (13) is approximated by using the Gauss–Legendre quadrature on the interval $[0, t_{\max}]$ ⁠. The function q is analytic, and therefore, there is no singularity at t² = −z. Since we can compute derivatives of q, the error introduced by this quadrature can be estimated using results in Sec. 5.2 of Ref. 19. For example, we obtain

|\frac{1}{\sqrt{π}} \int_{0}^{t_{\max}} \frac{e^{z} - e^{- t^{2}}}{t^{2} + z} d t - \frac{1}{\sqrt{π}} \sum_{m = 1}^{M^{g}} w_{m}^{g} \frac{e^{z} (1 - e^{- (t_{m}^{2} + z)})}{t_{m}^{2} + z}| \leq ε^{g}

with M^g = 16 and ɛ^g ≈ 10⁻¹⁴, where $t_{m,} w_{m}^{g}$ are the standard Gauss–Legendre nodes and weights on the interval $[0, t_{\max}]$ ⁠. In the second term in (13), we drop $e^{- t^{2}}$ (since its contribution is less than $e^{- t_{\max}^{2}} \approx 4.2 \cdot 1 0^{- 15}$ ⁠) and obtain

\frac{1}{\sqrt{π}} \int_{t_{\max}}^{\infty} \frac{e^{z}}{t^{2} + z} d t = \frac{e^{z}}{\sqrt{π z}} Arctan (\frac{\sqrt{z}}{t_{\max}}) .

(14)

While we obtain an explicit expression, computing arctangent of a complex argument is relatively expensive. For a complex argument, we evaluate arctangent using

Arctan (z) = \frac{1}{2} i \log \frac{1 - i z}{1 + i z} .

As a result of dropping $e^{- t^{2}}$ in the second term of (13), our approximation in (14) has a singularity at $z = - t_{\max}^{2}$ ⁠. In order to avoid using (14) in the vicinity of singularity, we use two different parameters t_max and t_max,1 and switch to the version with t_max,1 if $|z + t_{\max}^{2}| \leq 1 / 2$ ⁠, where $t_{\max, 1} = \sqrt{t_{\max}^{2} + 1}$ ⁠.

We note that it is possible to increase the number of terms in the quadrature in order to avoid evaluating arctangent. This might be of interest on a parallel [Graphics Processing Unit (GPU) or multi-core] computer since computation of quadrature terms is trivially parallel. As a result, we obtain approximation

\begin{align} |e^{z} F (0, z) - \frac{e^{z}}{\sqrt{π z}} Arctan (\frac{\sqrt{z}}{t_{\max}}) - \frac{1}{\sqrt{π}} \sum_{m = 1}^{M^{g}} w_{m}^{g} \frac{e^{z} (1 - e^{- (t_{m}^{2} + z)})}{t_{m}^{2} + z}| \\ \leq \tilde{ε}, |z + t_{\max}^{2}| > 1 / 2, \end{align}

(15)

where $\tilde{ε} \approx 1 0^{- 14}$ ⁠. For $|z| > t_{\max}^{2}$ ⁠, we have a converging series for the second integral in (9) as follows:

\begin{matrix} \int_{0}^{t_{\max}} \frac{e^{- t^{2}}}{t^{2} + z} d t & = & \frac{1}{z} \int_{0}^{t_{\max}} \frac{e^{- t^{2}}}{t^{2} / z + 1} d t \\ = & \frac{1}{z} \sum_{j = 0}^{\infty} {(- 1)}^{j} z^{- j} \int_{0}^{t_{\max}} e^{- t^{2}} t^{2 j} d t \\ = & \frac{1}{2 z} \sum_{j = 0}^{\infty} {(- 1)}^{j} z^{- j} (Γ (j + 1 / 2) - Γ (j + 1 / 2, t_{\max}^{2})) \\ = & \frac{1}{z} \sum_{j = 0}^{\infty} {(- 1)}^{j} z^{- j} t_{\max}^{2 j + 1} F (j, t_{\max}^{2}) \end{matrix}

(16)

so that we can use

\begin{align} |F (0, z) - [\frac{\sqrt{π}}{2 \sqrt{z}} - \frac{e^{- z}}{2 \sqrt{π} z} \sum_{j = 0}^{J} {(- 1)}^{j} z^{- j} (Γ (j + 1 / 2)))) \\ ((- Γ (j + 1 / 2, t_{\max}^{2}) ()]| \leq ϵ_{t_{\max}}, |z| > t_{\max}^{2}, \end{align}

(17)

instead of (12) and

\begin{align} |e^{z} F (0, z) - [\frac{e^{z} \sqrt{π}}{2 \sqrt{z}} - \frac{1}{2 \sqrt{π} z} \sum_{j = 0}^{J} {(- 1)}^{j} z^{- j} (Γ (j + 1 / 2)))) \\ ((- Γ (j + 1 / 2, t_{\max}^{2}) ()]| \leq ϵ_{t_{\max}}, |z| > t_{\max}^{2}, \end{align}

(18)

instead of (15). Since the parameter t_max is fixed, the coefficients of the series are computed offline.

TABLE II.

Weights and exponents of the approximation of $g_{12} (s)$ on $[0, 1]$ in (22). With these parameters, the absolute error in (23) is ɛ ≈ 2.5 · 10⁻¹³.

m	η_m	w_m
1	0.707 194 313 205 700 10 · 10¹ + 0.164 872 912 507 521 15 · 10²i	0.364 436 324 028 985 01 · 10⁻¹⁰ + 0.264 117 510 721 075 04 · 10⁻¹⁰i
2	0.707 194 313 205 700 10 · 10¹ − 0.164 872 912 507 521 15 · 10²i	0.364 436 324 028 985 01 · 10⁻¹⁰ − 0.264 117 510 721 075 04 · 10⁻¹⁰i
3	−0.571 432 717 151 916 35 + 0.132 785 794 532 336 33 · 10²i	0.181 852 503 467 536 33 · 10⁻⁶ − 0.218 604 589 713 993 52 · 10⁻⁵i
4	−0.571 432 717 151 916 35 − 0.132 785 794 532 336 33 · 10²i	0.181 852 503 467 536 33 · 10⁻⁶ + 0.218 604 589 713 993 52 · 10⁻⁵i
5	−0.471 930 213 303 925 06 · 10¹ + 0.998 352 571 123 710 32 · 10¹i	−0.994 891 692 720 557 48 · 10⁻³ − 0.230 490 791 052 030 73 · 10⁻³i
6	−0.471 930 213 303 925 06 · 10¹ − 0.998 352 571 123 710 32 · 10¹i	−0.994 891 692 720 557 48 · 10⁻³ + 0.230 490 791 052 030 73 · 10⁻³i
7	−0.717 046 627 728 950 89 · 10¹ + 0.667 123 608 398 207 68 · 10¹i	−0.256 252 169 858 790 06 · 10⁻¹ + 0.358 183 352 748 769 82 · 10⁻¹i
8	−0.717 046 627 728 950 89 · 10¹ − 0.667 123 608 398 207 68 · 10¹i	−0.256 252 169 858 790 06 · 10⁻¹ − 0.358 183 352 748 769 82 · 10⁻¹i
9	−0.848 997 470 547 246 99 · 10¹ + 0.334 348 041 684 674 91 · 10¹i	0.165 068 015 448 807 23 + 0.322 739 644 717 760 45i
10	−0.848 997 470 547 246 99 · 10¹ − 0.334 348 041 684 674 91 · 10¹i	0.165 068 015 448 807 23 − 0.32 273 964 471 776 045i
11	0.365 644 143 631 509 73 · 10²	−0.201 046 416 615 651 64 · 10⁻²⁵
12	−0.324 242 392 559 219 54 · 10¹	−0.395 635 369 550 420 78 · 10⁻³
13	−0.890 660 477 331 007 53 · 10¹	0.723 499 458 050 852 92

m	η_m	w_m
1	0.707 194 313 205 700 10 · 10¹ + 0.164 872 912 507 521 15 · 10²i	0.364 436 324 028 985 01 · 10⁻¹⁰ + 0.264 117 510 721 075 04 · 10⁻¹⁰i
2	0.707 194 313 205 700 10 · 10¹ − 0.164 872 912 507 521 15 · 10²i	0.364 436 324 028 985 01 · 10⁻¹⁰ − 0.264 117 510 721 075 04 · 10⁻¹⁰i
3	−0.571 432 717 151 916 35 + 0.132 785 794 532 336 33 · 10²i	0.181 852 503 467 536 33 · 10⁻⁶ − 0.218 604 589 713 993 52 · 10⁻⁵i
4	−0.571 432 717 151 916 35 − 0.132 785 794 532 336 33 · 10²i	0.181 852 503 467 536 33 · 10⁻⁶ + 0.218 604 589 713 993 52 · 10⁻⁵i
5	−0.471 930 213 303 925 06 · 10¹ + 0.998 352 571 123 710 32 · 10¹i	−0.994 891 692 720 557 48 · 10⁻³ − 0.230 490 791 052 030 73 · 10⁻³i
6	−0.471 930 213 303 925 06 · 10¹ − 0.998 352 571 123 710 32 · 10¹i	−0.994 891 692 720 557 48 · 10⁻³ + 0.230 490 791 052 030 73 · 10⁻³i
7	−0.717 046 627 728 950 89 · 10¹ + 0.667 123 608 398 207 68 · 10¹i	−0.256 252 169 858 790 06 · 10⁻¹ + 0.358 183 352 748 769 82 · 10⁻¹i
8	−0.717 046 627 728 950 89 · 10¹ − 0.667 123 608 398 207 68 · 10¹i	−0.256 252 169 858 790 06 · 10⁻¹ − 0.358 183 352 748 769 82 · 10⁻¹i
9	−0.848 997 470 547 246 99 · 10¹ + 0.334 348 041 684 674 91 · 10¹i	0.165 068 015 448 807 23 + 0.322 739 644 717 760 45i
10	−0.848 997 470 547 246 99 · 10¹ − 0.334 348 041 684 674 91 · 10¹i	0.165 068 015 448 807 23 − 0.32 273 964 471 776 045i
11	0.365 644 143 631 509 73 · 10²	−0.201 046 416 615 651 64 · 10⁻²⁵
12	−0.324 242 392 559 219 54 · 10¹	−0.395 635 369 550 420 78 · 10⁻³
13	−0.890 660 477 331 007 53 · 10¹	0.723 499 458 050 852 92

Note that the series in (16) and (17) is related to the asymptotic expansion of $F (0, z)$ (see, e.g., Ref. 4),

\begin{align} F (0, z) & \sim \frac{1}{2} \frac{\sqrt{π}}{\sqrt{z}} - \frac{\sqrt{π}}{2} \frac{e^{- z}}{z} \sum_{j = 0}^{J} \frac{z^{- j}}{Γ (\frac{1}{2} - j)} \\ = \frac{1}{2} \frac{\sqrt{π}}{\sqrt{z}} - \frac{e^{- z}}{2 \sqrt{π} z} \sum_{j = 0}^{J} {(- 1)}^{j} z^{- j} Γ (j + 1 / 2) . \end{align}

(19)

We use (17) and (18) for $|z| \geq 100$ so that it is sufficient to keep only seven terms, yielding an error of less than 10⁻¹³.

For $|z| \leq r_{0}$ ⁠, we use the Taylor expansion of (6),

F (0, z) = \sum_{j = 0}^{\infty} \frac{{(- 1)}^{j} z^{j}}{j! (2 j + 1)},

(20)

and we need 10 terms to maintain an accuracy of about 13 digits. While selecting parameters as above leads to algorithms with a reasonable speed, we did not optimize these choices as they may depend on several factors, e.g., computer architecture.

Since the Boys function $F (0, z)$ is related to the error function (and can be used to compute it), we compared the speed of our algorithm with that of the well-known algorithm by Gautschi²⁰ for computing the error function with a complex argument using a rational approximation of the closely related Faddeeva function. The speed of that algorithm was measured in comparison with the speed of computing $\exp (z)$ ⁠. In Ref. 20, it is stated that with an accuracy of ∼10 digits, the code is 7–15 times slower than the speed of computing $\exp (z)$ ⁠. Using the same comparison for our algorithm, this ratio is ∼12 for an accuracy of about 13 digits. Our algorithm is implemented using Fortran 90 compiled by Intel’s ifort with compiler flags -O₃ -ipo -static and running on a laptop with $\approx 2.3$ GHz chipset. We timed our code by performing 10⁶ evaluations, yielding $\approx 0.92 \cdot 1 0^{- 7}$ s/evaluation in comparison with $\approx 0.79 \cdot 1 0^{- 8}$ s/evaluation for $\exp (z)$ with a complex argument.

While algorithms for computing the Fresnel integrals appear to be somewhat faster than using the Boys function in (5) (see, e.g., Ref. 21), we note that the generalized Fresnel integrals, e.g., $\int_{0}^{x} e^{i t^{n}} d t$ ⁠, n ≥ 2, can be evaluated using our approach and plan to consider algorithms for these oscillatory special functions elsewhere.

III. APPROXIMATION OF $F (n_{\max}, z)$ FOR REAL AND COMPLEX ARGUMENTS

The function

g_{n} (s) = {(1 - s)}^{n - 1 / 2}

(21)

decays monotonically on $[0, 1]$ ⁠, and we use Algorithm 1 in Ref. 22 to construct its near optimal approximation via exponentials. We refer to Ref. 22 and references therein for the details of the algorithm that we use to obtain the necessary parameters (for an example, see Table II).

We obtain approximation

|g_{n} (s) - \sum_{m = 1}^{M} w_{m} e^{η_{m} s}| \leq ε for s \in [0, 1],

(22)

where $w_{m}, η_{m} \in C$ ⁠. We note that n should be sufficiently large (e.g., n ≥ 7) to avoid the impact on the approximation of the singularity of the nth derivative of g_n. Its numerical effect makes the accuracy of the current double precision implementation of Algorithm 1 in Ref. 22 insufficient to reliably produce approximation (22) for 1 ≤ n ≤ 6.

FIG. 1.

View large Download slide

The function $g_{12} (s)$ in (21) and the error of its near optimal approximation via exponentials in (22) with parameters described in Table II.

Substituting the approximation of $g_{n} (1 - s)$ into the integral (1), we arrive at

\begin{array}{l} F (n, z) - \frac{1}{2} \sum_{m = 1}^{M} w_{m} \int_{0}^{1} e^{- z s} e^{η_{m} (1 - s)} d s \\ = \frac{1}{2} \int_{0}^{1} e^{- z s} (s^{n - 1 / 2} - \sum_{m = 1}^{M} w_{m} e^{η_{m} (1 - s)}) d s \end{array}

and estimate

\begin{array}{l} |F (n, z) - \frac{1}{2} \sum_{m = 1}^{M} w_{m} \int_{0}^{1} e^{- z s} e^{η_{m} (1 - s)} d s| \leq \frac{ε}{2} \int_{0}^{1} e^{- R e (z) s} d s \\ = \frac{ε}{2} \frac{1 - e^{- R e (z)}}{R e (z)} . \end{array}

Since

\frac{1}{2} \sum_{m = 1}^{M} w_{m} \int_{0}^{1} e^{- z s} e^{η_{m} (1 - s)} d s = \frac{1}{2} \sum_{m = 1}^{M} w_{m} e^{η_{m}} \frac{1 - e^{- (z + η_{m})}}{z + η_{m}},

we have

|F (n, z) - \frac{1}{2} \sum_{m = 1}^{M} w_{m} e^{η_{m}} \frac{1 - e^{- (z + η_{m})}}{z + η_{m}}| \leq \frac{ε}{2} \frac{1 - e^{- R e (z)}}{R e (z)} \leq \frac{ε}{2} .

(23)

Indeed, denoting the factor on the right-hand side of (23), $q (z) = (1 - e^{- z}) / z$ ⁠, we have

q (2 z) = \frac{1}{2} (e^{- z} q (z) + q (z)),

and therefore, for $R e (z) \geq 0$ ⁠,

|q (2 z)| \leq |q (z)| .

This implies that $|q (z)|$ reaches it maximum at z = 0, where $q (0) = 1$ ⁠.

If $R e (z) < 0$ ⁠, we compute $e^{z} F (0, z)$ instead of $F (0, z)$ ⁠,

e^{z} F (n, z) = \frac{1}{2} \int_{0}^{1} e^{z (1 - s)} s^{n - 1 / 2} d s = \frac{1}{2} \int_{0}^{1} e^{z s} {(1 - s)}^{n - 1 / 2} d s .

Using (22), we obtain

\begin{array}{l} e^{z} F (n, z) - \frac{1}{2} \sum_{m = 1}^{M} w_{m} \int_{0}^{1} e^{z s} e^{η_{m} s} d s = \frac{1}{2} \int_{0}^{1} e^{z s} [g_{n} (s) - \sum_{m = 1}^{M} w_{m} e^{η_{m} s}] d s \end{array}

and the estimate

|e^{z} F (n, z) - \frac{1}{2} \sum_{m = 1}^{M} w_{m} \frac{e^{z + η_{m}} - 1}{z + η_{m}}| \leq \frac{ε}{2} \frac{e^{R e (z)} - 1}{R e (z)} \leq \frac{ε}{2} .

For computing values of $e^{z} F (n, z)$ for 0 ≤ n ≤ n_max for $R e (z) < 0$ ⁠, we use recursions

e^{z} F (n, z) = \frac{n - 1 / 2}{z} e^{z} F (n - 1, z) - \frac{1}{2 z}

(24)

instead of (2) and

e^{z} F (n - 1, z) = \frac{2 x}{2 n - 1} e^{z} F (n, z) + \frac{1}{2 n - 1}

(25)

instead of (3).

IV. IMPLEMENTATION DETAILS

The speed of computation of values of $F (n_{\max}, z)$ for n_max ≥ 7 depends on the number of terms M in approximation (22). We demonstrate the results of approximating $F (12, z)$ and display function $g_{12} (s)$ in Fig. 1. Using only 13 terms (see Table II), we achieve accuracy for $F (12, z) ε \approx 2 \cdot 1 0^{- 14}$ [e.g., accuracy of evaluation of $F (12, 0)$ is 2.08 · 10⁻¹⁴].

In implementing this approximation, we need to isolate cases where z is close to −η_m by using the Taylor expansion for $\frac{1 - e^{- (z + η_{m})}}{z + η_{m}}$ ⁠. Since most of η_m have an imaginary part, it is a minimal effort if z is real since η_m is real in only three terms in our example in Table II. In addition, for the real argument z, we need to use only five terms with complex valued parameters as they come in complex conjugate pairs.

We implemented these algorithms using Fortran 90 on a laptop described in Sec. II. Computing the Boys functions $F (n, z)$ for n = 0, …, 12 for the real argument takes $\approx 0.34 \cdot 1 0^{- 7}$ s. The subroutine for the complex valued argument is slower and takes $\approx 0.21 \cdot 1 0^{- 6}$ s.

V. CONCLUSION

Since their introduction in Ref. 1, the Boys functions with a real argument have widely been used for computing Gaussian integrals. When using mixed Gaussian/exponential bases, one needs to evaluate the Boys functions with complex argument. Such mixed bases are appropriate for scattering problems and for bound state problems where using only plane waves becomes too expensive near singularities. Consequently, mixed Gaussian/exponential bases provide a greater flexibility in formulation and solving problems of quantum chemistry, and we present our results, in part, to facilitate their use.

While for a real argument the Boys functions can be easily tabulated in regions where their asymptotic is not accurate, it is more difficult to apply such a straightforward implementation for a complex argument. A careful reading of Refs. 3, 4, and 8 reveals shortcomings of the existing approaches (relying mostly on expansions) to computing the Boys functions of a complex argument (see, e.g., conclusion in Ref. 8). For our approach, a better comparison is offered by Gautschi’s algorithm²⁰ for the error function of complex argument since it is related to $F (0, z)$ ⁠, as in (4); see Sec. II. Our approach of approximating a part of the integrand so that the resulting integral can be evaluated explicitly is simpler and yields tight accuracy estimates. Note that the part of the integrand we are approximating is real, while the Boys functions we are computing are complex valued. As a side remark, we note that the Boys function $F (0, z)$ remains bounded for a complex argument with $R e (z) \geq 0$ [and $e^{z} F (0, z)$ for $R e (z) < 0$ ] and, for this reason, provides a good alternative approach for computing the error function of a complex argument.

We avoid the direct timing comparisons with existing algorithms since such comparisons are generally misleading. Given different hardware (single core, multi-core, GPU, etc.), different compilers and compiler flags, and different implementations, it is hard to compare algorithms by simply running them. Instead, one can look at algorithmic possibilities an approach offers. Our code is compact, and it is easy to simply count the total number of operations. We note that computation of each term in sums (12) and (23) is trivially parallel and only recursions in (2) and (3) require a sequential implementation (with just three multiplications and one addition per function). Thus, timing of our algorithms implemented on a multi-core or GPU computer will be much faster than the quoted timings of our implementation on a single central processing unit (CPU).

SUPPLEMENTARY MATERIAL

See the supplementary material for five Fortran 90 subroutines implementing, as an example, algorithms for computing the Boys function with indices n = 0, …, 12. The subroutine dboysfun12.f90 evaluates the Boys functions $F (n, z)$ for real non-negative argument z. The subroutines zboysfun12.f90 and zboysfun00.f90 evaluate the Boys functions $F (n, z)$ for complex argument z with a non-negative real part. Finally, the subroutines zboysfun00nrp.f90 and zboysfun12nrp.f90 evaluate the functions $e^{z} F (n, z)$ for complex argument z with a negative real part.

ACKNOWLEDGMENTS

S.S. was supported by the NSF under Grant No. CHE-1800584. S.S. was also partly supported through the Sloan research fellowship.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

APPENDIX: CONSTRUCTION OF QUADRATURE IN (11)

Changing variables in (7) t = e^τ/2, we rewrite it as

s^{- 1 / 2} = \frac{1}{\sqrt{π}} \int_{- \infty}^{\infty} e^{- s e^{τ} + τ / 2} d τ, 0 \leq s \leq 1,

(A1)

and discretize it (following Ref. 23), yielding

|s^{- 1 / 2} - \sum_{m = 1}^{M} w_{m} e^{- η_{m} s}| \leq ε s^{- 1 / 2}, δ \leq s \leq 1,

(A2)

where η_m, w_m > 0 are arranged in an ascending order, and we estimate that

\begin{array}{l} |F (0, z) - \frac{1}{2} \int_{0}^{1} e^{- z s} (\sum_{m = 1}^{M} w_{m} e^{- η_{m} s}) d s| \\ \leq \frac{1}{2} \int_{0}^{1} e^{- R e (z) s} |s^{- 1 / 2} - \sum_{m = 1}^{M} w_{m} e^{- η_{m} s}| d s \\ \leq \frac{ε}{2} \int_{0}^{1} e^{- z R e (z)} s^{- 1 / 2} d s \\ = ε F (0, R e (z)) \leq ε . \end{array}

Using δ = ɛ = 10⁻¹³ in (A2) results in approximation with M = 210. We also need this approximation to satisfy

\begin{align} |\frac{1}{2 \sqrt{π}} \int_{- \infty}^{\infty} \frac{1}{e^{τ} + z} e^{τ / 2} d τ - \sum_{m = 1}^{M} \frac{w_{m}}{η_{m} + z}| \\ = |\frac{1}{2} \sqrt{\frac{π}{z}} - \sum_{m = 1}^{M} \frac{w_{m}}{η_{m} + z}| \leq ε, |z| \geq r_{0}, \end{align}

(A3)

in order to obtain

|F (0, z) - (\frac{1}{2} \sqrt{\frac{π}{z}} - \frac{e^{- z}}{2 \sqrt{π}} \sum_{m = 1}^{M} \frac{w_{m} e^{- η_{m}}}{η_{m} + z})| \leq 2 ε .

(A4)

The exponents and the weights in (A2) grow as $η_{m} \approx e^{τ_{m}}$ and $w_{m} \approx e^{τ_{m} / 2}$ (see Ref. 23) so that in (A4), it is sufficient to use a subset of terms with $η_{m} \leq e^{τ_{\max}}$ ⁠. Selecting τ_max = 7/2 so that $t_{\max} = e^{τ_{\max} / 2}$ in (10), the error $ϵ_{t_{\max}} \approx 5.9 \cdot 1 0^{- 18}$ ⁠. Consequently, we only need the 22 terms displayed in Table I and obtain the approximation of (8) in (12).

REFERENCES

1.

S. F.

Boys

, “

Electronic wave functions-I. A general method of calculation for the stationary states of any molecular system

,”

Proc. R. Soc. London, Ser. A

200

(

1063

),

542

–

554

(

1950

).

https://doi.org/10.1098/rspa.1950.0036

Google Scholar

Crossref

2.

O.

Matsuoka

, “

Field and field gradient integrals based on Gaussian type orbitals

,”

Comput. Phys. Commun.

3

(

2

),

130

–

135

(

1972

).

https://doi.org/10.1016/0010-4655(72)90060-4

Google Scholar

Crossref

3.

J. A. C.

Weideman

, “

Computation of the complex error function

,”

SIAM J. Numer. Anal.

31

(

5

),

1497

–

1518

(

1994

).

https://doi.org/10.1137/0731077

Google Scholar

Crossref

4.

P.

Cársky

and

M.

Polávsek

, “

Incomplete Gamma F_m(x) functions for real negative and complex arguments

,”

J. Comput. Phys.

143

(

1

),

259

–

265

(

1998

).

https://doi.org/10.1006/jcph.1998.5975

Google Scholar

Crossref

5.

M.

Peels

and

G.

Knizia

, “

Fast evaluation of two-center integrals over Gaussian charge distributions and Gaussian orbitals with general interaction kernels

,”

J. Chem. Theory Comput.

16

(

4

),

2570

–

2583

(

2020

).

https://doi.org/10.1021/acs.jctc.9b01296

Google Scholar

Crossref

PubMed

6.

M.

Polávsek

and

P.

Cársky

, “

Efficient evaluation of the matrix elements of the Coulomb potential between plane waves and Gaussians

,”

J. Comput. Phys.

181

(

1

),

1

–

8

(

2002

).

https://doi.org/10.1006/jcph.2002.7124

Google Scholar

Crossref

7.

B. A.

Mamedov

, “

On the evaluation of Boys functions using downward recursion relation

,”

J. Math. Chem.

36

(

3

),

301

–

306

(

2004

).

https://doi.org/10.1023/b:jomc.0000044226.49921.f5

Google Scholar

Crossref

8.

R. J.

Mathar

, “

Numerical representations of the incomplete gamma function of complex-valued argument

,”

Numer. Algorithms

36

(

3

),

247

–

264

(

2004

).

https://doi.org/10.1023/b:numa.0000040063.91709.58

Google Scholar

Crossref

9.

A. K. H.

Weiss

and

C.

Ochsenfeld

, “

A rigorous and optimized strategy for the evaluation of the Boys function kernel in molecular electronic structure theory

,”

J. Comput. Chem.

36

(

18

),

1390

–

1398

(

2015

).

https://doi.org/10.1002/jcc.23935

Google Scholar

Crossref

PubMed

10.

G.

Mazur

,

M.

Makowski

, and

R.

Łazarski

, “

Boys function evaluation on graphical processing units

,”

J. Math. Chem.

54

(

10

),

2022

–

2047

(

2016

).

https://doi.org/10.1007/s10910-016-0668-x

Google Scholar

Crossref

11.

G.

Beylkin

and

L.

Monzón

, “

Efficient representation and accurate evaluation of oscillatory integrals and functions

,”

Discrete Contin. Dyn. Syst.

36

(

8

),

4077

–

4100

(

2016

).

https://doi.org/10.3934/dcds.2016.36.4077

Google Scholar

Crossref

12.

T. N.

Rescigno

,

C. W.

McCurdy

, Jr., and

V.

McKoy

, “

Low-energy e⁻-H₂ elastic cross sections using discrete basis functions

,”

Phys. Rev. A

11

(

3

),

825

(

1975

).

https://doi.org/10.1103/physreva.11.825

Google Scholar

Crossref

13.

N. S.

Ostlund

, “

Polyatomic scattering integrals with Gaussian orbitals

,”

Chem. Phys. Lett.

34

(

3

),

419

–

422

(

1975

).

https://doi.org/10.1016/0009-2614(75)85528-x

Google Scholar

Crossref

14.

R.

Colle

,

A.

Fortunelli

, and

S.

Simonucci

, “

A mixed basis set of plane waves and Hermite Gaussian functions. Analytic expressions of prototype integrals

,”

Nuovo Cimento D

9

(

8

),

969

–

977

(

1987

).

https://doi.org/10.1007/bf02464849

Google Scholar

Crossref

15.

R.

Colle

,

A.

Fortunelli

, and

S.

Simonucci

, “

Hermite Gaussian functions modulated by plane waves: A general basis set for bound and continuum states

,”

Nuovo Cimento D

10

(

7

),

805

–

818

(

1988

).

https://doi.org/10.1007/bf02450141

Google Scholar

Crossref

16.

L.

Füsti-Molnar

and

P.

Pulay

, “

Accurate molecular integrals and energies using combined plane wave and Gaussian basis sets in molecular electronic structure theory

,”

J. Chem. Phys.

116

(

18

),

7795

–

7805

(

2002

).

https://doi.org/10.1063/1.1467901

Google Scholar

Crossref

17.

W. P.

Reinhardt

, “

Complex coordinates in the theory of atomic and molecular structure and dynamics

,”

Annu. Rev. Phys. Chem.

33

,

223

(

1982

).

https://doi.org/10.1146/annurev.pc.33.100182.001255

Google Scholar

Crossref

18.

K.

Wolinski

,

J. F.

Hinton

, and

P.

Pulay

, “

Efficient implementation of the gauge-independent atomic orbital method for NMR chemical shift calculations

,”

J. Am. Chem. Soc.

112

(

23

),

8251

–

8260

(

1990

).

https://doi.org/10.1021/ja00179a005

Google Scholar

Crossref

19.

D.

Kahaner

,

C.

Moler

, and

S.

Nash

,

Numerical Methods and Software

(

Prentice-Hall, Inc.

,

1989

).

Google Scholar

20.

W.

Gautschi

, “

Efficient computation of the complex error function

,”

SIAM J. Numer. Anal.

7

(

1

),

187

–

198

(

1970

).

https://doi.org/10.1137/0707012

Google Scholar

Crossref

21.

J.

Boersma

, “

Computation of Fresnel integrals

,”

Math. Comput.

14

(

69-72

),

380

(

1960

).

https://doi.org/10.1090/s0025-5718-1960-0121973-3

Google Scholar

Crossref

22.

G.

Beylkin

,

L.

Monzón

, and

I.

Satkauskas

, “

On computing distributions of products of non-negative independent random variables

,”

Appl. Comput. Harmon Anal.

46

(

2

),

400

–

416

(

2018

); arXiv:1707.07762.

https://doi.org/10.1016/j.acha.2018.01.002

Google Scholar

Crossref

23.

G.

Beylkin

and

L.

Monzón

, “

Approximation of functions by exponential sums revisited

,”

Appl. Comput. Harmon Anal.

28

(

2

),

131

–

149

(

2010

).

https://doi.org/10.1016/j.acha.2009.08.011

Google Scholar

Crossref

2021

Author(s)

A fast algorithm for computing the Boys function

I. INTRODUCTION

II. APPROXIMATION OF $F (0, z)$ FOR COMPLEX VALUED ARGUMENT

A. The case $R e (z) \geq 0$ ⁠.

B. The case $R e (z) < 0$

III. APPROXIMATION OF $F (n_{\max}, z)$ FOR REAL AND COMPLEX ARGUMENTS

IV. IMPLEMENTATION DETAILS

V. CONCLUSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX: CONSTRUCTION OF QUADRATURE IN (11)

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A fast algorithm for computing the Boys function

I. INTRODUCTION

II. APPROXIMATION OF F0,z FOR COMPLEX VALUED ARGUMENT

A. The case Rez≥0⁠.

B. The case Rez<0

III. APPROXIMATION OF Fnmax,z FOR REAL AND COMPLEX ARGUMENTS

IV. IMPLEMENTATION DETAILS

V. CONCLUSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX: CONSTRUCTION OF QUADRATURE IN (11)

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

II. APPROXIMATION OF $F (0, z)$ FOR COMPLEX VALUED ARGUMENT

A. The case $R e (z) \geq 0$ ⁠.

B. The case $R e (z) < 0$

III. APPROXIMATION OF $F (n_{\max}, z)$ FOR REAL AND COMPLEX ARGUMENTS