Calculation of Thermodynamic Properties of Liquid Argon fromLennard-Jones Parameters by a Monte Carlo MethodBY I. R. MCDONALD AND K. SINGERDept. of Chemistry, Royal Holloway College (University of London),Englefield Green, Surrey.Received 16th January 1967Extensions of the Monte Carlo method of Metropolis et al. are described which permit bothisochoric and isothermal extrapolations of Monte Carlo data. Thermodynamic properties arecalculated for liquid argon from Lennard-Jones parameters for a wide V, T-range. The agreementbetween calculated and experimental values is, on the whole, satisfactory.The Monte Carlo (MC) method devised by Metropolis et aZ.1 has been applied tothe calculation of equilibrium properties of a variety of systems.1-12 Of particularinterest in the present context is the work of Wood et al.who calculated the radialdistribution function, internal energy, pressure and specific heat of argon over a largerange of densities 4 at 55”C, and for one density 5 of the liquid at - 147°C fromLennard-Jones (LJ)-parameters.This paper explores the possibility of calculating the principal thermodynamicfunctions of a simple classical liquid from a LJ-potential over an appreciable range ofvolume and temperature. The project has become practicable (a) because the greaterpower of computers now available permits the execution of computations of the typereported in ref. (4) without excessive cost in computer-time ; and (b) because the MCmethod has been modified so that properties for a certain range of V, T-values can beobtained from MC data for one V, T-point.Argon is chosen for this work becausereliable LJ-parameters 13 and excellent P- V-T data 14-16 covering the entire liquidrange are available.MONTE CARL0 METHODThe essential features of the method of ref. (l), which has been discussed,lp 49 17are these. A fluid “sample” is represented by N coordinate triplets confined withina cube of volume v. The potential energy @, assumed to be a sum of pair energies,is calculated for the initial configuration. One of the N “particles” receives a randomdisplacement and the resultant change A@ of the potential energy is determined. IfA@ is negative, the new configuration is accepted ; if it is positive, the move is acceptedonly with the probability p = exp(-PA@), l/kT.If the move is accepted, thenew configuration and the new potential energy become the “current” properties ;if it is rejected the system returns to the previous state. The sequence of configurationsgenerated by repetition of this procedure has been shown to tend towards a Boltzmann-weighted ensemble. The mean value of any function of the 3N coordinates over asufficiently long MC sequence is therefore an approximation to the ensemble averageof this quantity. Periodic boundary conditions are imposed to minimize the error4I . R . MCDONALD AND K . SINGER 41resulting from the small size of the sample, which, for reasons of economy, usuallycannot exceed N - 102.The mean potential energy (@> so obtained determines the molar configurationalinternal energy Ut ; the mean virial,( Y ) = <xii.via>idetermines the pressure according toIsothermal changes of the Helmholtz free energy can be calculated from (">-valuesand isochoric changes of the configurational entropy from (@)-values byPV = R T -(No/3N)(Y>AA = -(NOIN) Pdv (1) and ASt = CLdT/T, CJ = (No/N)(a(@}/aT),. (2)In these and subsequent formulae, N is the number and v the volume of the MCsample ; molar quantities are denoted by NO and V ; the dagger 7 means "configura-tional" and the angular brackets (> indicate ensemble averages.The present extension of this method is designed to evaluate the classical con-figurational partition function in the formQt(N,v,T) = (VN/N. l)rmaxg(0) exp ( - pO)d@,f:: 1::Q, minwhere g(CD)dCD is the fraction of all accessible configurations for which the potentialenergy has a value between @ and @+do.Boltzmann-weighted sequences ofconfigurations and averages of @ and Y are calculated as in ref. (4) ; but in addition,a histogram is accumulated of the numbers of configurations with potential energieswithin small specified ranges of width 6. The entry for the n-th interval is the numberof configurations with potential energies between CD, - 6/2 and Qn + 612, denoted byFor sufficiently small 6, Fn is assumed to converge statistically toFn F(Qn - 6/2,mn + 612).where K depends on v but not on T. Sincethe histogram determines the partition function apart from the normalizing constantK .Qt at another temperature, T' can be obtained by reweighting the histogram:where p' = l/kT'.In practice, the range of temperatures which can be reached by thisprocedure is limited by the width of the parent histogram.Since K cannot be determined by the present method, isothermal changes of Qfare calculated indirectly, by means of (1). To compute changes along a neighbouringisotherm T', a reweighting procedure is applied. Let {a), be the set of G, con-figurations for which the potential energy lies between @,-a12 and @,,+6/2, Yl.n,orthe virial of one of these configurations, and Tn the average for the set (a)n, i.e.42 THERMODYNAMIC PROPERTIES OF LIQUID ARGONThe ensemble average of the virial at Tis(yy>T = CF,(T)qn/(xFn(T))n nand the value (Y)Tt is obtained by substituting Fn(T) = Fn(T)exp[-(D’-Q)@,J forFn(T) in (6).Tables of qyr,-values are calculated during the MC runs.Another extension of the range of application of the MC method arises from thepossibility of calculating from MC data obtained for the LJ-parameters ( E , a) propertiescorresponding to different parameters (d, a’). This can be done if the differencesbetween the “new” and the “old” parameters are small.The LJ-potential is defined by$(r) = 4~[(a/r)’~ - ( ~ / r ) ~ ] (7)modified by a hard core, i.e., $(r) = 00 for r < rc, introduced to increase the speed ofcomputation. Since the potential function and the corresponding virial function,$(r) = rd4/dr, are linear combinations of r-12 and r-6, the potential function 4’(r)and the virial function $’(r) for the parameters (E’, a’) can be written as linear com-binations of $ and $.The same relationships hold for the total potential energy@n,a and the total virial Yn,a of any configuration a belonging to the set ( c l I ncalculated with different LJ-functions, i.e.,= a@,,,+ b v n , a , (8)(9) y ; , a = c@n,a + d y n , a ,wherewith z = E’/E and A = d/a.The ensemble average for the “new” potential energy is<af> = CC x (a@n+byn,J ~ X P [-D(aQn+byn,a)I)/(C C ~ X P [-D(aan+byn,a)J>,n taln n {ahwhere @n,= is replaced by an since @n,ag@n. The further substitutionsN -Yn,a f Yn+&, a@,+ b!Pn f @n,lead to<Q‘) = {C c ( G n + bAn,a) ~ X P [ - B ( G n + b ~ n , a ) l ~ / ( ~ c ~ X P C - B ( G n + b ~ n , a > l ) * (12)n Ialn n {ahEqn.(lOa) shows that b -g 1 if 10-0‘1 $a ; in this case exp ( -#lbAn,a) may be expandedand terms of higher than the second order in pbAn,a neglected. If the substitutionsare made in the resulting expressions, one obtains(1 3) (a’) = { E F ‘ n [ 6 n ( l + p 2 b 2 ~ / 2 ) - B b 2 ~ l ) l ( C ~ n ( 1 +B 2 b 2 2 An/2))n nThe analogous formula for the virial iswhere yn=c@n+dYn. Eqn. (13) and (14) are exact if the distributions of A%,a arI . R. MCDONALD AND K. SINGER 43Gaussian. If j P b 2 z < 1 for all n, one may use the simpler expressions,(Y‘) E (CFn[‘u,-pbdii3)/CFn. n (16)n -The values of YYn may be calculated during MC runswithout loss of computing speed.These formulae can be used (a) to investigate the effect of changes of the LJ-parameters on the equilibrium properties, and (b) to extrapolate data calculated foru,T to u‘,T’. Eqn.(1 3)-( 16) are applicable to isothermal changes because substitutionof the LJ-parameter 0 by 0’ = a/A is equivalent to the scaling of all linear dimensionsof a sample configuration by A, i.e., r i p d r i j for all i, j and u-u’ = A321. Changes ofE to E’ are equivalent to changes of Tand can be dealt with by the reweighting procedureaccording to (4) and (6). The isothermal extrapolation is limited to small changes ofvolume because, apart from the approximation made in the derivation of eqn. (13)and (14), the reweighting of the histogram Fn to becomes inaccurate if the twohistograms do not overlap sufficiently.MC data for the calculation of thermodynamic properties may also be obtainedby a method not based on ref.(1). In this, a sequence is generated by a random walkin 3N-dimensional configuration space without Boltzmann-weighting. A move isrejected only if it raises the potential energy above a preset upper bound. Thehistogram G(Qn - 6/2, Qn + 6/2) = Gn so obtained is assumed to tend to g(Q) in (3), andthe partition function is calculated fromi.e., the variances ofi.e., by Boltzmann-weighting of the histogram according to the temperature required.The limitation of this method of “unweighted configurations” arises from the veryrapid increase of g(@) with a. A considerable number of MC runs with differentupper energy bounds is therefore required to cover a @-range sufficient for the calcula-tion of properties over a wide T-range.Isothermal changes and extrapolations canbe dealt with as described above. The disadvantages of the “unweighted” methodincrease with density and with N ; but they are not too serious in the gaseous range,even at high pressures. Results obtained by this method for gaseous argon over awide P-V-T range will be reported elsewhere.18 The two methods have not beencompared systematically .RESULTS AND DISCUSSIONCalculated values of selected thermodynamic properties of argon at five tempera-tures between the triple point (8343°K) and the critical temperature (150-7°K) and arange of densities are shown in table 1. The experimental results for the same V, T-points are given where these are available ; properties of the orthobaric liquidobtained by graphical interpolation are shown in table 2.Experimental P-V-T dataat 86.3 and 97.0 OK are derived from the empirical equation of state given in ref.(15) ; those for 108.2, 136.0 and 14802°K by interpolation of the data of ref. (16).All calculations are based on samples with N = 108 and on the LJ-parametersElk = 119*76”K, 0 = 3-405&13 with rc = 2-7A (cf. eqn. (7)). Not less than 168,000configurations were used to obtain the data for a given u,T-point. The mutualconsistency, the broad agreement with experiment and the magnitude of the statisticalfluctuations of the calculated mean values indicate that MC sequences of this lengt44 THERMODYNAMIC PROPERTIES OF LIQUID ARGONusually yield reliable statistical data.The standard deviations (calculated fromsub-averages over 4000 configurations) are m0-3 x 10-15 ergs for (<DIN) and =2 to4 x lO-.15 ergs for ('FIN). The relative error in the virial is therefore largest whenT"K86.397.0a97*Ob108.2c136*Od148-2av <;>cm3 ergx 101424-78 11-9025.60 11.512644 11.1526.90 10.9327.50 10.5627.50 10.3028.03 10.1228.48 9.9528.95 9.7629.68 9.5030.92 9-1428.48 9.7529.68 9.3627-50 10.1428-03 9.9728.48 9.8128-95 9.6229.68 9.3828.48 9.5929.68 9-2433.51 7.9736.58 7.3340.00 6.7233.51 7.8735-00 7.5836.58 7.2538.19 6.9540.00 6.64TABLE 1 .-THERMODYNAMIC PROPERTIES OF ARGON-1 avcm2 dyne-1 dyne cni-2 deg-1F(B)T (.E)V=(E)TY - c;x 10-7-ui 12 R $P('> dyne cm-2 x 10-8 xefop14 MC a p t .MC expt. MC expt. ~ 1 0 1 0MC expt. MC expt.1.64 1.57 9.98 0.91 1.60 1 *744.12 -0.43 9.66 0.96 1.66 1.82668 -2.36 9.36 1.27 1.71 2.447.32 -2.79 9.17 1 *236.36 -2.03 8.85 1 *280.66 2-13 2.10 864* 1*35* 1.87 1.73 2.30 1.942.51 0.76 0.92 8.49 1.03 2.17 1.91 1.73 2.073-34 0.17 0.09 8.35 1.09 2.34 2-06 1-77 2.154.18 -0.42 -0.67 8.19 1.02 2.50 2.19 1.72 2.225.55 - 1.34 7-97 0-85 2-76 1 *407.58 -2.60 0.68 3.10 1.05-0.14 2.93 2.16 7-28* 1*00*2.91 0.81 0.35 6.99 0.9 1-2.21 4.55 4.24 7.57* 0*88* 1.88 1.94 1.52 1.210.01 2.87 3.01 7.44 0.93 2.18 2.14 1.47 1.390.85 2.23 2.16 7.32 0-83 2.31 2.29 1.42 1-501.60 1.68 1-39 7.18 0.91 2.45 2.43 1.50 1.593.43 0.40 0.35 7.01 0.71 2-71 2.64 1-18 1.69-2.87 5.61 6-41* 1*13*0.78 2.54 229 6.18 0.781.42 2.52 2.37 4.24 4.15 0.56 0.58 4.73 4.53 0.83 1.214.20 0.79 0.95 3.90 3.82 0.47 0.57 6.70 8.30 0.60 0.975.88 -0.12 0.28 3.58 0.46 0.58 17-31 23.50 0.56 0.61-0.35 3.89 3.85 3.75 0-61 0.61 2.89 0.881.68 2.49 2.79 3.69 3.60 4.27 4.722.80 1.93 2.01 3.54 3.46 0.50 0.57 6-50 6.58 0.64 0.893.38 1.45 1.52 3.39 3.32 9.15 9.664.54 0.80 1.13 3.25 3.18 0.44 0.57 15-42 16.73 0.54 0.65MC = present work ; expt = experimental results ; * see text.deg-1 x 103MC expt.2.783.024.17430 3.362.75 3.954-14 4.434.30 4863.863.322.86 2.353-20 2.973.28 3.443-68 3.863.20 4.463.91 5.264.03 8.089.76 14.292.554-16 5.868.25 10.87direct MC calculation ; reweighted from 86.3"K data ; reweighted from 97.O"K data ; reweighted from 148~2°K data.( Y / N ) is small, i.e., when PVMRT and the relative error in P is largest when P%O( ( Y / N ) M 3kT).Estimated probable statistical errors are : 6U+* 5 cal/mole and6P M 2 x 107 dyne cm-2. The isochoric reweighting gives reliable results for variationsof up to 15 %in T.PRESSUREThe variation of PV/RT is shown in the figure ; with the exception of one point(28.48 cm3, 97-0"K) the agreement with experiment may be regarded as satisfactorythroughout the range of stability of the liquid. At 86-3"K the calculated pressuresform distinct "solid" and "liquid" branches. Although the pressure calculated forV = 27.50 cm3 from an initial face-centred cubic arrangement remained much belowthat calculated from an initial disordered liquid-type configuration, the pressure ofthe "solid" is anomalous, being higher than that calculated for Y = 26.90 cm3.This arises almost certainly from the presence of liquid-type configurations or regionsin the MC sample, indicating a fluctuation between two pha~es.4~10 At 148.2"Kand Y = 40.00 cm3 the calculated pressure falls well below the experimental valueI .R . MCDONALD AND K . SINGER 45This is interpreted as the onset of a van der Waals’ loop, a view supported by moreextensive computations 18 in this region with N = 32. The values of P at 9700°Kobtained by isochoric reweighting of the 86.3”K-data are in particularly good agree-ment with experiment.The reason for the discrepancy between these results and1.51.00 5FY w s4 0-- 0 5- 1.0-----35I 1 1 1 1 1 1 I I I24 26 28 30 32 34 36 38 4 0molar volume cm3FIG. 1 .-D- V isotherms for argon. 0, data from direct MC calculation ; 0 , data from isochoricreweighting. Full curves : expt. results ; dashed curve : sketched through MC points. Curve 1 :86~3°K (“solid” branch) ; curve 2 : 86.3”K (“liquid” branch, expt. equation of state extrapolated tonegative pressures) ; curve 3 : 97.O”K ; curve 4 : 136-0°K ; curve 5 : 148.2”K.the value obtained directly from MC runs for 28.48 cm3 is not known. Inspectionof the statistical data reveals signs of long-term fluctuations which suggests thatsubstantially longer sequences than 170,000 configurations are occasionally required.INTERNAL ENERGY AND LATENT HEATSThe agreement between observed 14 and calculated values of Ut at 148.2 andComparison with experimental data for the lower temperaturesThe latent heat AH,,, of vaporization at 86~3°K can be136~0°K is good.is made in table 2.TABLE 2.-pROPERTIES OF ORTHOBARIC LIQUID ARGONT P V - Uf/RT C p OK dyne cm-2x 10-6 cm3calc.expt. calc. expt. calc. expt.86.3 0.92 28.62 28.52 8.34 8-19 1 a07 0.8597.0 2-59 29.91 29-95 6.97 6.86 0.65 0.78108.2 5.97 31.13 31-71 5-92 5.69 0.35(?) 0.65calc = graphical interpolation of MC data of table 1 ; expt = linear interpolation of experimentalresults.21calculated from the experimental vapour pressure 20 (0.91 x 106 dyne cm-2) and theMC data.This gives a value of 1543 cal/mole compared with the measured value 20of 1564 at 857°K. The discontinuity of the pressure at 86.3”K and I/ = 27.50 cm46 THERMODYNAMIC PROPERTIES OF LIQUID ARGONindicates a phase transition in the MC model ; in a real system the phase equilibriumis characterized by the equality of the Gibbs free energy and of P at different densities.Properties relating to the liquid-solid transition calculated from MC data and theTABLE 3.-LATENT HEAT OF FUSION AND VOLUMES OF MELTING SOLIDAND FREEZING LIQUID FOR ARGON AT T = 86~3°Kmelting pressure 19 1.03 x 108 dyne cm-2Monte Carlo expt.Vsol (cm3) 25.0 24,619Kiq (c&) 27.9 27.919AHfu, (cal/mole) 246 284 (T = 84"K)zoexperimental melting pressure 19 are shown in table 3.excellent agreement with experiment but VsOl is markedly too high.AHf,, is-by about 14 %-lower than the observed value.20AHf,, are insensitive to the choice of pressure.The calculated VIis is inAs a consequenceBoth AHvap andHEAT CAPACITYThe computed and observed 14 values of Ct are in good agreement at 136.0 and148.2"K.The agreement, however, becomes significantly poorer as the volumeincreases, i.e., as the critical density is approached. The MC method is expectedto underestimate near the critical point because the periodic boundary conditionsuppresses all large density fluctuations. At 86~3°K the MC values for Ci/R showa clear maximum in the neighbourhood of the pressure discontinuity (V = 27.50 ~1113).Experimental and calculated values for the orthobaric liquid are compared in table 2.Ref.(16) tabulates Cy as a function of P at 90, 100 and 110°K but the data show asmuch scatter as those given for 97.0 and 108.2"K in table 1 and interpolation isdifficult. The values of Ct/R measured at 100°K between P = 4 x 106 and 300 x 106dyne cm-2 lie between 1-09 and 1.78, suggesting that the MC results are too low.OTHER PROPERTIESThe compressibilities (column 10 of table 1) are calculated by least-squares fittingof a parabola to the MC P-V data for each temperature. The thermal pressurecoefficients (column 1 1) are determined by isochoric reweighting (eqn. (6)), AT = 1°K).The expansivity (column 12) is the negative product of the two quantities. The goodagreement at 97-0°K is probably fortuitous.Elsewhere the agreement is mostly fair.ISOTHERMAL EXTRAPOLATIONValues of 2 have not yet been reliably determined over a range of densities. Theresults obtained for V = 28.03 and 28- 95 cm3 from the MC data for V = 28.48 cm3and T = 86.3"K by means of eqn. (15) and (16) are shown in table 4. The differencesbetween the extrapolated and the independently calculated values of (@/iV) and(Y/N} are smaller than the statistical error of the MC calculation. For moderatelysinall extrapolations of (<DIN) (AV/Y-2 %) values of ," are not required. This isalso true for small extrapolations of (!PIN) (AV/V-0-5 %) when the second term ineqn. (16) can be neglected. Extrapolation from 28.48 to 28-38 cm3 yielded a com-pressibility of 2.24 x 10-10 cm2 dyne-1 compared with 2.34 x 10-10 cm2 dyne-1obtained from the MC isotherm. These results are regarded as encouragingI .R . MCDONALD A N D K . SlNGER 47TABLE d.-ISOTHERMAL EXPRAPOLATION OF MONTE CARL0 DATA FORV = 28-48 cm3, T = 86~3°Kdirect calculation extrapolationV <;> <;> @ (;)6cm3) (erg x 1015) (ergx 1015)28.03 - 10.12 2.5 1 - 10.08 2-6428.95 - 9.76 4-18 - 9.81 4.34CONCLUSIONSCalculations based on an MC method and an LJ-potential have been shown toyield quantitatively significant values of thermodynamic propel ties of a simple liquid(argon) over a considerable V,T range, including phase transitions. Satisfactoryresults have been obtained for the internal energy. For the pressure, the agreementwith experimental values is less good but still, on the whole, satisfactory.For thederivatives of UP and P the agreement with experiment is better than qualitative.In a number of cases the discrepancies between calculated and experimentalresults are larger than the statistical sampling errors inherent in the MC method.It is hoped that future work will show to what extent the discrepancies are due to theapproximations of the MC method, e.g., small sample size and the periodic boundarycondition, or to inadequacy of the potential function.By calculating not merely the mean values but also the frequency distributions forthe potential energy and the virial, it has been possible greatly to extend the T-rangecovered by a given number of MC calculations. Further economies in computationbecome possible by means of an isothermal extrapolation procedure for which onlylimited results have so far been obtained.This method will also permit the systematicinvestigation of the effects of changes of the LJ-parameters on the thermodynamicproperties.The computing time required to obtain results of the type reported in this paper islarge but not prohibitive. MC calculations for more complex systems, e.g., mixturesof simple fluids or systems governed by non-central forces, need no longer be regardedas impractical.APPENDIXAPPROXIMATIONSEFFECT OF DISTANT PARTICLES. TO reduce computing labour the sum Of all pairenergies is replaced by one in which the potential due to a uniform density is substitutedfor the effect of the distant particles,4 i.e.,@ = 4(rij)+(A7/2u)[m +(r)4nr2dr = @*(r,).In this work r , = 6.5 A.MC runs yield ensemble averages of the “effective”potential energy ( W ( r , = 6.5 A)). These values are further corrected by computingthe values of @*(r, = 6-5 A) and 4i)*(rz = 8-4 A) for a small number of configurationsat each density and adding the mean difference of (@*(rm = 8.4 A) - @*(rm = 6.5 A)]to (@*(rm = 6.5A)) in all cases. The same procedure is applied to the virials. Forthe potential energy the “effective” long-range potential amounts to 15-20 % of thetotal; for the virial, the long-range contribution is of opposite sign and of the sameorder of magnitude as the exactly calculated short range contribution. The “effective”rij <rm r48 THERMODYNAMIC PROPERTIES OF LIQUID ARGONlong range contributions and their corrections vary smoothly with density. It isestimated that the error resulting from this procedure does not exceed 0.5 % for 4and 5 % for Y.THE SAMPLE SIZE.In all calculations N = 108. The mean values of all macro-scopic equilibrium properties are assumed to be determined by the mean values ofthe MC sample. This involves errors because the periodic boundary conditionexcludes configurations which in reality occur, and because density fluctuations aresuppressed. Density fluctuations in small volume elements can be calculated fromthe isothermal compressibility. For N = 108 such fluctuations would not contributesignificantly to the macroscopic entropy if the microscopic subsystems are treated asindependent.This assumption is, I: ever, unjustified because adjacent subsystemsshare a large number of pair interaciiwns which must be expected to lead to a correla-tion of their fluctuations. No analysis of these errors has been carried out so far.Properties calculated for some u,T-points with samples of size N = 32 and N = 108are in good agreement.18THE COMPUTER PROGRAMMEA large part of the computation is done by means of integer arithmetic (i.e.,24-bit words). This reduces the execution times of transfers, additions and sub-tractions and makes an extensive use of tables possible. The sample cube has thelength L = 1000; the 3 x 108 coordinates are integers between 0 and 999. At theprevailing densities a hard core radius re = 2.7 A corresponds to > 100 units. Noserious distortions are likely to result from the discreteness of the variables. Changesof volume are taken into account by scaling of the LJ parameters. The periodicboundary condition requires that in the calculation of a pair distance rfj2, the com-ponents x ~ j , yij, zfj must be “modified”, i.e., if xi-.q>L/2, it is replaced byxg +L- xj, etc.To determine the interaction potential between particles i and j ,the computer calculates xi - xj, yi -J)], zi - zj from the listed particle coordinates ;it then obtains the squares of the “modified” differences from a table, adds these andlooks up the values of the potential and of the virial functions for the rig2 from tables.Lists are kept of all coordinates, of all pair potentials in the “present” configuration(Wfj) and of the contribution of each particle i to the “present” value of iD (Ai =+cW& ‘&4g = a).When a move is attempted by a displacement of particle k,a list is formed of the provisionally altered pair potentials (Vkj, all ji:k), and theprovisional new contribution of k to @ ( A = Eva,). The provisional new @ isaccepted or rejected according to the criterion of ref. (1) (the exponential functionsrequired for exp (-PA@) are tabulated). If the move is accepted, the lists of co-ordinates, pair interactions, particle contributions and the value of @ are updated(e.g., V k p Wkj, A+&) and an entry corresponding to the new value of @ is made inthe histogram. If the move is rejected an entry corresponding to the “present” valueof @ is made (again) in the histogram.The treatment of the virials by means oftables is analogous.The programme generates about 3000 configurations per minute on the Mark 1Atlas computer of the University of London. Calculations for a v,T-point werebased on not less than six ten-minute runs (168,000 configurations) of which the firstone was discarded to allow for “equilibration”. The formation of tables for thepotential and the virial functions at the beginning of each run does not add signifi-cantly to the computing time; calculations with any other type of pair potentialj4i ij 4 I . R . MCDONALD AND K . SINGER 49instead of the LJ-potential would not present any problem (though eqn.(13)-(16) forthe isothermal extrapolation would not be applicable).We thank the Institute of Computer Science (University of London) for theallocation of sufficient computer time, Mr. K. R. Knight for his contribution to thecomputer programme, and the Science Research Council for financial support. Weare grateful to Prof. J. S. Rowlinson for bringing to our attention some importantexperimental data.1 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J, Chem.2 M. N. Rosenbluth and A. W. Rosenbluth, J. Chem. Physics, 1954, 22, 881.3 W. W. Wood and J. D. Jacobson, J. Chem. Physics, 1957,27,1207.4 W. W. Wood and F. R. Parker, J. Chem. Physics, 1957, 27, 720.5 W. W. Wood, F. R. Parker and J. D. Jacobson, Nuuuo Cimento, supp. 1, vol. 9, series X, 1958,6 Z . L. Salsburg, J. D. Jacobson, W. Fickett and W. W. Wood, J. Chem. Physics, 1959, 30, 65.7 E. B. Smith and K. R. Lea, Trans. Faraudy Soc., 1963,59, 153.8 A. Rotenberg, J. Chem. Physics, 1965, 43, 1198.9 A. Rotenberg, J. Chem. Physics, 1965,43,4377.1oB. J. Alder, Physic. Rev. Letters, 1966, 16, 88.11 A. A. Barker, Austral. J. Physics, 1965, 18, 119.12 S. G. Brush, H. L. Sahlin and E. Teller, J. Chem. Physics, 1966, 45, 2102.13 A. Michels, Hub. Wijker and Hk. Wijker, Physica, 1949, 15, 627.14 A. Michels, J. M. Levelt and G. J. Wolkers, Physica, 1958, 24, 769.15 A. van Itterbeek and 0. Verbeke, Physicu, 1960,26,931.16 A. van Itterbeek, 0. Verbeke and K. Staes, Physica, 1963, 29, 742.17 W. W. Wood and J. D. Jacobson, Proc. Western Joint Computer Conference, (San Francisco,18 I. R. McDonald and K. Singer, to be published.19 F. Din, Thermodynamic Functions of Gases, vol. 2, (Butterworths, London, 1956), p. 146-201.20 P. Flubacher, A. J. Leadbetter and J. A. Morrison, Proc. Physic. SOC., 1961, 78, 1449.21 J. S. Rowlinson, Liquids and Liquid Mixtures, (Butterworths, London, 1959), p. 64.Physics, 1953, 21, 1087.133.1959), 261