Szabo and Ostlund's calculation of repulsion integral on page 172 via A.41 Returning to Szabo &amp; Ostlund after a hiatus I tried to recreate his STO-3G calculation for HeH+ which begins on p. 168. It all went fine except for the two-electron, two-center integrals. I employed Szabo's values for these on p. 172 and was able to get his result for the energy and so on.
What I was not able to do is duplicate his p. 172 values for the 2-e integrals from his expression on page 416 (A.41). For example, on page 172, he gives $(\phi_{1}\phi_{2}|\phi_{1}\phi_{2}) = .1773.$
There are other approaches to this problem and I have looked at a few, but I would like to understand it as Szabo presents it. Because I was able to get the nuclear repulsion terms and so on I think my expressions for the contracted gaussians are fine but am misinterpreting some aspect of (AB|CD) in A.41.
If anyone is familiar with this particular calculation, and could show the proper calculation of (say) (12|12)...This is a very picayune question in the sense that I am interested in the particular calculation in Szabo: how to get 0.1773 on p. 172? Thanks for any illumination.
Edit
The Fortran program alerted me to a confusion. The &quot;contracted&quot; gaussians, using (11/11) as an example, look like:
$$G_a= .444635(2a_1/\pi)^{3/4}e^{-a_1|r-R_A|^2}+.535328(2a_2/\pi)^{3/4}e^{a_2|r-R_A|^2}+.154239(2a_3/\pi)^{3/4}e^{a_3|r-R_A|^2}$$
Regardless of how we do the integral it is:
$$\int G_a^2 \frac{1}{r_{1,2}}G_a^2$$ and we are raising the coefficients associated with $G_a$ to the 4th power and they come out of the integral. If we could do this integral in the usual naive way using linearity we would have 81 terms (some repeated) involving constants like:
$$(d_id_jd_kd_l)(2a_i/\pi)^{3/4}(2a_j/\pi)^{3/4}(2a_k/\pi)^{3/4}(2a_l/\pi)^{3/4} $$ where $i,j,k,l$ range from 1 to 3 respectively, where the $d_i$ are the coefficients $0.444365,$ etc.
That is why on p. 420 in the fortran program the quantities D1 include weights as factors, and a product of four D1's are pre-multiplied on p. 421 by the constants calculated in &quot;TWOE&quot; on p. 423. Using Fourier transforms to do the integral does not change the constants.
In short, the program is doing a lot of work here. It should still be possible to exhibit (11/11) easily but no wonder A.41 is not used in practice--it represents no economy over the fortran.
Restricting attention to the integrals representing &quot;self-interaction &quot; of helium in HeH+ we have to find coefficients of the integral,
$$\int A\cdot A\frac{1}{r_{12}}\cdot A  \cdot A$$
in which for an STO-3G calculation $A=a\cdot g_1+b\cdot g_2+c\cdot g_3$ with $g_1$ a gaussian and $a$ a coefficient formed as on page 153 sec. 3.203. The coefficients will fall out as a sum of 81 terms, some of which are identical:
$$(a+b+c)^4 = a^4+4a^3b+6a^2b^2+4ab^3+b^4+4a^3c+12a^2bc+12ab^2c+4b^3c+6a^2c^2+12abc^2+6b^2c^2+4ac^3+4bc^3+c^4.$$
That there are 81 terms may be checked by viewing the fortran at page 420 and 421. The expression for V1111 is formed by pre-multiplying by D(i)D(j)D(k)D(l) in which i,j,k,l range from 1 to 3. The coefficients DDDD do not appear explicitly in A.41.
The 12 terms in $ab^2c$ and so on may represent different quantities, and care must be taken to keep the accounting correct because (using notation defined below)
$$12ab^2c\cdot 2\pi^{5/2}/((a_1+a_2)(a_2+a_3)\sqrt{a_1+a_2+a_2+a_3})$$
$$\neq 12ab^2c\cdot 2\pi^{5/2}/((a_2+a_2)(a_1+a_3)\sqrt{a_1+a_2+a_2+a_3}) $$
Letting $\zeta_1 = 2.0925$ (page 170) we have:
$a_1=0.109818\cdot\zeta_1^2,~~ a_2=0.405771\cdot\zeta_1^2,~~ a_3=2.22766\cdot\zeta_1^2; $
$a =0.444635\cdot (2a_1/\pi)^{3/4},~b=0.535328\cdot(2a_2/\pi)^{3/4},~ c= 0.154329\cdot(2a_3/\pi)^{3/4}$
The sum represented in the appendix A.41 for the 2-center 2-electron integral
$$ \int He \cdot He \frac{1}{r_{12}}\cdot He\cdot He$$ is given below, and in Mathematica I get 1.307238. Substituting in the constants for hydrogen I get about 0.77466. These are Szabo and Ostlund's values on page 172. They are a good test because they avoid calculation of the error function used in cross terms and they are close to results for the integration of Slater orbitals using spherical coordinates, which are $(5/8)\cdot \zeta. $
The calculation exhibits the coefficients in A.41 explicitly and is in essence the fortran code for V1111.
$ a^4\cdot 2 \pi^{5/2}/((a_1 + a_1)\cdot(a_1 + a_1)\cdot\sqrt{4\cdot a_1})+$
$4 a^3b\cdot 2 \pi^{5/2}/((a_1 + a_1) \cdot(a_1 + a_2)\cdot\sqrt{3 a_1 + a_2})+ $
$ 2 a^2 b^2\cdot 2 \pi^{5/2}/((a_1 + a_1) (a_2 + a_2)\cdot\sqrt{2 a_1 + 2 a_2}) +$
$ 4 a^2 b^2\cdot2 \pi^{5/2}/((a_1 + a_2)\cdot(a_1 + a_2)\cdot \sqrt{2 a_1 + 2 a_2}) +$
$4 a b^3\cdot 2 \pi^{5/2}/((a_1 + a_2) (a_2 + a_2)\cdot \sqrt{a_1 + 3 a_2}) +$
$ b^4\cdot 2 \pi^{5/2}/((a_2 + a_2) (a_2 + a_2)\cdot \sqrt{2 a_2 + 2 a_2}) +$
$ 4 a^3 c\cdot 2 \pi^{5/2}/((a_1 + a_1)*(a_1 + a_3)\cdot \sqrt{3 a_1 + a_3}) +$
$ 4 a^2 b c\cdot 2 \pi^{5/2}/((a_1 + a_1)\cdot(a_2 + a_3)\cdot \sqrt{2 a_1 + a_2 + a_3}) +$
$ 8 a^2 b c\cdot 2 \pi^{5/2}/((a_1 + a_2)*(a_1 + a_3)\cdot \sqrt{2 a_1 + a_2 + a_3}) +$
$  4 a b^2 c\cdot 2 \pi^{5/2}/((a_1 + a_2) (a_2 + a_3)\cdot \sqrt{a_1 + 2 a_2 + a_3}) +$
$  8 a b^2 c\cdot 2 \pi^{5/2}/((a_2 + a_2) (a_1 + a_3)\cdot \sqrt{a_1 + 2 a_2 + a_3}) +$
$ 4 b^3 c\cdot 2 \pi^{5/2}/((a_2 + a_2)\cdot(a_2 + a_3)\cdot\sqrt{3 a_2 + a_3}) +$
$ 2 a^2 c^2\cdot 2 \pi^{5/2}/((a_1 + a_1)\cdot(a_3 + a_3)\cdot \sqrt{2 a_1 + 2 a_3}) +$
$ 4 a^2 c^2\cdot 2 \pi^{5/2}/((a_1 + a_3)\cdot(a_1 + a_3)\cdot\sqrt{2 a_1 + 2 a_3}) +$
$ 8 a b c^2\cdot 2 \pi^{5/2}/((a_1 + a_2)\cdot(a_3 + a_3)\cdot \sqrt{a_1 + a2 + 2 a3}) +$
$ 4 a b c^2\cdot 2 \pi^{5/2}/((a_1 + a_3)\cdot(a_2 + a_3)\cdot \sqrt{a_1 + a_2 + 2 a_3}) +$
$ 2 b^2 c^2\cdot 2 \pi^{5/2}/((a_2 + a_2)\cdot(a_3 + a_3)\cdot \sqrt{2 a_2 + 2 a_3}) +$
$ 4 b^2 c^2\cdot 2 \pi^{5/2}/((a_2 + a_3)\cdot(a_2 + a_3)\cdot \sqrt{2 a_2 + 2 a_3}) +$
$  4 a c^3\cdot 2 \pi^{5/2}/((a_1 + a_3)\cdot(a_3 + a_3)\cdot \sqrt{a_1 + 3 a_3}) +$
$ 4 b c^3\cdot 2 \pi^{5/2}/((a_2 + a_3)\cdot(a_3 + a_3)\cdot \sqrt{a_2 + 3 a_3}) +$
$ c^4\cdot2 \pi^{5/2}/((a_3 + a_3)\cdot(a_3 + a_3)\cdot \sqrt{4 a_3})$
Some copy-and-paste Mathematica code --don't know if it works in Wolfram Alpha.
z11 = 2.0925; z21 = 1.24;

