New methods of computing radial distribution functions of fluids

Material Information

New methods of computing radial distribution functions of fluids
Added title page title:
Radial distribution functions of fluids
Lado, Fred F., 1938-
Publication Date:
Copyright Date:
Physical Description:
xi, 121 leaves. : ill. ; 28 cm.


Subjects / Keywords:
Approximation ( jstor )
Coulomb potential ( jstor )
Diameters ( jstor )
Differential equations ( jstor )
Distribution functions ( jstor )
Equations of state ( jstor )
Fourier transformations ( jstor )
Radial distribution function ( jstor )
Sine function ( jstor )
Thermodynamics ( jstor )
Dissertations, Academic -- Physics -- UF
Kinetic theory of liquids ( lcsh )
Physics thesis Ph. D
Statistical mechanics ( lcsh )
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis -- University of Florida.
Bibliography: leaves 117-120.
Additional Physical Form:
Also available on World Wide Web
General Note:
Manuscript copy.
General Note:

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
000543389 ( AlephBibNum )
13143435 ( OCLC )
ACW7097 ( NOTIS )


This item has the following downloads:

Full Text





August, 1964

A mis padres.


This dissertation is concerned with two problems,

relatively independent of one another. The first is the

question of the effect of long range forces on the radial

distribution function g and the second deals with

attempts to improve the integral equations from which g

may be computed for any potential.

Chapter I presents a brief review of the major

results obtainable from the cluster diagram expansion of g

pioneered by Mayer and Montroll. It provides the background

for both Parts One and Two. The special case of long range

forces is taken up in the second chapter, where two expres-

sions are obtained, based on different approximations,

which relate a given "short range" g and a long range

potential to the total g arising from the sum of the long

and short range potentials. These expressions are then

solved numerically for the Gaussian model and comparisons

made to previous solutions. The solutions are reported in

Chapter III. Chapters IV and V are taken up by the generali-

zation to two-component systems of the expressions derived

in Chapter II and their numerical solution for a model of

an electrolyte.

In principle, the radial distribution function for

any fluid may be found "exactly" by the Monte Carlo (MC)

technique, given the pair potential of the fluid. In

practice, however, this technique is a voracious consumer

of computer time, so that it is usually used as a means of

producing "touchstones" against which other, more frugal

approximate solutions may be tested. These approximation

methods usually result in integral equations for g which

have the added advantages over MC that they provide a means

of computing experimental potentials from scattering data

and may indeed have analytic solutions, as witnessed by

Wertheim's solution of the Percus-Yevick (PY) equation for

hard spheres.

The diagram development of Chapter I suggests a

possible avenue of approach in the improvement of the

integral equations for g derived there. This matter is

taken up again in Chapter VI, which begins Part Two. The

upshot of the arguments presented here is a new integral

equation for g which contains a parameter, m, determined

by the peculiar potential and temperature of the system

studied. This equation reduces to the well-known Percus-

Tevick (PY) or convolution-hypernetted chain (CHNC) integral

equations upon the proper choice of the value of m Solu-

tions of the equation have been obtained numerically for

the Lennard-Jones, Coulomb, and hard sphere potentials and

the results compared to the PY and CHNC solutions, as well

as to the "standard" solutions available. These numerical

results constitute Chapter VII.

The numbering of equations throughout the disserta-

tion is by chapter and serial order therein. Thus Eq.

(5.11) refers to the eleventh equation in Chapter V. The

same rule applies to tables and figures.


This is the fourth dissertation produced by the

statistical mechanics group under the direction of Professor

Arthur A. Broyles. As such, its author has enjoyed the

benefits of the accumulated experience, both published and

private, of his predecessors, Drs. H. L. Sahlin, A. A. Khan,

and D. D. Carley. Wherever possible, specific acknowledg-

ments have been made in the text, but in many cases the

debts are of a more subtle nature.

It is not possible, by a simple statement of

gratitude, to repay the large debt that has accrued to

Professor Broyles in the course of his guidance of the

author's research. Beyond the professional aspect,

Professor Broyles shares with the author's wife his heart-

felt gratitude for their always cheerful counsel and

encouragement, especially during the lower points of the

cyclic process known as doctoral research.

The author was able to enjoy the computing stages

of this research thanks to the presence of Mr. R. A. Smith,

who stood always ready to pinpoint those more insidious

programming errors which tend to mar the relationship be-

tween man and computing machine. All calculations were

carried out on the IBM 709 computer of the University of

Florida Computing Center.

The physical appearance of this dissertation is due

to the efforts of Mrs. Thyra Johnston, who typed the

manuscript, and Miss D. Smoleny, who drew the figures.

This is, finally, a welcome opportunity to thank the

members of the supervisory committee, who, especially in

their role as teachers, have given much of their time and

efforts towards the author's benefit.




PREFACE . . . . . . . . . . iii

ACKNOWLEDGMENTS . . . . . . . . ... vi

LIST OF TABLES. . . . . . . . . . x

LIST OF FIGURES . . . . . . . . ... xi




FUNCTION . . . . . . . . . 2

II. LONG RANGE FORCES. . . . . . . ... 15



SOLUTIONS. . . . . . . . ... 55




VII. NUMERICAL SOLUTIONS. . . . . . ... 72

The Lennard-Jones Potential. . . . ... 77

The Coulomb Potential. . . . . ... 88

The Hard Sphere Potential. . . . ... 89



A. FOURIER TRANSFORMS . . . . . .. 96


STATE. . . . . . . . . ... .105


LIST OF REFERENCES . . . . . . . . 117

BIOGRAPHICAL SKETCH. . . . . . . . ... .121



6.1 Representative Values of the Parameter m .

6.2 Effect of Temperature on m for the LJ
Potential . . . . . . . .

7.1 Comparison of PT, CHNC, and "m" Integral
Equation Solutions to Monte Carlo for the
LJ Potential: T* = 2.74, n* 2/5 . .

7.2 Comparison of PY, CHNC, and "m" Integral
Equation Solutions to Monte Carlo for the
LJ Potential: T* = 2.74, n* = 5/6 . .

7.3 Comparison of PY, CHNC, and "m" Integral
Equation Solutions to Monte Carlo for the
LJ Potential: T* 2.74, n* = 1.0 .

7.4 Thermodynamic Quantities for the LJ
Potential . . . . . . . .

D.1 The Computed Functions g
Potential . . . .

D.2 The Computed Functions g
Coulomb Potential . .

D.3 The Computed Functions g
Hard Sphere Potential.

D.4 The Computed Functions g
Hard Sphere Potential.

D.5 The Computed Functions g
Hard Sphere Potential.

D.6 The Computed Functions g
Hard Sphere Potential.

D.7 The Computed Functions g
Hard Sphere Potential.

and H for the

and H for the

and H for the

and H for the

and H for the

and H for the

and H for the










. 113

. 114


. 116


Figure Page

3.1 Radial distribution functions for the
Gaussian model: n 0.4, m = 0.9 . . .. 35

3.2 Radial distribution functions for the
Gaussian model: n = 0.4, m 0.8 . . .. 36

3.3 Radial distribution functions for the
Gaussian model: n = 0.4, m = 0.7 . . .. 37

3.4 Radial distribution functions for the
Gaussian model: n 0.4, m 0.6 . . .. 38

3.5 Gaussian model potentials . . . . .. 39

3.6 RMS deviations. . . . . . . . ... 40

3.7 RMS deviations. . . . . . . . ... 41

3.8 Density dependence of the RMS deviations. . 42

5.1 Radial distribution functions of a 0.00911
molar electrolyte . . . . . . ... 60

6.1 Mayer f bonds for various potentials. ... .70

7.1 Radial distribution function for the LJ
potential: T' = 2.74, n* = 2/5 . . ... 80

7.2 Radial distribution function for the LJ
potential: T* 2.74, n* 5/6 . . ... 82

7.3 Radial distribution function for the LJ
potential: T* = 2.74, n* 1.0 . . . .. .84

7.4 Radial distribution function for the Coulomb
potential: 9 w 0.4 . . . . . ... 90

7.5 Radial distribution function for the Coulomb
potential: 9 = 1.0 . . . . . . .. 91

7.6 The equation of state of a fluid of hard
spheres . . . . . . . .... ... 94






Statistical mechanics seeks to deduce the macroscopic

thermodynamic behavior of a physical system from the laws

governing the microscopic behavior of its constituent

particles. This goal can be embodied in the calculation of

a single quantity from which all thermodynamic information

can be obtained. Thus, for a classical canonical ensemble,

this end is achieved by a knowledge of the partition

function Q defined by

Q- 1 -...f&H N d3p, N d
Q e ~d pi N d3ri (1.1)
N J J i=l i=l


S 1/kT (1.2)

Here T is the absolute temperature and k and h are

Boltzmann's and Planck's constants, respectively. The

N-particle Hamiltonian H describing the system is a function

of the 3N components of momentum and the 3N components of


The free energy F of the system is then given by

-BF = In Q (1.3)

and all other thermodynamic quantities can be obtained from
the free energy. For example, the energy E of the system
is given by

E b(BP)/b (1.4)

Another quantity that may be used to compute the
thermodynamic quantities is the pair distribution function

