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 = r1 r2
(0.5)
Define the single variable g to be the average of (0.3) over r2 for r1 = r + r2
(0.6)
Most inter-particle potentials depend only on the magnitude of r, thus define
(0.7)
This is crudely related to the scattering function in ..\definitions\gofr\IbachLuth.htm. A prescription for finding g(r) using
Finally
(0.8)
The Ornstein Zernike equation
(Proc. Acad. Sci. Amsterdam, 17, 793 (1914))
Define the function h(r) to be
(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)=hMC(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 N6 function evaluations, the evaluation is now reduced to 2 3d transforms to find H and C each of which require N3 log2(N) operations, and a back transform to find d(r) involving N3 log2(N) operations. To put numbers to it. A smart calculation involving only 100 pts per dimension requires 1012 operations to do the straightforward integrals required to find d(r). To go the transform route involves 106x20x3 or < 108 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πr2 E(r)=-(4/3)πr3 for r<R
and 4πr2
E(r)=-(4/3)πR3 for r>R
so that
v(r)=1000/(3r3) r>R
v(r)=50-r2/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