# 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.}$$

$$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.

• Can you put your code on github? Difficult to answer the question without seeing your formulae. Your result could be due to a floating point error ... I've found this when reproducing the computations of other scientists. – user1945827 Nov 19 '15 at 6:37
• @user1945827 I was thinking about the same thing, but many calculations (see Geoff's answer) gives a value of -74.96(...). My code is already on GitHub (github.com/RMeli/Hartree-Fock) but I didn't gave the link because I didn't wanted too much help! ;) The routine computing the integrals is the ELECTRONIC in INTEGRALS.PY (the program is in Python because I can code fast, but I will migrate it to fortran when it works fine). – user23061 Nov 19 '15 at 9:04
• consider scholar.google.ru/… – permeakra Nov 19 '15 at 14:32
• @permeakra I spent much time using Google and Google Scholar. However, Google Scholar results are normally about very efficient algorithms to compute two-electron integrals. What I would like to implement at first is an exact formula (like the one proposed in Cook's book or in Szabo's book for 1s orbitals). – user23061 Nov 19 '15 at 14:41
• @R.M. This article is on the first page of the google scholar link from above.... I dunno, seems suspiciously like what your want. rspa.royalsocietypublishing.org/content/royprsa/258/1294/… – permeakra Nov 19 '15 at 15:33

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.

• This DO solve my problem!!! Thank you!!! I will certainly look to Obara-Saika or other efficient schemes, but for the time being the goal is to implement BO-MD and not to have an efficient code! ; ) – user23061 Nov 25 '15 at 23:49

My advice is to implement the Obara-Saika recurrence formulae that are outlined in "Molecular Electronic-Structure Theory" by Helgakar, et al. I would stick with Cartesian functions, since a) they are easier and b) spherical harmonics don't matter for molecules anyway.

I did this years ago (in Mathematica) when I was in a similar place -- having completed Szabo and Ostlund and wanting an arbitrary "angular momentum" code.

The book will cost you a small fortune, but if you can borrow it through Interlibrary Loan, it is a phenomenal text.

• I will certainly look at this book, thank you! – user23061 Nov 24 '15 at 23:09
• Unfortunately your suggestions on Boys function does not solve this problem (as I hoped). The error should be directly in the computation of two-electron integrals. I will look to Obara-Saija recurrence. – user23061 Nov 25 '15 at 0:26
• Check normalizations of you basis functions – Eric Brown Nov 25 '15 at 0:29
• Why you suspect about normalization? – user23061 Nov 25 '15 at 1:00
• Did Sij evaluate to 1? – Eric Brown Nov 25 '15 at 1:01

My suggestion would be to use another existing code and run the calculation.

For example, if I do an HF/STO-3G calculation on $\ce{H2O}$, I get:

$$E_\mathrm{H_2O}=-74.9659011\:\mathrm{a.u.}$$

I don't have Cook's book on hand, so I can't look for the error, but I'd suspect some error as you do.

• It will help if I write down explicitly Cook's formulas? I already checked on internet for other HF/STO-3G results and all of them start with -74.96(...) a.u., this is why I am concerned. – user23061 Nov 19 '15 at 9:06
• @R.M. Just because other more modern techniques are efficient does not mean they are inexact. Feel free to post the formulas, but I'd adapt a different method. – Geoff Hutchison Nov 19 '15 at 14:52
• @Goeff I know there are very good (and efficient) techniques, but this will be the next step in my project. At the beginning I just want to code the simplest algorithm possible, even if it is inefficient. Cook's analytical solution only imply a series of loop, which is extremely easy and straightforward to implement (in addition, this is how I implemented the other integrals). However, if I will not find an analytical solution quickly, I will certainly spend some time to implement a more difficult technique. – user23061 Nov 19 '15 at 15:13