# Boys function for Gaussian integrals in ab initio calculations

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?

• I recommend asking questions about code and/or libraries at the computational science SE. I'm not saying is off topic, I just think it's better suited there. Commented Nov 27, 2015 at 1:52
• I just found a library called "LibInt" for computing molecular integrals: sourceforge.net/projects/libint It contains a file 'Boys.h', so what you want seems to be implemented. I havn't tested it yet (I'm currently working on writing a Full CI code, so I keep running into problems similar to yours ^^), but it's used in several quantum chemistry packages, so I'm sure it'll give good results. Commented Dec 4, 2015 at 23:53
• Thank you for the hint. I managed to implement this computation using the incomplete gamma function (as suggested in the answer) contained in SLATEC, but I will certainly give a look at the library!
– user23061
Commented Dec 5, 2015 at 20:23
• As an insider to ORCA, I can state that the Boys function is computed by various methods depending on the argument and that higher orders are computed recursively. I can also recommend looking at open-source implementations such as NWChem for inspiration on this. Commented Mar 12, 2017 at 12:08

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


The simplest option is a quadrature of

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

using QuadGK

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.

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

• What a great response. Thank you!
– Rick
Commented Jul 6, 2023 at 22:28
• @Rick glad it helped! Commented Jul 7, 2023 at 23:17

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  *)


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.

• I already found this analytical formula in one publication, but the incomplete gamma function is not part of Fortran standard library. I will look for Python, though. The values of n I have at the moment range from 0 to 3, therefore I don't think an overflow will be a problem.
– user23061
Commented Nov 24, 2015 at 23:04
• Commented Nov 24, 2015 at 23:12
• I don't know how it is possible I did not saw it! Thank you very much!
– user23061
Commented Nov 25, 2015 at 0:09
• Unfortunately this does not solve the mentioned total energy problem... =(. I will have to search elsewhere.
– user23061
Commented Nov 25, 2015 at 0:24

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.

• the leading subscript 1 in front of F appears to be a typo in the WIREs paper. StackExchange does not allow me to correct this because the edit is less than 6 characters. Commented May 9, 2017 at 18:38
• @Jeff Do you have a citation? I don't see any typo here. The relationship to Kummer's function 1F1 and the Boys function Fn is derived in Sec (9.8.6) of Helgaker, Jorgensen, and Olsen; Molecular electronic-structure theory; (2000). Commented May 10, 2017 at 16:42
• I see now. ${}_{1}F_{1}$ is a completely different $F$ from $F_{1}$. How nice. Commented May 11, 2017 at 18:40
• @Jeff Correct, $_1F_1$ is the Kummer confluent hypergeometric function, not the Boys function. The Boys function $F_n$ is a weighted version of $_1F_1$. I agree, the notation can be very confusing here. Commented May 11, 2017 at 18:41
• I think that the CHGM have troubles calculating the Boys function for what is needed for Molecular Integrals. It performs well enough for H2O/6-31++G**. But for bigger molecules it seem to fail p-nitrophenolate/4-31G. Do you know any other implementations of Boys function in Fortran? Commented Aug 21, 2017 at 18:35

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 .

• While the asker may find this useful, note that the OP requests implementations of python and/or fortran specifically. Commented Nov 23, 2017 at 10:25
• It is worth stating that it is roughly as easy to call C from Python as it is to call Fortran from Python, so this is useful in that regard. Less useful is source code embedded in a PDF. Commented Nov 23, 2017 at 15:53