Boys function for Gaussian integrals in ab initio calculations

Question

A few days ago I mentioned a problem with an Hartree-Fock program I am writing (HF using cartesian Pople's STO-3G basis set). I can reproduce overlap and kinetic integrals of some references for ($\ce{H2}$, $\ce{HeH+}$ and $\ce{H_2O}$). Things start to change (at the third/fourth decimal digits) in the nucleus-electron attraction integral. At the beginning I thought the error I get on the total energy was mainly caused by electron-electron repulsion integrals (as errors on nucleus-electron attraction seem minimal). But now I am quite convinced the problem comes from my calculation of the Boys function $$ F_n(x) = \int_0^1 t^{2n}e^{-xt^2}\, dt. $$ A clue of this comes from the fact that for $\ce{H2}$ and $\ce{HeH+}$ everything works fine: in this case Boys function is only computed for $n = 0$.

In order to compute this function I used numerical integration routines of SciPy (Python library) or directly QUADPACK (I traduced my Python program into Fortran90). However, looking at the literature it seems that nobody compute Boys function via numerical integration. I tried to implement different versions (recursive, asymptotic values, ...) of Boys function but all of them eventually had problems (depending on the range of $x$ or on the value of $n$).

Since I am loosing too much time on this problem, which seems to be numerical, I was wondering if there is some standard implementation of Boys function that I can simply use in my code.

Do you know any Python or Fortran library containing a good (safe and efficient) implementation of Boys function? If not, there is a simple algorithm that works well in any range of $x$ and is not sensible to numerical errors?

Answer 1

I know this is an old question, but I would like to give a small comparison regarding efficiency when evaluating the Boys function $F_n(x)$. Below are some implementations (in Julia) with simple benchmarks.

In the end, I also give some useful approximations to the Boys function, both for small and large values of $x$.

---

For benchmarking, I'll use BenchmarkTools to measure the average time it takes for a particular implementation to calculate $F_1(x)$ $\forall x \in \hat{x}$, where $\tilde{x}$ is the 10000-vector below:

using BenchmarkTools

n = 1
x̃ = range(1e-15, 15, length=10000)

Implementations are in order, from fastest to slowest. (Of course, benchmarks should be compared with each other; different computers will produce different results.)

Evaluations

Using the incomplete gamma function

As mentioned by Eric above, $F_n(x)$ can be written in terms of the incomplete gamma function as

$$\frac{1}{2 x^{\frac{1}{2} + n}} \Gamma \left( \frac{1}{2} + n \right) \gamma \left( \frac{1}{2} + n, x \right)$$

using SpecialFunctions

# Incomplete gamma function implementation.
boys2(n, x) = gamma(0.5 + n) * gamma_inc(0.5 + n, x, 0)[1] / (2x^(0.5 + n))
@benchmark boys2.(n, x̃)

BenchmarkTools.Trial: 
  memory estimate:  406.52 KiB
  allocs estimate:  1005
  --------------
  minimum time:     1.317 ms (0.00% GC)
  median time:      1.345 ms (0.00% GC)
  mean time:        1.380 ms (0.33% GC)
  maximum time:     2.880 ms (0.00% GC)
  --------------
  samples:          3614
  evals/sample:     1

Using quadratures

The simplest option is a quadrature of

$$\int_0^1 t^{2n} e^{-x t^2} dt$$

using QuadGK

# Adaptive Gauss-Kronrod quadrature implementation.
boys1(n, x) = quadgk(t -> t^(2n) * exp(-x * t^2), 0.0, 1.0)[1]
@benchmark boys1.(n, x̃)

BenchmarkTools.Trial: 
  memory estimate:  8.32 MiB
  allocs estimate:  326405
  --------------
  minimum time:     23.747 ms (0.00% GC)
  median time:      24.629 ms (0.00% GC)
  mean time:        25.839 ms (3.78% GC)
  maximum time:     45.680 ms (0.00% GC)
  --------------
  samples:          194
  evals/sample:     1

Using the confluent hypergeometric function

Yet another relation, mentioned by Joshua above, uses the Kummer confluent hypergeometric function in

$$F_n(x) = \frac{{_1}F_1 \left( \frac{1}{2} + n, \frac{3}{2} + n, -x \right)}{2n + 1}$$

using using HypergeometricFunctions

# Confluent hypergeometric function implementation.
boys3(n, x) = pFq([0.5 + n], [1.5 + n], -x) / (2n + 1)
@benchmark boys3.(n, x̃)

BenchmarkTools.Trial: 
  memory estimate:  18.03 MiB
  allocs estimate:  162806
  --------------
  minimum time:     58.511 ms (0.00% GC)
  median time:      59.524 ms (0.00% GC)
  mean time:        60.007 ms (1.19% GC)
  maximum time:     65.223 ms (0.00% GC)
  --------------
  samples:          84
  evals/sample:     1

---

Approximations

Small values of $x$

We can expand the integrand $t^{2n} e^{-x t^2}$ as a series around $x = 0$ and obtain successive polynomial approximations to the Boys function. Below is a systematic way of doing that (using SymPy in Julia):

using SymPy

I = t^(2n) * exp(-x * t^2)

order = 2  # The order of the final polynomial.
I_taylor = series(I, x, 0, order)
B = SymPy.integrate(I_taylor, (t, 0, 1)) |> simplify

$$B^2_n(x) = \frac{3 + 2n - (2n + 1) x}{(2n + 1) (2n + 3)} + O(x^2)$$

Such approximations can be used for small values of $x$ (they are safer for $x \approx 0$ than some implementations above). An approximation of order 10 is compared below with a quadrature of the Boys function.

Large values $x$

Since the Boys function is decreasing, we could, for large $x$, approximate

$$\int_0^1 t^{2n} e^{-x t^2} \approx \int_0^\infty t^{2n} e^{-x t^2} = \frac{1}{2x^{\frac{1}{2} + n}} \Gamma \left( \frac{1}{2} + n \right) = \frac{(2n - 1)!!}{2^{n + 1}} \sqrt{\frac{\pi}{x^{2n + 1}}}$$

(I learned this trick here.) Compare this to the Boys function below.

---

Further reading

The above are obviously toy examples. Here are some nice reads:

• Journal of Mathematical Chemistry 54, 2022-2047 (2016)
• Journal of Computational Chemistry 36, 1390-1398 (2015)
• Journal of Mathematical Chemistry 40, 179-183 (2006)

Answer 2

I am not aware of any existing Fortran code for direct numerical quadrature of this problem, but it is worth pointing out that Mathematica can perform this integral symbolically:

Integrate[t^(2 n) Exp[-x t^2], {t, 0, 1}]

(* 1/2 x^(-(1/2) - n) (Gamma[1/2 + n] - Gamma[1/2 + n, x]) *)

where the $\Gamma$ function can be computed by exponentiating easy-to-find $\ln \Gamma$ functions, often called lngamma(). (One of them is the Euler function $\Gamma(z)$ and the other incomplete $\Gamma(a,z)$ function).

As a check:

With[
  {x = 1, n = 2},
  1/2 x^(-(1/2) - n) (Gamma[1/2 + n] - Gamma[1/2 + n, x])
  ] // N

(* 0.100269  *)

and numerical quadrature yields:

With[
 {x = 1, n = 2},
 NIntegrate[t^(2 n) Exp[-x t^2], {t, 0, 1}]
 ]

(* 0.100269  *)

One has to be a little bit careful, as the $\Gamma$ function is closely related to the factorial, and it can explode as the arguments get too large. I would experiment and see whether you can get by with this analytic solution, which essentially takes you from numerical quadrature into evaluation of special functions.

The $\Gamma$ function approach is faster by over two orders of magnitude than the quadrature approach. This is critical when evaluating molecular integrals. I have seen implementations where the Boys function is actually interpolated rather than evaluated, and there is discussion about the error in interpolation being almost insignificant if suitable interpolation basis functions are used.

Answer 3

The Boys function is $F_n(x)$ a special case of the Kummer confluent hypergeometric function $M(a,b,x) = {_1}F_1(a,b,x)$, which can be found in many special function libraries, such as scipy.special. According to equation (9) of this paper, the relationship is

$F_n(x) = \frac{{_1}F_1(n+\frac{1}{2},n+\frac{3}{2},-x)}{2n +1}$

If you use scipy you can get ${_1}F_1$ here.

Otherwise, you can grab a FORTRAN implementation here. I believe the routine you want is CHGM. There are probably other implementations.

Answer 4

I published an implementation in C/C++ in http://vixra.org/abs/1709.0304 which uses Gauss-Jacobi quadratures. It is limited to quantum numbers of 130 and to real arguments; a relative accuracy near 1.e-15 is achieved. The source code is available in https://zenodo.org/record/3603249 .