Szabo and Ostlund's calculation of repulsion integral on page 172 via A.41

Restricting attention to the integrals representing "self-interaction " of helium in HeH+ we have to find coefficients of the integral,

$$int Acdot Afrac{1}{r_{12}}cdot A cdot A$$

in which for an STO-3G calculation $A=acdot g_1+bcdot g_2+ccdot 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^2ccdot 2pi^{5/2}/((a_1+a_2)(a_2+a_3)sqrt{a_1+a_2+a_2+a_3})$$

$$neq 12ab^2ccdot 2pi^{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.109818cdotzeta_1^2,~~ a_2=0.405771cdotzeta_1^2,~~ a_3=2.22766cdotzeta_1^2; $

$a =0.444635cdot (2a_1/pi)^{3/4},~b=0.535328cdot(2a_2/pi)^{3/4},~ c= 0.154329cdot(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 Hecdot 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^4cdot 2 pi^{5/2}/((a_1 + a_1)cdot(a_1 + a_1)cdotsqrt{4cdot a_1})+$

$4 a^3bcdot 2 pi^{5/2}/((a_1 + a_1) cdot(a_1 + a_2)cdotsqrt{3 a_1 + a_2})+ $

$ 2 a^2 b^2cdot 2 pi^{5/2}/((a_1 + a_1) (a_2 + a_2)cdotsqrt{2 a_1 + 2 a_2}) +$

$ 4 a^2 b^2cdot2 pi^{5/2}/((a_1 + a_2)cdot(a_1 + a_2)cdot sqrt{2 a_1 + 2 a_2}) +$

$4 a b^3cdot 2 pi^{5/2}/((a_1 + a_2) (a_2 + a_2)cdot sqrt{a_1 + 3 a_2}) +$

$ b^4cdot 2 pi^{5/2}/((a_2 + a_2) (a_2 + a_2)cdot sqrt{2 a_2 + 2 a_2}) +$

$ 4 a^3 ccdot 2 pi^{5/2}/((a_1 + a_1)*(a_1 + a_3)cdot sqrt{3 a_1 + a_3}) +$

$ 4 a^2 b ccdot 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 ccdot 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 ccdot 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 ccdot 2 pi^{5/2}/((a_2 + a_2) (a_1 + a_3)cdot sqrt{a_1 + 2 a_2 + a_3}) +$

$ 4 b^3 ccdot 2 pi^{5/2}/((a_2 + a_2)cdot(a_2 + a_3)cdotsqrt{3 a_2 + a_3}) +$

$ 2 a^2 c^2cdot 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^2cdot 2 pi^{5/2}/((a_1 + a_3)cdot(a_1 + a_3)cdotsqrt{2 a_1 + 2 a_3}) +$

$ 8 a b c^2cdot 2 pi^{5/2}/((a_1 + a_2)cdot(a_3 + a_3)cdot sqrt{a_1 + a2 + 2 a3}) +$

$ 4 a b c^2cdot 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^2cdot 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^2cdot 2 pi^{5/2}/((a_2 + a_3)cdot(a_2 + a_3)cdot sqrt{2 a_2 + 2 a_3}) +$

$ 4 a c^3cdot 2 pi^{5/2}/((a_1 + a_3)cdot(a_3 + a_3)cdot sqrt{a_1 + 3 a_3}) +$

$ 4 b c^3cdot 2 pi^{5/2}/((a_2 + a_3)cdot(a_3 + a_3)cdot sqrt{a_2 + 3 a_3}) +$

$ c^4cdot2 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])