GENERAL DISCUSSIONDr. J. A. Barker and Dr. D. Henderson * (C.S.I.R.O., Melbourne) said: At present,reliable values for the radial distribution function and the equation of state are avail-able only for the hard-sphere fluid. Thus, a promising method for dealing withsystems of molecules with attractive forces is to consider the attractive potential asa perturbation on the hard-sphere potential. We have recently used perturbationtheory to calculate the equation of state of a fluid of molecules interacting accordingto the square-well potential. This is a particularly good potential because the effectof the attractive forces is not complicated by the “ softness ” of the repulsive partof the potential. In addition, Monte Carlo and molecular dynamics calculationshave provided quasi-experimental data with which our calculations can be comparedwithout the uncertainty due to presence of three-body forces and the lack of know-ledge of the intermolecular potential which is inevitable in applications to realfluids.We divide the intermolecular potential into the “ unperturbed ” hard-spherepotential uO(R) and the “ perturbation ” u,(R).In our calculations, u,(R) = --8for 0 < R < 1-50 and is zero otherwise. If N , is the number of intermolecular distancesin the interval a<R< I.%, then the Helmholtz free energy can be written :where Fo is the free energy of the unperturbed system and the angular brackets mean“ average over the configurations of the unperturbed system ”.The average ( N , ) is given by an integral over the radial distribution function ofthc unperturbed system, go(R).Thus, the first-order term is identical with that givenby Z ~ a n z i g . ~ The second-order term is equivalent to that given by Zwanzig but issimpler and more suggestive. For example, it is intuitively plausible that the devia-tions of N1 from its mean value should be least important at high densities so thatthe perturbation expansion (1) converges best at high densities. Also N1 can beregarded as representing the number of molecules in a spherical shell surrounding acentral molecule. If this shell were a large macroscopic volume, then<m - <w2 = <N,>kT(aP/aP), (2)where (dpldp) is the macroscopic compressibility. An even better approximationis obtained if we use the local compressibility instead of the macroscopic compressi-bility, i.e., if we replace (dp/dp)ogo(R) by (d/dp)[pg(R)],.Our results are based onthis local compressibility approximation. However, similar results are obtained ifthe macroscopic compressibility approximation is used. We are presently makingMonte Carlo calculations of ( N f ) and ( N , ) which will eventually supersede (2).At present we have only a few scattered results. However, these results do indicatethat (2) is a reasonable approximation.* permanent address : University of Waterloo, Canada.’ B. J. Alder, unpublished results.A. Rotenberg, J. Chern. Physics, 1965, 43, 1148.R. W. Zwanzig, J. Chern. Physics, 1954, 22, 1420.5GENERAL DISCUSSlONUsing (2), the free energy is51where q = npo3 /6.In obtaining (3) we have used the Percus-Yevick (PY) compressi-bility isotherm to obtain (aplap),. In evaluating F we used the PY expressionfor go(R). To order T-2, (3) gives the exact second virial coefficient.8-40 0.2 0.4 0 0Vol vFIG. 1 .--Equation of state for the square-well potential. The points are the Monte Carlo values ofRotenberg and the curves are isotherms calculated from eqn. (3). The points given by 0 0 , $ and 0were calculated using 256 molecules at E/kT = 0,0-3, 1 and 2, respectively, whife the points given by0 were calculated using 864 molecules at E/kT = 1.In fig. 1 we have compared the equation of state which results from numericaldifferentiation of (3) with the Monte Carlo calculations of Rotenberg.' The agree-ment is good-even at the lowest temperature, kT/& = 0.5, which is far below thecritical temperature.The agreement with Alder's molecular dynamics calculationsis even better.We have generalized this procedure to include potentials with a soft core suchas the 6 : 12 potential. For an arbitrary potential u(R) we define a modified potentialM. S. Wertheim, Physic. Rev. Letters, 1963,10, 321 ; J. Math. Physics, 1964,5, 643.J. Chem. Physics, 1963, 39, 474.A. Rotenberg, J . Chem. Physics, 1965, 43, 1148.B. J. Alder, unpublished results.E. Thiele52 GENERAL DISCUSSIONby the equations, ( " i d ) , d + - - a , R - dt'(a,y,d,Q; R) = u d+-UR - d Q - i (= O , a<d+-<d+-,U &= yu(R), R > Q. (4)When cc = y = 0, v becomes the hard-sphere potential with diameter d, while foru = y = 1 the original potential u(R) is recovered ; u is an inverse steepness parameter1 I I I I I0.0 0.2 0.4 0.6 0,B 1.0No '1 VFIG.2.-Equation of state for the 6 : 12 potential, kT/c = 2-74. The solid line gives the resultsobtained from eqn. (5) and the broken line gives the results. of McQuarrie and Katz. The pointsgiven by 0 and 0 were calculated by Wood and Parker for 32 and 108 molxules, respectively.for the repulsive region and y a depth parameter for the attractive region. By amethod based on the ideas of Rowlinson and Zwanzig the first derivatives of theconfiguration integral with respect to u and y can be calculated in the limit u, y+O.Thus, we find for the Helmholtz free energyF, -?go( y, ,I)[ JI exp (- u( R)/kT)dR -(a - d ) +z7Jom go ( Nd3 7'2 ' r m R 2 d R kT + higher order terms, ( 5 )1 F -- -NkT - NkTJ.S. Rowlinson, Mol. Physics, 1964, 7, 349; 1964, 8, 107.R. W. Zwanzig, J. Chem. Physics, 1954, 22, 1420GENERAL DISCUSSION 53where Fo and go are the free energy and radial distribution function respectively forthe hard-sphere system. For a = y = 1 we obtain an expression for the free energycorresponding to the potential u(R). We choose CT so that u(a) = 0 and determine dso that the second term in (5) is zero. These choices should minimize the contributionof the higher-order terms and are a generalization of Rowlinson's pr0cedure.l Inaddition, these choices ensure that (5) reduces to (3) if the square-well potential is used.Thus, d is an effective core diameter which may depend on temperature but not ondensity. Pressures calculated by differentiation of (5), using the Percus-Yevickradial distribution function, the 6 : 12 potential, and neglecting the higher-orderterms, are compared with the Monte Carlo results of Wood and Parker for~ T / E = 2.74 in fig.2. The contribution of the term in y 2 calculated by our " localcompressibility " approximation is very small at this temperature. Also shown areresults given by the equation of McQuarrie and K a t ~ . ~ The results given by eqn. (5)appear to be satisfactory.Dr. P. Hutchinson (A.E.R.E., Harwell) said: I report the results of some calculationsin which the Percus-Yevick (P.Y.) and hypernetted chain (H.N.C.) approximationsare tested against the molecular dynamics calculations of Rahman at densitiescorresponding to the liquid phase of argon.From thepair correlation function h(r), computed by Rahman, the direct correlation functionC(r) was calculated. It is then possible from the P.Y. and H.N.C. equations tocalculate the potentials 4p.y. and #H.N.C. which would give such forms of C(r) andh(r) in these approximations. These can then be compared with the potential 4(r>which was the starting point of Rahman's calculation. This is a 6-exp Buckinghampotential truncated at 7.65 A ; the temperature was T = 85.5"K, and the densitywas 0.02103 atm/A3 (equivalent to an argon density of 1.407 g/cm3). The temperatureis thus significantly below the critical temperature for argon.The numerical procedure is to calculate C(r) via a double Fourier transform.Fourier transformation of the Ornstein Zernike equation giveswhereThe procedure is that suggested by Johnson, Hutchinson and March.5a k ) = L(k)/(l +pJi(k)), (1)(14 h{k) = J exp (-ik .r)h(r)d3r.Hence to calculate C(r) we first evaluate h(k) by Fourier transformation. Fromthis we may obtain C(k) from (I), and c ( k ) may be transformed to obtain C(r).Great care is necessary to obtain accurate results from a numerical Fourier transformsmethod. The overall accuracy of the method may be assessed by evaluatingh,(r) = - [exp (ik . r)h"(k)d3k.Formally, h,(r) = lz(r), so that the difference between the two functions is an indicationof the error in the method.Comparison of these two functions showed that thedifferences were about 1 % except at the origin where they may rise to -5 %. How-ever, errors near the origin are unimportant as this region contributes little to the(W3J. S. Rowlinson, Mol. Physics, 1964, 7, 349 ; 1964, 8, 107.W. W. Wood and F. R. Parker, J. Chern. Physics, 1957, 27, 720.D. A. McQuarrie and J. L. Katz, J. Chem. Physics, 1966, 44, 2393.A. Rahman, Physic. Rev. Letters, 1964, 12, 575.M. D. Johnson, P. Hutchinson and N. H. March, Proc. Roy. Soc. A , 1964,282,28354 GENERAL DISCUSSIONtransformation (la). An additional check is available by direct solution of theOrnstein Zernike equation for C(r). Given h(r) the 0-Z equation is simply aninhomogeneous linear integral equation which may be solved by direct inversion.While this last method is subject to greater numerical error it produces the samegeneral result for C(r) as that described below.The results of the Fourier transform method are shown in fig.1 where C(r) isplotted with h(r) and - 4(r)/kT. The values for C(r) are reliable for r <4 A, in thatthe values obtained there are insensitive to the truncation point of h(r) and the choice1.7 2.0 3.0 4.0 5 . 0 6 . 0 7.0 &O 8.5r, (A)- , &I; - - -Y C(4; 0, cp(r)/kT.Fro. i.-Comparison of the direct and total correlation functions and the potential.of mesh for the integral over k. Beyond 6A the results are strongly sensitive to theabove variations and decrease in magnitude when C(k) is smoothed.This indicatesthat the last part of the curve is due to cancellation error (noise) and that C(r) isvanishingly small in this region. The most notable feature of C(r) is that it has thesame general shape as the potential but dies away rather more rapidly. This iscontrary to the expectation of both the P.Y. and H.N.C. approximations thatC(r)- -4(r)/kT for r large. In fig. 2 are shown the two approximate potentials4p.y.(r) and $H.N.C.(r) in comparison with &), given in dimensionless units. The fitobtained by 4H.N.C.(r) is much the better of the two. The well depth of q5p.y.(r) is muchtoo large and displaced too far to the right. The comparison between well depths is+p,ye(4.6)/kT = -4.01, #H.N.C.(4*0)/kT = - 1.28, +(3*8)/kT = - 1-40GENERAL DISCUSSION 55It has been shown by Gaskell that when the compressibility is small the H.N.C.approximation gives values of the pressure which are several orders of magnitude toolarge the virial equation is used. This situation would hold here as calculation giveskT(ap/ap), N 0.05.2.0I .51.0k5BI 0 ' 5nW0.0- 0 5- 1 .01 :0i ' iI Ir, (4FIG. 2.Xomparison of exact and approximate potentials.-? -rg(r)/kT; - - -$ -WXNc(r)/kT; 0, -cpF.Y,(r)/kTThe virial equation is2n O0 P = pkT - ,p2 [ r34'(r)g(r)dr.J oThis minimum value of $H,N.C,(~) is shifted just far enough to bring the large positivevalues of -@(Y) into coincidence with the peak in g(r), and this produces the grossover-estimate of the pressure.The virial equation is thus a very sensitive test ofany theory in the low-pressure, high-density state of a liquid. Finally we may con-clude that despite the poor results obtained from the virial equation the H.N.C.approximation is superior to the P.Y. equation for a liquid well below its criticaltemperature. Grateful thanks are due to Dr. A. Rahman for making his resultsavailable for this calculation.Prof. J. S . Rowlinson (Imperial College, London) said: In his lecture Rushbrookementioned that he and Silbert have been able to show that the inclusion of a tripletT. Gaskell, Proc. Physic. SOC., 1966, 89, 231.G. S. Rushbrooke and M. Silbert, Mol. Physics, 1967, 12, 50556 GENERAL DISCUSSIONpotential in the hyper-netted chain theory is " particularly simple " since it leads onlyto the introduction of a new elementary graph, x1(1,2) in his notation, which is oflower order in the density than those arising from the two-body potential.There isan analogous result for the Percus-Yevick the0ry.l The new graph is here alwaysof the form [e(l,2)xl(l,2)], which leads at once to the physically reasonable resultthat the triplet forces are of importance only outside the repulsive regions of the two-body potential, i.e., outside the regions in which e(1,2) = 0. Moreover, the simplegraphical expansion that is so characteristic of the Percus-Yevick theory is stillpreserved when the triplet potential is included. It can be shown that g(1,2) is stilla sum of all graphs in which the interiorf(i,j) bonds do not cross.The only changeneeded is the inclusion of those graphs in which one or more of the peripheral bondsare notf(i,j), but [e(i,j)xl(i,j)].Dr. S. Levine (Manchester University) (partly communicated) : The Percus-Yevicktheory has been described in terms of " turning on " an external field which enablesone to vary arbitrarily the local number density of molecules in the fluid in the grandcanonical ensemble. This density is regarded as an independent variable and thena suitable dependent variable must be so chosen as to give a rapidly convergentfunctional Taylor expansion. Are there not other and indeed simpler and moremeaningful explanations giving the physical basis of this theory ?Prof. J. S. Rowlinson (Imperial College, London) said: I shall attempt to answerthis question in my summarizing remarks.Mr.A. Moreton and Prof. G. H. A. Cole (University of Hull) (communicated):We too would stress the importance of the Born-Green-Yvon equation as an exactequation which shows a fluid-solid transition. The need to improve upon theKirkwood closure procedure is then of central importance in the theory of liquids.We have recently used the procedure of functional Taylor expansion to express g ( 3 )in terms of gt2) and the pair potential function u(r). The general method of functionalexpansion has been set down by Verlet amongst others.For our purposes we expand the functional p(l)(3/u12) [exp (Pu,(3)) - 11 as a seriesin terms of p(')(4/u1)[exp (Pu1(4)) - 11. Here p ( l ) is the single particle distribution,uI2 is an external potential due to the presence of two extra particles at positions rland r2.The expansion is made about the condition u1 = 0; u12 = u1 +u2 reducesto u12 = u2 in this limit. Taken to the linear term, our expansion gives the Kirkwoodsuperposition statement, which is (1) of the Rice-Young paper but with z = 0.The inclusion of the quadratic term in the expansion gives9'3'(1,2,3) = g'2'(1,2)g'2'(l,3)g'2'(2,3) 1 + p dr4y(1,4)f(1,4) x 1 1Here we useY(r> = g'2'(r) exp p 4 9 ,f ( r ) = exp (- pu(r)) - 1.In the law-density limit the various functions in (1) can be expressed as a simpledensity expansion using an iterative method and (1) is then exact up to the linear termJ. S. Rowlinson, MoZ.Physics, 1967, 12, 513.L. Verlet, Physics, 1966, 32, 304GENERAL DISCUSSION 57in the density. When used in conjunction with the Born-Green-Yvon equation,(1) then provides as an exact expression for the fourth virial coefficient. The tripletterm in the functional Taylor expansion has also been obtained. Although it is toolengthy an expression to set down here we have found that, in the low density limit,it leads to an exact expression for g ( 3 ) up to the quadratic density term, and so givesan exact expression for the fifth virial coefficient.For higher densities than those where g ( 2 ) can be expressed as a simple densityexpansion, the above expansion procedure is still applicable and the expression (1)for g ( 3 ) can be used as a closure expression for the Born-Green-Yvon equation.The solution of this more general situation, for given u(r), could apply to higher fluiddensities. We are at present calculating g(2) by using (1) and the Born-Green-Yvon equation, for a function u of the Lennard-Jones type but with a hard core atvanishing interparticle separation.Our method could be extended to includecorrelations higher than the triplet. In this way it can be considered with the methodof Fisher.'Prof. Stuart A. Rice and Dr. David A. Young (University of Chicago) said: Themediocre quality of the thermodynamic functions which are obtained from thesolutions of the Yvon-Born-Green equation seem to exclude this equation as a usefultool in the study of fluids. However, by modifying the YBG equation we mayobtain useful and interesting results.The YBG equation arises from the introductionof the Kirkwood superposition approximation into the second member of an infinitesequence of integral equations. This latter equation, which involves the three-particle distribution g( 129, is exact, and we may construct increasingly good approxi-mations to this exact equation by multiplying the superposition product by a seriesof triplet correlation terms. By introducing the first two correlation terms, we getgood results for rigid spheres and discs, and for 6-12 and square-well potentials.We may even expect good results from the simple YBG equation for long-rangepotentials of the Coulomb type.2However,Kirkwood has shown that periodic, solid-like solutions may exist at densities abovethis singular point.Also, by introducing an approximation to the long-rangebehaviour of g(123) into the modified equations we may be able to correlate theYBG singularity with the fluid-solid phase transitions predicted by molecular dynamics.Despite the limitations imposed by the speed of present-day computers, we hope thatthese studies will add a little to our understandirlg of the statistical theory of fluidsand phase transitions.Prof. J. Walkley and Dr. Wing Y. Ng (Simon Fraser University) said: In twopapers use has been made of the Lennard-Jones 12-6 potential. In each case anagreement between theory and experiment is found, despite the inadequacy of thispotential function to express the two-body second virial coefficient behaviour.Inan N-body situation the total potential experienced by any one particle cannot beregarded as the simple pair-additive sum of the interaction with all other particles.For the 3-body case the perturbation of the pair-additive s u n coming from the mutualdisposition of the three particles has been calculated for the dispersion r6 attractionterm.4 The use of the 12 : 6 pair potential in a strictly pair-additive manner appearsto generate an (apparently) acceptable N-body potential. In calculations pertainingI . Z . Fisher, Statistical Theory ofLiquids (Univ. of Chicago Press, 1964), p. 152.C. W. Hirt, Physics Fluids, 1967, 10, 565.J. G. Kirkwood and E. Monroe, J. Chern. Physics, 1941,9, 514.R. J. Bell and A.E. Kingston, Proc. Physic. Soc., 1966, 88,901.The exact meaning of the singularity in the YBG equation is unclear58 GENERAL DISCUSSIONto the solid phase of the inert gases, it was similarly observed that the simple bi-reciprocal 12 : 6 potential allowed theoretical volume-temperature and heat capacity-temperature dataare in good agreement.l For arange of two parameterm : 6 potentials(viz., 4(rl1) = &f(rij/u), with usual rotation) the characteristic parameters E and Q canbe obtained by fitting the theory (a harmonic Einstein theory) to certain propertiesextrapolated to OK. The total lattice energy U, and the molar volume V , wereoriginally used. The " best '' set of parameters were then obtained by a plot of thetheoretical reduced zero-point energy (A:) and the reduced " experimental " zero-point energy (9R0, /8N&, with 0, the limiting Debye frequency) simultaneously asfunctions of the de Boer parameter (A = h/(m&)%o).The intersection of the twocurves give the " best " value of E and Q without specifying the " best " potentialfunction. However, this set of values allowed both second virial and viscosity datafor the inert gases to be reduced to a corresponding states pattern of behaviour.2TABLE 1.-nz : n KIHARA CORE POTENTIAL PARAMETERS FOR ARGONV 20 .(A1 B (cal mole-1) m : n E/kCK)9 : 6 123.2 3.3932 0.157 0 18710 : 6 123.0 3.3948 0.1056 0 18711 : 6 122.0 3.3960 0.0567 0 18712 : 6 120.5 3.3987 0.0098 0 187Experimental data : Uo = - 1846 cal mole-' ; Ro = 3.7549 (A) ;KT = 0-375 cm dyne-2 x lolo ; Lo = = (9/80,) = 187 cal mole-'.TABLE 2.-m : n KIHARA CORE POTENTIAL PARAMETERS FOR XENONB (cal mole-1) V A0 m : n ~lk("K) d A )9 : 6 240 3.9483 0.157 0 13310 : 6 238 3.9510 0.106 0 13311 : e 236 3.9521 0.0567 0 13312 : 6 234 3.9599 0-0098 0 133Experimental data : Uo = - 3856 cal mole-' ; Ro = 4.3356 (A) ;KT = 0-28 cm dyne-2 x 10'' ; = (9180,) = 123 cal mole-'.We have extended these calculations to include a non-additivity term in the staticlattice potential sum.In keeping with theory we consider only the largest termcoming from the 3-body interaction and written as Av/a$, where A is the appropriatelattice constant, a i j the nearest neighbour distance and v an unknown non-additivecoefficient. Calculations were also made more flexible by considering a series ofrn : 6 Kihara core pair potentials, viz.,where r = rlj/a tc.The calculation involves four unknown coefficients, the reductionparameters e and a, the coefficient v and the reduced hard-core cut off 6. Fourzero-point properties are required for their evaluation, and besides the three usedabove, U,, Yo and Ro, an accurate value of the zero-point isothermal compressibilityKT could be used for both argon and xenon. The parameters obtained by fittingthe theoretical model (again a harmonic Einstein model) to this data are given intables 1 and 2.I. H. Hillier and J. Walkley, J. Chem. Physics, 1965, 43, 3713.J. Walkley, J . Chem. Physics, 1966,44, 2417GENERAL DISCUSSION 59From these tables the following observations may be made: (i) for all m : 6potentials examined the non-additivity coefficient is effectively zero ; (ii) for bothargon and zenon, for the same value of rn the resulting value of 6 is identical, inti-mating an essentially corresponding states behaviour pattern, and (iii) for both argonand zenon each of the rn : 6 Kihara potentials gives, through the different 6 values,curves which are virtually identical and which are exceedingly close to that of a 12 : 6Lennard-Jones potential.The calculations are, then, consistent with the accepteduse of the Lennard-Jones 12-6 potential as an acceptable '' many body " potential.Prof. U. v. Weber (Rostock University) said: We studied the deviations fromideality of vapours by an isochoric method to get the second virial coefficients andby means of them the intermolecular potentials.Density was varied to take careof higher coefficients. To make the adjustment of parameters and powers in runambigously, the range of temperature must be as large as possible. The bestfitting powers of attraction obtained by the method of Buckingham are r9 forbenzene and n-octane and r-l0 for hexane and n-octane with reasonable diameters.I restrict myself to quasispherical molecules.Only r6 is theoretically justified and is indeed a fair potential for rare gases, butnot satisfying for large molecules. In our Institute, W. Lichtenstein investigatedvarious potential models. Kihara allows for extension to hard cores, but his assump-tion of effect only between the two nearest points is not justified. Laeuger consideredr6 potentials between the corners of two cubes. Hoover and Rocco regardedspherical shell potentials with centres of attraction and repulsion on the shell, butthey integrated only numerically. Lichtenstein considers two spherical models :attraction in -6 power between every volume element of the first particle to everyvolume element of the second and between every surface element of the first particleto the second. Repulsion was assumed in either model between the two centres inpowers -24, to obtain an analytically integrable expression. The shell model fittedbetter. The argon, tetrachlormethane and benzene data were fitted as well asdesirable by this potential with reasonable diameters. To avoid singularity in the caseof contact, recession of the points of attraction up to a distinct amount was supposed.This seems to me a reasonable step towards a dynamic model. The potential well ismuch deeper and more narrow than for a Lennard-Jones potential, a fact whichmight be taken into account in the more complicated theory of liquids