OZeqn2.doc1/20/20191/1

The Ornstein Zernike equation

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

Define the function h(r) to be

(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 defined in gofr.doc and gofr2.doc. This probability in a liquid becomes one for large enough r. The Ornstein Zernike equation defines the direct correlation function c(r)

(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. The value of (2) comes from the fact that it relates a short range function in r to a long range function in r. Names such as Percus Yevick and Hyper-Netted-Chain or BBKGY are associated with approximations to c(r). A technique due in part at least to Hanson and Levesque is to define c(r) by

h(r)=hMC(r) r<R/2

c(r)=0 r>R/2

The integral in (2) is a 3-d convoution integral. Direct evaluation of this is detailed in Convolution in 3d.doc. Usually c and h are functions of |r| which reduces the integrals to one dimensional integrals as detailed in Sine transform.doc. In transform space ConvDetail.doc

(3)

or

(4)

or

(5)

Monte Carlo

Assume that as a result of the Monte Carlo calculation described in gofr2.doc, the values of g(r) or equivalently hMC(r) are known for vales of ri i = 1 to M. Define

(6)

where c is the set of c values at each point ri. Since there are M constants and M points in the sum, it should be possible to make 2 become zero. This will not always result in a reasonable function c(r). Three possible approaches are detailed below.

  1. A simple way to minimize 2 is to make it the ftbm of the AMOEBA.
  2. Extremal.htm discusses non-linear minimization of 2. In general this useses nlfit\Welcome.htm. It requires inputting the M values of hMC(ri) as data and the writing of a routine Poly which calculates h(ri,c). To efficiently calculate these values using Sine transform.doc, they need to all be calculated on a first call to Poly, stored and then used later on subsequent calls. The routine Poly also needs the partial of h with respect to c(ri). This can be done numerically, but the code is greatly simplified if these derivatives can also be calculated on a first call and stored (Derivatives of h.doc).
  3. If the sum in 2 is extended to all N values of r between r1 = 0 and r1 = N where h is assumed small, and a wi is defined such that wi = 1 for i < M and wi = 0 for i > M, then (6) becomes similar to that in WeightedLeast SquaresbyIterationLambda.htm.doc. The difference is that in the fit discussed there, H is fitted rather than C. A version of this designed for this fit is WeightedLeast SquaresbyIterationLambda.doc. Each step in this method produces a new estimate for c(r). These can be combined using the technique discussed in Mixing.htm.doc

A single constant c.

It is not necessary to set up a Monte Carlo integration to deduce some things. In particular one can assume that for a system of particles interacting by means of a Lennard Jones Potential that there is no chance of finding two particles close (at r = 0) and that this determines most of the two body correlation function. Assume that c(r) =  for r < r0 and is zero elsewhere. Then determine the single constant  by requiring h(0) = -1.

Figure 1 Expected values of h(r) and c(r) in the single constant (alpha) model.

The integral of g(r) over a box of radius R, where 4R3/3 =  is equal to the volume.

(7)

So that

(8)

  1. Set c(r) =  for r<r0
  2. Sine transform.doc derives the fact that (9)
    . Change to standard form by letting x=2fr . Then the integral is
  3. In the limit that f0
    This is the volume of the short range function c(r).
  4. Sine transform.doc derives the fact that .
    In the limit that r0 this becomes
    This is not a convergent integral, H(f) for large f goes as cos(2fr0)/(f2) so that this integral becomes
  5. . Note that the non-convergent part is small and oscillatory. There is some advantage to making Fr0 such that the sin(2Fr0)=0, but in general, this is simply the abiguity of a sharp cutoff in r.
  6. Define Then select a reasonable F and r0 and solve this for  using Bracketing.htm.doc.
  7. Finally now that  has been found and F and ro assumed. The two body correlation function h(r) is found by making an FFT of H using code contained in Sine transform.doc

Simulated Monte-Carlo 2

At low densities, we expect g(r) to be approximately 1 until the potential becomes repulsive. Then we expect g(r) to be exp(-VLJ(r)/(kT)). Thus we have as our first guess

Transform this using Sine transform.doc

Form

Set I =1


  1. If 2 exit
  2. Set I=I+1
  3. Form
  4. Goto 1

The above procedure may not converge. There are, however, M values of h for r <  and M values of c, so in principle it is possible to form these.

Monte-Carlo

It is only slightly harder than the above to use the Markov chain with a box of side L, to find a realistic estimate of g(r) for r<r0 gofr2.doc ..\integration\MonteCarlo\gofr2.htm

R

+ +

+ + R + +

using Monte Carlo techniques. Obviously there is an enormous premium for making small N give the large N result. The close in parts for r<R/2 rather rapidly reach the same result that one would find in an infinite system, while the parts for r>R are obviously very biased by the fact that they are averages over periodic repeats.

This is only one of many possible ways to define c(r), others have 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 (10)
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
(11)
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.