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.


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

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

  • $\begingroup$ 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. $\endgroup$ Nov 19 '15 at 6:37
  • $\begingroup$ @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). $\endgroup$
    – user23061
    Nov 19 '15 at 9:04
  • $\begingroup$ consider scholar.google.ru/… $\endgroup$
    – permeakra
    Nov 19 '15 at 14:32
  • $\begingroup$ @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). $\endgroup$
    – user23061
    Nov 19 '15 at 14:41
  • $\begingroup$ @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/… $\endgroup$
    – 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


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.

  • 1
    $\begingroup$ 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! ; ) $\endgroup$
    – 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.

  • $\begingroup$ I will certainly look at this book, thank you! $\endgroup$
    – user23061
    Nov 24 '15 at 23:09
  • $\begingroup$ 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. $\endgroup$
    – user23061
    Nov 25 '15 at 0:26
  • $\begingroup$ Check normalizations of you basis functions $\endgroup$
    – Eric Brown
    Nov 25 '15 at 0:29
  • $\begingroup$ Why you suspect about normalization? $\endgroup$
    – user23061
    Nov 25 '15 at 1:00
  • $\begingroup$ Did Sij evaluate to 1? $\endgroup$
    – 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:


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.

  • $\begingroup$ 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. $\endgroup$
    – user23061
    Nov 19 '15 at 9:06
  • $\begingroup$ @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. $\endgroup$ Nov 19 '15 at 14:52
  • $\begingroup$ @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. $\endgroup$
    – user23061
    Nov 19 '15 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy