# 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. – Martin - マーチン Nov 27 '15 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. – Giogina Dec 4 '15 at 23:53
• While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review – ringo Dec 5 '15 at 1:07
• That is true. As I said, I'm still in the process of figuring this out myself, so I will probably expand on my answer later. – Giogina Dec 5 '15 at 6:41
• 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 Dec 5 '15 at 20:23

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 Nov 24 '15 at 23:04
• – Eric Brown Nov 24 '15 at 23:12
• I don't know how it is possible I did not saw it! Thank you very much! – user23061 Nov 25 '15 at 0:09
• Unfortunately this does not solve the mentioned total energy problem... =(. I will have to search elsewhere. – user23061 Nov 25 '15 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. – Jeff May 9 '17 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). – jjgoings May 10 '17 at 16:42
• I see now. ${}_{1}F_{1}$ is a completely different $F$ from $F_{1}$. How nice. – Jeff May 11 '17 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. – jjgoings May 11 '17 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? – Erik Kjellgren Aug 21 '17 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.

• While the asker may find this useful, note that the OP requests implementations of python and/or fortran specifically. – airhuff Nov 23 '17 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. – pentavalentcarbon Nov 23 '17 at 15:53