
Citation 
 Permanent Link:
 https://ufdc.ufl.edu/UF00097932/00001
Material Information
 Title:
 New methods of computing radial distribution functions of fluids
 Added title page title:
 Radial distribution functions of fluids
 Creator:
 Lado, Fred F., 1938
 Publication Date:
 1964
 Copyright Date:
 1964
 Language:
 English
 Physical Description:
 xi, 121 leaves. : ill. ; 28 cm.
Subjects
 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 )
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis  University of Florida.
 Bibliography:
 Bibliography: leaves 117120.
 Additional Physical Form:
 Also available on World Wide Web
 General Note:
 Manuscript copy.
 General Note:
 Vita.
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 nonprofit 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:
 021508206 ( AlephBibNum )
13143435 ( OCLC ) ACW7097 ( NOTIS )

Downloads 
This item has the following downloads:

Full Text 
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

Full Text 
PAGE 1
NEW METHODS OF COMPUTING RADIAL DISTRIBUTION FUNCTIONS OF FLUIDS By FRED LADO A DISSERTATION PRESENTED TO THE GRADUATE COUNQL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THB DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA August, 1964
PAGE 2
UNIVERSITY OF FLORIDA 3 1262 08552 5565
PAGE 3
A mis padres. 11
PAGE 4
PREPACB 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 taien up in the second chapter, where two expressions 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 generalization to twocomponent systems of the expressions derived in Chapter II and their numerical solution for a model of an electrolyte. ill
PAGE 5
In prittclple, the radial distribution function for any fluid may be found "exactly" by the Monte Carlo (MO) 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 advemtages over MC that they provide a means of computing experimental potentials from scattering data and may indeed have analytic solutions, as witnessed by Werthelm'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 PercusYevick (PY) or convolutlonhypemetted chain (CHNC) Integral equations upon the proper choice of the value of m . Solutions of the equation have been obtained numerically for the LennardJones, Coulomb, and hard sphere potentials and It
PAGE 6
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 dissertation is bj chapter and serial order therein. Thus Bq. (5.11) refers to the eleventh equation in Chapter V. The same rule applies to tables and figures.
PAGE 7
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. Sahlln, A. A. Khan, and D. D. Carlej, Wherever possible, specific acknowledgmenta 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 heartfelt 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 reseaxch 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 between man and computing machine. All calculations were i
PAGE 8
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 meinuscript, 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
PAGE 9
TABL3 OF CONTENTS Page PREPACS Hi ACKNOWLEDGMENTS vl LIST 0? TABLES x LIST 0? FIGURES xi PART ONE PEHTUHBATION 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. ... ^3 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 10
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 ix
PAGE 11
LIST OP TABLES Table P^SÂ® 6.1 Representative Values of the Parameter m . . . 71 6.2 Effect of Temperature on m for the LJ Potential 71 7.1 Comparison of PT, CHNC, and "m" Integral Equation Solutions to Monte Carlo for the LJ Potential! T* = 2.7^, n* 2/5 81 7.2 Comparison of PY, CHNC, and "m" Integral Equation Solutions to Monte Carlo for the LJ Potential: T* = 2.7^, n* = 5/6 83 7.5 Comparison of PY, CHNC, and "m" Integral Equation Solutions to Monte Carlo for the LJ Potential: T* 2.7'^, n* = 1.0 85 7.^ ThermodyTiamic Quantities for the LJ Potential 86 D.l The Computed Functions g aind H for the LJ Potential 109 D.2 The Computed Functions g and H for the Coulomb Potential Ill D.3 The Computed Functions g and H for the Hard Sphere Potential 112 D.4The Computed Functions g and H for the HsLTd Sphere Potential 115 D.5 The Computed Functions g ajad H for the Hard Sphere Potential 11^ D.6 The Computed Functions g and H for the Hard Sphere Potential 115 D.7 The Computed Functions g and H for the Hard Sphere Potential 116
PAGE 12
LIST OF FIGURES Figure Page 3.1 Radial distribution functions for the Gaussian model: n 0.^, m = 0.9 55 5.2 Radial distribution functions for the Gaussian model: nÂ»0.4,mÂ«0.8 36 3.5 Radial distribution functions for the Gaussian model: n = 0.'t,mÂ»0.7 57 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 ^ 3.7 RMS deviations ^1 3.8 Density dependence of the RMS deviations. . . . A2 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.7^, n* = 2/5 80 7.2 Radial distribution function for the LJ potential: T* = 2.7^, nÂ» Â» 5/6 82 7.3 Radial distribution function for the LJ potential: T* > 2.7^, n* 1.0 8* 7.4 Radial distribution function for the Coulomb potential: Q 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 zl
PAGE 13
PART ONE PERTURBATION CORRECTION OF THE RADIAL DISTRIBUTION FUNCTION FOR LONG RANGE FORCES
PAGE 14
CHAPTEH I CLUSTBH EXPANSION 0? 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 J J i=i 1=1 where 3 lAT . (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 $N components of position. The free energy P of the system is then given by PF . In Q (1.3)
PAGE 15
3 and all other thermodTnamic quantities can be obtained from the free energy. For example, the energy E of the system is given by E d(3P)/d0 . (1.^) Another quantity that may be used to compute the thermodynamic quantities is the pair distribution function g(r, tPo) (defined below) if the particles are assumed to interact only through twobody forces. The integrations over the momenta in Eq. (1.1) may be carried out immediately, yielding Q Z/NI A^^ (1.5) where and 3U 3^ ^3 e d^r^^ ... d^rjj (1.6) A h(27nnkT)"^/^ . (1.7) Here U(r, , . . . ,rj,) is the potential energy of the system and V the volume. Wo have assumed that all particles have the same mass m so that the kinetic energy is T Â» 5^ p. /2m . Then we may write 0F Â» ln(NIA^^) + In Z (1.8) and, from Eq. (1.4),
PAGE 16
4^ ENkT + /.../e U d^r^^ . . . d^v^ , (1.9) Implementing the assumption of padrwise interaction, we write U = L (r.,) , (1.10) i
PAGE 17
5 ?or 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., eiv^tv^) 6(^12^ Â• (1.15) This has given rise to the use of the name "radial distribution function." The normalization of the pair, or radial, distribution function is \ sMd^T 1 I . (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 psirticles. It has the additional property that its Fourier transform may be directly measured by xray and neutron scattering experiments [2,5]. Ignoring terms of order 1/N , g(r) is defined as g(r) ^ r... fe^^ d^r^ ... d\ . (1.1?) 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
PAGE 18
of Eq. (1,17) is expanded in Mayer cluster integrals. This results in the wellknown expansion of 6(r) in powers of the density [^9] : G(r) Â« g(r) 1 (1.18) f(r) + [1 + f(r)]C(r) . (1.19) Here C(r) is defined by oe ^ C(r) ^ n*i; ^^ n f^. d^r^ ... d^r,^^ t1 J (1.20) where fij fir^.) = 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 L . The integrals may be represented by diagrams according to the following prescription. For each f. , 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 s 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
PAGE 19
integrands having identical products of fj ^ . 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 CD + /A*(^*/^ ^  > Â• Â• . (1.22) where the density and symmetry factors have been explicitly written. Any diagram can easily be written out explicitly in terms of the integrations it represents by applying the prescription given above in reverse order. Thus we have <^ / ^15^5^^1^^^2^S^H Â• (123) 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 emd 2 must pass. (b) A parallel diagram contains at least two paths between 1 and 2 which are connected only at 1 and 2.
PAGE 20
(c) A bridge diagram is one which is neither series nor parallel. The sum of all series, parallel, and hridge diagrams will "be denoted by S(r) , P(r), and B(r), respectively. To the second power in density these sets are S(r) Â» ncr\)^ n P(r) I n2 qAo /~\ *,^*/^ + Â• Â• n2 ^ + o and or B(r) 5 In terms of these sets, we have C(r) Â» S(r) + P(r) + B(r) G(r) S(r) + P(r) + B(r) + f(r) [1 + S(r) + P(r) + B(r)] g(r)e^^^^^ 1 + S(r) ^ P(r) * B(r) . + . . . (1.24) (1.25) (1.26) (1.27) (1.28) (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)Cl + S(r) + P(r) ^ B(r)] (1.50) g(r)f(r)e^^^^^ * P(r) ^ B(r) (1.31)
PAGE 21
and thus G(r) S(r) + T(r) . (1.32) It is important to note that the sets S, P, and B do not tncludÂ© 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 Bq. (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 18 S(r) n / nT^^)(iir^2^d}r^ , (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 subdia grams. The same applies for the integral over r, ,
PAGE 22
10 with Eq. (1.52) we have thus an integral equation for G(r) , Q(r) T(r) + n / T(r')G(r r'  )dV . (1.5^) This is the integral equation of Ornstein and Zernioke [10] which gives G(r) when the direct correlation function T(r) is known. Unfortunately, we can see from Bq. (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 Bq. (1.5^) 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 Bqs. (1.32) and (1.53) we have that S(k) nT^(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
PAGE 23
d^l [ c/N, 11 2 (1.37) The branch diagrams may be members of the sets S and B , but not of P . We gire Interested in the contribution of all pareillel diagrams with t branches, which we will denote by P^Cr) . Consider a typical diagram with t branches and let the contribution from the iU) branch be a^(r) , so that the contribution of the whole diagram is Ca3_(r)][a2(r)] . Ca^(r)] . We can generate a subset of p^(r) by holding a2(r) , . . . ,a^(r) constant and allowing ou(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)] ... [a^(r)] . A still larger subset of P^(r) , which includes the one above as a subset of itself, may be found by allowing the second brsinch to be composed of any member of S(r) or B(r) as well as the first, giving the contribution I CS(r) + B(r)]CS(r) + B(r)][a^(r)] ... [a^(r)] . Diagrams with nonequal subdiagrams are counted twice in the multiplication of the first two bracketed terms and the
PAGE 24
12 factor 1/2 removes this redundancy. Diagrams vrith 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 P^(P) . CS(r)^ B(r)3^ ^ (1.58) For the parallel set P(r) , t may be any number from 2 to infinity, so that CO tÂ»2 (1.39) e^^""^ ^ ^^""^ S(r) B(r) 1 . (1.^) Thus, from Bq. (1.29), g(r)eP^(^) . eS(^) * ^(r) ^^^^^^ In g(r)e 30(r) Â» S(r) + B(r) . (1.^2) Then substituting Eqs. (1.^1) and (l.'<2) into (1.^0) givea a useful relation for P , viz.. P(r) g(r)e 3lZ)(r) 1 In g(r)e 30(r) (1.^3) which can be used to eliminate P(r) from previous equations,
PAGE 25
13 Thus the direct correlation function T "becomes, from Eq. (1.51). T(r) G(r) ln[g(r)e^^'^^^ ] > B(r) . d.'^^) 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) ln[g(r)e^^^^^ ] . (1.^5) Inserting this into Eq, (1,5^) gives the convolutionhypemetted chain (CHNC) integral equation [1116] G(r') In g(r')e^^^^'^ ln[g(r)e0^(^)] n X G( r r'l )d^r' . (1.^) Another approximation to B which has proved useful is the assumption that B(r) ^ P(r) (1.^7) or B(r) P(r) Â«= . (1.^8) Neglect of these two sets gives T(r) g(r)f(r)e^^^''^ (1.^9) as can be seen easily from Bq. (1.51). Combined with
PAGE 26
1* Eq. (1.5^) Â» this approximation leads to the PercusYevick (PY) integral equation [1720] B(r)e^^(r) = 1 + n TgCr' )f (r' )e^^(^' ^ X G(r r')d^r' . (1.50) A third and older integral equation for g which is often found in the literature is the equation of BomGreenYvon (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 [2^,25]. Ve 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 . [25] 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.
PAGE 27
CHAPTBH 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,2?]. 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 , p' Â« IAt' , to be a perturbation of 30 at T and applying the corresponding 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. 15
PAGE 28
16 Thus the change in the potential could be the addition of a wealc long range interaction to a short range repulsive core. The derivation of a van der Vaalslike 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 C 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 Htlckel [55] and reads g(r) exp where pe^ rA (2.1) \ [^ t; nPe^]"^/^ (2.2) and Â• is the electron charge. This result can easily be obtained from the development of the last chapter. Prom Eq. (1.^1) we get sMe^^^'^ . eS(^) (2.5) by neglecting B(r) , as in the CHNC approximation. The corresponding direct correlation function T(r) is then given by Bq. (1.^5) t which can be written
PAGE 29
17 T(r) G(r) ln[l ^ G(r)] 00(r) . (2.4) ?or a weak long range force, G(r) will be small over a large region so that the logarithm in Eq. (2,^) may be expanded. Thus, T(r) 00(r) + I G^(r) + . . . . (2.5) If only the first term is retained we obtain a noniterative solution for g by combining Eq. (2.5) with S(k) n[00(k)]2/Cl ^ n30(k)] , (2.6) obtained from Sq. (1Â»56) and the above approximation to T Â• Then S(k) 00(k) 30(k)/[l 4n00(k)] (2.7) and g(r) exp CO 2nrr Jq 1 Â» n00(lc) k sin krdk (2.8) where the angiilar integrations in the Fourier transform have been carried out. For the Coulomb potential, P0(k) ^ n3e?^^ . (2.9) (See Appendix A). By inserting this expression into Bq, (2.8) and performing the integration C5^], we obtain the DebyeHttckel (DH)g . Bq. (2.1).
PAGE 30
18 The two essential approximations of the DH expression are thus B(r) ^ (2.10) T(r) ^ P0(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) 0^^(r) + 0^^(r) . (2.12) We will assume that g corresponding to the short range potential is a known function, which we call g . From Eq. (1.^1) we may write g^^(r) expCS^^(r) + B^^(r)] , (2.13) where the label "sr" on S and B indicates they are the sets correspondinLg to , Simileirly, for the complete potential , we may write another equation identical to Eq. (1.^1). Then dividing this equation by Eq. (2.15) we have formally
PAGE 31
19 g(r) g'*^(r)exp[p0'^'(r) + AS(r) + AB(r)] (2.1^) where AS(r) S(r) S^'^(r) (2.15) and AB(r) is similarly defined. Prom Bqs. (1.55) and (1.56) we have further that 2 ^srp Ai(k) '''^ W ^T W (2.16) 1 nT(k) 1 nT^^(k) 1 n[l+nGÂ®^(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 30 . This leads then to the final expression g(r) g^^(r)e^^^) (2.18) with H(k) . [1 ^ ^"wAi'^C^^ . 1 * n[l * nGÂ«'(k)]00^=^(k) ^ ^^'
PAGE 32
20 Here g(r) is given completely in terms of the long range potential suid the short range g , which is the desired form C55]. 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) f^^(r) + e^^^^^'^^^f^^Cr) (2.20) f^^(r) + Af(r) . (2.21) If we replace each f "bond in a typical diagram of the expansion of T(r) , for example, o^ = /f3_2fi5f32^^^5 ' ^^^^^ S3? by the sum of f and Af and expand the products, we obtain the following, ^ ,Â•> ,^. A, (2.23) Here a wavy line connecting two points indicates an ar Q0^^ Ir f bond and a dashed line represents Af e ^ f
PAGE 33
21 Since we have assumed that ^ 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 AP bond in each diagram, and so on. The sum of the first set gives immediately T^^(r) . We write ,sr T(r) T'''(r) + 00 + n CTÂ—to ^ o*^'^~0 . A> + ... . (2.2^) The sum of all diagrams having one AP bond is the desired AT . The previous approximation to AT may be recovered by taking only the single bond diagram after T as AT , writing it in two parts, 0.0 eP<^''"(^)fl^(r) (2.25) f^^(r) + f^^(r)f^^(r) (2.26) and further approximating by neglecting the second part and taking f^^(r) P0^''(r) . (2.27)
PAGE 34
22 That is, AT(r) ^ P0^^(r) . (2.28) We can, however, retain the whole first diagram, giving instead AT(r) e"^^ f^^(r) . (2.29) Further, all diagreoas with a direct Af bond between 1 and 2 may be summed to give AT(r) =Af(r)g^^(r)e^^ ^^^ (2.50) g^^(r)f^^(r) . (2.31) It will turn out, however, that these higher approximations, Eqs. (2.29) and (2.51)i are not as satisfactory as the original approximation, Eq. (2.28). This will be discussed further in the next chapter. Consider now the logarithm of both sides of Bq,, (2.18), Ing(r) Ing^'^(r) H(r) . (2.52) In the region where g(r) is close to one, we may approximate Ing(r) ^ g(r) 1 (2.53) Ing^'^(r) =g^'^(r) 1 . (2.3^)
PAGE 35
Then Eq. (2.52) becomes , sr G(r) G'=*^(r) H(r) or 23 (2.35) 1 tnG(k) 1 h iiG^^(k) 1 + nCl + nG^^(k)]3?^^(k) (2.36) G(k) i 1 K QG^^(k) T^rsv, :^rF, 1 1 + nCl ^ nG^'(k)]00^^(k) Â• (2.37) This is precisely the equation obtained by Broyles, Sahlin, and Carley (BSO) by a method based on collective coordinates C36,373. The BSC equation, Eq. (2.57), 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 Ing(r) is not valid in this region. If we now further restrict the region of r under sr consideration to that where g is practically one, then G in Eq. (2,35) may be neglected, giving G(r) H(r) . (2.38) This equation was obtained by Eemmer as the first order expression for g(r) in the large r region [38], We note also that in the limit (r) Â» bo that
PAGE 36
24 0"^^(r) Â— 0(r) , Bqs. (2.18) and (2.19) give the DH g when is specified as the Coulomb potential. An equation for g(r) employing a PYtype approximation, where AP(r) + AB(r) is neglected, is easily obtained. From Eq. (1,29), we write g^^(r)e^^^''<^) 1 ^ S^^(r) * P^^(r) ^ B^^(r) . (2.39) A similar equation holds for the total g . Subtracting one from the other then gives g(r)eP<^(^) Â» g^^(r)eP^'''(^) ^ AS(r) > AP(r) + AB(r) (2.^0) g(r) g^^(r)eP<^^''(^) . e'^^^^^ t^M (2.^1) 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 00 ^ we have finally g(r) g^^(r)e^<^^''(^> * e'^^^^^ [00^^(r)H(r)] ,(2./^3) where H(r) is the same as in Eq. (2.19). Once again higher approximations to AT may be incorporated into thÂ«
PAGE 37
25 final result, but it la found that the results are less satisfactory, Ve have found two expressions which give g(r) for the total potential 0(r) iZl^^(r) + 0^^(r) (2.4^) sr when g Is known. Using a CHNCtype approximation, AB = , we were led to Eq. (2,18), A PYtype approximation, AB + AP =*= , gave rise to Sq, (2 A3). 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.
PAGE 38
CHAPTBH III A mJMEHICAL TEST: THE GAUSSIAN MODEL Tiie 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 sr an exact g 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 negative Gaussian function. This model was used to carry out numerical computations of Eqs. (2.18) and (2.^3), as well as of the equations resulting from the higher approximations to AT . The BSO 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 [593. 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. 26
PAGE 39
27 classified the series and bridge diagrams having five or less field points and evaluated these for the Gaussian model. This data wsis then used to obtain the coefficients ga.(r) t in g(r) 30(r) oo . Z n^ t1 'g^(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 n ScC^) should have an almost negligible effect on the total g(r) . The Mayer f function is given by f(r) e^^/*)^ (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., ^ d ^oo / f(r)i^dr / r^ dr / e'^^^^^ r^ dr Jo Jo Jo so that or d^ T a5n^/2 ui 1/3 Ti a 1.100 a (3.3) (3. A) (5.5)
PAGE 40
28 The reduced units are then given bj r/d , n Nd^A . (3.6) The potential corresponding to a solution g(x) is given by .2 00(x) ln(l e'^*^^ ^ ) . (3.7) We can obtain a g(x) for a slightly different potential in the following way. Perform a coordinate transformation sr, x' Â» x/m , where t1 (nm5)*g^(x/m) . (3.10) The long range potential is given by the equation Ir sr> (x) 30(x) 30*'^(x) . (3cll)
PAGE 41
29 The two quantities on the right are computed directly from Eqs. (5.7) and (5.9). This long range potential and the short range g computed from Helfand's data by means of Eq. (5.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 Bq. (5.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 , x Â„ , must be chosen. The range ' max ' Â° between and x is then divided into N intervals of equal size A . A function f(x) is then considered determined if it is known at the N + 1 points x. = jA , J J = 0,IT . The calculations of 00(jA) , 00Â®^(jA), 00^^(JA) ancL s(jA) are straightforward and require no comment. To compute g (jA) we first recompute g(jA) for the density 3 nm and then treinsform to a new coordinate i = j/n to obtain g (iA) g(jA/m) . Depending on the size of m , i will reach its maximum value at some point before i does, leaving g (lA) undetermined for the remaining points. In this region g (iA) was assigned the value one, since
PAGE 42
30 SI* in every case g had already settled to its asymptotic veilue 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) ^r / F(k)e^^'^ d\ (5.12) (2n)5J F(k) /F(r)e"^^*^ d^r . (3.13) Carrying out the angular integrations, F(r) ~L> / lEp(k) sin kr dk (3.1'*) rF(r) sin kr dr . (5.15) For numerical calculation F(r) and F(k) must be made discrete functions. Put r iAr , k = jAk . (3.16) Then we have
PAGE 43
31 P(i) ^2^^ L dP(J) sin (ijATAk) (5.17) 2n lAr J0 N F(J) ^^^7^ ^ iJ'Ci) sin (ijATAk) . (3.18) Â»0 ^S^ioC It is shown in Appendix A that in order to conserve the orthogonality of the sine functions in this discrete representation, and thus the basic relation of a function to its Fourier transform, we must put ArAk 2k/(2N^1) (3.19) where N is, as above, the nvimber of points. Then, putting A 2n/(2N>l) , (3.20) we have finally P(i) Â— T^^ r Z dF(d) sin (iJA) (3.21) 2n''i(Ar)^ j0 ^i p(j) . i^JlUu) ij.(i) sin (ijA) . (5.22) Eq. (3.22) is used to compute the Fourier transforms of 30(i) and GÂ®^(i) which then determine H(J) . The discrete version of H(x) is then found with Eq. (3.21), giving the final answer S(dA) g^^(jA)e^^J^^ . (3.23)
PAGE 44
52 Hemmer's small r solution for g(x) , Eq. (35) of Reference 58 was solved in the following form: g(x) g^^(x) D[g^^(x) ^ A(x)] (3.2^) where 2it^ J 1+ nTl0^ (k) k^dk , (3.25) oo and 7/ 1 + /4n n / G^^(x) x^dx , (3.26) .sr. 5 A(x) eP^ <) Y. ^^ V^Cx) . (3.27) t=l Calculations were made with Eqa. (2,18) and (2,^3), the BSC equation, and Hemmer's two asymptotic equations, for densities of 0.1, 0.2, 0.3, and 0.^. 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 (x) . The results for n o 0.^ and m = 0.9, 0.8, 0.7, and 0.6 sire shown in Figures 3.13.^. 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
PAGE 45
53 two aaymptotlc equations yield good agreement in the small and large x limits for which they are applicable hut 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 expledned on the basis of cancellations between the sets P and B [20]. The potentials for these cases are shown in Figure 5.5. The results of other cases for Bqs. (2,18) and (2.^3) 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 xrms a N i I Cg^^jA) 6^Â°Â°^P(dA)]' 1/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.^3) 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
PAGE 46
3^ by 30 . The basis of this assertion lies in the comparison of results for the Gaussian model of the three approximations for AT , Eqs. (2.28), (2.29), and (2.51). The situation is exemplified by Figure 5.8, where the rms values for the equations resulting from the PYtype approximation 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.
PAGE 47
4Â»
PAGE 48
CO 36 o B H O o w n a p o a d o Â•rl P O o H Â•H Pi +3 m rH Â•d Â•% M I rvj ba
PAGE 49
o
PAGE 50
O 58
PAGE 51
39 ?ig. 5. 5. Gaussian model potentials. These ootentials were used in obtaining the g's of Figures 3.15.^,
PAGE 52
^ 0.01 1.0 0.8 0.6 0.02,rms 0.01 _ 1.0 0.8 0.6 Pig. 3.6. RMS deviations. The rms deviations from the exact Gaussian model g(x) for g computed from Eq. (2.18) (dashed line) and Eq. (2.45) (solid line).
PAGE 53
41 0.02 ,rms 0.01 0.6 u.u<^
PAGE 54
42 0.03rrms 0.02 0.01 Eq. (2.28) Bq. (2.29) Sq. (2.31) Fig. 3. 8. Density dependence of the RMS 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 a is 0.6.
PAGE 55
CHAPTER IV GENERALIZATION TO A TWOCOMPONENT FLUID The equations of Chapter II may be readily generalized to a twocomponent fluid, such as an electrolytic solution or a fused salt. For such a system, three radial distribution functions axe in general necessary. If the fluid components are labelled 1 and 2, the required distributions are Bii% S22' ^^^ ^12 Â° ^21* Ve will assume, as in Chapter I, that the potential U of the system may be written as a sum of pair interactions, as follows: ^^^1 V Z. %C^ij) ^ L <^22^^iJ^ lÂ£i
PAGE 56
i^4 2 from N, + 1 to N, the total number, where N Nj_ + N2 . ('<.2) Again we assume that we may write where J25^g is the short range interaction between components a and & , and so on. In terms of diagrams, the pair correlation fxmction G(r) is again described by an equation like Bq. (1.28), where now three different G's must be distinguished. Ve write Â»a3<^> Â• ^aP^^) * ^^ * ^aP^^^^^^aP^^) ^ ^aP^^) + B^p(r)] , a,P 1,2 {^A) where S n 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 P . Similar definitiona 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 P will always be understood to take on the values 1 and 2.
PAGE 57
*5 Direct correlation functions '^ncS^'^^ ^^Â® defined in complete amalogy with the onecomponent fluid, viz., but now the expressions for the series diagrams become ^c^&M ^ n^ / T^(r')G^3(r?')dV . (4.6) \l,2 J Here n, and np 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 S g are necessarily members of the a and 3 components, one from each, the intervening field points may be of either component 1 or 2, The field point r' in Bq. (4,6) is occupied by a member of component \ . It is treated as a fixed point by the diagrams in the sets T , 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, gTi(r) . Here we have o^'^o + no Â©^''^o ni \j \j T ixp
PAGE 58
46 where the circles denote component 1 and the squares component 2. In the one component fluid there was only one diagram for this coefficient. [See Bq. (1.22).] Equation (^.^4) may be rewritten to read 30^B(r) The parallel set may be eliminated as before to yield ga3(r)e ^^ e ""^ ^*^ . (4.8) That this is so can be made plausible by considering the first few terms of a specific g, say giT . Then we have n, o'^'^ >o + nÂ« cf\ + n, cr >o S^j^(r) Â» n^ o^ >o+ np 1 12 "*" 12 x> + np o^ td + n, o^^^^'^o + ^1^2 ^^^^^^ + ^1^2 ^^^^^^^ + n2 o^^^'>o+ Hi cr^^^D + n, np o"''^^"^^ + nj^n2 0^'*''^+ n2^c/*^^+ . . . (4.9) and
PAGE 59
^7 (4.10) where the circles and squares again refer to components 1 and 2, respectively. [Of. Bqs. (1.24) and (1.26)]. Then the equation oo Pa0<) Z irrfSapC^) * B^g(r)]Â» (4.11) mÂ«2 [see Bq. (1.59)] becomes, for ^^^^^(r) , P^3^(r) I ^x^M + S^^(r)B^3_(r) + \ \^W + ... (^.12) I ^1^ <^ ^ ^1^ <^ * I ^2^ <^ * Â• Â• Â• (^.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 Bq. (4.8) may be written, %ri.th all quantities labelled with an "sr" superscript. Then the ratio of the two equations gives, as in Chapter I, (4.14)
PAGE 60
where and a similar definition holds for AB Â« . ocp Before proceeding to the evaluation of the quantities AS Q and AB q , we need to obtain explicitly the ap ap relations "between the Fourier transforms of the three G Â« ap and the three T g implied by Eqs. (^.5) and (^.6). By taking the Fourier transforms of these equations, with the aid of the convolution theorem, we get the following: \l,2 These equations may be solved simultaneously for either they give Go in terms of T , or vice versa. In the first case ap ap [lnpTp2(lc)]T.,(k) + npT,p^(k) G,(k) :r^^ ^: ^ ^^ ^ M C^.l?) ^^ [ln3^T^^(k)][ln2T22(k)3 n3^n2T^2^^) [ln,T,,(k)]Tpp(k) + n,T,p^(k) Qpp(k) =iii ^1 ^ ^^ ^ '^ (^.18) ^"^ [ln3_T^3^(k)][ln2T22(k)] n^n2T^2^^^ T,p(k) 0,p(k) Â» ^ ^^:z :rÂ— 2 .(^.19) ^'^ [ln3^T^3^(k)][ln2T22(k)] n3^n2T^2^^)
PAGE 61
*9 For the inverse solution we obtain [l+npGpp(k)]G,,(k) npG,p^(k) T,(k) zr^^^ ~ ^ ^ '3 ^^.20) ^^ [l+n3^G3^j^(k)][l+n2G22(k)] n3^n2G3^2^^) [l+n,G,,(k)]Gpp(k) n,G,p^(k) J (k) ^ ^= (^,21) 22 [l+n;_G3_3_(k)][l+n2G22(lc)] n^n2G^2 (^^ G,p(k) rp (y\ , iÂ£ .(^1.22) ^2 [l^.n3_G^^(k)][l+n2G22(k)] ^^^2^12^(1^)' Prom the above it is easy to show that the following are true: 1 + n^^G^^^Ck) . Cln2T22(k)] / D^(k) (^.23) 1 + n2G22(k) [ln^T3^3^(k)]/D^(k) (4.2^) 1 n^^T^^Ck) . [l^.n2G22(k)]/Dg(k) (^.25) 1 n2T22(k) [l+n^G^^^Ck) ] / Dg(k) (^.26) Dg(k) l/D^(k) (^.27) where Dg(k) [l+n^G]^^(k)][l+n2G22(k)] n3^n2G\2^(l'^) (^.28) 7; 2, D^(k) [ln3_T3_3^(k)][ln2T22(^)^ n^n2T3_2^(k) .(^.29)
PAGE 62
50 We now return to the problem of finding expressions for the AS^g similar to Bq. (2.17). We will do the case of ^^11 Â• The others are similar. By definition, we have Aa^l(k) G^3_(k) \^W mW + f^i(k) (^.50) \^W mW ^\^W . (^.31) Prom Eq. (^.25), we write 1 + ih.^ii(k) Â» [ln2T^(k) n2AT22(k)] /D^(k) (^.32) where 0+;^^'^ ^^ written in an expanded form by replacing D^(k) Cln3^T^J^(k)][l.n2T^(k)] n^n2T^^(k) n3^AT^3^(k)[ln2T(k)] n2AT22(k) [ln^TjJ(k)] 2n^n2T^2^k)AT^2^^^ ' nj^a2[AT3^3^(k)AT22(k) ATj^g^Ck)] . (^.33) We divide top aind bottom of the right hand side of Eq. (^.52) by D^^(k) , where all the sets of Eq. (^.29) acquire an "sr" label, and use Eqs. (^.23), (4. 2^), and (^.2?) to obtain
PAGE 63
51 1 + n^G^^^Ck) [l+n^^G^JCk) n^^T^^WD^'^W} / C(k) where C(k) 1 n^AT^3^(k)[l+n3^Gj[(k)] n2AT22(lc) C1+ii2G2(1c)] 2nj^H2^12^^^^^12^^^ * n^n2[AT^j^(k)AT22(lc) Putting Gj^.j^(k) from Bq. (^^.3^) into Eq. (^.51) gives, finally, after some straightforward algebra, Â£^S^^W [l+n^G^^(k)]^AT3^l(k) + n2^G^2^(k)AT22(k) + 2ii2G^2^k)AT3^2^^^'^^'^^l^^^^ n^n2^^^11^^)^^22(k)AT^2^(k)][l4.ii^G^^(k)] X D^^(k) / C(k) . (^.56) Note that if 1I2 Â» or AT,2 ^^^22 * ^ *^Â® above reduces to Eq. (2.17). We will take AT^pCk) . B0llW (^.37) as in Chapter II. For purposes of simplification, we will
PAGE 64
52 restrict the following dtacussion to those long range potentials that satisfy the relation AT^3^(k)AT22(k) AT^2^(k) , (4.58) 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 0llW 0llM 0llM 0^^(r) (4.39) which hold for the Coulomb potential. Then, with Eqs. (4.57) and (4.59), Eq. (4.36) becomes [lHw,(k)]23?^^(k) ASTT(k) ::=i :=T^ :^y^^ 1 + n^[l+Wj^(k)]3r^(k) + n2Cl+W2(k)]30^^(k) a\iW C^.^o) where w^(k) n3^GjJ(k) n2Gj2^k) (4.41) W2(k) n2G2(k) nj^G^2(^) Â• (4.42) We have finally, then, by neglecting AB^n(r) in Eq. (4.14) ,
PAGE 65
53 fir ^ll^^^ gll(r) sllMe ^^ (^.^3) where [l+w.(k)]230^^(k) In a completel7 similar manner, we obtain gl2(^) 812^^)Â® (^Â•''^^ where Hp2 is obtained by interchanging 1 and 2 in Eq. (^.^6), and H,2 ^s gotten from the following equation, [l4ir (k) ] [l^Wp(k) ] 0^w H, o(k) ^ ;=r=^ = ^^rz: Â• ^2 1 + n^[l+vij_(k)]30^^(k) + n2[l4W2(k)]P0^^(k) It is important to note that the above equations for H g apply with the restriction on the long range potentials imposed by Eq. (^.59). The general H^g may be obtained from Eq. (^.36) and the two other equations similar to it for AS22 and AS, 2 Â• The Pltype approximation may be applied as in Chapter II , giving
PAGE 66
54 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 , sr '^^^""M 30(r) ,_ gll(r) gii(r)e + e C30^''(r) E^^M:\ , (4.49) sr, s 0^^''(^) 30(^) ir g22(r) g2(r)e + e [30^''(r) H22(r)] , (^.50) and (^.51) where H q is defined as above. cxp
PAGE 67
CHAPTER V A NUMERICAL APPLICATION: ELBCTROLTTIC 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: (5.1) where and 0llM
PAGE 68
56 0^^(r) z^e^/Dr 0lf(r) 0^f(r) , (5.5) (5.6) where D is the dielectric constant of the solvent. The last two equations show that this model satisfies Bq. (^.^1), so that we may use the resulting simplified equations of the last chapter. The short range distribution functions maj be most easily obtained from the PercusTevick integral equation. Lebowitz [^0] has recently solved this equation analytically for binary mixtures of hard spheres, generalizing a method due to Wertheim [^1,^2], In Lebowitz' s solution, the direct correlation functions of the hard sphere mixture are given by the following expressions: (5.7) T (r) < [a b r dr^] , r ^ R Â± Â± Â± , r > R^ T_^_(r) a Â«a^ Â—[b+^Xdx+dx^^] Â» X ^ r Â£ R^_ (5.8) , r > R^^ where X r \ . The constants a,b,a,b__,b,d, and \ are defined in Appendix B. We have assumed, for def initeness, that
PAGE 69
57 R ^ H in vn^itiag the equations above. The equations for R > R_ can be easily found by interchanging all + and labels. With the direct correlation functions on hand, thÂ« corresponding radial distribution functions are found from Eqs. (^.19)(^.21), which relate the Fourier transforms. Then the total g's may be computed by applying the results of the last chapter, Eqs. (^.^5), (^.^7), (^.A8), and (^.51)(^.53)Â» vdth the H's determined from Bqs. (^.46) and (^.^9). This procedure gives a noniterative scheme for computing the radial distribution functions describing an electrolytic solution, in so far as the model chosen represents such a physical system, A Monte Carlo solution of this model has recently been carried out by Shaw [^^3]. 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 concluslonfl that the ensemble of 100,000 configurations 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.
PAGE 70
58 Shaw considered an electrically neutral collection o of cations and anions of equal size (diameter ^.5 A) and opposite charge (> one electronic charge) in a dielectric continuum with the pemnittivity of water at 25" C. The density was that of a 0.00911 molar electrolyte, that is, 6 2 a number density of 2.7^ x 10 ions/A^ 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 ft 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 concentrations 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 approximations 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 haird core diameters.
PAGE 71
59 The results of the calculation are shovm in Figure 5.1. The curves labelled DH were computed by Shaw [^33 using a somewhat modified DebyeHtickel 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 "Pltype" and a "CHNCtype" approximation, as discussed previously, both gave essentially the same answers and a single curve has been drawn to represent both.
PAGE 72
60
PAGE 73
PART TWO inPROVMENT OF THE INTEGRAL EQUATIONS FOR THE RADIAL DISTRIBUTION FUNCTION
PAGE 74
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 OrnsteinZemicke integral equation, Eq. (1.5A), when the direct correlation function T is known. In addition, T could be written as in Eq. (1.51), 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 PT 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) (CHNC) , (6.2) where P is the set of parallel diagrams, 62
PAGE 75
65 Numerical solutions of the PT and CHNC equations have shown that 1) the PY equation Is superior to the CHNC equation for the hard sphere potential [253 2) the PT equation Is superior to the CHNC equation for the LennardJones (612) potential at relatively high temperatures [25] 5) the CHNC equation is probably superior to the PY equation for the LennardJones potential at relatively low temperatures [20] ^) the CHNC equation is superior to the PY equation for the Gaussian model [^^]. Khan [20] has shown how these results may be qualitatively 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 Bp ; 1 . e . , At sufficiently low densities, Pp and Bp will describe the approximate behavior of their respective total sets. We will suppose in what follows that this is the case.
PAGE 76
6^ The difficulty with the set B is already evident in Bp . 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 PT and CJHNC equations circumvent this difficulty by, in effect, substituting a constant for the cross bond in Bp > equal to 1 and , 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 haird sphere potential is precisely 1 up to the hard sphere diameter, as approximated by the PT equation. Beyond this point it is zero, but the PT approximation 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 CHNO 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.
PAGE 77
65 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 Bp 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 and 1 , respectively.
PAGE 78
66 Tlie 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 t f r30'(r)g(r)d5r (6.6) and V ll Â» [1 ^ 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 C3_/C2 (6.8) where C^^ /[f(r)S(r) + I rf'(r)S(r) *  nrf '(r)^i^]d5r (6.9) and C2 Â» / [P(r) + f(r)P(r) ^  rf '(r)P(r) ^ 1 arf'(r)^2iÂ£l]d5r (6.10) b fjn Â•
PAGE 79
67 The sets S(r) and P(r) in the definitions of C, and Gp may then be obtained approximately from, say, the PT integral equation. The above procedure, however, has serious computational 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 PT equation, as a first approximation. Thus the advantage of generality entails a corresponding computational disadvantage. A more limited condition for m is derived from the virial expansion for the pressure and has been used by Hutchinson and Rushbrooke [^5] in computing virial coefficients of a hard sphere gas. We consider the virial diagram (where all variables are integrated out) derived from Bp and replace the cross bond by a constant; i.e., we write Â•" <> (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
PAGE 80
68 / ^13^52^1^^2^5^'^^^1^V^^3'^^^^ m Â» Â— r Iw. , (6.12) T 13 52" 14" ^2^ "1^ "2" "3 ^ 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 einy 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 PT value of 1 than to the CHNC value of . 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 DebyeHtickel , 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
PAGE 81
69 g(r)Â«^^^^^ 1 (l4.in)P(r) ^ n / [g(r')f(r')e^^^'''^+ (l4m)P(r ') ]G(trr'  )d5r' , (6.13) where P is given in terms of g by Eq. (1.^3) and m is defined bj Bq. (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.
PAGE 82
70 f(r) LJ (Low Temp.) Gaussian Model LJ (High Temp.) Hard Sphere 1.0 Fig. 6.1. Mayer f bonds for various potentials.
PAGE 83
71 TABLE 6.1 REPRESENTATIVE VALUES OF THE PARAMETER m Hard Spheres
PAGE 84
CHAPTBH 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 approximation 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 PT and CHNC equations, as reported in the literature cited, and to a solution taken to be the correct 72
PAGE 85
73 answer. The latter was usually a Monte Carlo solution, with the exception of the hard sphere equation of state, as noted helow. 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 consider 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 ^(^12> U(.r^^)fir^2^d^T^ (7.1) where f is the Player function. The Fourier transform of this equation reads F(k) f2(k) . (7.2) In terms of P , Eq, (6.12) becomes oo ^(r)f(r)r^dr P^(r)r^dr m ^B Â— Â» (7.3) i 2 which defines the average of f over the distribution F 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
PAGE 86
7^ Appendix A. Th.e transform of P Is then computed from Bq. (7.2) and the Fourier inversion is carried out. Finally Bq. (7.3) is evaluated using Weddle's rule [^73. An interesting check on the accuracy of the numerical procedure is affoorded 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^<^(^)l (7.4) Q(r) g(r)eP^^^^g(r) (7.5) and P'(r) (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)eP<^(r)[l eWr)^ (^.^^ f(r)[l + H(r)] (7.8) p'(r) (1+m) [H(r) ln[l + H(r)]] (79) and G(r) H(r) Q(r) . (7.10)
PAGE 87
75 Then Eq. (6.13) may "be written In terms of H, Q, and P , giving H(r) P'(r) + n / [P'(rO Q(r')][H(? r'\ ) q(r t'\ )]d5r' . (7.11) Using the convolution theorem we may write the Fourier transform representation of this equation as H(k) P'(k) + n[P'(k) Q (k)][H (k) Q (k)] (7.12) whence we get finally that ^ ^ [Q (k) P'(k)] H(k) Q(k) ^ ~ . (7.13) 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, H , is made for H . This is used to compute Q and P from Eqs. (7.8) and (7.9). The Fourier transforms of Q and P are then found and used in Eq. (7.13) to obtain a new H , denoted H, . This quantity H, is then an input quantity for the computation of Q, and P , , which in turn are used to calculate a new E^ Â» 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
PAGE 88
76 of selfconsistency Is given by the rootmeansquare (rms) difference between successive H's , where the rms is defined as follows: rms ^ 2 T 1/2 (7.1^) Here N is the number of points in the nximerical solution, A is the interval between points, and H, is the value of H after the kti> iteration. In practice, a solution was considered converged, i.e., selfconsistent, if the rms between successive iterations was of the order of 10 ^ or less. An additional criterion was that the input and output 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 [A8]. When H is sufficiently close to selfconsistency, it is found to approach its final value exponentially. Then from any three successive iterations, k1, k, and k+1 , we can compute the value that H would have after an infinite number of iterations, viz.,
PAGE 89
77 H^l^CdA) Hj^^^(jA)H^^^(dA) 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 ktt iteration, Hjj.''Â®*(jA) a Hj^(dA) + (la)Hj^^3^(dA) . (7.16) Normally the constant a is set equal to one. VHien 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 formulation. It is given by the expression Hr) ^e[(a/r)^2 (a/r)^] (7.17) where the constants Â£ , having the dimensions of energy, and
PAGE 90
78 a f with the dimensions of length, are to be determined from experiment. Expressions involving the LJ potential may be put in dimensionless form by taking C as the unit of energy and a as the unit of length. Then we write 30(r) '<(l/r^l)/T*r^ (7.18) where the (dimensionless) reduced temperature T* is defined by T* kT/e . (7.19) The corresponding dimensionless density is n* Â• na^ . 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.7^ isotherm. Broyles et al . [23 'Â•9] have computed the solutions to the PT and GHNO 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.7^ and n* Â» 2/5, 5/6, and 1.0 in Figures 7.17.5, 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 OHNC equation is in poorer agreement.
PAGE 91
79 All foxLT plotted curves would lie confusingly close and they have been compared in tabular form in Tables 7Â»17.3. These tables give a measure of the deviation from the MO solution of the g's computed from the integral equations in the form (g^^^/g)! Â» as in Reference ^9. 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 PT and then the CHNC equations. Once the radial distribution function is known the pressure p and energy E may be easily found from the expressions [1] 00 p/nkT 1 (27ni/5) T 30'(r)g(r)r^dr (7.20) and B/KkT I E'/NkT (7.21) 00 2Tm r 30(r)g(r)r^dr . (7.22) Values of pressure and temperature for the three oases mentioned above are shown in Table 7.^. The MC values were chosen from the various computations reported by Wood and Peirker [27] on the basis of the run having the largest
PAGE 92
80 lA \ Â« Â— Â• a OJ H P Â« P O Pi Â•3 ^ (D P r^ O a o H +J O o =1 ,Q H Pi P <0 H I
PAGE 93
81 TABLE 7.1 COMPARISON OP PY, CHNC, AND "m" INTEGRAL EQUATION SOLUTIONS TO MONTE CARLO FOR THE LJ POTENTIAL: T* Â» 2.7^, n* = 2/5^ r
PAGE 94
82 1.0 Fig. 7. 2fÂ— Radial distribution function for the LJ potential: T* 2.7^, n* = 5/6.
PAGE 95
83 TABLE 7.2 OCMPASISON OF PTÂ» CHNC, AND m" INTEGRAL EQUATION SOLUTIONS TO MONTE OAHLO FOR THE LJ POTENTIAL: T* 2.?^, n* 5/6* r
PAGE 96
84 2.2
PAGE 97
85 TABLE 7.5 COnPARISON OP PY, CHNC, AND "m" INTEGRAL EQUATION SOLUTIONS TO MONTE CAHLO POH THE LJ POTENTIAL: T* = 2.?^, n* 1.0^ r (Sfic/8m>^ ^WSPY^^ (Smc/^Pj'^^ 0.90
PAGE 98
86 TABLE 7.^ THERMODYITAMIG QUANTITIES FOR THE LJ POTENTIAL n*
PAGE 99
87 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 PY are small. The results for the "m" equation in Table 7.^ 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 EVnkT 2iin / 30(r)g(r)r^dr Jo CO 30(r)g(r)r2dr (7.23) H + 2Tm J (27iikT)^^^p * Ag , (7.2^) Here the first term on the right is the value obtained on the computer. The correction, Ag 2Tm I 00(r)g(r)r''dr , (7.25) may be approximated by setting g equal to one everywhere In the region. If R is sufficiently large this approximation will be quite good. For the LJ potential, then, we easily obtain
PAGE 100
^B
PAGE 101
89 P0(r) l/(Â©r) (7.51) where r is in units of a . Note that Â© is a function, of the density as well as of the temperature. Equation (6.15) has been solved for the Coulomb potential for the two cases Â© O,'*and Â© 1.0. The computed values of m for these cases were 0.566 and O.25I, respectively. It was pointed out in the last chapter that due to the divergent nature of integrals containing Coulomb f functions, the corresponding DebyeHtlckel, or exponentially shielded Coulomb, potentials were used to compute m . The computed results are given in Figures 7.^ and 7.5Â» where comparison is made to solutions of the PY said OHNC equations obtained by Carley [50,51] and to Monte Carlo calculations by the same author [26], A clear improvement over the PT and CHNC results is apparent in both cases by use of the selfadjusting mapproximation of the last chapter. The results of the Â© 1.0 case axe especially good, the computed g diverging from the Monte Carlo average only in the region where g approaches one. The Hard Sphere Potential As a physical model, a fluid of hard spheres is adequate only for describing systems at very high temperatures. Nonetheless, the simple analytical form of the hard
PAGE 102
90 1 1 I'
PAGE 104
92 sphere potential lias invited ei^austive studies of this model, in that otherwise intractable expressions take on manageable forms. For hard spheres of diameter a , the potential is given by I oo Â» r < a 30(r) { (7.32) I , r > a and is temperature independent. We will use the hard sphere diameter a as the natural unit of length, so that the reduced density is n* na^ . Prom Bq. (7.32) it is easy to see that the Mayer f function is here a step function, viz., 1 , r < 1 f(r) \ (7.33) , r >1 where r is in units of a . The derivative of f is thus the Dirac delta function, f'(r) 6(rl) (7.5^) and 80, from Eq. (C.7) we get for the (virial) equation of state j 1 ^ ^ n*g(l) . (7.35) We have solved Eq. (6.15) for the hard sphere model for nine densities, from n* 0.1 to 0.9 , and have obtained the pressure in each case using Eq. (7.55). The computed g's and H's are given in Tables D.5 to D.7. The
PAGE 105
93 computed equation of state is shown in Pig. 7.6, along with the PT and CHNC equations of state obtained "by Klein [24, 25], also using Bq. (7.55). The reference isotherm was put together by Klein by taking the equation of state obtained from the five known virial coefficients for n* less than 0.5 and the molecular dynamics calculations for n* greater than 0.7. For 0.5 ^ n* ^ 0.7, the reference isotherm is a graphical extrapolation of the former onto the latter. Numerical solutions of the integral equations for g for the hard sphere potential are extremely sensitive to the size of the interval, as has been pointed out by Sahlin C56]. It is found in practice that g at the point 1 will decrease as the interval size is decreased. Thus the solution shown in Pig. 7.6, computed for an interval size of 0.025 (in reduced units), probably lies slightly above what would be the correct solution for an infinitesimal interval size. Nonetheless, it is possible to conclude, on the basis of the available numbers, that the Â•'m" equation yields improved solutions over the PY and CHNC equations for the hajxi sphere potential.
PAGE 106
9^ 0.^ 0.8 ^, 1.2 Fig. 7. 6. The equation of state of a fluid of hard spheres.
PAGE 107
APPENDICES
PAGE 108
APPENDIX A FOURIER TRANSFORMS There is a certain amount of ambiguity in the definition of Fourier transforms in that factors of 2n may be distributed in several ways. This is harmless so long as a given definition is applied consistently. We have used throughout this dissertation the following definitions: F(r) P(k) (2ti) ^ rP(k)e^^*^d^k (A.l) jF(r)e"^^*Vr . (A. 2) If F is a function of the magnitude of r alone, then F is a function only of the magnitude of k , and the integrations over angles above may be carried out, giving oo p(r) . Â— i/ F(k) k sin krdk (A. 3) P(1j:) Â» ^ / F(r) r sin krdr . (A.A) i The Fourier transform of the Coulomb l/r potential may be found most easily by means of a wellknown trick. Write 96
PAGE 109
97 l/r 11m e'^/r . (A. 5) a Â» o Then the transform of l/r may be written lim ^ / e'^ sin krdr lim ^ ^^ (A. 6) aooj aoo a+k ^hA^ . (A. 7) The Fourier transform of 30(r) for the Coulomb potential is thus 30(k) ^n3e^/k^ , (A. 8) as stated in Chapter II. For numerical work the infinite sum over infinitesimal intervals, represented by the symbol / dx , must be replaced by a finite sum over finite intervals. So long as only the value of an integral is desired, the machine calculation may be performed by the best convenient rule, such as Weddle's rule [^73. For Fourier integrals, however, more care must be taken. The basic property of the Fourier transform that makes it eminently useful is the fact that it is an expansion in orthogonal functions, making possible the reciprocal relation of a function to its transform indicated by Eqs. (A. 5) and (A.^). If, in a given numerical computation, this reciprocal property is to be used, the
PAGE 110
98 discrete representation, of the expansions must maintain orthogonality. We will now consider how this is done. The following technique is due to Professor A. A. Broyles. Choose a range (0,r ) for the r variable and ^ ' max (0,k ) for the k variable. Placing N equallyspaced points in the range of x divides it into N intervals of equal size Ar Â» r /N . A continuous function F(r) is then represented by the discrete function F(nAr) , n = 0,N Similarly, the range of k is divided into M equal intervals dJs. o ^max^^ * Equations (A. 5) and (A.^) then become, respectively, by the trapezoidal rule [A7] , 2 " P(nAr) ^^^ Â— y F(mAk)m sin (mnArAk) (A. 9) 2 ^ F'(niAk) ^^^^^ ^ F(nAr)n sin (mnArAk) . (A. 10) nÂ»l (The m Â» and n Â» terms are not written in the sums since they are identically zero.) Thus we have M p(n^) . 2^Ar ^ F(n'Ar)n' Y sin (mn' ArAk n' Â»1 m^ X sin (mnArAk) . (A. 11) Consider now the second sum in Eq. (A. 11),
PAGE 111
J Â• > sin (mn'ArAk) sin (uinArAk) m1 99 (A. 12) M i V cos[(n' n)mArAk] cos[(n' +n)mArAk] m=l With the identity [52] M 2 ^ cos mÂ© Â» 1 + mÂ«l we may write sin (M + ^)Â© T f sin I Â© (A. 15) (A. 14) J sin I (2M+l)(n'n)ArAk sin H (n'n)ArAk sin I (2M+l)(n'+n)ArAk sin ^ (n' +n)ArAk (A.15) We now want to impose orthogonality. Write J c5Â„Â„' (A. 16) nn where U aa ^^ "^^^ Kronecker delta, 5u(A. 17) 1 , i " J , i / d and C is a constant to be determined. Then for n / n% we require
PAGE 112
100 sin i (2M+l)ArAk(n'ii) sin i (2M+l)ArAk(n' +a) 2 . ^~^r ^ .(A. 18) sin 4 ArAk(n'n) sin ^ ArAk(n' +n) This is satisfied for arbitrary integers n and n' , n ^ n' , if (2M+l)ArAk 2n (A. 19) since both sides then vanish identically. For n' = n, the second term in Eq. (A. 15) vanishes also, using the result Bq. (A. 19). The first term gives n .n' sin ^ li^ si^ ; i^^'A , lim n cos n (^O (A.20) 1' sin ^^^ n^n' ^ cos . i^ 2M + 1 . (A.21) Thus, and Eq. (A. 11) becomes nn^) . ^^ ^ ^ n Fin ^)6^. (A.25) n' 1 P(nAr) . (A. 2^) If we had started from P instead of F , we would have found that we needed the condition (2N+l)ArAk 2n ' (A.25)
PAGE 113
101 for orthogonality. Since the intervals are the same, this implies ve must have N Â» M . Summing up, then, we find that the reciprocal set of Fourier transforms for continuous functions, Bqs. (A. 3) ancL (A.^), is adequately represented in discrete form by the following set: 2 ^ P(n) Â— A Â— y m P (m) sin (nmA) (A. 26) m>l N ^ ^"'^ ^ ^ n P (n) sin (nmA) (A.27) mal where and A 27i/(2N+l) (A. 28) D = (Ar)^ . (A. 29) Note that once the range and interval size of the r variable are chosen, the orthogonality condition fixes the range and interval size of the k variable. Furthermore, the development above can only be carried out for the trapezoidal rule of numerical integration, which gives equal weights to all the points [^7].
PAGE 114
APPENDIX B PARAMETERS IN THE LEBOWITZ SOLUTION FOR A MIXTURE OP HARD SPHERES Lebowitz's exact solution of the PY equation for a mixture of hard spheres [^0] gives the three direct correlation functions of the system, Eq.s. (5.7) and (5.8), as a function of r and seven parameters, called here a,b,a,b,b,d, and \ . These parameter are given below. Let n and n be the number densities of the + positive and negative hard sphere ions, respectively, and R and R the corresponding hard sphere diameters. Ve will assume that R ^ R . We then define the following Â— + terms: and w Ttn /6 (B.l) + + \ (R_ R^)/2 (B.2) R^_(R_ + R_^)/2 (B.3) zwR^+wR^ . (B.^) + + 102
PAGE 115
103 Lebowitz found that the pressure p of the mixture la given, in the PY approximation, by 0P [ ^% * ^)(1 + z + z'^) iw w (R R )^[R +R (RR(wR^ '1/ ^. w_R 2)] I / (1 _ 2)3 (B.5) The hard core radial distribution functions, labelled here with a "shortrange" superscript, were found to have the following values at the hard core diameters: 6^^(R^) [1 + ^ z ^I w_R_2(R^ . R_)] /(I 2)2 (B.6) gffCRj [1 + I z + I w^R^2(g^ . R^)] / (1 _ 2)2 (B,7) s!!(V) fRs+H.(R+) + VÂ— ^^^^/^V Â• (B.8) With these defined quantities, we complete the set of required parameters by listing the following quantities 1 + z + z^ + g \^(^+ * n_)(l + 2z) 3w (R R )2[R + R + R R (2w R 2 Â•Â• w R ^)] / (1 z)5 ^. ^ R^^0p/(1 z) (B.9)
PAGE 116
10* b^ 6Cw R 2g 2^^j ^ ^ ^ 2 2^ j^ ^3^^Q^ b 6[w R g_(R ) + w R g (R )] (B.ll) d [w^a_^ + w_a_] . (B.12) Tbe constants a and b axe obtained from a and b by interchanging w ,R with w_,R_ .
PAGE 117
APPENDIX C THE VIRIAL AND COMPRESSIBILITY EQUATIONS OP STATE The equation of state of a fluid composed of molecules which interact through spherically symmetric pair forces may be obtained in two ways, one proceeding from the virial theorem, the other based on the fluctuations in density. Derivations of both expressions may be found in, e.g., Hirschfelder, Curtiss, and Bird [53]. They are (a) from the virial theorem, ^ . , . 2^ f 3 dlrl g(^)^3^^ (,^,) and (b), from the compressibility, kT(dp/^n)^^ 1 + ^ Tin / [s(r)l]r^dr . (0.2) Equation (C.2) may be put in a more convenient form by using the fact that the Fourier transforms of G and T are related by the expression [see Chap. I] 1 + nG(k) 1/[1 nT(k)] , (C.3) so that kTCdp/^n)"'1 + nG(0) (C.^) 1/[1 nT(0)] (C.5) 105
PAGE 118
106 and ^ ^ 1 ^ Ttn / T(r)r^(ir . (0.6) J A Equation (C.l) may also be conveniently rewritten as follows, ^ = 1.2^/ g(r)eP^^^)f'(r)r5dr (C.7) where f ' is the derivative of the Mayer f function. It can, of course, be shown [5^] that if the exact diagram expansions of T and ge''^ are inserted in the equations (0.6) and (0.7) Â» both give the correct virial expansion. When approximations to the bridge set are introduced, however, the two equations give inconsistent answers. This situation could be exploited to determine the parameter m of Chapter V by imposing consistency between Eqs. (0.6) and (0.7) even with the approximated B set. That is, we write / T(r)r^dr ' \\ g(r)e^'^^^^f '(r)r^dr ^ n ^ rg(r)e^^^^)f'(r)r5dr . (0.8) Using the diagram expansions of Eqs. (I.30) and (1.29), we find that this leads to the condition
PAGE 119
L 107 oo [p(r) + B(r) + f(r)[P(r)+B(r)] ^  f '(r)CP(r)+B(r)] ^[P(r)H.B(r)]l Ar / [f(r)S(r) > f'(r) n ^ + I f (r)S(r) + I f '(r) n ^isilr^dr . (C.9) If B is now approximated by mP , we get for m the equation oo (1+m) / [P(r) + f(r)P(r) +  f (r)P(r) oo * I ^'(^^ ^T^^""^^ " " / t^^^^^S^^) * f f'(r)S(r) ^ I f'(r) n ^^Ir^dr . (CIO) The sets P and 3 may then be related to g obtained from, say, the PY equation and the integrals evaluated to give m .
PAGE 120
APPENDIX D THE CMIPUTBD FUNCTIONS g AND H FOR THB LJ, COULOMB, AND HARD SPHERE POTENTIALS In this appendix we present the results of computations for the radial distribution function g and the function H defined in Bq. C?.'*) from the integral equation presented in Chapter VI, Eq,. (6.13). Table D.l lists the three LJ cases studied. The Coulomb solutions are given in Table D.2, and the hsurd sphere solutions take up the remaining tables, D.3 to D.7. 108
PAGE 121
109 TABLE D.l THE COMPUTED FUNCTIONS g AND H FOR THE LJ POTENTIAL " ""
PAGE 122
110 Table D.l (continued)
PAGE 123
Ill TABLE D.2 THE COMPUTED FUNCTIONS g AND H FOR THE COULOMB POTENTIAL
PAGE 124
112 TABLE D.3 THE COMPUTED FUNCTIONS g AND H FOR THE HARD SPHERE POTENTIAL
PAGE 125
115 TABLE D.^ THE COMPUTED FUNCTIONS g AND H FOR THE HARD SPHERE POTENTIAL
PAGE 126
11* TABLE D.5 THE COMPUTED FDUCTIONS g AND H FOR THE HARD SPHERE POTENTIAL
PAGE 127
115 TABLE D.6 THE COMPUTED FUNCTIONS g AND H FOR THE HARD SPHERE POTENTIAL
PAGE 128
116 TABLE D.7 THE COMPUTED FUNCTIONS g AND H FOR THE HARD SPHERE POTENTIAL n* 0.9 r/a g H ^__ 0.2 0.5 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.5 2.4 2.5 2.6 2.8 3.0 5.2 5.^ 5.6 5.8 4.0 0.000
PAGE 129
LIST OF REPBRENCBS 1. T. L. Hill, Statistical Mechanics , McGrawHill, New Tork, 1956. 2. N. S. Gingrich, Rev. Mod. Phys. r^, 90 (19^3). 3. H. S. Green, The Molecular Theory of Fluids , NorthHolland, Amsterdam (1951). 4. J. de Boer and A. Michels, Physica 6, ^09 (1939). 5. J. E. Mayer and E. Montroll, J. Ghem. Phys. 2i 2 (19^1) 6. J. de Boer, Physica 15, 680 (19^9). 7. E. E. Salpeter, Ann. Phys. ^, 183 (1958). 8. A different approach to the derivation of the cluster expansion which is. in many respects, more lucid is presented in the first chapter of R. Brout and P. Carruthers, Lectures on the Many Bl ec tron Problem , Interscience, New TorF7 1963. 9. Here and in the following we have adopted, with some minor variations, the notation of M. Klein and M. S. Green, J. Chem. Phys. ^, 136? (1963). 10. L. 3. Ornstein and P. Zemicke, Proc. Acad. Sci. Amsterdam 12Â» 793 (191^). 11. J. M. J. Van Leeuwen, J. Groeneveld, and J. de Boer, Physica 2^, 792 (1958). 12. M. S. Green, J. Chem. Phys. 2^, 1403 (I960). 13. E. Meeron, J. Math. Phys. 1, 192 (I960). 14. L. Verlet, Nuovo Cimento 18, 77 (I960). 15. G. S. Rushbrooke, Physica 26, 259 (I960). 16. T. Morita and K. Hiroike, Prog. Theor. Phys. 23, 1003 (I960). ^ 117
PAGE 130
118 17. J. K. Percus and G. J. Tevick, Pliys. Rev. 110, 1 (1958), 18. J. K. Percus, Phys, Rev. Letters 8, ^2 (1962). 19. G. Stall, Physica 22, 517 (1963). 20. A. A. Khan, Phys. Rev. 1^, A367 (1964). 21. M. Born and H. S. Green, A General Kinetic Theory of Liquids , Camhridge , London, 19^^. 22. J. Tvon, Actualltj^s scientifiiues et industriellee , Hermann , Paris, 1935. 23. A. A. Broyles, S. U. Ohung, and H. L. Sahlin, J. Chem, Phys. YL^ 2462 (1962). 24. M. Klein, J. Chem. Phys. ^, 1388 (1963). 25. M. Klein, Phys. Fluids 2. 591 (1964). No computations for the BGT equation are reported in this paper. However, it is shown that the PY equation leads to better results than the CHNC equation, which in turn was found to be superior to the BGT equation in the previous reference. 26. D. D. Carley, "Recent Studies of the Classical Electron Gas," (to be published). 27. The same problem arises from the LennardJones (612) potential. W. W. Wood and P. R. Parker, J. Chem. Phys. 22, 720 (1957). 28. J. D. van der Waals, Dissertation, Leiden, 1873. 29. L. 8. Ornstein, Dissertation, Leiden, 1908. 50. M. Kac, G. B. Uhlenbeck, and P. C. Hemmer, J. Math. Phys. 4, 216 (1963). 51. G. B. Uhlenbeck, P. C. Hemmer, and M. Kac, J. Math. Phys. 4, 229 (1963). 32. P. C. Hemmer, M. Kac, and G. E. Uhlenbeck, J. Math. Phys. 2, 60 (1964). 33. P. Debye and E. Hflckel, Physik. Z. 24, 185 (1923). 34. D. Bierens de Haan, Nouvelles Tables , Hafner, New York, 1957. Bq. (4), p. 22T:
PAGE 131
119 35Â» Bqs. (3.18) eind (3.19) were originally obtained by a collective coordinate method. The derivation presented, which has the advantage of greater clarity, is based on a suggestion by A, A. Khan, who proposed the viewpoint embodied in Eq. (5.1'<). 36. H. L, Sahlin, Dissertation, University of Florida, 1963. 37. A. A. Broyles, H. L, Sahlin, and D. D. Carley, Phys. Rev. Letters 10, 319 (1963). 38. P. C. Hemmer, J. Math. Phys. ^, 75 (1964). 39. E. Helfand and R. Kornegay, "Study of the Pair Correlation Function Diagrams," Physica (to be published). 40. J. L. Lebowitz, Phys. Rev. 1^, A895 (1964). 41. M. Wertheim, Phys. Rev. Letters 10, 321 (1965). 42. M. Wertheim, J. Math. Phys. ^, 643 (1964). 43. J. N. Shaw, Dissertation, Duke University, 1963. 44. D. D. Carley, "Radial Distribution Function for the Gaussian Model," (to be published). 45. P. Hutchinson and G. S. Rushbrooke, Physica 29, 675 (1963). "~ 46. J. A. Barker and J. J. Monaghan, J. Chem. Phys. 36, 2564 (1962). 47. H. Margenau and G. M. Murphy, The Mathematics of Physics and Chemistry , Van Nostrand, New Tork,"~r961. Chapter 13. 48. A. A. Broyles, J. Chem. Phys. ^, 456 (I960). 49. A. A. Broyles, J. Chem. Phys. ^, 493 (1963). 50. D. D. Carley, Dissertation, University of Florida, 1963. 51. D. D. Carley, Phys. Rev. 1^1 , 1406 (1963).
PAGE 132
120 52. L. B. V. Jolley, Summation of Series , Dover, New York, 1961. Equation (^18;, p. 7^ 55. J. 0. Hirschfelder, 0. P. Curtiss, and R. B. Bird, Molecular Theory of Gases and Liquids , Wiley, New Tork, 1954. See pages 13^ and 1277 54. G. S. Rushbrooke and H. I. Scoins, Proc. Roy. Soc. A216, 205 (1953).
PAGE 133
BIOGRAPHICAL SKETCH Fred Lado was 'born June 5Â» 1938, in La Coruna, Spain. In June, 1956 Â» bÂ© graduated from Jesuit High School in Tampa, Florida. He attended the University of Florida and received the degree of Bachelor of Science with High Honors in February, I960. In the same month he was admitted to the Graduate School of the University of Florida. From September, I960, to February, 1961, he attended Syracuse University with a Woodrow Wilson Fellowship. He returned to the University of Florida in February, 1961. While pursuing his graduate studies at this institution, he has held fellowships from the National Science Foundation and has worked as a teaching assistant and a research assistant. Fred Lado is meirried to the former Maria Dolores Pastor and is the father of one child. 121
PAGE 134
This dissertation was prepared under the direction of the chairman of the candidate's supervisory committee and has been approved by all members of that committee. It was submitted to the Dean of the College of Arts and Sciences and to the Graduate Council, and was approved as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August 8, 1964 Dean ts and Sciences Dean, Graduate School Supervisory Committee: Chairman ^ ^.H /â€¢ciM.^i
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EC91YUHNJ_VW33LR INGEST_TIME 20170719T20:35:49Z PACKAGE UF00097932_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES

