Reference value for two-electron repulsion integral over GTO's

I am currently trying to implement a Full CI program from scratch. The energies I get are a bit too high, so I'm looking for the mistake.

One possibility is my implementation of the two-electron repulsion integrals, $$\int_{-\infty}^\infty \int_{-\infty}^\infty \frac{e^{-\alpha_1 \boldsymbol{r}_1^2} e^{-\alpha_2 \boldsymbol{r}_1^2} e^{-\beta_1 \boldsymbol{r}_2^2} e^{-\beta_2 \boldsymbol{r}_2^2}} {\left| \boldsymbol{r}_1-\boldsymbol{r}_2 \right|\ } \mathrm{d}\boldsymbol{r}_1^3 \mathrm{d}\boldsymbol{r}_2^3$$

with $\boldsymbol{r}_1=(x_1,y_1,z_1), \boldsymbol{r}_2=(x_2,y_2,z_2).$

However, I cannot find any values I could compare my results to. I tried (numerical) integration in Maple, but that's way too slow and/or numerically unstable because of the singularity at $\boldsymbol{r}_1=\boldsymbol{r}_2$. I tried installing the libraries libint and Libcint, and the python module pyscf, which all should be able to do this kind of computations, but I horribly fail at installing things not made for Windows (MINGW is only half working, and I don't have a proper Linux installation available right now...).

So, could someone who has such a program installed please give me a number this integral evaluates to, for whatever numbers $\alpha_1 \ldots \beta_2$ ?

• I successfully installed libint some years ago using cygwin (free posix emulation layer for windows.) – permeakra Dec 11 '15 at 8:20
• Thank you, but the thing is, I already wasted about three days on trying to get some library to work (and still didn't, I know, that's pathetic ^^), and I'm exhausted. All I really need is a single number to compare... Here I think it's more likely someone has a quantum chemistry program installed that's willing to easily output such a number. It would really help a lot. Well, might try the python thing again tomorrow on a friends laptop, if I can borrow it... – Giogina Dec 11 '15 at 11:22
• Try a live DVD/USB if you want to use Linux. – Greg Dec 12 '15 at 17:31

You don't state which method you implemented in order to compute two-electron integrals, therefore I will list all the main references at first.

Cook's book [1] contains analytical formulas for overlap, kinetic, electron-nucleus attraction and electron-electron repulsion integrals. The analytical formula for electron-electron repulsion integrals is wrong in the book, but you can look at this discussion for the errata.

However, a Full CI calculation is computationally costly. For this reason I believe a better approach would be to use a more efficient scheme to compute integrals. You can look at the following:

• Obara-Saika scheme
• McMurchie-Davidson scheme

You can find a good explanation of these methods in Ref. [2].

Now, if you want to check what you already implemented, you can find a list of the two-electron integrals in the STO-3G basis set for $\ce{HeH+}$ in Szabo's book appendix [3]. Here I can give you a set of two-electron integrals I obtain for $\ce{H2}$ (always in the STO-3G basis set) for a bond length of $1.4$ Bohr:

 (           1           1           1           1 )   0.77460834925515787
(           1           1           1           2 )   0.44410904384277344
(           1           1           2           1 )   0.44410904384277350
(           1           1           2           2 )   0.56967771733030592
(           1           2           1           1 )   0.44410904384277361
(           1           2           1           2 )   0.29702946944511982
(           1           2           2           1 )   0.29702946944511982
(           1           2           2           2 )   0.44410904384277333
(           2           1           1           1 )   0.44410904384277333
(           2           1           1           2 )   0.29702946944511982
(           2           1           2           1 )   0.29702946944511982
(           2           1           2           2 )   0.44410904384277361
(           2           2           1           1 )   0.56967771733030592
(           2           2           1           2 )   0.44410904384277350
(           2           2           2           1 )   0.44410904384277344
(           2           2           2           2 )   0.77460834925515787


Note that these results came from a program I wrote myself, but it usually match very well Szabo's [3] and Gaussian09 values for the total energy. If these results match your calculations, then the problem might be in integrals with higher angular momentum.

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

[2] T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-Structure Theory, Wiley, 2000.

[3] A. Szabo and N. Ostlund, Modern Quantum Chemistry, Dover, 1996.

• Awesome, thank you very much! My integrals are indeed all about 0.02 too high, enough to explain the difference in my final result. I implemented the analytical solution, and you are right, the computational cost is incredible. I will look into the integration schemes you mentioned! – Giogina Dec 11 '15 at 14:35
• @Giogina You can stick with the analytical implementation in order to have a working program at first. Make sure you implemented well your calculation (Cook's correct analytical formula is quite straightforward) and make sure your problem does not comes from Boys function! – user23061 Dec 12 '15 at 12:41
• Hm, on the first glimpse my re-derivation of the formula does look very similar to Cook's one, but with 15 summations it's so easy to have a small mistake somewhere, I'll probably just have to rewrite the whole thing. I trust my boys function, since I'm still in Maple, it should do that well. – Giogina Dec 12 '15 at 14:37