Citation
New methods of computing radial distribution functions of fluids

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:
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 )
non-fiction ( marcgt )

Notes

Thesis:
Thesis -- University of Florida.
Bibliography:
Bibliography: leaves 117-120.
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 non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
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 two-component systems of the expressions derived

in Chapter II and their numerical solution for a model of

an electrolyte.







In principle, the radial distribution function for

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

technique, given the pair potential of the fluid. In

practice, however, this technique is a voracious consumer

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

producing "touchstones" against which other, more frugal

approximate solutions may be tested. These approximation

methods usually result in integral equations for g which

have the added advantages over MC that they provide a means

of computing experimental potentials from scattering data

and may indeed have analytic solutions, as witnessed by

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

hard spheres.

The diagram development of Chapter I suggests a

possible avenue of approach in the improvement of the

integral equations for g derived there. This matter is

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

upshot of the arguments presented here is a new integral

equation for g which contains a parameter, m, determined

by the peculiar potential and temperature of the system

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

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

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

tions of the equation have been obtained numerically for

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








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

as to the "standard" solutions available. These numerical

results constitute Chapter VII.

The numbering of equations throughout the disserta-

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

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

same rule applies to tables and figures.












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 TWO-COMPONENT 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 Lennard-Jones 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

N-particle 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 two-body 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 3-particle, 4-particle, etc., distribution functions
can be similarly defined. In general, the n-particle
distribution function is given by

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

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


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

All thermodynamic quantities can be expressed in terms of

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









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

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

Here. C(r) is defined by


C(r) nt f dr ... drt+2
'77707 fij d3 "' d rt+2
t-1 J
(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 non-nodal diagrams in the
expansion of G(r), Eq. (1.28), by T(r); i.e.,

T(r) P(r) + B(r) + f(r)[l + S(r) + P(r) + B(r)]
(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 1-2 bond. Diagrams
of this type are put in explicitly by the definitions, as
in the last term of Eq. (1.28). Thus also the non-nodal set
T defined by Eq. (1.30) includes not only the sets P and

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

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

point. Let the first nodal point encountered along a path

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

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

be some member of the set T The subdiagram between 3

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

restriction. Hence in general it will be some member of

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

members of T and between 3 and 2 all members of G That
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 non-equal subdiagrams are counted twice in
the multiplication of the first two bracketed terms and the









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


(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
t-2


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 [11-16]

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

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

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

B(r) P(r) (1.47)
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 Percus-Yevick
(PY) integral equation [17-20]


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

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

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

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

dissertation.

For a realistic potential made up of a repulsive core

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












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 Waals-like equation of state is
based upon this type of potential [28,29]. A one-dimensional
model of this type has been studied by Kac, Uhlenbeck, and
Hemmer in a succession of papers [30-32] and was found to
display a phase transition even in lowest order.
The classic example of a long range potential is
the Coulomb potential which is, in principle, of infinite
range. An explicit approximation to the radial distribution
function for a system of particles interacting with Coulomb
forces was found by Debye and Huckel [33] and reads

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

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 non-iterative
solution for g by combining Eq. (2.3) with

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

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

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


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

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

P0(k) 4 ne2/2 (2.9)


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









The two essential approximations of the DH
expression are thus

B(r) 0 (2.10)

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

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

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

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)expC-00r(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 + nG-Sr(k)2AT(k) AT(k). (2.17)
1 n[l+nGsr(k)]AT(k)

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

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

with


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









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

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) + e-BSr(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) + 0----o + n --o + o' %


+ 6 o + ** (2.24)

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


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


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

and further approximating by neglecting the second part and
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) e-Bsrflr(r) (2.29)

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


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


gSr(r)flr(r) (2.31)

It will turn out, however, that these higher
approximations, Eqs. (2.29) and (2.31), are not as satis-
factory as the original approximation, Eq. (2.28). This
will be discussed further in the next chapter.
Consider now the logarithm of both sides of Eq,
(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 + n--r(k) (2.36)
1 + n[1 + nGsr(k)] r(k)



(k) 1 1 + nGr(k) 1
G(k) n + n[ + nGsr(k)101r(k) l
(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 PY-type approxi-
mation, where AP(r) + AB(r) is neglected, is easily
obtained. From Eq. (1.29), we write


gsr(r)e08(r)e ) 1 + Sr(r) + Psr(r) + BSr(r)
(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) + e-0(r)AS(r) (2.41)

where we have used the approximation

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

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


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

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










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

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 CHNC-type approximation,

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

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

In the next chapter these equations are tested on a

model for which near exact answers are available for

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 so-called Gaussian

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

tive Gaussian function. This model was used to carry out

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

as of the equations resulting from the higher approximations

to AT The BSC equation was applied to the Gaussian model

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

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

results reported.

The Gaussian model has been recently investigated by

Helfand and Kornegay [39]. These authors generated and


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









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

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

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 e-1 21 x (3.7)

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


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

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

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


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

The long range potential is given by the equation


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


(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)e-ikr 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 J-0
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) j-0


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


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

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









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

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

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.1-3.4. We see that, as
expected, g computed from the BSC equation compares
favorably with the exact g in the large x limit, but
becomes far too small as it approaches the origin. Hemmer's










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

no means of interpolation for intermediate values of x .

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

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

diagrams has been neglected in its derivation. As with the

PY integral equation itself, this success is explained on

the basis of cancellations between the sets P and B

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

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 PY-type approxi-

mation and the above-mentioned approximations for AT are

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

that at small densities the higher approximations do indeed

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

far worse answers at larger densities.















0 .4



a 0

M ca 0 e 'o
60 bo wM M

I I I
Id I "


0
I




0

Sr


r-l












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

m-0.8


1.5
.5 m=0.6





1.0 m=0.6


m=0.7


0.5 m-0.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.1-3.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 TWO-COMPONENT FLUID

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


U(~9.....) ll(rij) + 022( )
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 one-component fluid, viz.,

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

but now the expressions for the series diagrams become


S B(r) .- Z n Ta,(r')G(15-'1)dr' (4.6)
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)
m-2
[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 two-component system.
For the short range potentials alone, an equation
identical to Eq. (4.8) may be written, with all quantities
labelled with an "sr" superscript. Then the ratio of the
two equations gives, as in Chapter I,


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


(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

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


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


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








For the inverse solution we obtain

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


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


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) [l-n2T22(k)] /Dt(k) (4.23)

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

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

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

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

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

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









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


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

From Eq. (4.25), we write


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

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


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


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


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


T12(k)] (4.53)

We divide top and bottom of the right hand side of Eq.
(4.32) by Dsr(k) where all the sets of Eq. (4.29) acquire
an "sr" label, and use Eqs. (4.23), (4.24), and (4.27) to
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- 2-sr2
l+nl%11(k)) ATll(k) + n2 G12 (k)AT22(k)


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


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


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


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

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

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


Tll(k) (4.40)

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)e-22(r) (4.45)
22(r) g2222


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

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


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

(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 PY-type approximation may be applied as in
Chapter II, giving









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


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


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


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



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





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


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

(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 Percus-Yevick integral equation.

Lebowitz [40] has recently solved this equation analytically

for binary mixtures of hard spheres, generalizing a method

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

direct correlation functions of the hard sphere mixture

are given by the following expressions:

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


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

where

x r-X.

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

Appendix B. We have assumed, for definiteness, that









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

R > R can be easily found by interchanging all + and -
+
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 non-iterative scheme for

computing the radial distribution functions describing an

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

sents such a physical system. A Monte Carlo solution of

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

have computed, for purposes of comparison, the ionic radial

distribution functions for the same input conditions as the

Monte Carlo solution, although it is pointed out by Shaw

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

rations generated for the Monte Carlo computation "yields

radial distribution functions still showing fluctuations

large enough to preclude strictly quantitative comparison

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

afford more than a good qualitative test, due to insufficient

computer runs.










Shaw considered an electrically neutral collection
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 10-6 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 cut-off

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

being nearly one everywhere outside the hard core diameters.









The results of the calculation are shown in Figure

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

using a somewhat modified Debye-Hickel expression obtained

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

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

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

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

The comparison with the Monte Carlo solution is, as

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

present solution is at least qualitatively correct. The

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

approximation, as discussed previously, both gave essentially

the same answers and a single curve has been drawn to


represent both.









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 two-body forces alone,

may be obtained with the aid of the Ornstein-Zernicke

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

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

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

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

ability of the diagram technique, and a closed set of

equations for g was not achieved.

This impasse has necessitated the introduction of

approximate integral equations, foremost among which are

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

viewed as two different approximations to the set B ,

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 Lennard-Jones (6-12) potential at rela-
tively high temperatures [23]
3) the CHNC equation is probably superior to the PY
equation for the Lennard-Jones potential at
relatively low temperatures [20]
4) the CHNC equation is superior to the PY equation
for the Gaussian model [44].

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

P2 and B2; i.e.,


P2(r) <>o (6.3)


B2(r) (6.4)

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









The difficulty with the set B is already evident

in B2 The presence of the cross bond connecting the

two field points prevents the factorization of the diagram

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

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

