Equation of State of a Monatomic Fluid with 6-12 PotentialBY STUART A. RICE AND DAVID A. YOUNG*Institute for the Study of Metals and Dept. of Chemistry,The University of Chicago, Chicago, Illinois 60637Received 24th October, 1967The modified Yvon-Born-Green equation of Rice and Lekner is solved for a monatomic fluidwith 6-12 potential at kT/c = 2.74. The equation of state and excess internal energy are calculatedand compared with Monte Carlo calculations ; the agreement is good. Calculations with hard-spheremodifications and with the Kihara potential are also discussed. The Kirkwood limit of stability forthe fluid phase is calculated.A study,’ by Rice and Lekner, of the equation of state of a rigid-sphere fluid, hasshown that the Yvon-Born-Green equation may be so improved as to yield anequation of state in close agreement with the molecular dynamics data of Alder andWainwright.2 In this paper we present an analogous study of a fluid with moleculeswhich interact through a 6-12 Lennard-Jones potential.Monte Carlo calculations 3of the equation of state corresponding to this potential have been reported for asupercritical isotherm at kT/E = 2-74, where E is the depth of the potential well.Our calculations, based on the improved YBG equation, will be at this temperatureso as to permit comparison of theory and “experiment”.Recent work 4 has provided the solutions of the Yvon-Born-Green (YBG),Percus-Yevick (PY), and the Convolution Hypernetted Chain (CHNC) equationsfor the 6-12 potential at kT/E = 2.74. As for the hard sphere potential, the PYtheory yields the best agreement with experiment over a wide density range, but givesno explicit evidence for a liquid-solid phase transition.Our motive for choosing toextend the YBG theory is that it predicts a transition to occur at high density, whilethe PY and CHNC theories do n0t.t Hence we regard the YBG equation to containthe qualitative features necessary to describe the limits of stability of the fluid phase.TRIPLET DISTRIBUTION FUNCTIONThe exact, formal, expression for the triplet distribution function of a fluid isg(3)(123) = g(2)(12)g(2)(13)g(2)(23) exp ~(123, p),~(123, P ) = C pngfl+3(123).(1)(2)wherecofl=lHere p is the fluid number density. The 6, are expansion terms developed and dis-cussed by Meeron 5 and Salpeter.6 The superposition approximation may be obtainedby placing z = 0.This approximation yields the correct third virial coefficient fora fluid, but incorrect higher virial coefficients. Inclusion of terms up to 6, yieldscorrect virial coefficients up to Bm.*Fannie and John Hertz Fellow.?Note added in proof: The PY equation admits solid-like solutions above a critical density butdoes not prohibit simultaneous fluid-like solutions.1S . A . RICE AND D . A . YOUNG 17For the hard sphere potential, Rice and Lekner (referred to as RL henceforth)found that 6 4 and 8 5 could be calculated with a reasonably small amount of computertime, but that 8 6 could not. For the 6-12 potential we find that 84 can be calculatedaccurately, that 6s can be estimated fairly well, and that 156 cannot be calculated at all,again owing to the computer time involved.Starting with 64 and 65, we shall usethe technique of the Pad6 approximant to achieve an estimate of the entire 6, series.The derivation of the diagrams corresponding to the 6 is fully explained in RL andhere we only point out those aspects of the calculation which are special to the 6-12potential.The calculation of 64(123) was carried out numerically by means of Gaussianquadrature. Now0'= p f 1 4 f 2 4 f 3 4 ,I6,( 123) = a4/ \(3)2 0 ' 0 3where hj = exp( - u(ij)/kT) - 1 and d4 is the element of volume swept out by particlenumber 4. There are no boundary problems such as arise in the hard sphere casebecause the 6-12 potential is continuous.Since for moderately large distancesr,j&(r) is small, the integration in (3) may be carried out over a small region. Takingadvantage of the symmetry of the plane of the triplet (123) with respect to the per-pendicular z axis, we find that a hemisphere of radius 3.0 molecular diameters servesthe purpose adequately. All distances will be given in units of the molecular diameter,the diameter being defined as the distance CT for which the potential u(a) = 0.The integration (3) was divided into three parts. In the interior region(O-O<r<l-O), where the integrand is large, the integral was evaluated by using Carte-sian coordinates and 16 point Gaussian quadrature. Then the contributions fromthe two spherical shells (1.0 GrG1-5) and (1.5 <r G3.0) were evaluated in sphericalpolar coordinates using 12 point Gaussian quadrature.The results comparefavourably with those obtained by using larger numbers of Gaussian points, viz., 20,16, and 16. The average error is about 1 %. The integrations required in (3) werecarried out for all triplet configurations (R,S,T) such that 0.0 <R,S,T<2*6 andAR = A S = AT = 0.20. R = 2.6 was found to be satisfactory cutoff since64(2.6,2.6,2.6) = 0.0015, which is close enough to zero to be negligible.The 65( 123) diagrams were grouped and evaluated using the same techniques as inRL. There are 13 diagrams to evaluate. By using the functionspij = d5fi5 fj5; tijk = S,(ijk); eij = fij+ 1 = exp [I -u(ij)/kT],the calculation may be reduced to three 3-fold integrations and one 6-fold integration ;-5.J d4f 1 4 f 2 4 f 3 4 I d5f 1 5 f 2 5 f 3 5 f 4 5 * (4)The first three integrals may be evaluated by using values of 64 tabulated with aninterval of 0.1. Since the storage of 64 in this state in the computer required approxi-mately 19,700 locations, it was not feasible to increase the state of subdivision of 84.Hence the tijk required in the integrand had to be truncated to the nearest tabulatedvalue, at the expense of accuracy.This difficulty was avoided in the work with hard spheres because 84 is known interms of elementary functions and can be rapidly computed for any desired tripletconfiguration. Exactly how much accuracy is lost by the truncation is not known18 EQUATiON OF STATE OF A MONATOMIC FLUIDHowever, very accurate values of pij were calculated and these were often much largerthan the tijk which appeared with them in the integrands. Hence the inaccuraciesin tijk are masked to some extent by this discrepancy in size.Again the individual integrals were divided into three parts, this time using 12, 12,and 8-point Gaussian quadrature in a sphere of radius 3.0.Comparison with moreaccurate calculations shows an average error of 5 % due to imprecise integration.The total error is not known because of the truncation of the ?i& terms. Once again65(R,S,T) was obtained for all configurations up to R = S = T = 2-6 andAR = A S = AT = 0.2.64 and 8 5 are thus obtained as three-dimensional arrays with a spacing of 0-2 diam.units.These values were interpolated by using a linear Taylor expansion. In theexamples below, R, S, and Tare the tabulated triplet lengths (multiples of 0.2) and x, y ,.3.2. I0 -. I-.2-.3-.4-. 5-. 6-.7FIG. 1 . 4 4 and 85 for the 6-12 potential at ~ T / E = 2.74 in the equilateral configuration (R,R,R).and z are the interpolated lengths (odd multiples of 0.1).in this way. There are three cases to be considered :Both d4 and 8s are tabulatedI.11.111.G(R,S,z) = .5(6[R, S, ( Z + el)] + S[R, S, ( Z - .1)]>,G(R,y,z) = .5(6[R, ( y + *1), ( z - .1)] + 6[R, (y - .l), (2 + el)]},G(x,y,z) = .5(6[(x + -1), (y - *l), ( z - d)] + 6[(x - -1), (y + *1), ( z - .1)] +6[(x - -1), (y - .1), (2 + -1)] - 8[(x - .1), ( y - .1), ( z - ml)]}.( 5 4(5b)(54In each case the 6 in the parentheses on the right side are tabulated values, or they arenonexistent owing to the violation of the restriction that R, S, and Tmust be the sidesof a triangle.In the event that one of the 6 necessary for the interpolation does notexist, a simpler approximation is used. The result is some 2700 pieces of data foreach of 6 4 and 8s.Several special values of 64 were calculated and checked against interpolatedresults in order to test the interpolation scheme. The error ranged from 3 to 20 %(for case 111), with an overall average error of 5 %. Hence the overall picture of thetabulated 84 and 6s values is one of fair but not high accuracy.Both 84 and 8 5 take on positive and negative values, and both are large andnegative near (O,O,O).In this sense they differ strongly from the hard sphere terms inwhich 64<0 and 85>0 over their respective ranges. 8 4 and 8 5 for the equilateraltriplet configuration are shown in fig. 1S. A . RICE AND D . A . YOUNG 19YVON-BORN-GREEN-EQUATlONThe YBG equation has been reduced to a form which can be used for numericalcomputation. Details of the reduction are described by Hill.7 We write, followingthe same arguments,kT In g(R) = - zl(R) + npIW .'(y)g(y)dy"y y 2 -(R-X)2[g(x) - Ilxdx +0 R-yu(R) = 4&( - A).Here R, x, y and z are reduced distances with units of molecular diameters. Eqn. (6)is exact, but since only the first two terms of z(x,y,z) are known, we must resort toapproximation if a solution is to be found.A good approximation to the entireseries is the Pad6 approximantwhich we henceforth adopt.NUMERICAL INTEGRATION OF THE MODIFIED YBG EQUATIONWe have obtained numerical solutions over a range of densities for(i) z = 0 (superposition) ; (ii) z = p64 ;(iii) z = pS4+p265 ; (iv) z = p64/(1 -p&/64) = z(Pad6).(v) z = p64(HS) ; (vi) z = p64(HS) +p265(HS) ;(vii) z = Pad6 (HS) ;Also we have tried mixing terms for hard spheres and the 6-12 potential:(viii) z = p64(6-12)+p265(HS).Of the terms (v)-(viii) only (v) turns out to be of real interest.In each of the above calculations the parameter p was selected in the range0.1 < p G1.2, and kT/& was fixed at 2.74. Solutions of the YBG equation wereobtained by the iteration procedure commonly used in solving nonlinear integralequations.This is discussed in RL. The input for the (k + 1)th iteration is taken to beFor low densities, we found that the optimum CI can be as high as 0.75, and that forhigh densities CL must be reduced to 0.2, where convergence is correspondingly slower.If a is too large for a given value of p the iterative solution undergoes a divergentoscillation. For a given z and p the first input used was the convergent solution of thepreviously chosen density in the series of calculations.Convergence was observed by calculating the reduced pressure for each iteration :g r ( r ) = d n ( r ) + aiX"t(r) - gt(r)l- (8)u'(r)g(r)r4nr2dr. (9)The approach to a stable solution could be either monotonic or oscillatory.Thelarger the value of a the greater the tendency for oscillatory behaviour. Convergenc20 EQUATION OF STATE OF A MONATOMIC FLUIDwas usually taken to mean that the input and output g(r) differed by not more than0.5 % at the highest density. The accuracy of the solutions depends on : (i) theinterval Ar used in the numerical integration ; (ii) the maximum value of R to whichthe integration is carried ; (iii) the number of iterations performed ; (iv) the accuracyFor the superposition approximation, some solutions were obtained for Ar = 0.05.The pressures calculated from these results differed from results interpolated fromAr = 0.1 by about 3 % for p = 0.8. The superposition pressures at high densitycoincide with those of Broyles 4 to within 2-3 %.Hence it may be assumed that theresults obtained for Ar = 0.1 are reasonably accurate.Interpolations of g ( r ) were quadratic (i.e. second order in the Taylor expansion)and the pressures obtained from interpolation (Ar = 0.05) differ from those of theoriginal results (Ar = 0-1) by about + 5 %. Most of the contribution to the pressureoccurs near r = 1.0 where the function g(r) is steep. A special interpolation techniquewas used to achieve an increment of 0.05 in the integration of the triple integral ineqn. (6). For z = p&+pZ&, for densities above p = 0.6, divergent results wereobtained with the increment 0.1. These divergences disappeared when the smallerincrement was used, showing that they were numerical in origin.Comparison of thepressures obtained for the two increments showed that atp = 0-6 the accuracy obtainedis approximately 1 %.Condition (ii) was significant only for the high density solutions in the super-position approximation. All solutions were obtained with R(max) = 4.0, and inaddition R(max) for p = 1.1 and 1.2 was extended to 5.0. The discrepancy inpressures between the two limits in the latter cases are 5 % for p = 1-2 and 2 % forp = 1.1. Condition (iii) was met by iterating until the desired convergence wasobtained, viz., a self-consistency of better than 0-5 % in the reduced pressure.of ‘C(x,r,d.TABLE 1.-YBG REDUCED PRESSURES PplpkT FOR A SERIES OF ‘C FUNCTIONS.PARENTHESES ARE NOT EXACT OWING TO SLOW CONVERGENCE.VALUES INTHE SUPERPOSITION PRESSUREAT p = 1-2 IS PROBABLY DIVERGENT.*1-4-6-7-8-91.01.11.2superposition 6 4-975 -9791.154 1-2151.537 2.013 - 2.9462.259 (5.06)27243.2293-6834,074obs.P64+ P265-979 -979 -979 - -9761.212 1.206 1 -226 1.30 1.211.973 2.079 1.912 2.03 1-892.803 - 2.55(4.44) 3.3604.3725.626Pad6 p&(HS) Kihara (Levelt)The accuracy of the z functions may be checked by evaluating virial coefficients.According to the theory, inclusion of terms in 1: up to S, should result in correctvirial coefficients up to Bm.For the superposition approximation B3 was calculatedto be 0.365 B2. The exact value is 0.365 2.32. B4 was found to be 0.104 B3. Thecorrect value is found, by interpolation of Barker’s results,s to be 0.116 B3.Similarlyfor BS we calculate 0.056 B4 and Barker’s value is 0.049 234. Here B = 2m3/3. Thediscrepancies cited reflect the limited accuracy of 84 and 8s.Table 1 contains values ofp/pkTfor 7 = 0 (superposition), ‘C 5 pS4, z = p&+p26~and z = PadC. AIso z = p&(HS) is included. These results are shown in fig. 2 anS . A . RICE AND D . A . YOUNG 213. The excess internal energy of the fluid may be calculated from the pair distributionfunction, i.e.,-!- NkT = LJrn 2kT u(r)g(r)4nr2dr. (10)The calculated excess internal energy U/NkT is given in table 2 for each of the casesIi' 8.06.0 '''1 DATAMONTE CARL0h45.0\4.03.02.01 .o -0.0 --uPO3.I .2 .3 .4 .5 .6 .7 .8 .9 1.0 1.1 1.2 1.3FIG.2.-Equation of state for the 6-12 potential fluid at kT/c = 274. The Monte Carlo data aretaken from ref. (3).4.0 -3.5 -3.0 -bl 2.5- %.52 . 0 -1.5 -1.0 -I I I I-I.4 .5 .6 .;I .8PO3FIG. 3.Detail of the fluid branch of the equation of state at kT/E = 2-7422 EQUATION OF STATE OF A MONATOMIC FLUIDTABLE 2.-YBG EXCESS INTERNAL ENERGIES U/NkT FOR A SERIES OF z FUNCTIONS. VALUESIN PARENTHESES ARE NOT EXACT OWING TO SLOW CONVERGENCE. THE S~PERPOS~~ION ENERGYATP = 1.2 IS PROBABLY DIVERGENT.obs.P superposition p84 p6,+,28, Yade p&dHS) Kihara (Levelt)--218 *1 -*223 -*222 --222 --222 --222 -*4 --877 -a859 --865 -*858 -a889 --744 --.829.6 -1.310 -1.267 -1.286 -1.269 -1.357 -1.118 -1.205-7'8 - 1.758 (- 1.625) -(1.630) - 1.863*9 -1.989 -2.1691.0 -2.234 -2.6161.1 -2.5171.2 (-2.84)- - 1 -455 - 1-477 -TABLE TH THE PAIR DISTRIBUTION FUNCTION Q(R) FOR Z =pdq+p2&.R0.10.20.30.40.50.60.70.80.91.01.11-21.31.41-51.61-71.81.92.02 12.22.32.42-52.62.72.82-93.03.13-23.33.43.53.63.73.83-94-0p = 0.40.00.00.00.00.00.00.00.00.1231.2361.6011 -4221.2191 -0780.99 10-9440.9260.9330.9600.9951 -0201 -0271.0231.0151 -0061 -0000.9970.9950.9960.9981 -0001 -0021 -0021 -0021 -0021.0011 -0001 *ooo1 -0001.000p = 0.60.00.00.00.00.00.00.00.00-1731 -5691.8541.5111 -2051,0080.8940.8440.8400.8800.9551 -0341 -0761 so731 -0481 ,0180.9940.98 10,9750.9790.9850.9991 a0071.0101 -0091 *OOG1 -0020.9980.9970.9970.9981 a000p = 0.70.00.00.00.00.00.00.00.00.2161.8352.0321 -5601.1830.9530.8270.7790-7880.8530.9641 -0721.1221.1071 -0621.0150.9800-9650.9590.9670 9841 -0021-01 51.0191.0161 -0091 -0020.9960-9940.9940.9961.000p = 0.80.00.00.00.00.00.00.00.00-2902.3372.3801.6401.1 120.8150-6680.6280.6660.79 11.0001 -2001 ~2721 -2091.0930-9890-9200.8940.8940-9260.9781 -0281 -0581 -0601 -0421.0150.99 10.9770.9730.98 10.9951.01S.A. RICE AND D. A. YOUNG 23cited, and is plotted in figs.4 and 5. For 0 <p G0.8, the superposition excess internalenergy is nearly linear in density, with an equation U,/NkT = -2.1975~03. Pairdistribution functions are listed in table 3.t20 MONTE CARLO DATA-. 2-.4-.69 8-1.0$ -1.4-1.6-1.8-2.0-2.2-2.4-2.6-2.8i-( -1.2.I .2 -3 .4 .5 .6 3 .8 .9 1.0 1.1 1.2 1.3p03FTG. excess internal energy of the 6-12 potential fluid at kT/E = 2-74. The Monte Carlo dataare taken from ref. (3).I 1 1 1 1 1 1 1oMONTE CARLO-1.1-1.2 ---1.6---1.7 SUPERPOSITION\-.8-\-.9-1.0-1.3- 2 -1.4-----1.5-1.7\-Is8 1.9 ir -*.+ , , , , , , , , , _I.4 .5 .6 .? .6p03FIG. 5 . D e t a i l of the excess internal energy of the fluid at /cT/E = 2-74.COMPARISON WITH EXPERIMENTThe 6-12 potential, although qualitatively sound, is inadequate to describequantitatively potential functions in real gases. One of the best pair potential func-tions is that due to Kihara.In particular, we have chosen the form used b24 EQUATION OF STATE OF A MONATOMIC FLUIDBarker :9U(R) = 4 e [ ( _ 9 > L 2 - ( z r ] R-‘1 R-‘1 .This potential function was inserted into the YBG equation with z = pS4(6-12)+p2a5(6-12) and the pressures and energies were calculated from the solutions of thisequation. We believe this to be a reasonable approximation since the main effects ofthe change in potential should be in the term zi(R), and not in the several 6,.Comparison with the argon data of Levelt 10 and the results for the pure 6-12potential show that the Kihara pressures are too large, and the energies too small (inmagnitude).In fact, the 6-12 results are better than the Kihara results. This is anillustration of the importance of the repulsive part of the potential in determining theequation of state of a dense fluid. The Kihara potential has a steeper repulsive partthan the 6-12 and this is reflected in the higher pressures. Evidently in a dense fluidthe Kihara potential is a less accurate description of the real effective potential thanthe 6-12.DISCUSSIONFrom the entries in the tables it is seen that the approximations z = 0, z = pS4,and z = pS4+p2& yield increasingly accurate pressures. For the internal energy,it appears that z = p64 is closer to the Monte Carlo curve than z = pS4+p265.Inboth the calculated pressure and energy, the Pad6 results lie nearer to the z = pS4curve than to the z = p64+p265 curve. Theoretically, the Pad6 approximant shouldlead to better pressures and energies than either of the other z functions. We believethat it is the inaccuracies in the 64 and S5 functions that are responsible for the failureof the Pad6 approximant to provide good agreement with experiment.The results with 64(HS) represent a “mixture” of the 6-12 and hard-spherepotentials. These calculations were motivated by the conjecture that the very steeprepulsive part of the 6-12 potential causes it to behave much like a hard core. Hencethe first-order correction using the hard sphere 6 4 should yield an improved equationof state.The results are indeed better than the superposition approximation, andsolutions may be obtained for the full range of densities. However, the excessinternal energies for 64(HS) are worse than for the superposition approximation. Thisemphasizes the qualitative difference between the 6- 12 and hard-sphere potentials.Also, the higher approximations z = p64(HS)+p265(HS) and z = Pad6 (HS) yieldpressures which are poorer than for p64(HS).From studies of the hard sphere and disc fluids,l, 11 the Kirkwood criterion for thestability of the fluid phase is not exact,l2 but coincides with the limit of stability forthe fluid in the superposition approximation (and for short-range modifications ofthis approximation). The limit of stability is given by the first real root of1 +?I: [kR cos kR - sin kR]u’(R)g(R,Q,p)dR = 0.By evaluating the left-hand side with a series of g(R), we obtain a function of kwith damped oscillations.This function resembles the Bessel functions obtainedfrom this equation for hard spheres and disks. At some density p,, the first minimumof the function will become zero, defining the upper limit of stability of the fluid.Since we do not have g(R,p) for a continuous range of densities, we must use inter-polation to locate the density cut-off. For the 6-12 potential, the superpositionapproximation leads to k , = 6.0 and pc = 1.12. For z = p64+p265, k, = 6.0 anS . A . RICE A N D D. A . YOUNG 25pc = 0433. Although these results are approximate, they are supported by the factthat we have found the solutions of the superposition approximation for p> 1.2 to bedefinitely divergent.At the highest fluid densities the YBG equation of state diverges markedly from theMonte Carlo equation of state.Although some of this discrepancy is probably dueto inaccuracies in the calculations, it seems likely that 8 4 and 85 are not adequate toaccount for all of the important long-range interactions in the 6-12 fluid. Clearly, adifferent approach is needed.We conclude that an adequate theory of the fluid must : (a) treat the short-rangerepulsions accurately, i.e., more accurately than by using the hard-sphere potential ;(b) take into account relatively long-range correlations, i.e., of the order of 3-5 atomicdiameters. A mean field treatment, superimposed on the 64, 8 5 fluctuation termsdescribed herein, is indicated.We thank the Directorate of Chemical Sciences, AFOSR for financial support.We have also benefited from the use of facilities provided by ARPA for materialsresearch at the University of Chicago. We are indebted to John Lekner for helpfuldiscussions. One of us (D. A. Y.) thanks the Hertz Foundation for generous financialassistance .1 S. A. Rice and J. Lekner, J . Chem. Physics, 1965, 42, 3559.2 B. J. Alder and T. E. Wainwright, J. Chem. Physics, 1960, 33, 1439.3 W. W. Wood and F. R. Parker, J. Chem. Physics, 1957, 27, 720.4A. A. Broyles, S. U. Chung, and H. L. Sahlin, J. Chem. Physics, 1962, 37, 2462.5 E. Meeron, J. Chem. Physics, 1957, 27, 1238.6 E. E. Salpeter, Ann. Physics, 1958, 5, 183.7 T. L. Hill, Statistical Mechanics (McGraw-Hill Book Co., Inc., New York, 1956). * J. A. Barker, P. J. Leonard, and A. Pomp, J. Chem. Physics, 1966,44,4206.9 J . A. Barker, W. Fock, and F. Smith, Physics Fluids, 1964,7, 897.10 J. M. H. Levelt, Physica, 1960, 26, 361.1 1 D. A. Young and S. A. Rice, to be published.12 J. G. Kirkwood, in Phase Transformations in SoZids, ed. R. Smoluchowski, J. E. Mayer, andW. A. Weyl (John Wiley and Sons, Inc., New York, 1951), chap. 3