quantum chemistry - Computing two-electron integrals with an STO-3G basis set

I am trying to implement a restricted Hartree-Fock calculation using an STO-3G basis set, for fun. I managed to perform this calculation where only $\mathrm{1s}$ orbitals are present ($\ce{H2}$ and $\ce{HeH+}$) as explained in Szabo and Ostlund's book. In this book, authors give explicit formulas for overlap, kinetic, nuclear-electron, and electron-electron integrals for $\mathrm{1s}$ orbitals and they work correctly.

In order to generalize my calculation to systems containing $\mathrm{2s}$ and $\mathrm{2p}$ orbitals (for $\ce{H2O}$ and $\ce{N2}$), I used the general formulas I found in Cook's book for the electron-nuclear and electron-electron integrals. In this case, I obtain results that are slightly different from Szabo's book:

$$E_\text{tot}(\ce{H2O}) = -74.4442002765 \text{ a.u.}$$

instead of

$$E_\text{tot}^\text{Szabo}(\ce{H2O}) = -74.963 \text{ a.u.}$$

This is obviously problematic since orbital energies suffer from the same error and this leads to an erroneous ionization potential (0.49289045 a.u. instead of 0.391 a.u., a difference of approximately 63 kcal$\cdot$mol$^{-1}$).

Since I checked my code multiple times and I wrote the two-electron computation code from scratch twice, I was wondering if there is a typo in Cook's book. Is there is a good reference where I can find the (correct) formula to compute two-electron integrals of gaussian functions (in Cartesian coordinates) with arbitrary angular momenta? At the moment I am not looking for a very efficient (recursive) algorithm to perform this task, I only need an exact formula like the one proposed in Cook's book.

Sources:

[1] Szabo and Ostlund, Modern Quantum Chemistry, Dover, 1989.

[2] Cook, Handbook of Computational Chemistry, Oxford University Press, 1998.

Answer

Actually there is a mistake in the analtical expression in Cook's Book. On his web page he has a pdf with the corrected verison

http://spider.shef.ac.uk/

Maybe this solves your problem, but I would also recommend to implement the Obara-Saika Scheme or rys-Quadrature since they are really much more efficient. If your are programming in Python, you might have a look at the PyQuante project, which implements all this stuff. Concerning Obara-Saika you might also read about the Head-Gordon Pople Scheme. It is in principial an adapted version of Obara-Saika which reduces the FLOP count.