quantum chemistry - 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?

Answer

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.