a1 = .109818*z11^2; a2 = .405771*z11^2; a3 = 2.22766*z11^2; 

a = .444365 (2 a1/Pi)^(3/4); b = .535328 (2 a2/Pi)^(3/
     4); c = .154329 (2 a3/Pi)^(3/4);

a^4*2 Pi^(5/2)/((a1 + a1)*(a1 + a1)*Sqrt[4 a1]) +
 4 a^3 b*2 Pi^(5/2)/((a1 + a1) (a1 + a2)*Sqrt[3 a1 + a2]) +
 2 a^2 b^2*2 Pi^(5/2)/((a1 + a1) (a2 + a2)*Sqrt[2 a1 + 2 a2]) +
 4 a^2 b^2*2 Pi^(5/2)/((a1 + a2)*(a1 + a2)*Sqrt[2 a1 + 2 a2]) +
 4 a b^3*2 Pi^(5/2)/((a1 + a2) (a2 + a2)*Sqrt[a1 + 3 a2]) +
 b^4*2 Pi^(5/2)/((a2 + a2) (a2 + a2)*Sqrt[2 a2 + 2 a2]) +
 4 a^3 c*2 Pi^(5/2)/((a1 + a1)*(a1 + a3)*Sqrt[3 a1 + a3]) +
 4 a^2 b c*2 Pi^(5/2)/((a1 + a1)*(a2 + a3)*Sqrt[2 a1 + a2 + a3]) +
 8 a^2 b c*2 Pi^(5/2)/((a1 + a2)*(a1 + a3)*Sqrt[2 a1 + a2 + a3]) +
  4 a b^2 c*2 Pi^(5/2)/((a1 + a2) (a2 + a3)*Sqrt[a1 + 2 a2 + a3]) +
  8 a b^2 c*2 Pi^(5/2)/((a2 + a2) (a1 + a3)*Sqrt[a1 + 2 a2 + a3]) +
 4 b^3 c*2 Pi^(5/2)/((a2 + a2)*(a2 + a3)*Sqrt[3 a2 + a3]) +
 2 a^2 c^2*2 Pi^(5/2)/((a1 + a1)*(a3 + a3)*Sqrt[2 a1 + 2 a3]) +
 4 a^2 c^2*2 Pi^(5/2)/((a1 + a3)*(a1 + a3)*Sqrt[2 a1 + 2 a3]) +
 8 a b c^2*2 Pi^(5/2)/((a1 + a2)*(a3 + a3)*Sqrt[a1 + a2 + 2 a3]) +
 4 a b c^2*2 Pi^(5/2)/((a1 + a3)*(a2 + a3)*Sqrt[a1 + a2 + 2 a3]) +
 2 b^2 c^2*2 Pi^(5/2)/((a2 + a2)*(a3 + a3)*Sqrt[2 a2 + 2 a3]) +
 4 b^2 c^2*2 Pi^(5/2)/((a2 + a3)*(a2 + a3)*Sqrt[2 a2 + 2 a3]) +
  4 a c^3*2 Pi^(5/2)/((a1 + a3)*(a3 + a3)*Sqrt[a1 + 3 a3]) +
 4 b c^3*2 Pi^(5/2)/((a2 + a3)*(a3 + a3)*Sqrt[a2 + 3 a3]) +
 c^4*2 Pi^(5/2)/((a3 + a3)*(a3 + a3)*Sqrt[4 a3])