By completing squares Gaussian integrals over all space can be reduced to the same basic form. This enables the analytic evaluation of forms containing rij. These will then be compared with Monte-Carlo estimates of the same quantities in order to test the accuracy of various Monte-Carlo estimates.
The integral form considered here is
The basic procedure of completing the squares to integrate
(2)
is in Intemp.doc .htm.
Expanding the double sum in equation (1) yields
(3)
so that
(4)
or defining
(5)
(6)
Sorting out the last integration variable
(7)
Then using the integration template the last integral becomes
(8)
Putting the exponential with the similar term in the integral over r1 yields
(9)
Simplify the second exponential
(10)
And now recognize the integral as
(11)
so that we have
(12)
FUNCTION H(A,BETA0,N)
BETA=BETA0 -- INPUT VALUE
FACT=1
DO NM=N,3,-1
FACT=FACT*
BETA=A*BETA/(A-2*BETA)
ENDDO
For N>2, the sequence is independent of m and n.
H=FACT*H(A,beta,2,m,n)
The bottom layer for m and n = 0 is found to be
in H200.doc .htm. The Fortran code is in testing\hfuns.for.
The worry comes when evaluating terms that could be singular. The first of these is the nuclear electron potential 1/R1 and its square (part of p4) which goes as 1/R12. The change needed in equation (1)is only to the r1 integral
(14)
The r2 integration is exactly the same as in H200.doc .htm leading to equation 6. Thus
(15)
Change to spherical coordinates to take advantage of the lack of angular dependency in the integrand.
(16)
The cases of interest are m=0,1 and 2
(17)
In practice, both of these will be calculated on each test run and they will be summed over all N electrons.
NMC= 1048576
THE NUMBER OF ELECTRONS IS 2
INTEGRAL IS 4.448057117553909E-001 +- 1.676313541383681E-003
AINTEGRAL IS 4.423603938062348E-001 RES= 1.458747357691631
RI= 1.376489907610462 +- 3.932147188571423E-003
ARI= 1.372032810009793 RES= 1.133502228406804
RI2= 3.345312960514719 +- 7.052165793538759E-003
ARI2= 3.342278530980440 RES= 4.302833516844589E-001
The change to r1, r2 and r12 in the exponent is detailed in Htor12.doc .htm. The result is
where
(19)
The variable change from to is detailed in r12b.doc .htm.
(20)
The case n=0, the normalization which is equal to the results derived without this variable change is detailed in r1r2x.doc .htm to arrive at
with equation (22) in explicit agreement with equation (13).
At this point note that
(23)
This means that
(24)
Using the value of Ω(A,B,0) from (21)
(25)
This integral is given by Jefferey[1] as
(26)
In our integral the value at the lower limit is zero, while the infinity at the upper limit gives
(27)
so that
(28)
For n=2, we use the relationship
(29)
This leads to
(30)
Again using the value of Ω(A,B,0) from (21)
(31)
Equation (18) can be used to put the n=1 and n=2 results into the H form needed by the code.
(32)
and
(33)
The numerical integrals are in testing\testing.mdp (This assumes Power Station Fortran)
The list below is of the main codes. These will run in watfor or watcom. They have includes to the necessary subroutines.