# Maximum required order for the Boys function

I am planning to pretabulate a lookup table in order to evaluate the Boys function. We can expand the Boys function as:

$$F_n (x_t + \Delta x) = \sum _{k=0}^\infty \frac{F_{n+k}(x_t)(-\Delta x)^{k}}{k!} \hspace{16pt} \tag{9.8.12}$$

We can then use downward recursion (upward is numerically unstable) in order to evaluate the function for any value of $n$ (order):

$$F_n (x) = \frac{2xF_{n+1}(x)+\exp(-x)}{2n+1} \tag{9.8.14}$$

I would like somebody to explain how one chooses which value of $n$ the lookup table will take? My working/logic so far as:

• The Hermite integral has the following recursion (McMurchie-Davidson scheme, $x$ direction):

$$R^{n}_{t+1,u,v} = tR^{n+1}_{t-1,u,v} + X_\mathrm{PC}R^{n+1}_{t,u,v} \tag{9.9.18}$$

$\hspace{16pt}$ with:

$$\hspace{16pt} R^{n}_{0,0,0}=(-2p)^{n} F_{n}(p R_\text{PC}^{2}) \tag{9.9.14}$$

$\hspace{16pt}$ and so the recursion should limit $n$ to $n\leq t+u+v$.

• In 2-electron integrals, the Hermite integrals are bound by $t \leq l_\mu + l_\lambda + l_\nu + l_\sigma$. Therefore $n$ is bound by the sum of the total angular momenta of the basis functions, i.e. $n \leq L_\mu + L_\lambda + L_\nu + L_\sigma$ where $L_\mu = l_\mu + m_\mu + n_\mu$.

• Therefore, considering the extreme case of four $L=3$ basis functions, $n$ should be bound by $n \leq 12$.

Does this mean that a lookup table should start at $n=12$ and always recur down from this, or should I make different lookup tables and assess which one to use depending on the largest angular momentum functions in each individual computation? Are there any flaws or holes in my above logic?

All equations taken from Molecular Electronic Structure Theory, Helgaker et al. (2000).

I think your logic seems fine.

When you make your look up tables, you want to make a table for each $n$ up to some value of $n_\text{max}$. Consider you want to calculate an integral that only requires $n=2$, then it would be very inefficient to calculate $F_{12}$ and then down recur.

In your equation (9.8.12) we can notice that in the sum we have the Boys function as $F_{n+k}$. According to "Molecular Electronic Structure Theory", we might need the first $6$ terms of the sum. I.e. when doing your tabulation you will need to calculate all $n$ up to $n_\text{max}+6$.

Also when you do your tabulation you need a maximum range of your tabulated Boys function, i.e. a maximum value for where you can start using the long range approximation afterwards (else your tables will be HUGE). The long range approximation is given as:

$$F_n(x)\leq\frac{(2n-1)!!}{2^{n+1}}\sqrt{\frac{\pi}{x^{2n+1}}}\tag{9.8.10}$$

For which values of $x$ to start using this approximation for different $n$ you can take a look at Table II in Efficient recursive computation of molecular integrals over Cartesian Gaussian functions. I am sure newer litterature on this topic exist, but I haven't encountered it.

As a side note, you can also calculate the Boys function with out tabulation (this might be less effective). By using the following formula:

$$F_n(x)=\exp(-x)\sum_{i=0}^\infty\frac{(2n-1)!!(2x)^i}{(2n+2i+1)!!}$$

This formula makes an implementation of the Boys functions much easier. If you are making an integral code from scratch, I would suggest to implement the above formula. You can always speed up the code later. Note with the above formula you can also just calculate the Boys function for the maximum $n$ you need and then use downwards recursion. If your goal is to implement the Boys function via. tabulation, then you can also use the above formula to generate the values for your look-up tables.

• Hi Erik - SlowQuant was one of my original inspirations to write my own program. I am currently calculating the Boys function with the simple implementation you mention but I read that a lookup table offers significant optimisation. Your explanation answers my question. Helgaker says that $\Delta x = 0.1$ is sufficient - should this be made lower at very small values of $x$? Jun 20, 2018 at 20:04
• I think $\Delta x = 0.1$ is sufficient for the entire interval (I am not 100 % sure). Very glad that it could be an inspiration! For an actual implementation of the Boys function (that is also readable) you can take a look at github.com/theochem/horton if you need inspiration. They have a Python generator script github.com/theochem/horton/blob/master/tools/boys.py And then the actual Boys function in C++: github.com/theochem/horton/blob/master/horton/gbasis Look for files called boys_(something) Jun 20, 2018 at 20:15
• By the way, I think you 9.8.10 is wrong - the double factorial is $(2n-1)$. I can't edit because it's not more than 6 characters. Jun 21, 2018 at 16:26
• Yes, you are right, changed it Jun 21, 2018 at 20:02
• It just means to the power of $i$ @HemanthHaridas Feb 10, 2022 at 13:17

I know this question is pretty old, but for the benefit of people who would come here looking to evaluate boys function, I attach an implementation of the formula mentioned above. The reference implementation is in Fortran, and requires no additional libraries.

  SUBROUTINE BOYS(N,T,RES)
INTEGER, INTENT(IN)           ::  N
DOUBLE PRECISION, INTENT(IN)  ::  T
DOUBLE PRECISION, INTENT(OUT) ::  RES
INTEGER                       ::  CNT1
DOUBLE PRECISION              ::  ERROR, AUX1, AUX2, AUX3, PRES

CNT1  = 0
ERROR = 1.0
CALL DOUBLEFACTORIALD(2*N-1, AUX1)
CALL DOUBLEFACTORIALD(2*N+1, AUX2)
ERROR = (AUX1/AUX2)
RES   = (AUX1/AUX2)
PRES  = RES
DO WHILE (ABS(ERROR) .GT. 1.0D-15)
CNT1  = CNT1+1
CALL DOUBLEFACTORIALD(2*N+2*CNT1+1, AUX2)
AUX3  = (AUX1/AUX2)*((2*T)**CNT1)
RES   = RES+AUX3
ERROR = RES-PRES
WRITE(*,*)  RES, PRES, ERROR, CNT1
PRES  = RES
END DO
RES = EXP(-1*T)*RES
! WRITE(*,*)  ERROR
END SUBROUTINE BOYS


This code reproduced the results provide here for Electron-Nuclear Interaction integrals in water at STO-3G level. The calls to DOUBLEFACTORIALD is a subroutine that calculates the double factorial of a number and returns it as a double.