Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UF00097932/00001
## Material Information- Title:
- New methods of computing radial distribution functions of fluids
- Added title page title:
- Radial distribution functions of fluids
- Creator:
- Lado, Fred F., 1938-
- Publication Date:
- 1964
- Copyright Date:
- 1964
- Language:
- English
- Physical Description:
- xi, 121 leaves. : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Approximation ( jstor )
Coulomb potential ( jstor ) Diameters ( jstor ) Differential equations ( jstor ) Distribution functions ( jstor ) Equations of state ( jstor ) Fourier transformations ( jstor ) Radial distribution function ( jstor ) Sine function ( jstor ) Thermodynamics ( jstor ) Dissertations, Academic -- Physics -- UF Kinetic theory of liquids ( lcsh ) Physics thesis Ph. D Statistical mechanics ( lcsh ) - Genre:
- bibliography ( marcgt )
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:
- 000543389 ( AlephBibNum )
13143435 ( OCLC ) ACW7097 ( NOTIS )
## UFDC Membership |

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 |