NEW METHODS OF COMPUTING RADIAL
DISTRIBUTION FUNCTIONS OF FLUIDS
By
FRED LADO
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
August, 1964
A mis padres.
PREFACE
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 twocomponent 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 PercusYevick (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 wellknown Percus
Tevick (PY) or convolutionhypernetted chain (CHNC) integral
equations upon the proper choice of the value of m Solu
tions of the equation have been obtained numerically for
the LennardJones, 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.
ACKNOWLEDGMENTS
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.
vii
TABLE OF CONTENTS
Page
PREFACE . . . . . . . . . . iii
ACKNOWLEDGMENTS . . . . . . . . ... vi
LIST OF TABLES. . . . . . . . . . x
LIST OF FIGURES . . . . . . . . ... xi
PART ONE
PERTURBATION CORRECTION OF THE
RADIAL DISTRIBUTION FUNCTION
FOR LONG RANGE FORCES
Chapter
I. CLUSTER EXPANSION OF THE RADIAL DISTRIBUTION
FUNCTION . . . . . . . . . 2
II. LONG RANGE FORCES. . . . . . . ... 15
III. A NUMERICAL TEST: THE GAUSSIAN MODEL. ... 26
IV. GENERALIZATION TO A TWOCOMPONENT FLUID. . 43
V. A NUMERICAL APPLICATION: ELECTROLYTIC
SOLUTIONS. . . . . . . . ... 55
PART TWO
IMPROVEMENT OF THE INTEGRAL EQUATIONS FOR
THE RADIAL DISTRIBUTION FUNCTION
VI. A NEW INTEGRAL EQUATION. . . . . ... 62
VII. NUMERICAL SOLUTIONS. . . . . . ... 72
The LennardJones Potential. . . . ... 77
The Coulomb Potential. . . . . ... 88
The Hard Sphere Potential. . . . ... 89
viii
Page
APPENDICES
A. FOURIER TRANSFORMS . . . . . .. 96
B. PARAMETERS IN THE LEBOWITZ SOLUTION FOR A
MIXTURE OF HARD SPHERES. . . . . .102
C. THE VIRIAL AND COMPRESSIBILITY EQUATIONS OF
STATE. . . . . . . . . ... .105
D. THE COMPUTED FUNCTIONS g AND H FOR THE LJ,
COULOMB, AND HARD SPHERE POTENTIALS ... .108
LIST OF REFERENCES . . . . . . . . 117
BIOGRAPHICAL SKETCH. . . . . . . . ... .121
LIST OF TABLES
Table
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
Page
71
71
81
83
85
86
SI111
112
. 113
. 114
S115
. 116
LIST OF FIGURES
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
PART ONE
PERTURBATION CORRECTION OF THE RADIAL DISTRIBUTION
FUNCTION FOR LONG RANGE FORCES
CHAPTER I
CLUSTER EXPANSION OF THE RADIAL DISTRIBUTION FUNCTION
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
where
S 1/kT (1.2)
Here T is the absolute temperature and k and h are
Boltzmann's and Planck's constants, respectively. The
Nparticle Hamiltonian H describing the system is a function
of the 3N components of momentum and the 3N components of
position.
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 twobody forces. The integration
over the moment in Eq. (1.1) may be carried out immediately,
yielding
Q Z/NIA35 (1.5)
where
Z = .. ee d3 r ... d3rN (1.6)
and
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)
i
where
rij i j (1.11)
and 0 (r) is the pair potential. With this, Eq. (1.9)
becomes
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
32")
(1.13)
The 3particle, 4particle, etc., distribution functions
can be similarly defined. In general, the nparticle
distribution function is given by
g.(n) n f... d rn+ ... d3rN
(1.14)
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 N2 particles.
It has the additional property that its Fourier transform
may be directly measured by xray 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 wellknown expansion of g(r) in powers of
the density [49]:
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
t1 J
(1.20)
where
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
2.
(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)
and
G(r) S(r) + P(r) + B(r)
+ f(r) [1 + S(r) + P(r) + B(r)] (1.28)
or
g(r)eB0(r) 1 + S(r) + P(r) + B(r) (1.29)
We will denote the set of all nonnodal 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)]
(1.30)
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 12 bond. Diagrams
of this type are put in explicitly by the definitions, as
in the last term of Eq. (1.28). Thus also the nonnodal set
T defined by Eq. (1.30) includes not only the sets P and
B but also all diagrams with a direct 12 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
is
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
that
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 nonequal 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
(1.58)
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
oo
P(r) = Sr) B(r
t2
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)
(1.43)
which can be used to eliminate P(r) from previous equations.
(1.39)
(1.40)
(1.41)
(1.42)
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 [1116]
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)
or
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 PercusYevick
(PY) integral equation [1720]
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 BornGreen
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 LennardJones (612) 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
dissertation.
For a realistic potential made up of a repulsive core
and an attractive part such as the LennardJones (612)
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.
CHAPTER II
LONG RANGE FORCES
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 Waalslike equation of state is
based upon this type of potential [28,29]. A onedimensional
model of this type has been studied by Kac, Uhlenbeck, and
Hemmer in a succession of papers [3032] 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 er/ (2.1)
where
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 noniterative
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
Then
S(k) 8p(k) 00(k)/[l + n00(k)] (2.7)
and
g(r) exp (k k sin krdk (2.8)
2rr 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
DebyeHtckel (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;
i.e.,
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)expC00r(r) + AS(r) + AB(r)] (2.14)
where
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 + nGSr(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)eH(r) (2.18)
with
H(k) + nr(k)2 r(k (2.19)
1 + nCl + nGr(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
C353.
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) + eBSr(r)flr(r)
= fSr(r) + Af(r) .
(2.20)
(2.21)
If we replace each f bond in a typical diagram of the
expansion of T(r) for example,
S = f12f13f32d3r
(2.22)
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
represents
0^;0 + 'o *
(2.23)
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) + 0o + 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
taking
flr(r) ~0lr(r) .
(2.27)
That is,
AT(r) 0lr(r) (2.28)
We can, however, retain the whole first diagram, giving
instead
AT(r) eBsrflr(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,
(2.18),
lng(r) Ingsr(r) H(r) (2.32)
In the region where g(r) is close to one, we may approxi
mate
lng(r) g(r) 1 (2.33)
lng r(r) gsr(r) 1 .
(2.34)
Then Eq. (2.32) becomes
G(r) Gsr(r) H(r) (2.55)
or
1 + nG(k) 1 + nr(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
(2.37)
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 PYtype 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)
(2.39)
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) + e0(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)e0lr(r) + eBO(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
satisfactory.
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 CHNCtype approximation,
AB = 0 we were led to Eq. (2.18). A PYtype 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
comparison.
CHAPTER III
A NUMERICAL TEST: THE GAUSSIAN MODEL
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 socalled 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)
in
g(r) ee(r) 1 + E t(r) (3.1)
tl
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. )
or
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
2
00(x)  In(l e1 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 el121(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) eO(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) .
(3o11)
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)eikr d3r (3.13)
Carrying out the angular integration,
00
F(r) l f kF(k) sin kr dk (5.1i)
oo
F(k) rF(r) sin kr dr (3.15)
0
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 J0
N
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) j0
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)eH(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)
where
oo
D = r(k)  k2dk (3.25)
2n 1+ n~B r(k)
ad 1 + 4 n GSr(x) x2dx (3.26)
and
A(x) e sr(x t+ ngtr(x) (5.27)
t=l
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.13.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 PYtype 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
3.5.
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
34
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 PYtype approxi
mation and the abovementioned 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
I I I
Id I "
0
I
0
Sr
rl
o
*d
O
ue
"4
o
O 0
o
.4
0
'I
o
0
a3
O
0
U\
'4
o
o* 0
000
4
0 0 0
O 0
O o
,H
.4
0 4,
Kd
aa
^ 0 I I I
N 00
0
4 .
I 1 0 I
O *
\
\
ie
o
a
o
cd
0
4,
.10
O
.I
rl
0 T
4,
4.
r< O
uI '*
I I
R 38
o
lr; a
x C
o
0
0
4
o
P
0
4
4,
o
rt
Sd
* 0
0
'4
P4
a!
0 K
o b
fci
2.5
001r
~ 0sr
2.0
m=0.9
m0.8
1.5
.5 m=0.6
1.0 m=0.6
m=0.7
0.5 m0.8
m10.8
m=0.9
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.13.4.
0.02 
rms
0.01 
0 "
1.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).
0.02
rms
.01 
n = 0.3
,
1.0
0.02
rms
0.01
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
line).
___
0.03
rms
rs Eq. (2.28)
0.02 Eq. (2.29)
 Eq. (2.31) /
0.01
0
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.
CHAPTER IV
GENERALIZATION TO A TWOCOMPONENT FLUID
The equations of Chapter II may be readily
generalized to a twocomponent 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( )
li
+ Z 12(rj) (4.1)
i=l,N1
jN1+l1,N
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 onecomponent 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)
X.1,2
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
(4.8)
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+ . .
(4.9)
B,(r) nl2 + nn2 ) + n2 c + ..
(4.10)
where the circles and squares again refer to components 1
and 2, respectively. [Cf. Eqs. (1.24) and (1.26)]. Then
the equation
oo
71
Pa(r) misa(r) + Ba (r))m (4.11)
m2
[see Eq. (1.39)] becomes, for P1(r) ,
P11(r) 81 S2(r) + Sll(r)B11(r)
+ B11(r) + ... (4.12)
nl2 + n1n2 + n2 + ...
(4.13)
from Eq. (4.9). These terms satisfy the definition of the
parallel diagrams for the twocomponent 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
(4.14)
where
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)
1=l,2
These equations may be solved simultaneously for either
G in terms of TB or vice versa. In the first case
they give
Eln2T22(kl)T(k) + n2T12 (k)
G1(k) . (4.17)
) [lnlT11(k)]ln2T22(k) nln2T12 (k)
[lnlTll(k)]T22(k) + nlT12 (k)
G22(k) 2 (4.18)
[ln1T11(k)][ln2T22(k)] nln2T12 (k)
T12(k)
12(k) .T12(k) .(4.19)
1 n1 11(k)][ln2T22(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) = (lG2 12 (4.21)
S [l+nG11(k)][l+n2G22(k)] n1n2G12 (k)
G12(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) [ln2T22(k)] /Dt(k) (4.23)
1 + n2G22(k) [1nTll(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)
where
Dg(k) = [l+nlGll(k)][l+n2G22(k)] n1n2G122(k) (4.28)
Dt(k) [1nlT11(k)][ln2T22(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)
.sr
Gl(k) Gl(k) AT11(k) (4.31)
From Eq. (4.25), we write
1 + nlG11(k) [Cln2 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) [Cnlf Y(k)][ln2 Tr(k) nn2sr2(k)
n1All (k)[ln5T2(k)] n2022(k)[ln1~ (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
obtain
1 + nlG1(k) [l+nlG (k) n2T22(k)D(k)) /C(k)
(4.34)
where
C(k) 1 nGAT11(k)[l+nlG11(k)] n2 A22l(k)[l1 G22l(k)]
 2n1n2Gil(k)a 12(k) + n1n2CT11(k)AT22(k) 
AT12(k)]Dgr(k)
(4.35)
Putting G11(k) from Eq. (4.34) into Eq. (4.31) gives,
finally, after some straightforward algebra,
AS11(k) I
s+ r 2 2sr2
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)
(4.36)
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
(4.37)
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 twocomponent 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)
where
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)
where
Cl+wl(k) ]20r(k)
11(k) 1 + nl[+wL(k) (k) + n2l+w2(k) (k)
(4.44)
In a completely similar manner, we obtain
(rg22(r) r)e22(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)
(4.47)
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 PYtype approximation may be applied as in
Chapter II, giving
g (r) gr(r)e a()
+ eO0 (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)] ,
(4.49)
P0(r)
+ e [0(r) 0r(r) H22(r)] ,
(4.50)
+ e 00rH12(r) 0r(r)) ,
(4.51)
where Ha% is defined as above.
(4.48)
CHAPTER V
A NUMERICAL APPLICATION: ELECTROLYTIC SOLUTIONS
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
+
where
 (R+ + R_)/2 (5.4)
and
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 PercusYevick 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
where
x rX.
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 
+
labels.
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
(4.49).
This procedure gives a noniterative 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
0
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 106 ions/A3 for each ionic
species.
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 cutoff
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 DebyeHickel 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 "PYtype" and a "CHNCtype"
approximation, as discussed previously, both gave essentially
the same answers and a single curve has been drawn to
represent both.
60
0
0.0
o
o Io o
I 4. I Z o,
0 4>
44>
A
S I o o
 o
o o
4
I Io I
+ + + I ++ 0
WW bD bO El
14
4
0 0
0
0 L
00
0' N '
0 \4
O 0
o 0 a I
PART TWO
IMPROVEMENT OF THE INTEGRAL EQUATIONS FOR THE RADIAL
DISTRIBUTION FUNCTION
CHAPTER VI
A NEW INTEGRAL EQUATION
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 twobody forces alone,
may be obtained with the aid of the OrnsteinZernicke
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 ,
namely
B(r)  P(r) (PY) (6.1)
and
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 LennardJones (612) potential at rela
tively high temperatures [23]
3) the CHNC equation is probably superior to the PY
equation for the LennardJones 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
everywhere.
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
two.
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
considered.
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)
and
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)
where
C1 = f[f(r)S(r) + 1 rf'(r)S(r) + 1 nrf'(r)nr r]d3r
(6.9)
and
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 Jif3f32flf2d3rldr2d3r3dr (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 DebyeHtlckel, 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
then
69
g(r)e0(r) 1 (l+m)P(r)
+ n [g(r')f(r')eS(r')+ (l+m)P(r')]G(Vf5' j)d r' ,
(6.13)
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
potentials.
2.0
f(r)
1.0
0
1.0
m = 0.35
Fig. 6.1.Mayer f bonds for various potentials.
 LJ (Low Temp.)
Gaussian Model
.. LJ (High Temp.)
 Hard Sphere
TABLE 6.1
REPRESENTATIVE VALUES OF THE PARAMETER m
LJ LJ
Hard High Low
Spheres Temp. Temp. Coulomb Gaussian
0.73a 0.81b 0.45 0.37 0.35
aFrom
bFrom
Hutchinson and Rushbrooke, Reference 45.
Barker and Monaghan, Reference 46.
TABLE 6.2
EFFECT OF TEMPERATURE ON m FOR THE LJ POTENTIAL
kT/C
0.6
0.8
1.0
1.2
1.3
1.333
2.0
2.74
4.0
8.0
20.0
m
U.049
0.253
0.454
0.589a
0.625a
0.634a
0.679a
0.732
0.797a
0.813a
0.768a
aFrom Reference 46.
CHAPTER VII
NUMERICAL SOLUTIONS
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
oo
2(r)f(r)r2dr
m  (7.3)
SF2(r)r2dr
0
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)
and
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)
and
(7.10)
G(r) H(r) Q(r)
Then Eq. (6.13) may be written in terms of H, Q, and P ,
giving
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
selfconsistent. That is, when the new H calculated is
identical to the H used as input. A measure of the degree
of selfconsistency is given by the rootmeansquare (rms)
difference between successive H's where the rms is
defined as follows:
N 2 ] 1/2
rms N [ k+ ) H(A) (7.14)
jo
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., selfconsistent, if the rms
between successive iterations was of the order of 103 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 selfconsistent 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,
kl, 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)Hkl(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 socalled mixing process, where
we write, after the kt iteration,
knew(JA) a Hk(JA) + (la)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
rate.
The LennardJones 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/r61)/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.17.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)
o
and
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
0
I
80
LIn
a
f\
U
E
P
P
0
4'
a
0
o
0
to
4D
o
.i
d
.,
4
0
VP
r(
Ln 
TABLE 7.1
COMPARISON OF PT, CHNC, AND "m" INTEGRAL EQUATION SOLUTIONS
TO MONTE CARLO FOR THE LJ POTENTIAL: T* 2.74, n* 2/5a
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"
MC
Fig. 7.2,Radial distribution function for the LJ potential:
T* = 2.74, n* = 5/6.
TABLE 7.2
COMPARISON OF PT, CHNC, AND "m" INTEGRAL EQUATION SOLUTIONS
TO MONTE CARLO FOR THE LJ POTENTIAL: T* 2.74, n* 5/6a
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).
2.2
2.8
2.0
2.5
1.5 2.2
1.0
1.0
"m"
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.
TABLE 7.3
COMPARISON OF PY, CHNC, AND "m" INTEGRAL EQUATION SOLUTIONS
TO MONTE CARLO FOR THE LJ POTENTIAL: T* = 2.74, n* 1.0a
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).
THERMODYNAMIC
TABLE 7.4
QUANTITIES FOR THE LJ
n* 2/5 5/6 1.0
(p/nkT)MCa 1.21.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.
POTENTIAL
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.,
R
E'/nkT = 27n B (r)g(r)r2dr
oo
00
+ 2mn f (r)g(r)r2dr (7.23)
R
(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)
R
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)
3TER
8 (7.27)
3TR5
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
potential
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
