# Chemistry - Boys function for Gaussian integrals in ab-initio calculations

## Solution 1:

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.

## Solution 2:

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.

## Solution 3:

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)