g(rl,r2) (defined below) if the particles are assumed to
interact only through two-body forces. The integration
over the moment in Eq. (1.1) may be carried out immediately,
Q Z/NIA35 (1.5)

Z = .. ee d3 r ... d3rN (1.6)

A h(2lmkT)-1/2 (1.7)

Here U(51,...,-N) is the potential energy of the system
and V the volume. We have assumed that all particles have
the same mass m so that the kinetic energy is
T = pi2/2m Then we may write

pF - In(NIA3N) + In Z (1.8)

and, from Eq. (1.4),

E NkT + .. e- U dr ... d3r (1.9)

Implementing the assumption of pairwise interaction,
we write
U 0 (r ) (1.10)
rij |i j (1.11)

and 0 (r) is the pair potential. With this, Eq. (1.9)
1 N2 3
E NkT + 0 (rl2)g( 1 2)d rld r2 (1.12)

where we have defined the pair distribution function

g(71, 2) by [1]

)1 v2 2 -U
g(1,r) = (1 ) ... e dr ... drN
The 3-particle, 4-particle, etc., distribution functions
can be similarly defined. In general, the n-particle
distribution function is given by

g.(n) n f... d rn+ ... d3rN

where terms of order 1/N have been ignored.

For a fluid, which shall be our exclusive concern

here, g is a function only of the magnitude of the
separation between particles 1 and 2; i.e.,

g(rlr2) g(rl2) (1.15)

This has given rise to the use of the name "radial distri-

bution function." The normalization of the pair, or radial,
distribution function is

ifg(r)dr 1 (1.16)

Physically, the distribution function g(r) is

proportional to the probability of finding any given pair

of particles in the fluid separated by a distance r ,

irrespective of the positions of the remaining N-2 particles.
It has the additional property that its Fourier transform
may be directly measured by x-ray and neutron scattering
experiments [2,3].
Ignoring terms of order 1/N g(r) is defined as

2 ..... drr r
g(r) = ... e- d r ... drN (1.17)

All thermodynamic quantities can be expressed in terms of

g(r) .
The computation of g(r) can be made more tractable
by the use of Mayer f functions, whereby the right hand side

of Eq. (1.17) is expanded in Mayer cluster integrals. This
results in the well-known expansion of g(r) in powers of
the density [4-9]:

G(r) g(r) 1 (1.18)
f(r) + l[ + f(r)]C(r) (1.19)

Here. C(r) is defined by

C(r) nt f dr ... drt+2
'77707 fij d3 "' d rt+2
t-1 J

fij f(rij) = e 1 (1.21)

and a 1,2 index has been suppressed on r Each power of
the number density n = N/V is associated with several
integrals, denoted by Z The integrals may be represented
by diagrams according to the following prescription. For
each fij in the integrand draw the points i and j and
connect them with a line. Points 1 and 2, which are not
variables of integration, are called the fixed points and
are drawn as white (i.e., unfilled) circles. The t field
points are drawn as black circles. The symmetry number a
depends on the number of field points (t) and the diagram
type (a). It is equal to the number of permutations of the
field points that lead to the same diagram, that is, to

integrands having identical products of fij Equation
(1.20) may thus be read as the sum over all graphs having
two distinguishable points fixed, the contribution of each
graph being divided by its symmetry number.
The first few terms in the expansion of C(r) in
this shorthand notation are

C(r) n A + n2+

+ + * (1.22)

where the density and symmetry factors have been explicitly
written. Any diagram can easily be written out explicitly
in terms of the integration it represents by applying the
prescription given above in reverse order. Thus we have

S ff3ff34f4f42d3r3d53 (1.23)

These diagrams can be classed into three groups
according to their topological characteristics:
(a) A series, or nodal, diagram contains at least one
field point (called a node) through which all paths
connecting 1 and 2 must pass.
(b) A parallel diagram contains at least two paths
between 1 and 2 which are connected only at 1 and

(c) A bridge diagram is one which is neither series
nor parallel.
The sum of all series, parallel, and bridge diagrams will
be denoted by S(r), P(r), and B(r), respectively. To the
second power in density these sets are

S(r) noo+ .n2 + b + + (1.24)

P(r) 7 n2 Q + .. (1.25)

B(r) !. n2 Q P+ * *^. (1.26)

In terms of these sets, we have

C(r) S(r) + P(r) + B(r) (1.27)
G(r) S(r) + P(r) + B(r)
+ f(r) [1 + S(r) + P(r) + B(r)] (1.28)
g(r)eB0(r) 1 + S(r) + P(r) + B(r) (1.29)

We will denote the set of all non-nodal diagrams in the
expansion of G(r), Eq. (1.28), by T(r); i.e.,

T(r) P(r) + B(r) + f(r)[l + S(r) + P(r) + B(r)]
g(r)f(r)e00(r) + P(r) + B(r) (1.31)

and thus

G(r) S(r) + T(r) .(1.32)

It is important to note that the sets S, P, and B

do not include diagrams with a direct 1-2 bond. Diagrams
of this type are put in explicitly by the definitions, as
in the last term of Eq. (1.28). Thus also the non-nodal set
T defined by Eq. (1.30) includes not only the sets P and

B but also all diagrams with a direct 1-2 bond.

Consider a typical diagram of the set S By
definition it has at least one field point which is a nodal

point. Let the first nodal point encountered along a path

from 1 to 2 be labelled 3. Then the subdiagram between 1

and 3 will have no nodes, by construction, and hence will

be some member of the set T The subdiagram between 3

and 2 may or may not have further nodes. There is no

restriction. Hence in general it will be some member of

the set G Thus we can generate all possible members of
the set S by substituting between 1 and 3 all possible

members of T and between 3 and 2 all members of G That

S(r) n / T(rl)G(r2)d r3 (1.33)

The factor n is introduced because the point 3 is a field
point for the whole diagram, although it plays the role of

a fixed point for the subdiagrams. The same applies for the

integral over .

With Eq. (1.32) we have thus an integral equation for

G(r) ,

G(r) T(r) +n T(r')G(j 5' )dr' (1.34)

This is the integral equation of Ornstein and Zernicke [10]
which gives G(r) when the direct correlation function

T(r) is known. Unfortunately, we can see from Eq. (1.31)
that T involves the unknowns P and B We will see
how to eliminate the set P in the following, but the
bridge set will remain to prevent the formation of an exact
integral equation with only one unknown.
The Fourier transform of Eq. (1.34) may be taken,
using the convolution theorem, with the result that

G(k) = T(k)/[1 nT(k)] (1.35)

where the tilde denotes the Fourier transform. (See
Appendix A.) Similarly, from Eqs. (1.52) and (1.35) we have
S(k) = nT2(k)/Cl nT(k)] (1.36)

Consider now a typical member of the set P. Since
there is no integration over particles 1 and 2 the diagram
may be "factored" by cutting at the two fixed points and
separating out the various independent branches. For
example, we may write

6 [ I [ 0 2 (1.37)

The branch diagrams may be members of the sets S and B ,
but not of P
We are interested in the contribution of all parallel
diagrams with t branches, which we will denote by pt(r)
Consider a typical diagram with t branches and let the
contribution from the ib branch be ai(r) so that the
contribution of the whole diagram is

Cal(rr)]a2(r)] .. E[at(r)]

We can generate a subset of pt(r) by holding
a2(r),...,at(r) constant and allowing al(r) to be any of
the diagrams of P(r) and B(r) The contribution of
this subset is then
[P(r) + B(r)][a2(r)] ... [at(r)] .

A still larger subset of pt(r) which includes the one
above as a subset of itself, may be found by allowing the
second branch to be composed of any member of S(r) or
B(r) as well as the first, giving the contribution

i ES(r) + B(r)][S(r) + B(r)][ 3(r)] .. [at(r)] .

Diagrams with non-equal subdiagrams are counted twice in
the multiplication of the first two bracketed terms and the

factor 1/2 removes this redundancy. Diagrams with equal
subdiagrams are only counted once and here the factor 1/2
serves to account for the correct symmetry factor. We may
go on in this way, allowing the first 3 branches to take
on arbitrary members, and so on, obtaining for the total
contribution of diagrams having t branches