CHNC equations circumvent this difficulty by, in effect,

substituting a constant for the cross bond in B2 equal

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

the Mayer f bonds for the various cases mentioned above.

These are shown in Figure 6.1.

The f bond for the hard sphere potential is precisely

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

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

mation is not ruined, because the presence of other bonds

connecting the two field points in B2 tends to cut off

the integrals anyway, taking over the task of the real

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

poorer one, since it approximates the cross bond by zero

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 Jif3f32flf-2d3rld-r2d3r3d-r (6.12)


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







69

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

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

(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
self-consistent. That is, when the new H calculated is
identical to the H used as input. A measure of the degree









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

difference between successive H's where the rms is

defined as follows:


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

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

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

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

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

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

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

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

at all points.

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

about 10 iterations. For larger densities, the slow rate

of convergence may be improved by the following device, due

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

consistency, it is found to approach its final value

exponentially. Then from any three successive iterations,

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

have after an infinite number of iterations, viz.,










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

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

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

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

The Lennard-Jones Potential

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

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

where the constants E having the dimensions of energy, and










a with the dimensions of length, are to be determined

from experiment. Expressions involving the LJ potential

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

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

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

where the dimensionlesss) reduced temperature T* is

defined by

T* kT/E (7.19)

The corresponding dimensionless density is n* na3 .

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

pressure, and the energy for a fluid interacting with a

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

to the PY and CHNC equations at these same densities and

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

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

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

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

PT equation also yields good solutions for these cases,

while the CHNC equation is in poorer agreement.










All four plotted curves would lie confusingly close

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

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

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

meaningful. More substantial deviations appear in Table 7.2

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

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

closest agreement in general, followed in turn by the PY

and then the CHNO equations.

Once the radial distribution function is known the

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

p/nkT = 1 (2zn/3) f 0(r)g(r)r3dr (7.20)
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.2-1.5 4.08 6.97
(p/nkT)m 1.24 4.09 6.99
(p/nkT)pyb 1.24 4.01 6.8
(p/nkT)CHNCb 1.28 5.11 9.1

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

aFrom Reference 27.
bFrom Reference 23.


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)
3TE-R

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 two-component 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 Percus-Yevick (PY) equation for hard spheres. The diagram development of Chapter I suggests a possible avenue of approach in the improvement of the integral equations for g derived there. This matter is taken up again in Chapter VI, which begins Part Two. The upshot of the arguments presented here is a new integral equation for g which contains a parameter, m, determined by the peculiar potential and temperature of the system studied. This equation reduces to the well-known PercusYevick (PY) or convolutlon-hypemetted 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 TWO-COMPONENT 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 Lennard-Jones 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. . . . A-2 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 N-particle 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 two-body 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^ E-|NkT + |/.../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 N-2 psirticles. It has the additional property that its Fourier transform may be directly measured by x-ray 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 well-known 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,^^ t-1 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 • (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 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 non-nodal diagrams in the expansion of G(r), Eq. (1.28), by T(r) ; i.e., T(r) P(r) + B(r) ^ f(r)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 1-2 bond. Diagrams of this type are put in explicitly by the definitions, as in the last term of Eq. (1.28). Thus also the non-nodal set T defined by Bq. (1,30) includes not only the sets P and B but also all diagrams with a direct 1-2 bond. Consider a typical diagram of the set S . By definition it has at least one field point which is a nodal point. Let the first nodal point encountered along a path from 1 to 2 be labelled 3. Then the subdiagram between 1 and 3 will have no nodes, by construction, and hence will be some member of the set T * The subdiagram between 3 and 2 may or may not have further nodes. There is no restriction. Hence in general it will be some member of the set G . Thus we can generate all possible members of the set S by substituting between 1 and 3 all possible members of T and between 3 and 2 all members of G , That 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 non-equal 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 [11-16] 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 Percus-Yevick (PY) integral equation [17-20] 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 Bom-GreenYvon (BGY) [21,22], based on the superposition approximation [1]. It has been found to be less reliable than the PY or CHNC integral equations for a Lennard-Jones (6-12) potential by Broyles et al . [23] and for hard spheres by Klein [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 Lennard-Jones (6-12) 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 Vaals-like equation of state is based upon this type of potential [28,29]. A one-dimensional model of this type has been studied by Kac, Uhlenbeck, and Hemmer in a succession of papers C 30-32] and was found to display a phase transition even in lowest order. The classic example of a long range potential is the Coulomb potential which is, in principle, of infinite range. An explicit approximation to the radial distribution function for a system of particles interacting with Coulomb forces was found by Debye and 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 non-iterative 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 2-nrr 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 Debye-Httckel (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) + 0----0 + 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 e-P<^''"(^)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 PY-type 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)e-P<^^''(^) . 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 CHNC-type approximation, AB = , we were led to Eq. (2,18), A PY-type 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 so-called 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^ t-1 '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 t-1 (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 J-0 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)^ j-0 ^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 + /4-n n / G^^(x) x^dx , (3.26) .sr. 5 A(x) e-P^ <-) 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.1-3.^. 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 PY-type approximation gives a better result for g , although an additional infinite set of diagrams has been neglected in its derivation. As with the PY integral equation itself, this success is 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 PY-type approximation and the above-mentioned approximations for AT are shown as functions of the density for m 0.6. It is seen that at small densities the higher approximations do indeed lead to successively more accurate g*s, but at the price of far worse answers at larger densities.

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.1-5.^,

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.4-5) (solid line).

