The Boys function, the integral behind Gaussian basis sets

Evaluate the Boys functions $F_{n} (T) = \int_{0}^{1} t^{2 n} e^{- T t^{2}} d t$ with vp.special.boys, which returns orders 0 through n at once, and measure their relative accuracy against mpmath at 50 digits over $T$ from 0 to 200. Every electron-repulsion and nuclear-attraction integral over Gaussian orbitals reduces to these functions, so a Hartree-Fock code calls them millions of times. This one is lib/hf/boys.viv, the routine Viv's Hartree-Fock library computes its integrals with.

The functions

$F_{n} (0) = 1/ (2 n + 1)$ , and for large $T$ ,

F_{n} (T) \approx \frac{( 2 n - 1 )!!}{2 ^{n + 1}} \frac{π}{T ^{2 n + 1}} .

Higher orders fall off faster.

import numpy as np
import vipyr as vp
from vipyr.figure import figure

Ts = np.concatenate([np.linspace(0.0, 2.0, 21), np.geomspace(2.2, 200.0, 60)])
F = np.array([vp.special.boys(n_max=16, T=float(T)).F for T in Ts])
print(
    "F_n(0):",
    np.round(F[0, :5], 6),
    "  1/(2n+1):",
    np.round(1 / (2 * np.arange(5) + 1.0), 6),
)
fig = figure("Boys functions $F_n(T)$", "$T$", "$F_n(T)$")
for n in (0, 2, 4, 8, 16):
    fig.plot(Ts[1:], F[1:, n], "n = %d" % n)
    fig.style("n = %d" % n, "nomarker")
fig.xscale("log")
fig.yscale("log")

F_n(0): [1.       0.333333 0.2      0.142857 0.111111]   1/(2n+1): [1.       0.333333 0.2      0.142857 0.111111]

Accuracy against mpmath

The reference is $F_{n} (T) = Γ (n + \frac{1}{2}) P (n + \frac{1}{2}, T) / (2 T^{n + 1/2})$ , with $P$ the regularized lower incomplete gamma function, at 50 digits. The program's tolerance is 5 × 10⁻¹⁵ relative, worst over every order.

import mpmath

mpmath.mp.dps = 50


def boys_mp(n, T):
    if T == 0:
        return mpmath.mpf(1) / (2 * n + 1)
    a = mpmath.mpf(n) + mpmath.mpf(1) / 2
    return mpmath.gammainc(a, 0, T) / (2 * mpmath.mpf(T) ** a)


orders = (0, 4, 8, 16)
rel = {
    n: [
        float(
            abs(F[i, n] - boys_mp(n, mpmath.mpf(float(T))))
            / boys_mp(n, mpmath.mpf(float(T)))
        )
        for i, T in enumerate(Ts)
    ]
    for n in orders
}
fig = figure("Relative error against mpmath, 50 digits", "$T$", "relative error")
for n in orders:
    fig.plot(Ts[1:], np.maximum(rel[n][1:], 1e-18), "n = %d" % n)
    fig.style("n = %d" % n, "noline circle")
fig.hline("5e-15", 5e-15)
fig.xscale("log")
fig.yscale("log")
for n in orders:
    print(
        "n = %2d   worst relative error over %d values of T: %.1e"
        % (n, len(Ts), max(rel[n]))
    )

n =  0   worst relative error over 81 values of T: 1.2e-15
n =  4   worst relative error over 81 values of T: 1.1e-15
n =  8   worst relative error over 81 values of T: 8.7e-16
n = 16   worst relative error over 81 values of T: 9.0e-16

How it's checked

check() reruns the recorded cases against the same mpmath formula; the tolerance comes from lib/hf/boys's own measurement, 1.7 × 10⁻¹⁵ on its table.

print("check():", vp.special.boys(n_max=16, T=2.0).check().passed)

check(): True

Source

gallery/science/plot_boys_function.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is lib/hf/gallery/boys.viv: its atlas, the defs and the calls between them, or open it in the workbench.