pt(r) S(r B(r)
pt(r) t!L~1

For the parallel set P(r) t
from 2 to infinity, so that

may be any number

P(r) = Sr) B(r

SeS(r) + B(r) S(r) B(r) 1

Thus, from Eq. (1.29),

g(r)e(r) = eS(r) + B(r)

in [g(r)e O(r)] =

S(r) + B(r)

Then substituting Eqs.
a useful relation for

(1.41) and (1.42) into
P viz.,

(1.40) gives

P(r) g(r)e(r) - n[g(r)eP0(r)


which can be used to eliminate P(r) from previous equations.





Thus the direct correlation function T becomes, from
Eq. (1.31),

T(r) G(r) ln[g(r)eO(r)] + B(r) (1.44)

Unfortunately the set B has no exploitable properties
which would effect its elimination and it is at this point
that approximations must be introduced. Perhaps the simplest
approximation which can be made is to neglect B altogether,
hoping that it is small. With this approximation,

T(r) G(r) In[g(r)eO0(r)] (1.45)

Inserting this into Eq. (1.34) gives the convolution-
hypernetted chain (CHNC) integral equation [11-16]

ln[g(r)e00(r) nf[G(r') In g(r')e 0(r')

x G( F'I )d3r' (1.46)

Another approximation to B which has proved useful
is the assumption that

B(r) P(r) (1.47)
B(r) + P(r) 0 (1.48)

Neglect of these two sets gives

T(r) g(r)f(r)e00(r) (1.49)

as can be seen easily from Eq. (1.31). Combined with

Eq. (1.54), this approximation leads to the Percus-Yevick
(PY) integral equation [17-20]

g(r)e(r) = 1 + n fg(r')f(r')e 0(r')

x G(IJ r'l)d3r (1.50)

A third and older integral equation for g which is
often found in the literature is the equation of Born-Green-
Yvon (BGY) [21,22], based on the superposition approximation

[1]. It has been found to be less reliable than the PY or
CHNC integral equations for a Lennard-Jones (6-12) potential
by Broyles et al. [23] and for hard spheres by Klein [24,25].
We will not have occasion to refer to it again in this


For a realistic potential made up of a repulsive core

and an attractive part such as the Lennard-Jones (6-12)
potential, Broyles et al. [23] have found that at relatively
high temperatures the PY equation is superior to the CHNC,
while the results obtained by Khan [20] suggest that at low
temperatures the CHNC equation may be the better of the two.
We will return to this temperature effect in Part
Two, where an attempt is made to obtain an integral equation
capable of adjusting itself to the particular potential and
temperature desired.



We will be concerned in this chapter with the effect

on g(r) and hence on the thermodynamic quantities, of a

small change in the pair potential, 0(r) describing the

system. A solution to this problem could be used in a

variety of applications. On a practical level, the need

for a method to correct g arises, for example, in Monte

Carlo calculations of the radial distribution function,

where the long range tail of a potential such as the

Coulomb potential must necessarily be truncated at some

finite distance. The effect of the neglected part of the

potential must be found for a complete solution [26,27].

Furthermore, if the function g is known for some

temperature T g for some slightly different temperature

T' may be easily found by considering 0'0 0' 1/kT',

to be a perturbation of 30 at T and applying the corres-

ponding correction. This obviates the need of another long,

iterative solution of the integral equations for g .

From a theoretical point of view, the solution of

this problem would allow the effects of the long and short

range parts of the potential to be studied independently.

Thus the change in the potential could be the addition of a
weak long range interaction to a short range repulsive core.
The derivation of a van der Waals-like equation of state is
based upon this type of potential [28,29]. A one-dimensional
model of this type has been studied by Kac, Uhlenbeck, and
Hemmer in a succession of papers [30-32] and was found to
display a phase transition even in lowest order.
The classic example of a long range potential is
the Coulomb potential which is, in principle, of infinite
range. An explicit approximation to the radial distribution
function for a system of particles interacting with Coulomb
forces was found by Debye and Huckel [33] and reads

g(r) exp [- 2 e-r/ (2.1)

S- [4 x npe2]-1/2 (2.2)

and e is the electron charge.
This result can easily be obtained from the develop-
ment of the last chapter. From Eq. (1.41) we get

g(r)e00(r) eS(r) (2.3)

by neglecting B(r) as in the CHNC approximation. The
corresponding direct correlation function T(r) is then
given by Eq. (1.45), which can be written

T(r) G(r) In[l + G(r)] B0(r) (2.4)

For a weak long range force, G(r) will be small over a
large region so that the logarithm in Eq. (2.4) may be
expanded. Thus,

T(r) - 00(r) + 1 G2(r) + (2.5)

If only the first term is retained we obtain a non-iterative
solution for g by combining Eq. (2.3) with

S(k) n[E0(k)]2/[l + nB0(k)3 (2.6)

obtained from Eq. (1.36) and the above approximation to T

S(k) 8p(k) 00(k)/[l + n00(k)] (2.7)

g(r) exp (k k sin krdk (2.8)
2r-r 1 + nB0(k)

where the angular integration in the Fourier transform have
been carried out. For the Coulomb potential,

P0(k) 4 ne2/2 (2.9)

(See Appendix A). By inserting this expression into Eq.
(2.8) and performing the integration C34], we obtain the
Debye-Htckel (DH)g Eq. (2.1).

The two essential approximations of the DH
expression are thus

B(r) 0 (2.10)

T(r) -00(r) (2.11)

Carley [26] has recently shown that the DH g is at least
as good as, if not better than, g's obtained from the PY
and CHNC integral equations, both of which evaluate many
more diagrams in the expansion of g .
Let us now consider a potential 0(r) which can be

written as the sum of a short range and a long range part;


0(r) 0Sr(r) + 0lr(r) (2.12)

We will assume that g corresponding to the short range

potential is a known function, which we call gr From
Eq. (1.41) we may write

gsr(r) expCSsr(r) + Bsr(r)] (2.13)

where the label "sr" on S and B indicates they are the

sets corresponding to sr Similarly, for the complete
potential 0 we may write another equation identical to
Eq. (1.41). Then dividing this equation by Eq. (2.13) we
have formally

g(r) gr(r)expC-00r(r) + AS(r) + AB(r)] (2.14)


AS(r) S(r) SSr(r) (2.15)

and AB(r) is similarly defined. From Eqs. (1.35) and
(1.36) we have further that

~2 -sr2
AS(k) nT nT (k) (2.16)
1 nT(k) 1 nTsr(k)

l + nG-Sr(k)2AT(k) AT(k). (2.17)
1 n[l+nGsr(k)]AT(k)

All we need now to complete the solution is a set of
suitable approximations for AT and AB Both of these
sets are due to the long range potential since, if it were
absent, they would be identically zero. Hence by analogy
with the long range problem solved by the DH g we neglect
AB and approximate AT by -00lr This leads then to the
final expression

g(r) gr(r)e-H(r) (2.18)


H(k) + nr(k)2 r(k (2.19)
1 + nCl + nG-r(k)]e0 r(k)

Here g(r) is given completely in terms of the long range
potential and the short range g which is the desired form

It is easy to obtain higher approximations to AT
Having written the potential as a sum of two parts, we may
write for the corresponding Mayer f function

f(r) fer(r) + e-BSr(r)flr(r)

= fSr(r) + Af(r) .



If we replace each f bond in a typical diagram of the
expansion of T(r) for example,

S = f12f13f32d3r


by the sum of fsr and A and expand the products, we
obtain the following,

2: 6-- :- %- + c=% + 6 +

+C. +~~~ +

Here a wavy line connecting
fsr bond and a dashed line

two points

0^-;0 + -'-o *

indicates an
Af e srflr
af W e- f

Since we have assumed that 01r is a weak potential we need
consider only terms up to first order in Af; that is,
diagrams having only one dashed line. Performing the above
expansion in each diagram of T gives rise to one set of
diagrams having nothing but short range bonds, another
having one AF bond in each diagram, and so on. The sum of
the first set gives immediately Tar(r) We write

T(r) Tar(r) + 0----o + n --o + o' %

+ 6 o + ** (2.24)

The sum of all diagrams having one AF bond is the
desired AT The previous approximation to AT may be
recovered by taking only the single bond diagram after Tsr
as AT writing it in two parts,

0--- e- r(r)flr(r) (2.25)

Sf1r(r) + fsr(r)flr(r) (2.26)

and further approximating by neglecting the second part and

flr(r) ~0lr(r) .


That is,

AT(r) 0lr(r) (2.28)

We can, however, retain the whole first diagram, giving

AT(r) e-Bsrflr(r) (2.29)

Further, all diagrams with a direct Af bond between 1 and
2 may be summed to give

AT(r) Af(r)gsr(r)e Bsr(r) (2.30)

gSr(r)flr(r) (2.31)

It will turn out, however, that these higher
approximations, Eqs. (2.29) and (2.31), are not as satis-
factory as the original approximation, Eq. (2.28). This
will be discussed further in the next chapter.
Consider now the logarithm of both sides of Eq,

lng(r) Ingsr(r) H(r) (2.32)

In the region where g(r) is close to one, we may approxi-

lng(r) g(r) 1 (2.33)

lng r(r) gsr(r) 1 .


Then Eq. (2.32) becomes

G(r) Gsr(r) H(r) (2.55)

1 + nG(k) 1 + n--r(k) (2.36)
1 + n[1 + nGsr(k)] r(k)

(k) 1 1 + nGr(k) 1
G(k) n + n[ + nGsr(k)101r(k) l

This is precisely the equation obtained by Broyles, Sahlin,
and Carley (BSC) by a method based on collective coordinates

[36,37]. The BSC equation, Eq. (2.37), was found to give,
in some cases, anomalous negative g's for small values of
r The derivation given above suggests that it is not
applicable at small r since the substitution of g(r) 1
for lng(r) is not valid in this region.
If we now further restrict the region of r under
consideration to that where gar is practically one, then
Gsr in Eq. (2.35) may be neglected, giving

G(r) - H(r) (2.38)

This equation was obtained by Hemmer as the first order
expression for g(r) in the large r region [38].
We note also that in the limit Osr(r) -- 0 so that

0ir(r) 0(r) Eqs. (2.18) and (2.19) give the DH g
when 0 is specified as the Coulomb potential.
An equation for g(r) employing a PY-type approxi-
mation, where AP(r) + AB(r) is neglected, is easily
obtained. From Eq. (1.29), we write

gsr(r)e08(r)e ) 1 + Sr(r) + Psr(r) + BSr(r)
A similar equation holds for the total g Subtracting
one from the other then gives

g(r)e0(r) s gS(r)e sr(r) + AS(r) + AP(r)
+ AB(r) (2.40)

g(r) gsr(r)e-~Ir(r) + e-0(r)AS(r) (2.41)

where we have used the approximation

AB(r) -AP(r) (2.42)

The set AS is determined in terms of AT by
Eq. (2.17), as before. If we now again approximate AT
by 0lr we have finally

g(r) gsr(r)e-0lr(r) + e-BO(r)[Olr(r)-H(r)] ,(2.43)

where H(r) is the same as in Eq. (2.19). Once again
higher approximations to AT may be incorporated into the

final result, but it is found that the results are less


We have found two expressions which give g(r) for

the total potential

0(r) Osr(r) + 01r(r) (2.44)

when gsr is known. Using a CHNC-type approximation,

AB = 0 we were led to Eq. (2.18). A PY-type approximation,

AB + AP 0 gave rise to Eq. (2.43). In both of these H
is determined by Eq. (2.19).

In the next chapter these equations are tested on a

model for which near exact answers are available for




The equations derived in the previous chapter may

be tested numerically if an exact solution of the system

being studied is available for comparison. We need also

an exact gSr for input, so that any deviations of the

computed g's from the exact g will reflect only the

approximations introduced in the previous chapter. Both

of these conditions can be met by the so-called Gaussian

model, where the Mayer f function is represented as a nega-

tive Gaussian function. This model was used to carry out

numerical computations of Eqs. (2.18) and (2.45), as well

as of the equations resulting from the higher approximations

to AT The BSC equation was applied to the Gaussian model

by D. D. Carley.* The two asymptotic expressions for g

obtained by Hemmer [38] have also been included in the

results reported.

The Gaussian model has been recently investigated by

Helfand and Kornegay [39]. These authors generated and

The author is indebted to Dr. Carley for permission
to include his results in the final comparisons, as well as
for providing an independent check on much of the Fortran
programs used in this part.

classified the series and bridge diagrams having five or less
field points and evaluated these for the Gaussian model.
This data was then used to obtain the coefficients gt(r)

g(r) -e-e(r) 1 + E t(r) (3.1)

up to t = 5. At small values of density the g's obtained
in this way are essentially exact. The criterion was that

n5 g(r) should have an almost negligible effect on the
total g(r) .
The Mayer f function is given by

f(r) -e (r/a)(3.2)

and the unit of distance is selected so as to make the
second virial coefficient of the pressure identical with
that of a gas of hard spheres of diameter d; i.e.,
00 d 00 2
f(r)r2dr r2 dr e-(r/a) r2 dr (3.3)

so that
d3 a3l/2
T5 -zi-- (3. )


d ( )/3/6 a 1.100 a (3.5)

The reduced units are then given by

x r/d n Nd3/V (3.6)

The potential corresponding to a solution g(x) is
given by
00(x) - In(l e-1 21 x (3.7)

We can obtain a g(x) for a slightly different potential in
the following way. Perform a coordinate transformation
x' x/m where 0 < m 1 Then the function 0sr(x) ,

0sr(x) 00(x/m) (3.8)

S-- n(l e-l121(x/m)2) (39)

will be of the same form as 0(x) but go to zero sooner.
That is, the transformation x' = x/m tends to compress
the potential in towards the origin. We will call Osr(x)
the short range potential. The corresponding gsr(x) is
obtained by replacing x by x/m and n by nm3 in
Eq. (3.1); that is,

gsr(x) e-O(x/m) 1 + Z (nm)g(x/m) (3.10)

The long range potential is given by the equation

BPlr(x) 00(x) p0sr(x) .


The two quantities on the right are computed directly from

Eqs. (3.7) and (3.9). This long range potential and the
short range g computed from Helfand's data by means of

Eq. (3.10) are then used as input quantities for the
equations obtained in the previous chapter, as well as the
BSC and Hemmer equations. The computed g's can then be
compared with the nearly exact g obtained also from the
Helfand and Kornegay data using Eq. (3.1).
The programs were written in the Fortran language.
They are straightforward programs and are not included here.
In order to illustrate the numerical procedure we
will consider a particular equation, say Eq. (2.18), with

H(x) determined by Eq. (2.19). For numerical calculation
a maximum value for x Xmax must be chosen. The range
between 0 and iXax is then divided into N intervals of
equal size A A function f(x) is then considered deter-
mined if it is known at the N + 1 points x = JA ,

j = 0,N The calculations of P0(jA) p0sr(jA), 00lr(jA)
and g(jA) are straightforward and require no comment. To
compute gsr(jA) we first recompute g(jA) for the density
nm3 and then transform to a new coordinate i = j/n to
obtain gsr(iA) g(j/m) Depending on the size of m ,
J will reach its maximum value at some point before i
does, leaving gsr(id) undetermined for the remaining points.
In this region g r(iA) was assigned the value one, since

in every case gsr had already settled to its asymptotic
value of one before reaching the undetermined region. The
short range g was then evaluated at the same points as the
other functions by linear interpolation of the results
computed above.
The Fourier transforms appearing in the equation for
H(k) are defined as follows: For any function of the
magnitude of r F(r) we define F(k) such that

F(r) 1 F(k)eik'r d3k (3.12)

F(k) F(r)e-ikr d3r (3.13)

Carrying out the angular integration,

F(r) l- f kF(k) sin kr dk (5.1i)

F(k)- rF(r) sin kr dr (3.15)
For numerical calculation F(r) and F(k) must be made
discrete functions. Put

r iAr k = jAk (3.16)

Then we have

2 2
p(i) j0 i(j) sin (ijArAk) (3.17)
2it iAr J-0
F(J) -. 4n )2L iF(i) sin (ijrk) (3.18)
A i=O0

It is shown in Appendix A that in order to conserve the
orthogonality of the sine functions in this discrete repre-
sentation, and thus the basic relation of a function to its
Fourier transform, we must put

Ark 2x/(2N+l) (3.19)

where N is, as above, the number of points. Then, putting

A 2/(2N+l) (3.20)

we have finally

F(i) A2 L j(j) sin (ijA) (3.21)
2n i(Ar) j-0

F(j) 4TEi ) r iF(i) sin (ijA) (3.22)

Eq. (3.22) is used to compute the Fourier transforms
of B0(i) and Gsr(i) which then determine H(j) The
discrete version of H(x) is then found with Eq. (3.21),
giving the final answer

g(JA) gsr(jA)e-H(JA) (3.23)

Hemmer's small r solution for g(x) Eq. (35) of
Reference 38, was solved in the following form:

g(x) gsr(x) DCgSr(x) + A(x)] (3.24)

D = r(k) -- k2dk (3.25)
2n 1+ n~B r(k)

ad 1 + 4 n GSr(x) x2dx (3.26)

A(x) e -sr(x t+ ngtr(x) (5.27)
Calculations were made with Eqa. (2.18) and (2.43),
the BSC equation, and Hemmer's two asymptotic equations, for
densities of 0.1, 0.2, 0.3, and 0.4. With these densities,
the truncated series for g(x) yields an essentially exact
solution. For each value of the density, four values of the
parameter m were taken, m 0.9, 0.8, 0.7, and 0.6,
representing successively larger perturbations of the
potential, or more significant values of 0lr(x) .
The results for n = 0.4 and m = 0.9, 0.8, 0.7,
and 0.6 are shown in Figures 3.1-3.4. We see that, as
expected, g computed from the BSC equation compares
favorably with the exact g in the large x limit, but
becomes far too small as it approaches the origin. Hemmer's

two asymptotic equations yield good agreement in the small
and large x limits for which they are applicable but allow

no means of interpolation for intermediate values of x .

Of the two equations given at the end of Chapter II, that

corresponding to the PY-type approximation gives a better
result for g although an additional infinite set of

diagrams has been neglected in its derivation. As with the

PY integral equation itself, this success is explained on

the basis of cancellations between the sets P and B

[20]. The potentials for these cases are shown in Figure

The results of other cases for Eqs. (2.18) and (2.45)

are given in Figures 3.6 and 3.7 in the form of rms deviation

from the exact g(x) where the rms is defined by

rms Cge(jA) gcomp(j)2] 2 (3.28)

Here N is again the number of points taken in the numerical
solution and A the interval. The superscripts on g(x)

refer to the exact and computed g's. We note from Figures

3.6 and 3.7 that the rms values obtained for Eqs. (2.18)
and (2.43) are relatively insensitive to changes in density.
It was mentioned in Chapter II that, although higher

approximations could be found for AT the results were

less satisfactory than those obtained with AT approximated


by -B0lr The basis of this assertion lies in the compari-
son of results for the Gaussian model of the three approxi-

mations for AT Eqs. (2.28), (2.29), and (2.31). The

situation is exemplified by Figure 3.8, where the rms

values for the equations resulting from the PY-type approxi-

mation and the above-mentioned approximations for AT are

shown as functions of the density for m 0.6. It is seen

that at small densities the higher approximations do indeed

lead to successively more accurate g's, but at the price of

far worse answers at larger densities.

0 .4

a 0

M ca 0 e 'o
60 bo wM M

Id I "








O 0






o-* 0


0 0 0

O 0
O o



0 4,

^ 0 I I I

N 00

4 .

I 1 0 I

O *









0 T


r< O

uI '*


R 38

lr; a

x C






* 0



0 K
o b



~ 0sr



.5 m=0.6

1.0 m=0.6


0.5 m-0.8

M.O. 9 9,

0 0.5 1.0 1.5 x 2.0

Pig. 3.5.-Gaussian model potentials. These potentials
were used in obtaining the g's of Figures 3.1-3.4.

0.02 -

0.01 -

0 "

n = 0.1

m 0.6

0.02 -

n = 0.2

0.01 L

m 0.6

Fig. 3.6.-RMS deviations. The rus deviations from
the exact Gaussian model g(x) for g computed from
Eq. (2.18) (dashed line) and Eq. (2.43) (solid line).



.01 -

n = 0.3






0.8 m

Fig. 5.7.-RMS deviations. The rms deviations
from the exact Gaussian model g(x) for g computed
from Eq. (2.18) (dashed line) and Eq. (2.43) (solid



rs- Eq. (2.28)

0.02---- Eq. (2.29)
--- Eq. (2.31) /


0 0.2 n 0.4

Fig. 3.8.-Density dependence of the EMS
deviations from the exact Gaussian model g(x) for
g's computed by neglecting AP + AB and approximating
AT by Eqs. (2.28), (2.29), and (2.31). The value of
m is 0.6.



The equations of Chapter II may be readily
generalized to a two-component fluid, such as an electro-
lytic solution or a fused salt. For such a system, three
radial distribution functions are in general necessary.
If the fluid components are labelled 1 and 2, the required
distributions are gll, g22, and gl2 = g21.
We will assume, as in Chapter I, that the potential
U of the system may be written as a sum of pair inter-
actions, as follows:

U(~9.....) ll(rij) + 022( )

+ Z 12(rj) (4.1)

Here 011 is an interaction between molecules of component
1, 022 between those of component 2, and 012 a "cross"
interaction between members of components 1 and 2. The
numbering of component 1 runs from 1 to N1 and of component

2 from N1 + 1 to N, the total number, where

N N1 + N2 (4.2)

Again we assume that we may write
0 sr +r lr ir

B(r) = s (r) + O~ (r) a, = 1,2 (4.3)

where 0sr is the short range interaction between components
a and 3 and so on. In terms of diagrams, the pair

correlation function G(r) is again described by an equation
like Eq. (1.28), where now three different G's must be
distinguished. We write

G (r) f (r) + El + f ( r)][S B(r) + P (r)

+ B a(r)] a,P 1,2 (4.4)

where Sa is the set of all nodal diagrams which begin
at a fixed point belonging to component a and end on a
fixed point belonging to component B Similar definitions
hold for the other sets. Note that these diagrams are not

the same as the diagrams of Chapter I. While the topological
shapes are the same for both, the addition of a second
component to the fluid allows many distinct diagrams to be
made up from one diagram of the single component type.
In the following a and 0 will always be under-
stood to take on the values 1 and 2.

Direct correlation functions TaO(r) are defined
in complete analogy with the one-component fluid, viz.,

T (r) G,(r) S (r) (4.5)

but now the expressions for the series diagrams become

S B(r) .- Z n Ta,(r')G(15-'1)dr' (4.6)

Here n1 and n2 are the number densities of components
1 and 2, respectively. Equation (4.6) expresses the fact
that though the two fixed points of the set Sa are

necessarily members of the a and 0 components, one from

each, the intervening field points may be of either
component 1 or 2. The field point ?' in Eq. (4.6) is
occupied by a member of component X It is treated as a
fixed point by the diagrams in the sets Ta and G \
and the fact that it is a field point of the whole diagram
accounts for the presence of the n. as well as the
integration over r' .

The differences in the diagrams may be simply
illustrated by considering the coefficient of the first
power of n in, say, gll(r) Here we have

n1 o +

where the circles denote component 1 and the squares com-
ponent 2. In the one component fluid there was only one
diagram for this coefficient. [See Eq. (1.22).]
Equation (4.4) may be rewritten to read

gB(r)e (r)

= 1 + Sap(r) + Pap(r) + BaB(r) (4.7)

The parallel set may be eliminated as before to yield

aB(r) (r) eSap(r) + BaB(r)
ga (r)e


That this is so can be made plausible by considering the
first few terms of a specific g, say gll. Then we have

Sl(r) = nI 01 + n2 oo + n12 2

+ n1n2C + +

+ na2 oI+

n1n2 0Al0 + n22 /"\a

n1n2 o o + n1n2

+ n22 o + n2 c n+ n1n2

+ n1n2 a"/ o+ n22 oAa o+ . .


B,(r) nl2 + nn2 ) + n2 c + ..

where the circles and squares again refer to components 1
and 2, respectively. [Cf. Eqs. (1.24) and (1.26)]. Then
the equation
Pa(r) misa(r) + Ba (r))m (4.11)
[see Eq. (1.39)] becomes, for P1(r) ,

P11(r) 81 S2(r) + Sll(r)B11(r)

+ B11(r) + ... (4.12)

nl2 + n1n2 + n2 + ...

from Eq. (4.9). These terms satisfy the definition of the
parallel diagrams for the two-component system.
For the short range potentials alone, an equation
identical to Eq. (4.8) may be written, with all quantities
labelled with an "sr" superscript. Then the ratio of the
two equations gives, as in Chapter I,

g (r) ga(r) exp [-BF30(r) + AS (r) + ABa(r)]
ao ao ao (X5



ASap(r) Sa(r) Sa(r) (4.15)

and a similar definition holds for B .
Before proceeding to the evaluation of the quanti-
ties ASa and ABa we need to obtain explicitly the
relations between the Fourier transforms of the three Ga
and the three TaB implied by Eqs. (4.5) and (4.6). By
taking the Fourier transforms of these equations, with the
aid of the convolution theorem, we get the following:

G p(k) T (k) + nT(k) (k) (4.16)

These equations may be solved simultaneously for either
G in terms of TB or vice versa. In the first case
they give

El-n2T22(kl)T(k) + n2T12 (k)
G1(k) .- (4.17)
) [l-nlT11(k)]l-n2T22(k) nln2T12 (k)

[l-nlTll(k)]T22(k) + nlT12 (k)
G22(k) 2 (4.18)
[l-n1T11(k)][l-n2T22(k)] nln2T12 (k)

12(k) .T12(k) .(4.19)
1 -n1 11(k)][l-n2T22(k)] nn2T12 (k)

For the inverse solution we obtain

[l+n2G22(k)]Gll(k) n2G12 (k)
Tll( [l+n1 Gk+nG11(k) +n222(k)] nn212 (0k) )

[l+nlGll(k)]G22(k) nlG12 (k)
T22(k) =- (l--G2 12 (4.21)
S [l+nG11(k)][l+n2G22(k)] n1n2G12 (k)

T12(k) .-(k 2 -(4.22)
2 [l+n1Gll(k)][l+n2G22(k)] nln2G12 (k)

From the above it is easy to show that the following
are true:

1 + nlGl1(k) [l-n2T22(k)] /Dt(k) (4.23)

1 + n2G22(k) [1-nTll(k)] / Dt(k) (4.24)

1 nlTl1(k) [l+n2G22(k)] / Dg(k) (4.25)

1 n2T22(k) [l+nlG11(k)] / Dg(k) (4.26)

Dg(k) 1/Dt(k) (4.27)

Dg(k) = [l+nlGll(k)][l+n2G22(k)] n1n2G122(k) (4.28)

Dt(k) [1-nlT11(k)][l-n2T22(k)] nln 2122(k) .(4.29)

We now return to the problem of finding expressions
for the ASaB similar to Eq. (2.17). We will do the case
of ASll The others are similar.
By definition, we have

AS11(k) 1(k) T11(k) Gl(k) + Tll(k) (4.30)
Gl(k) Gl(k) AT11(k) (4.31)

From Eq. (4.25), we write

1 + nlG11(k) [Cl-n2 2(k) n2T22(k)] /Dt(k) (4.52)

where Dt(k) is written in an expanded form by replacing
Tg(k) by TS(k) + &a%(k):

Dt(k) [C-nlf Y(k)][l-n2 Tr(k) nn2sr2(k)

n1All (k)[l-n5T2(k)] n2022(k)[l-n1~ (k)]

2n1n2 T1r(k)A2(k) n1n2[A11(k)22(k)

T12(k)] (4.53)

We divide top and bottom of the right hand side of Eq.
(4.32) by Dsr(k) where all the sets of Eq. (4.29) acquire
an "sr" label, and use Eqs. (4.23), (4.24), and (4.27) to

1 + nlG1(k) [l+nlG (k) n2T22(k)D(k)) /C(k)

C(k) 1 nGAT11(k)[l+nlG11(k)] n2 A22l(k)[l1 G22l(k)]

- 2n1n2Gil(k)a 12(k) + n1n2CT11(k)AT22(k) -



Putting G11(k) from Eq. (4.34) into Eq. (4.31) gives,
finally, after some straightforward algebra,

AS11(k) I

s+ r 2- 2-sr2
l+nl%11(k)) ATll(k) + n2 G12 (k)AT22(k)

+ 2n2Gs2(k)1AT2(k)Cl+n r(k)]

- n1n2 CrAT1(k)22(k)-T2(k)]Cl+nG r(k)


Note that if n2 0 or AT12 AT22 0 the above
reduces to Eq. (2.17).

We will take

ATa(k) - lr(k)

as in Chapter II. For purposes of simplification, we will


x Dsr(k) / C(k)
9pr] Ck

restrict the following discussion to those long range
potentials that satisfy the relation

AT1(k)a22(k) i22(k) 0 (4.38)

Since most of the two-component fluids which would be of
interest here are made up of charged particles, such as
fused salts and electrolytic solutions, we may further
restrict the equations above by considering only those
long range potentials that satisfy the relations

0r(r) 022(r) 0(r) 0lr(r) (4.39)

which hold for the Coulomb potential. Then, with Eqs.
(4.37) and (4.39), Eq. (4.36) becomes

[1+ w (k)]20 r(k)
ASII(k) - -
1 + nlCl+w(k)]0 r(k) + n21l+w2(k)]1r(k)

Tll(k) (4.40)


1(k) n1G (k) n12(k) (4.41)

w2(k) n2G22(k) nlG2(k) (4.42)

We have finally, then, by neglecting AB,1(r) in
Eq. (4.14) ,

gll(r) gsr(r)e (4.43)


Cl+wl(k) ]20r(k)
11(k) 1 + nl[+wL(k) (k) + n2l+w2(k) (k)

In a completely similar manner, we obtain

(rg22(r) r)e-22(r) (4.45)
22(r) g2222

gl2(r) g(r)e12(r) (4.46)

where H22 is obtained by interchanging 1 and 2 in Eq.
(4.46), and H12 is gotten from the following equation,

[l+wl(k) [l+w2(k) ]0r(k)
12(k) 1 + n1[l+w(k)]001r(k) + n2[l+W2(k)]~B0r(k)


It is important to note that the above equations for
Hap apply with the restriction on the long range potentials
imposed by Eq. (4.59). The general Hag may be obtained
from Eq. (4.36) and the two other equations similar to it
for AS22 and AS12
The PY-type approximation may be applied as in
Chapter II, giving

g (r) gr(r)e- a()

+ e-O0 (r)
+ e AS (r)

from Eq. (4.7) and its short range equivalent. Applying
the same restrictions on the long range potential as above,
we have the final expressions for the g's ,

S -Blr(r)
gll(r) gr(r)e

sr -00(r)
922(r) g22(r)e

gl2(r) g l2(r)e (r)

-P0(r) lr
+ e [~ (r) Hll(r)] ,


+ e [0(r) 0r(r) H22(r)] ,


+ e 00rH12(r) 0r(r)) ,

where Ha% is defined as above.




The distribution functions g g and g

describing an electrolytic solution may be found using the

results of the previous chapter. We will take as a model

of the electrolyte a system of nonpolarizable ions of charge

+ ze interacting with a hard core short range and a Coulomb

long range potential. Let R and R be the hard core
diameters of the positive and negative ions, respectively.

The interaction potentials are then the following:

oo r R
0sr(r) = (5.1)
0 r > R+

oo, r R
0sr(r) (5.2)
0 r >R

oo, r Z R
0sr(r) (5.5)
0 r >R


- (R+ + R_)/2 (5.4)


ir 22
0++(r) z2e2/Dr (5.5)

0lr(r) = 0r(r) (5.6)

where D is the dielectric constant of the solvent. The

last two equations show that this model satisfies Eq. (4.41),

so that we may use the resulting simplified equations of

the last chapter.

The short range distribution functions may be most

easily obtained from the Percus-Yevick integral equation.

Lebowitz [40] has recently solved this equation analytically

for binary mixtures of hard spheres, generalizing a method

due to Wertheim [41,42]. In Lebowitz's solution, the

direct correlation functions of the hard sphere mixture

are given by the following expressions:

+ + +
{[a+ + b r + dr3] r R+
T+(r) (5.7)
0 r > R

-a d
T (r) -a Cb+4dx+d x3 r Z R (5.8)
+- + r- +-
0 r > R


x r-X.

The constants a b a_, b_, b, d, and X are defined in

Appendix B. We have assumed, for definiteness, that

R R+ in writing the equations above. The equations for

R > R can be easily found by interchanging all + and -

With the direct correlation functions on hand, the

corresponding radial distribution functions are found from

Eqs. (4.19)-(4.21), which relate the Fourier transforms.

Then the total g's may be computed by applying the results

of the last chapter, Eqs. (4.45), (4.47), (4.48), and

(4.51)-(4.53), with the H's determined from Eqs. (4.46) and


This procedure gives a non-iterative scheme for

computing the radial distribution functions describing an

electrolytic solution, in so far as the model chosen repre-

sents such a physical system. A Monte Carlo solution of

this model has recently been carried out by Shaw [43). We

have computed, for purposes of comparison, the ionic radial

distribution functions for the same input conditions as the

Monte Carlo solution, although it is pointed out by Shaw

in his conclusions that the ensemble of 100,000 configu-

rations generated for the Monte Carlo computation "yields

radial distribution functions still showing fluctuations

large enough to preclude strictly quantitative comparison

with theory." That is, the Monte Carlo solution will not

afford more than a good qualitative test, due to insufficient

computer runs.

Shaw considered an electrically neutral collection
of cations and anions of equal size (diameter 4.5 A) and

opposite charge (+ one electronic charge) in a dielectric

continuum with the permittivity of water at 250 C. The

density was that of a 0.00911 molar electrolyte, that is,

a number density of 2.74 x 10-6 ions/A3 for each ionic


These parameters contain, from our point of view,

two unfortunate choices, tending to limit the model features

which may be tested. First, the distribution functions

++ and g__ will be identical, due to the assumption of
equal hard core diameters for the ions. The Monte Carlo

computation thus gives no idea of the effects of different

hard core ion sizes. Second, the very low ion concentra-

tions puts the problem in a region where practically any

theory will tend to give qualitatively correct solutions.

That is, the approximations usually introduced to obtain

tractable equations all tend to be better at small densities.

Solutions at higher densities would "strain" these approxi-

mations more, giving a better idea of their region and

degree of usefulness. At the present concentration, the

effect of the hard cores is limited to providing a cut-off

for the g's in the small r region, the short range g's

being nearly one everywhere outside the hard core diameters.

The results of the calculation are shown in Figure

5.1. The curves labelled DH were computed by Shaw [43]

using a somewhat modified Debye-Hickel expression obtained

by him, where an effort was made to avoid the usual DH

assumptions, viz., (1) that only the interaction between

oppositely charged particles is important, and (2) that the

total number of ions of either charge type is very large.

The comparison with the Monte Carlo solution is, as

expected, not conclusive, although it can be said that the

present solution is at least qualitatively correct. The

g's computed by using a "PY-type" and a "CHNC-type"

approximation, as discussed previously, both gave essentially

the same answers and a single curve has been drawn to

represent both.



o Io o-

I 4. I Z -o,
0 4>


S I o o
| o

o o
I Io I
+ + + I ++ 0
WW bD bO El

0 0

0 L
0' N '

0 \-4

O 0

o 0- a I






Let us recapitulate briefly the results of the

diagrammatic analysis of Chapter I. We found there that

the radial distribution function, g, for a fluid composed

of molecules interacting through two-body forces alone,

may be obtained with the aid of the Ornstein-Zernicke

integral equation, Eq. (1.34), when the direct correlation

function T is known. In addition, T could be written as

in Eq. (1.31), with two excess unknowns, or as in Eq. (1.44),

with only one. This last equation was the limit of the

ability of the diagram technique, and a closed set of

equations for g was not achieved.

This impasse has necessitated the introduction of

approximate integral equations, foremost among which are

the PY and the CHNC equations. We saw that they could be

viewed as two different approximations to the set B ,


B(r) - P(r) (PY) (6.1)


B(r) 0 (CHNC) (6.2)

where P is the set of parallel diagrams.

Numerical solutions of the PY and CHNC equations
have shown that
1) the PY equation is superior to the CHNC equation
for the hard sphere potential [25)
2) the PY equation is superior to the CHNC equation
for the Lennard-Jones (6-12) potential at rela-
tively high temperatures [23]
3) the CHNC equation is probably superior to the PY
equation for the Lennard-Jones potential at
relatively low temperatures [20]
4) the CHNC equation is superior to the PY equation
for the Gaussian model [44].