PAGE 53

41 0.02 ,rms 0.01 0.6 u.u<^

PAGE 54

4-2 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 TWO-COMPONENT FLUID The equations of Chapter II may be readily generalized to a two-component 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 one-component 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 n-i \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 g-iT . 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 two-component 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 [l-npTp2(lc)]T.,(k) + npT,p^(k) G-,(k) :r^-^ ^: ^ ^^ ^ M C^.l?) ^^ [l-n3^T^^(k)][l-n2T22(k)3 n3^n2T^2^^) [l-n,T,,(k)]Tpp(k) + n,T,p^(k) Qpp(k) =-i-ii ^1 ^ ^^ ^ '^ (^.18) ^"^ [l-n3_T^3^(k)][l-n2T22(k)] n^n2T^2^^^ T,p(k) 0,p(k) » ^ ^^-:z :r— 2 .(^.19) ^'^ [l-n3^T^3^(k)][l-n2T22(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) . Cl-n2T22(k)] / D^(k) (^.23) 1 + n2G22(k) [l-n^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) [l-n3_T3_3^(k)][l-n2T22(^)^ 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) » [l-n2T|^(k) n2AT22(k)] /D^(k) (^.32) where 0+;^^'^ ^^ written in an expanded form by replacing D^(k) Cl-n3^T^J^(k)][l-.n2T|^(k)] n^n2T^|^(k) n3^AT^3^(k)[l-n2T||(k)] n2AT22(k) [l-n^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+ii2G|2(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 two-component fluids which would be of interest here are made up of charged particles, such as fused salts and electrolytic solutions, we may further restrict the equations above by considering only those long range potentials that satisfy the relations 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 [lH-w,(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) n2G|2(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 Pl-type 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) g|2(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 Percus-Tevick 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), (^.A-8), and (^.51)-(^.53)» vd-th the H's determined from Bqs. (^.46) and (^.^9). This procedure gives a non-iterative 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 cut-off 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 Debye-Htickel 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 "Pl-type" and a "CHNC-type" approximation, as discussed previously, both gave essentially the same answers and a single curve has been drawn to represent both.

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 two-body forces alone, may be obtained with the aid of the Ornstein-Zemicke 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 Lennard-Jones (6-12) potential at relatively high temperatures [25] 5) the CHNC equation is probably superior to the PY equation for the Lennard-Jones 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 Debye-Htickel , 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^^^'''^+ (l4-m)P(r ') ]G(tr-r' | )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 e-Wr)^ (^.^^ f(r)[l + H(r)] (7.8) p'(r) (1+m) [H(r) ln[l + H(r)]] (7-9) 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 self-consistent. That is, when the new H calculated is identical to the H used as input. A measure of the degree

PAGE 88

76 of self-consistency Is given by the root-mean-square (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., self-consistent, 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 self-consistent solution usually being achieved in about 10 iterations. For larger densities, the slow rate of convergence may be improved by the following device, due to Broyles [A-8]. When H is sufficiently close to selfconsistency, it is found to approach its final value exponentially. Then from any three successive iterations, k-1, 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 so-called mixing process, where we write, after the ktt iteration, Hjj.''®*(jA) a Hj^(dA) + (l-a)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 Lennard-Jones 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.1-7.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 self-adjusting m-approximation 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(r-l) (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 well-known 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) a-ooj a-o-o 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 equally-spaced 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 [A-7] , 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) m-1 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) z-wR^+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 "short-range" 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) fR-s+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 d|lrl 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 , McGraw-Hill, 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 Lennard-Jones (6-12) 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 UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EC91YUHNJ_VW33LR INGEST_TIME 2017-07-19T20:35:49Z PACKAGE UF00097932_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES