The classical expectation value of the potential is defined to be

(0.1)

All coordinates are equivalent in this definition leading to

(0.2)

Define the two-body correlation function to be

(0.3)

So that

(0.4)

Let **r = r _{1} r_{2}**

_{
MPSetEqnAttrs('eq0005','',3,[[180,31,-3,-1,-1],[240,42,-3,-1,-1],[300,52,-4,-1,-1],[270,47,-4,-1,-1],[360,63,-5,-1,-1],[451,78,-6,-2,-2],[750,129,-11,-3,-3]])
MPEquation()
}(0.5)

Define the single variable g to be the average of (0.3)
over **r _{2} **for

_{
MPSetEqnAttrs('eq0006','',3,[[137,30,-3,-1,-1],[182,39,-3,-1,-1],[228,50,-4,-1,-1],[205,44,-4,-1,-1],[274,59,-5,-1,-1],[344,73,-6,-2,-2],[572,122,-11,-3,-3]])
MPEquation()
}(0.6)

Most inter-particle potentials depend only on the magnitude
of **r, **thus define

_{
MPSetEqnAttrs('eq0007','',3,[[152,32,-3,-1,-1],[205,42,-4,-1,-1],[257,52,-5,-1,-1],[230,47,-4,-1,-1],[307,63,-6,-1,-1],[385,78,-7,-2,-2],[637,129,-12,-3,-3]])
MPEquation()
}(0.7)

This is crudely related to the scattering function in ..\definitions\gofr\IbachLuth.htm. A prescription for finding g(r) using ** Contains HW 22 an interesting exercise of medium difficulty.**

Finally

_{
MPSetEqnAttrs('eq0008','',3,[[133,32,-3,-1,-1],[177,43,-4,-1,-1],[223,54,-5,-1,-1],[200,48,-5,-1,-1],[268,65,-6,-1,-1],[334,81,-8,-2,-2],[555,134,-13,-3,-3]])
MPEquation()
}(0.8)

**The
Ornstein Zernike equation**

(*Proc. Acad. Sci. Amsterdam, ***17, **793 (1914))

Define the function
h(**r**) to be

_{
MPSetEqnAttrs('eq0009','',3,[[74,14,-3,-1,-1],[100,18,-4,-1,-1],[124,23,-5,-1,-1],[111,21,-4,-1,-1],[148,28,-6,-1,-1],[186,35,-7,-2,-2],[310,60,-12,-3,-3]])
MPEquation()
}(1.1)

The function g(r) is the two-body correlation function, the probability of finding a particle a distance r from another particle at 0. This probability in a liquid becomes one for large enough r. The Ornstein Zernike equation defines the direct correlation function c(r)

(1.2)

This says that the probability of finding a second particle at a distance r from one at the origin differs from one by the direct correlation due to the particle at the origin plus the integral of the direct correlation from all of the other particles.

This would be just a defining relation for c(r), were it not for the fact that there are many reasonable definitions of c(r). My favorite way of defining c(r) is to calculate g(|r|) inside the usual periodically repeated box of side R, ..\definitions\gofr\gofr2.htm

R

+ +

+ + R + +

using

h(r)=h_{MC}(r) r<R/2

c(r)=0 r>R/2

This is only one of many possible ways to define c(r), others have names such as Percus Yevick and Hyper-Netted-Chain or BBKGY depending on the names of their originators or the shapes of the graphs used to arrive at the approximations.

As can be
seen to the right the 5 values of h computed by Monte-Carlo are the
result of 5 integrals over all r’ of the function c(r’)h(r-r’) with r taken as
each of the 5 values shown. This gives 5
equations for the 5 desired values of c(r) which in principle could be solved
by minimizing (1.3)

where the constant vector **c **in c can be taken simply as the values of c
at the 5 points where it is assumed to be non-zero. The integral uses h values between R/2 and R
which must be calculated from

(1.4)

for each assumed value of **c**. Note
that we have introduced a “fictitious” short ranged function c(**r**) which
produces the long range function h(**r**) and in addition given a scheme for
calculating this function. In essence the Ornstein Zernicke equation is
fundamental in liquid theory and many people without adequate thought about
Fourier transforms have tried to find efficient ways to do the 3 dimensional
integral for all 3d values of **r** amounting to a 6 dimensional set of
points to calculate.

(2.1)

Begin with the one-dimensional definition of the Fourier transforms.

(2.2)

Extended to 3 dimensions these become

(2.3)

The expression for a delta function is easily recovered by substituting the equation for h(r) into the equation for H(f)

(2.4)

from which we see that

(2.5)

Convolution

(2.6)

on allowing for the evaluation of the r’ integral first

(2.7)

The last integral is zero unless **f’=f **which reduces the set of integrals
to

(2.8)

which implies that

(2.9)

So that while the
original evaluation of the 3d integral for d(**r**) in terms of c(**r’**)
and h(**r-r’**), since this is needed for all **r **it involved 6d of
work or N^{6} function evaluations, the evaluation is now reduced
to 2 3d transforms to find H and C each
of which require N^{3} log_{2}(N) operations, and a back
transform to find d(**r)** involving N^{3} log_{2}(N)
operations. To put numbers to it. A
smart calculation involving only 100 pts per dimension requires 10^{12}
operations to do the straightforward integrals required to find d(**r). **To go the transform route involves 10^{6}x20x3
or < 10^{8 }operations (quite feasible on today’s computers).

** However, the date on the OZ equation is
1914. ** Obviously, the use of symmetries can reduce
this still further. For those cases
where h(**r**)=h(**|r|**) the transform H is given by

(2.10)

where we notice that
H(**f**) is also a function of |**f**| only. Explicitly making the back transform, we
notice that the - sign in the exponent cancels out of the integration and that
the back transform is the same as the forward transform.

(2.11)

So now, the only problem is how to make a sine transform. Note that this is similar to the symmetric region problem with the exception that the function H(f) now needs to be defined as f x H(f) and the function h(r) needs to be defined as r h(r). Sinetran.htm. This transform occurs so often, that Press gives us his version.

Of course it uses the FFT as its base,
but the integral is from 0 to infinity, not -+ infinity and involves the sine
not the exponential, so a factor of two savings in the number of data points
needed is in principle possible **Press (section 12.3 both editions)** gives
us the necessary code.

**Calculate by convolution methods for the following test
problem.**

Consider a charged sphere of radius 10 and charge density 1. Gauss’s law tells us that

4πr^{2} E(r)=-(4/3)πr^{3} for r<R

and 4πr^{2}
E(r)=-(4/3)πR^{3 } for r>R

so that

v(r)=1000/(3r^{3})
r>R

v(r)=50-r^{2}/6 r<R

We also have

which is a convolution integral. for\CONVOL.for

Interestingly the transform of 1/r is not trivial

Noting that the z axis of the r coordinate system can be located along f.

The integral over mu is easy

This nice integral does not converge. Introduce an epsilon as a convergence factor.

So that

For more detail than desired, see test_problem