Khan [20) has shown how these results may be quali-
tatively understood by examination of the first terms in
the diagram expansion of the sets P and B as done
below. [See eqs. (1.25) and (1.26)]. We shall call these

P2 and B2; i.e.,

P2(r) <>o (6.3)

B2(r) (6.4)

At sufficiently low densities, P2 and B2 will describe
the approximate behavior of their respective total sets.
We will suppose in what follows that this is the case.

The difficulty with the set B is already evident

in B2 The presence of the cross bond connecting the

two field points prevents the factorization of the diagram

in either direct space, as with the P set, or in Fourier

transform space, as with the series set S The PY and

CHNC equations circumvent this difficulty by, in effect,

substituting a constant for the cross bond in B2 equal

to -1 and 0 respectively. Consider now the forms of

the Mayer f bonds for the various cases mentioned above.

These are shown in Figure 6.1.

The f bond for the hard sphere potential is precisely

-1 up to the hard sphere diameter, as approximated by the PY

equation. Beyond this point it is zero, but the PY approxi-

mation is not ruined, because the presence of other bonds

connecting the two field points in B2 tends to cut off

the integrals anyway, taking over the task of the real

f bond. Clearly for this case the CHNC approximation is a

poorer one, since it approximates the cross bond by zero


The f bond for the LJ high temperature case is

similar to the hard sphere f the effect of the attractive

well of the LJ potential being minimized by the high

thermal energy of the particles of the system. The same

arguments apply as did for the hard sphere case, and the PT

approximation would appear to be more reasonable.

The LJ low temperature case is characterized by

the large peak in the f bond, resulting from the potential

well. In this case, as the diagram B2 is integrated over

the separation distance between the field points, the

negative contribution picked up at small r will tend to

be cancelled by the positive peak in f at larger r ,

resulting in a smaller net contribution than for the high

temperature case. The net effect will depend on the

particular temperature, but as the temperature is lowered

the change will be toward the CHNC approximation, as noted

by Khan [20).

That the effect of the cross bond is far from being

-1 everywhere is clear in the case of the Gaussian f and

here the CHNC approximation is the more reasonable of the


The above analysis is, of course, simply qualitative,

and it can say nothing of the expected behavior at high

densities. Nonetheless, the following procedure suggests

itself. Approximate the set B by

B(r) mP(r) (6.5)

where m is a constant and P is the parallel set, and let

m be chosen so as to reflect the peculiar circumstances and

potential of the system being studied. For the most part,

m will be expected to be between the CHNC and PY values of

0 and -1 respectively.

The problem is, of course, to find a consistent
method of calculating m for any given potential and
thermodynamic conditions. Several possibilities may be
An appealing and quite general approach would consist
of exploiting the existence of the two expressions for the
pressure. These are, from Appendix C,

= 1 fr~0 (r)g(r)d5r (6.6)


i- v = 1E + n G(r)) /nkT (6.7)

It is shown in Appendix C that these two equations, with
the approximation of Eq. (6.5), give the following condition
for m :
m - 1 /C2 (6.8)

C1 = f[f(r)S(r) + 1 rf'(r)S(r) + 1 nrf'(r)nr r]d3r


C2 f[P(r) + f(r)P(r) + 1 rf'(r)P(r) +

1 nrf (r)Pr) d3 (6.10)

The sets S(r) and P(r) in the definitions of C1 and

C2 may then be obtained approximately from, say, the PY

integral equation.

The above procedure, however, has serious computa-

tional drawbacks. In the first place, the derivatives with

respect to the density are not readily known and it has

proved difficult, at least for the LJ potential, to

obtain this quantity. Secondly, to obtain m by the above

procedure would entail solving another integral equation,

such as the PY equation, as a first approximation. Thus

the advantage of generality entails a corresponding compu-

tational disadvantage.

A more limited condition for m is derived from the

virial expansion for the pressure and has been used by

Hutchinson and Rushbrooke [45] in computing virial coeffi-

cients of a hard sphere gas. We consider the virial diagram

(where all variables are integrated out) derived from B2

and replace the cross bond by a constant; i.e., we write

<:> = m <> (6.11)

The constant m then corresponds, in some sense, to an

average value of the cross bond. It is a parameter dependent

on the particular potential and temperature used to evaluate

the virial diagrams of Eq. (6.11). Writing this equation

out explicitly, we have

l3f32 flf42f 34d3r d3r 2d3r 3d3r4
f 13 f32f 14 f42 34d r 2d r 14
m Jif3f32flf-2d3rld-r2d3r3d-r (6.12)

This is the procedure that has been actually followed,
with the results given in the next chapter. It has the
advantage of simplicity. Once a potential is chosen, m
may be immediately computed. It is the same for any given
isotherm and is independent of density, so that relatively
few calculations of m will be needed for any given problem.
Representative values of m corresponding to the f
functions illustrated in Figure 6.1 as well as the Coulomb
potential, are given in Table 6.1. They are also shown in
Figure 6.1 for the f functions given there. These values
of m tend to confirm the discussion made above. We see
that for the hard sphere and LJ high temperature cases,
m is closer to the PY value of -1 than to the CHNC value of
0 On the other hand, the m's for the LJ low temperature
case and the Gaussian model favor the CHNC approximation.
[The Coulomb m was computed here using an exponentially
shielded, or Debye-Htlckel, potential in order to avoid
divergencies.] The effect of temperature on m for the
LJ potential is shown in Table 6.2.
The new integral equation to be solved for g is


g(r)e0(r) 1 (l+m)P(r)

+ n [g(r')f(r')eS(r')+ (l+m)P(r')]G(Vf-5' j)d r' ,

where P is given in terms of g by Eq. (1.45) and m
is defined by Eq. (6.12).
In the next chapter we discuss the numerical
techniques used in solving Eq. (6.15) and display the
computed results for the hard sphere, LJ and Coulomb






m = -0.35

Fig. 6.1.-Mayer f bonds for various potentials.

-- LJ (Low Temp.)
Gaussian Model

-.. LJ (High Temp.)

-- Hard Sphere



Hard High Low
Spheres Temp. Temp. Coulomb Gaussian

-0.73a -0.81b -0.45 -0.37 -0.35


Hutchinson and Rushbrooke, Reference 45.
Barker and Monaghan, Reference 46.





aFrom Reference 46.



The discussion of the previous chapter has made

plausible the notion that the bridge diagram set B may

be approximated as being proportional to the parallel

diagram set P as written in Eq. (6.5). The PY and CHNC

equations are then but two particular cases of this approxi-

mation in which m is always taken to be -1 and 0,

respectively. We have seen how to allow more flexibility

in the approximation by permitting the conditions of the

problem to determine m rather than assigning it a value

once and for all. The arguments presented, however, merely

made the development reasonable. The resulting integral

equations are notable for their lack of "transparency" and

the ultimate test of any of them lies in the comparisons of

computed answers with some standard solution.

In this chapter we report the results of computations

with Eq. (6.13), which we will call the "m" equation, for

the hard sphere, LJ and Coulomb potentials. These

computed results are compared in each case to the answers

obtained from the PY and CHNC equations, as reported in the

literature cited, and to a solution taken to be the correct

answer. The latter was usually a Monte Carlo solution, with
the exception of the hard sphere equation of state, as noted

below. Comparisons with experiment are not reliable tests
of the accuracy of the computed results so long as the
experimental potential is not accurately known, as is the
present case. The computed functions g and H are given
in tabular form in Appendix D.
Before discussing the results, we will briefly con-

sider some of the numerical techniques used in the solutions.

Once a specific problem has been decided upon, the
parameter m must be found. This is done most easily in
the following manner. Define the function

F(rl2) ff(r13)f(r32)dr (7.1)

where f is the Mayer function. The Fourier transform of
this equation reads

F(k) f2(k) (7.2)

In terms of F Eq. (6.12) becomes
m- --- (7.3)
which defines the average of f over the distribution F2
The numerical procedure is then as follows. The function f
is computed for the given potential and temperature, and its

Fourier transform is obtained using the techniques of

Appendix A. The transform of F is then computed from
Eq. (7.2) and the Fourier inversion is carried out. Finally
Eq. (7.3) is evaluated using Weddle's rule [47]. An
interesting check on the accuracy of the numerical procedure
is afforded by the Gaussian f function, which permits m
to be computed analytically. The analytical and numerical
computations of m both gave the same number, m = -0.35355,
to five significant figures. This accuracy is not in general
to be expected, however, unless the interval is very small,
as the Gaussian is an especially smooth function.
The "m" equation was solved in its Fourier transform
representation by an iterative scheme. We define first
the quantities

H(r) g(r)e~0(r)- (7.4)

Q(r) g(r)e (r) g(r) (7.5)
P'(r) P (l+m)P(r) (7.6)

The functions Q, P', and G may all be written in terms
of H as follows:

Q(r) g(r)e (r)[1 e-(r)] (7.7)

f(r)[l + H(r)] (7.8)

P'(r) (l+m) [H(r) ln[l + H(r)]] (7.9)


G(r) H(r) Q(r)

Then Eq. (6.13) may be written in terms of H, Q, and P ,

H(r) P'(r) + n [P'(r) Q(r')][H(I i'I )

Q(Ir )]d3r' (7.11)

Using the convolution theorem we may write the Fourier
transform representation of this equation as

E(k) ''(k) + n[-'(k) Q (k)][5 (k) Q (k)] (7.12)

whence we get finally that

H(k) (7.1)
1 + n[Q (k) P'(k)]

Equations (7.13), (7.8), and (7.9) can be solved by
iteration. The procedure is as follows. An initial guess,
Ho is made for H This is used to compute Qo and Po
from Eqs. (7.8) and (7.9). The Fourier transforms of Qo
and P are then found and used in Eq. (7.13) to obtain
a new H denoted H1 This quantity H1 is then an
input quantity for the computation of Q1 and P1 which
in turn are used to calculate a new H2 and so on, for
any desired number of cycles.
The problem is considered solved when a solution is
self-consistent. That is, when the new H calculated is
identical to the H used as input. A measure of the degree

of self-consistency is given by the root-mean-square (rms)

difference between successive H's where the rms is

defined as follows:

N 2 ] 1/2
rms N [ k+ ) H(A) (7.14)

Here N is the number of points in the numerical solution,

A is the interval between points, and Hk is the value of

H after the kI iteration. In practice, a solution was

considered converged, i.e., self-consistent, if the rms

between successive iterations was of the order of 10-3 or

less. An additional criterion was that the input and out-

put H's should agree to at least two significant figures

at all points.

For small densities the rate of convergence is quite
fast, a self-consistent solution usually being achieved in

about 10 iterations. For larger densities, the slow rate

of convergence may be improved by the following device, due

to Broyles [48]. When H is sufficiently close to self-

consistency, it is found to approach its final value

exponentially. Then from any three successive iterations,

k-l, k, and k+l, we can compute the value that H would

have after an infinite number of iterations, viz.,

H2k(JA) Hk+l(jA)Hk-l(A)
ext( 2 Hk(JA) Hkl(JA) HkI(J) (7.15)

This extrapolated H is then used as the new H for the
input of the next iteration. Some care must be taken in
employing this trick, however. If it is used too soon, i.e.,
before the solution has settled down to an approximately
exponential rate of approach, it will cause the succeeding
iterations to "blow up."
An additional feature that is sometimes useful in
speeding convergence is the so-called mixing process, where
we write, after the kt iteration,

knew(JA) a Hk(JA) + (l-a)Hkl(jA) (7.16)

Normally the constant a is set equal to one. When a
solution tends to oscillate, however, a smaller a will
tend to dampen the oscillations and improve the convergence

The Lennard-Jones Potential

The LJ potential is a reasonably realistic potential
which still retains the advantages of an analytic formula-
tion. It is given by the expression

0(r) C [(a/r)12 (a/r)6] (7.17)

where the constants E having the dimensions of energy, and

a with the dimensions of length, are to be determined

from experiment. Expressions involving the LJ potential

may be put in dimensionless form by taking E as the unit

of energy and a as the unit of length. Then we write

PZ(r) 4(1/r6-1)/T*r6 (7.18)

where the dimensionlesss) reduced temperature T* is

defined by

T* kT/E (7.19)

The corresponding dimensionless density is n* na3 .

Wood and Parker [27] have used the Monte Carlo (MC)
method to calculate the radial distribution function, the

pressure, and the energy for a fluid interacting with a

LJ potential at several densities of the T* 2.74 iso-
therm. Broyles et al. [23,49] have computed the solutions

to the PY and CHNC equations at these same densities and

temperature. We have thus chosen these cases also for
purposes of comparison. The computed value of m was

-0.732 for this case.
The results of the "m" equation and the Monte Carlo
solution are shown for T* = 2.74 and n* 2/5, 5/6, and

1.0 in Figures 7.1-7.3, respectively. The agreement is in

general quite good, being poorest at the first peak. The

PT equation also yields good solutions for these cases,

while the CHNC equation is in poorer agreement.

All four plotted curves would lie confusingly close

and they have been compared in tabular form in Tables 7.1-

7.3. These tables give a measure of the deviation from the
MC 'solution of the g's computed from the integral equations
in the form (gC/g)-l as in Reference 49. All three

computed g's are in essential agreement at the reduced
density of 2/5 and the differences in Table 7.1 are not very

meaningful. More substantial deviations appear in Table 7.2

for the 5/6 density case, as well as for the 1.0 density in

Table 7.5. It is seen that the present equation gives the

closest agreement in general, followed in turn by the PY

and then the CHNO equations.

Once the radial distribution function is known the

pressure p and energy E may be easily found from the
expressions [1]

p/nkT = 1 (2zn/3) f 0(r)g(r)r3dr (7.20)

E/NkT- E '/NkT (7.21)

2m J 0(r)g(r)r2dr (7.22)

Values of pressure and temperature for the three cases

mentioned above are shown in Table 7.4. The MC values were
chosen from the various computations reported by Wood and

Parker [27] on the basis of the run having the largest















Ln -


r (gMC/gm)-1 (gMC/gpY)-1 (gMC/gCHNC)-1
0.90 -0.027 -0.016 -0.074
0.95 -0.020 -0.012 -0.049
1.00 -0.066 -0.060 -0.079
1.05 -0.043 -0.041 -0.048
1.10 -0.027 -0.029 -0.026
1.15 -0.011 -0.013 -0.006
1.20 0.003 0.002 -0.010
1.5 0.003 0.001 0.008
1.4 0.006 0.005 0.007
1.5 0.001 0.001 -0.003
1.6 -0.009 -0.006 -0.013
1.7 -0.020 -0.016 -0.024
1.8 -0,012 -0.013 -0.014
1.9 -0.019 -0.014 -0.019
2.0 -0.012 -0.012 -0.012
2.1 -0.003 -0.003 -0.002
2.2 0.002 0.000 0.003

aThe PY data are taken from Broyles, Reference 49.
The CHNC data is also due to Broyles (private communication).

* "m"


Fig. 7.2,-Radial distribution function for the LJ potential:
T* = 2.74, n* = 5/6.


r (g6O/gm)-l (gM/gp )-1 (g C/gCHNC)-1

0.90 0.714 0.772 -0.117
0.95 0.055 0.078 -0.141
1.00 -0.089 -0.083 -0.138
1.05 -0.037 -0.047 -0.001
1.10 -0.021 -0.041 0.065
1.15 0.004 -0.020 0.102
1.20 0.007 -0.016 0.089
1.3 0.020 0.020 0.040
1.4 0.007 0.060 0.007
1.5 0.011 0.036 -0.031
1.6 -0.013 0.003 -0.045
1.7 0.009 0.014 -0.010
1.8 0.010 0.006 -0.001
1.9 -0.004 -0.007 -0.009
2.0 -0.014 -0.021 -0.003
2.1 -0.002 -0.010 0.023
2.2 0.000 0.004 0.018

aThe PY data are taken from Broyles, Reference 49.
The CHNC data is also due to Broyles (private communication).





1.5 2.2



0 MG

0.5 II
1.0 1.5 2.0

Fig. 7.3.-Radial distribution function for the LJ
potential: T* = 2.74, n* = 1.0.


r (g
0.90 0.295 0.333 -0.304
0.95 -0.004 0.024 -0.214
1.00 -0.042 -0.042 -0.031
1.05 -0.068 -0.083 0.080
1.10 -0.050 -0.081 0.145
1.15 0.001 -0.033 0.169
1.20 0.047 0.021 0.157
1.3 0.102 0.145 0.031
1.4 0.055 0.116 -0.082
1.5 -0.009 0.023 -0.097
1.6 -0.037 -0.042 -0.080
1.7 -0.012 -0.030 -0.025
1.8 -0.013 -0.031 -0.017
1.9 -0.024 -0.033 -0.012
2.0 -0.007 -0.014 0.051
2.1 0.000 -0.002 0.053
2.2 0.045 0.063 0.044

aThe PY data are taken from Broyles, Reference 49.
The CHNC data is also due to Broyles (private communication).



n* 2/5 5/6 1.0

(p/nkT)MCa 1.2-1.5 4.08 6.97
(p/nkT)m 1.24 4.09 6.99
(p/nkT)pyb 1.24 4.01 6.8
(p/nkT)CHNCb 1.28 5.11 9.1

(E'NkT)MCa -0.86 -1.57 -1.60
(E2NkT)m -0.863 -1.58 -1.615
(E2NkT)pyb -0.865 -1.61 -1.67
(EBNkT)CHNCb -0.859 -1.40 -1.19

aFrom Reference 27.
bFrom Reference 23.


number of moves per particle, which might be expected to

give the best average. Here again the present equation

gives the closest answers, though the differences from PT

are small.

The results for the "m" equation in Table 7.4 include

a small correction for the truncation of the integrals in

Eqs. (7.20) and (7.22). This truncation is unavoidable in

machine calculations and was corrected for in the following

way. If the maximum value of r in the numerical solution

is R we write, e.g.,

E'/nkT = 27n B (r)g(r)r2dr
+ 2mn f (r)g(r)r2dr (7.23)

-(E'/nkT)comp + AE (7.24)

Here the first term on the right is the value obtained on

the computer. The correction,

AE 2Tn B (r)g(r)r2dr (7.25)
may be approximated by setting g equal to one everywhere

in the region. If R is sufficiently large this approxi-

mation will be quite good. For the LJ potential, then,

we easily obtain

AE 8T (1 1/R6) (7.26)

8- (7.27)
in reduced units. A similar correction was applied to the
machine solutions for the pressure. For the n* 1 case,
for example, we took R 4.5 (in units of a). This gave

a computed E'/nkT of -1.582 and a correction AE of

-0.033, or about 2 per cent of the total.

The Coulomb Potential

We consider here, following Carley [50,51], a system
of classical point charges embedded in a rigid uniform
background of charge of the opposite sign, such that the

net charge is zero. The pair potential is just the Coulomb

0(r) e2/r (7.28)

where e is the charge on an electron. For convenience we
will use the same reduced units employed in the above
references, viz., the unit of length is taken to be the
ion sphere radius a given by

a [3/(4m)]1/3 (7.29)

and the dimensionless reduced "temperature" is

9 kTa/e2 (7.30)

so that we may write