Predictive power of effective intermolecular pair potentials: MD simulation results for methane up to 1000 MPa

Fluid Phase Equilibria, 57 (1990) 35-46 Elsevier Science Publishers B.V., Amsterdam -

35 Printed in The Netherlands

PREDICTIVE POWER OF EFFECTIVE INTERMOLECULAR PAIR POTENTIALS: MD SIMULATION RESULTS FOR METHANE UP TO 1000 MPa BERTHOLD SAAGER and JOHANN FISCHER Institut ftir Thermo- und Fluiddynamik, Ruhr-Universitiit, D-4630 Bochum 1 (F.R. G.) (Received May 24, 1989)

ABSTRACT Saager, B. and Fischer, J., 1990. Predictive power of effective intermolecular pair potentials: MD simulation results for methane up to 1000 MPa. Fluid Phase Equilibria, 57: 35-46. Pressures, internal energies and specific heats of methane obtained from molecular dynamics simulations on a CYBER 205 using a Lennard-Jones potential are reported. The potential parameters were determined several years ago with a WCA-type perturbation theory by fitting to the vapour pressure and the saturated liquid density at one temperature. The agreement of the calculated pressures with the experimental pressures is on a gas isochore within rt 0.1 MPa, on two liquid isochores within f 1 MPa and on a supercritical isotherm within f2% up to 1000 MPa. The specific heats in the liquid agree with the experimental values within the accuracy of the latter. Further results given at state points where no experimental data exist may serve as support for the development of equations of state.

INTRODUCTION

With the increasing availability of computer facilities there is some challenge to predict the thermodynamic properties as well as the transport properties of fluids by computer simulations. Recent examples for the calculation of thermodynamic properties in the whole fluid region are the molecular dynamics (MD) simulations for propane (Lustig and Steele, 1988) and ethane (Lustig et al., 1989), and Monte Carlo (MC) simulations for carbon dioxide (Luckas and Lucas, 1989). Crucial points in such projects are the statistical uncertainties of the results, which will be dealt with later in this paper, and perhaps even more important, the underlying intermolecular potentials. Different routes have been suggested for the construction of effective intermolecular model potentials. One may use solid state properties, as done by MacRury et al. (1976) and by Murthy et al. (1983). Another method uses 0378-3812/90/$03.50

0 1990 - Elsevier Science Publishers B.V.

36

zero pressure densities and internal energies of the liquid (McDonald and Singer, 1972; Singer et al., 1977; Gupta et al., 1988). A further development along this line was to use saturated liquid densities and vapour pressures as determined from WCA-type perturbation theories (PT). This method was used to determine parameters for Lennard-Jones (LJ) and two-centre LennardJones (20 potentials (Fischer et al., 1984; Bohn et al., 1986b) as well as for n-centre Lennard-Jones (nCLJ) potentials (Lustig, 1986; Lustig, 1987) and Kihara potentials (Boublik, 1987). As the m-based method in particular yields LJ, 2CLJ and nCLJ potential parameters for about 20 substances it is interesting to see how accurately they predict thermodynamic properties if used in computer simulations. So far it has been found that the thermic and caloric properties of ethane can be described by a 2CLJ model in the whole fluid region with good accuracy (Fischer et al., 1987; Lustig et al., 1989; Saager et al., 1990) using the parameters obtained from PT. For propane, the original PT parameters gave pressures which are too high at triple point densities. However, reducing the size parameter u by only 0.4% led again to good predictions of all thermodynamic properties in the whole fluid region (Lustig and Steele, 1988). For oxygen and ethylene, simulations were performed with PT parameters for densities about two thirds of the critical; there the calculated pressures agree with the experimental pressures within the accuracy of the simulations (Saager et al., 1990). Another fluid for which the predictive power of PT parameters can be tested is methane, for which experimental high pressure puT data for the isotherm T = 298.15 K were published (Kortbeek et al., 1986; Kortbeek and Schouten, 1989) only after the presentation of our potential model. This lucky incident can be used for a fair test as it excludes the fact that the parameters have been tuned with respect to properties to be predicted. MD simulations for methane have already been performed by other groups. The first paper seems to be that of HanIey and Watts (1975b) who used an m-6-8 potential with coefficients determined from viscosity data for the dilute gas (Hanley and Watts; 1975a). The agreement between the calculated pressures and the experimental pressures was within 2 MPa close to the critical density and within 6 MPa at higher densities. Later Murad et al. (1979) used a five-centre expd potential suggested by Williams (1967). with that model the agreement between the calculated pressures and the experimental pressures was within 6 MPa close to the critical density and within 13 MPa at about 2.5 times the critical density. Murad and Gubbins (1978) also investigated a five-centre Lennard-Jones (SCLJ) potential. They determined the parameters by fitting the pressures, internal energies and specific heats in the density range 10-25.27 mol dme3 and the temperature range 123-276 K to experimental data. It transpired that the thermody-

37

namic quantities simulated with that potential showed equivalent accuracy to those found with the five-centre expd potential of Williams (1967). The 5CLJ potential then was used by van Waveren et al. (1986) to perform calculations on the isotherm T = 298.15 K at high densities. From Fig. 1 in the van Waveren paper one learns that at p - 32 mol dmm3 the calculated pressure is too high by about 60%; that density, however, amounts to only 90% of the highest experimental density. As a consequence, those authors resealed all size parameters by about 4% and then obtained good agreement with the experimental pressures at the high densities. After that resealing, however, the results at lower densities must become considerably worse. In the present paper we simply use a Lennard-Jones potential with the parameters (Fischer et al., 1984) c/k = 149.92 K, u = 3.7327 A, which were determined by fitting the saturated liquid density and the vapour pressure obtained from our WCA-type perturbation theory (Fischer, 1980) to the experimental values at one temperature (T = 126 K). Using the same set of parameters, several other thermodynamic properties have already been obtained using different tools of statistical thermodynamics: vapour pressures, bubble densities and dew densities from perturbation theory (Fischer et al., 1984), the second virial coefficient from the rigorous equation (Bohn et al., 1986a), isochoric heat capacities in the liquid from perturbation theory (Fischer et al., 1987) and pvT data in the gas phase from the HaarShenker-Kohler equation (Saager et al., 1989). At this point we should mention Jhat our parameters are rather similar to one set (e/k = 149.1 K, (T= 3.743 A) determined by McDonald and Singer (1972) from zero pressure densities and internal energies: One aim of the present paper is to test our LJ model by comparison with the experimental pressures on the isotherm T = 298.15 K from the critical density, pc, up to 3.5 p,. For densities lower than 2.5 p, comparison can be made with the data of Trappeniers et al. (1979); for the higher densities the measurements of Kortbeek et al. (1986) and their reevaluation (Kortbeek and Schouten, 1989) are available. However, in order to avoid the impression that this simple potential model works only for the pressure on this high temperature isotherm we also perform simulations for the gas as well as for the low temperature liquid. On the gas side we simulate close to 0.7 pc from 190 to 600 K and compare with the experimental pressure results of Kleinrahm et al. (1986) and of Douslin et al. (1964). On the liquid side we simulate close to the isochores 2.1 and 2.7 pc and compare with the experimental pressures of Goodwin and Prydz (1972). Moreover, we determine isochoric heat capacities close to these two isochores as well as for the density 2.8 p, which are compared to the experimental data of Younglove (1974). We should still mention that we do not compare internal energies, as these are not directly accessible experimental quantities and at

the moment it is not yet quite clear how well the results from the most recent correlation equations given by Friend et al. (1989) and by Setzmann (1989) agree. We do however, quote our results. In the next section we give details of the simulations and also discuss their statistical uncertainties. In the subsequent section the results will be presented. As the agreement between the calculated values and the experimental values is remarkably good, we also present results for other state points for which no experimental results are available. We believe that these data may serve as support for the development of equations of state. METHODOLOGY

OF SIMULATIONS

We have performed MD simulations in an NVT ensemble using a predictor-corrector algorithm of fifth order (Gear, 1971). The program is originally due to Haile (1980) and was vectorized by one of us (B.S.) to run on the CYBER 205. Periodic boundary conditions and the minimum image criterion were used. All runs were performed with 256 particles and the cut-off radius was taken to be half of the box length; long range corrections were taken into account. The temperature was kept constant by momentum scaling. The time step for the integration was chosen to be either 0.0030 or 0.0006 in units of a(m/c) ‘I2. The very short latter time step was chosen for the isotherm T = 298.15 K and the additional state points at high densities and high temperatures as a precautionary measure. The system was started either from an f.c.c. lattice or a previous configuration followed by an equilibration period over 50 000 time steps for the short time step and 10 000 time steps for the long time step. The length of the production runs varied from 60 000 to 150 000 time steps and is given together with the results. The quantities that were calculated directly in the simulations are the pressure and the residual part of the internal energy.

time steps. 10" -

Fig. 1. Running average of the pressure at the state point p = 15 mol drne3 ( po3 = O&9798), T = 298.15 K. This is an example for small pressure fluctuations.

39 26.5 Y26.4

-

26.3 26.2 26.1 2 6.0

0

20

40

60

80

100

time steps .111T3 -

Fig. 2. Running average of the pressure at the state point p = 35.563 mol dm- -3 1.113828), T = 298.15 K. This is an example for large pressure fluctuations.

- 3.1

0

20

40

60 time steps. 10e3 -

80

(pa3=

100

Fig. 3. Running averageof the residualinternal energy at the state point p =15 mol dme3 ( pu3 = O&9798),T = 298.15 K. This is an example for small energy fluctuations.

A major point in calculating thermodynamic properties from simulations is the statistical uncertainties which have to be considered carefully. Hence, we made for each run a plot of the running average of the pressure and the internal energy and estimated from this the statistical uncertainties which will be given, together with the results, in the next section. Examples for the

40

60 time steps. 11T3 -

80

100

Fig. 4. Running average of the residual internal energy at the state point p = 35.563 mol dmP3 ( po3 = 1.113828), T = 298.15 K. This is an example for large energy fluctuations.

fluctuations in the pressure as well as in the internal energy are shown in Figs. 1-4 for the highest density and a comparatively low density. Recording the running average also helped us to detect whether the system tended to freeze at a given state point. RESULTS

The results for the pressures and the residual internal energies, together with their statistical uncertainties, are given in Table 1 in reduced units for all simulation runs performed. Besides the runs from which the results are needed for the following comparison with the experimental data, some simulations were also made in the high density region where no experimental data exist. We should mention that we also tried to perform simulations at the state point T* = kT/c = 1.734258 and p* = pa3 = 1.142079. If starting from the f.c.c. lattice the system remained ordered. If the initial configuration was taken from a disordered state the internal energy remained constant for a certain period but then drifted towards lower values. Therefore, we conclude that this state point is definitely beyond the freezing line. Comparison of the simulated pressures with the experimental values is given in Table 2. At T = 298.15 K experimental results are available up to 25 mol dme3 from Trappeniers et al. (1979) and for higher densities from Kortbeek et al. (1986) which were revisited by Kortbeek and Schouten (1989). The revisited data show at a given pressure a slightly lower density which corresponds to an increase of the pressure of about 0.5% at a fixed density. It is seen that at the lowest density the simulated pressures are too low by about 1.5% whilst at the highest density they are too high by about 2%. We should remark that this agreement is within the estimated accuracy of the simulations given in Table 1 except at the three highest densities. For the gas, the isochore p = 7 mol dms3 was chosen because results from the Haar-Shenker-Kohler equation were given for that density (Saager et al., 1990) and experimental results are available from Douslin et al. (1964). An additional point closer to the dew line was taken from Kleinrahm et al. (1986). We see that the agreement between the simulated values and the experimental values is always within i-0.1 MPa and within the accuracy of the simulations. Finally, pressures are also compared for two liquid isochores, one at about 21 mol dmH3 and the other at about 27 mol dm-3 which is nearly the triple point density. We learn that the agreement is now within _+1 MPa and again within the accuracy of the simulations. The isochoric heat capacities presented in Table 3 were obtained from four to five results of the internal energy on one isochore which were fitted

41 TABLE 1 Results from MD simulations: residual internal energy, pressure and their estimated uncertainties for methane as LJ fluid. The last columns contain the time step At* in units of a( “/‘)r’2 and the number of production time steps T*

P*

U*

1.988727

2.001067 2.668090 3.335112 4.002134

1.113828 1.022406 0.933734 0.858164 0.812312 0.626397 0.469798 0.313199 0.839476 0.907806 0.985896 1.112795 1.073748 1.034702 0.995658

1.271045 1.822038 2.989574 3.990695 1.053895 1.213981 1.374066 0.680363 0.747065 0.840448 1.024840 1.049840 1.099840 1.124840 0.633191 0.658191 0.708191 0.733191 0.567609 0.592609 0.642609 0.642609 0.667609

0.214732 0.219239 0.219239 0.219239 0.669337 0.662666 0.658406 0.852902 0.847515 0.842348 0.667520 0.667520 0.667520 0.667520 0.857788 0.857788 0.857788 0.857788 0.882093 0.882093 0.882093 0.882093 0.882093

- 4.61 - 5.011 - 5.065 - 4.948 - 4.795 - 3.9013 - 2.961 - 2.0018 - 5.1354 -5.344 - 5.4246 - 4.605 - 3.962 - 3.451 - 3.038 - 3.026 - 1.644 - 1.456 - 1.2861 - 1.190 - 4.643 - 4.502 - 4.388 - 6.185 - 6.073 - 5.941 - 4.649 - 4.634 - 4.602 - 4.591 - 6.271 - 6.240 - 6.183 - 6.154 - 6.513 - 6.479 - 6.411 - 6.408 - 6.376

1.688567

P* 26.1 16.89 10.86 7.21 5.65 2.00 0.934 0.506 5.185 7.76 12.036 26.05 26.93 27.07 26.71 26.77 0.114 0.300 0.679 0.995 0.028 0.456 0.840 0.024 0.386 0.857 - 0.060 0.004 0.160 0.224 - 0.211 - 0.031 0.295 0.462 -0.311 -0.111 0.274 0.293 0.475

Au*

Ap*

At*

Production steps

0.02 0.010 0.015 0.015 0.015 0.0075 0.010 0.0075 0.0075 0.010 0.0075 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.0075 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.008 0.008 0.008 0.008 0.008

0.1 0.05 0.07 0.07 0.07 0.02 0.015 0.005 0.035 0.04 0.035 0.06 0.06 0.06 0.06 0.06 0.003 0.003 0.003 0.004 0.02 0.02 0.025 0.02 0.025 0.035 0.02 0.02 0.02 0.02 0.025 0.025 0.025 0.025 0.03 0.03 0.03 0.03 0.03

0.0006

100000 150000 100000 1cOOoo 100000 150000 100000 150000 150000 100000 150000 100000 100000 1OOooO lOOQO0 100000 60000 70000 70000 70000 7004xJ 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000 70000

0.0006

0.0030

0.0030

0.0030

0.0030

0.0030

0.0030

42 TABLE 2 Comparison of pressures from MD simulations with experimental results for methane

298.15

190.555 273.160 448.197 598.285 158 182 206 102 112 126

;)mol dm-‘)

PMD (MW

Pexp (MW

PCCp OfPa)

35.563 32.644 29.813 27.400 25.936 20.0 15.0 10.0 6.8561 7 7 7 21.371 21.158 21.022 27.232 27.060 26.895

1038.8 672.2 432.2 287.0 224.9 79.60 37.17 20.14 4.55 11.95 27.02 39.61 1.11 18.14 33.42 0.94 15.37 34.12

1013.04 a) 660.01 a) 423.55 a) 285.01 a) 222.27 a) 79.62 ‘) 37.44 =) 20.48 ‘) 4.54 d, 12.05 e, 27.12 e, 39.71@ 1.98 o 18.12 0 33.80 n 0.85 n 14.29 r, 33.34 rJ

1017.27 663.49 425.47 285.28 222.14

b, b, b, b, b,

a>Kortbeek et al., 1986; b, Kortbeek and Schouten, 1989; ‘) Trappeniers et al., 1979; d, Klein&m et al., 1986; e, Douslin et al., 1964; o Goodwin and Prydz, 1972.

to a least-square parabola. This table also contains previously obtained results from perturbation theory (Fischer et al., 1987). We see that the agreement between simulations and perturbation theory is quite satisfying and that the agreement of both with the experimental values of the total heat capacity is within 3% which is believed to be the experimental accuracy.

TABLE 3 Comparison of isochoric heat capacities from MD simulations and perturbation theory (PT) (Fischer et al., 1987) with the experimental results for methane (Younglove, 1974). The residual quantities cr and total values c, are given T (K) 161.140 102.424 92.592

P Cm01

C:es/R

cy/R

Cies/R

dm-‘)

MD

PT

exp

MD

21.313 27.388 28.164

0.65 1.16 1.37

0.54 1.18 1.32

0.66 1.14 1.25

3.66 4.16 4.37

c,/R

G/R

PT 3.55 4.18 4.32

C”/R

exp 3.67 4.14 4.25

43 CONCLUSIONS

We have found that the potential model for methane suggested by us some years ago (Fischer et al., 1984) predicts pressures which agree with the experimental values within the accuracy of the simulations except at very high densities. There the accuracy is within 2% up to pressures of 1000 MPa. Isochoric heat capacities of the dense liquid are predicted within the experimental accuracy. Comparing these results with those of the previous studies of Hanley and Watts (1975b) and of Murad et al. (1979) we learn that our potential model yields considerably better predictions of the thermodynamic properties. We remind that van Waveren et al. (1986) have shown the five-centre LennardJones potential of Murad and Gubbins (1978) to overpredict the pressure on the isotherm T = 298.15 K by about 60% at 90% of the highest density. The present study confirms the previous experience with ethane, oxygen, ethylene and propane (Fischer et al., 1987; Lustig and Steele, 1988; Lustig et al., 1989; Saager et al., 1990) that Lennard-Jones and n-centre Lennard-Jones potentials yield a good description of the thermodynamic properties of non-polar fluids if their parameters have been determined by fitting the vapour pressures and saturated liquid densities obtained from perturbation theory to experimental data. Concerning the m-6-8 potential of Hanley and Watts (1975b) we believe that it yields worse results for the thermodynamic properties than our LJ potential, mainly because its parameters were fitted to viscosity data of the low density gas. The problem with the five-centre potentials used by Murad and Gubbins (1978) and Murad et al. (1979) is presumably the fact that such potentials need seven parameters, and even if for the interaction between the C and H atoms the Lorentz-Berthelot rules are assumed, five parameters are left which are hard to determine properly. As it is now clear that methane can be described with good accuracy as a Lennard-Jones fluid and ethane as a two-centre Lennard-Jones fluid, it is also very likely that their mixture can be described on the basis of these potential models. As far as the unlike interaction is concerned it was already shown (Bohn et al., 1986a) that the Kohler new extended (KNE) combining rule yields good predictions for the liquid mixture excess properties without using any experimental mixture data. ACKNOWLEDGEMENTS

The authors thank Dr. P.J. Kortbeek and Prof. J.A. Schouten from the van der Waals Laboratory, Amsterdam, for having sent them their paper with the revised experimental data prior to publication.

44

Die Autoren danken Prof. W. Wagner (Bochum) fur die Anregung zu dieser Arbeit und Prof. F. Kohler (Bochum) ftir fruchtbare Diskussionen. Die Arbeit wurde durch die Deutsche Forschungsgemeinschaft im Rahmen des Schwerpunktprogramms Neue Arbeitsmedien in der Energie- und Verfahrenstechnik, AZ Fi 287/5, gefiirdert. LIST OF SYMBOLS AND CONSTANTS

Roman letters C” res C” m P P * =p d/c T T* =kT/c u U * = u/N<

isochoric heat capacity residual isochoric heat capacity molecular mass pressure reduced pressure temperature reduced temperature residual internal energy residual internal energy in reduced units

Greek letters AP AP*

simulation uncertainty for pressure simulation uncertainty for pressure in reduced units

At” = At/

bh/W21

AU

Au* f p p* = pa 3 0

time step in reduced units simulation uncertainty for residual internal energy simulation uncertainty for residual internal energy in reduced units energy parameter for LJ potential density reduced density size parameter for LJ potential

Constants k = 1.380658 X 1O-23 J K-‘Boltzmann’s constant R = 8.31434 J mol-’ K-’ universal gas constant

45

REFERENCES Bohn, M., Fischer, J. and Kohler, F., 1986a. Prediction of excess properties for liquid mixtures: results from perturbation theory for mixtures with linear molecules. Fluid Phase Equilibria, 31: 233-252. Bohn, M., Lustig, R. and Fischer, J., 1986b. Description of polyatomic real substances by two-center Leonard-Jones model fluids. Fluid Phase Equilibria, 25: 251-262. Boublik, T., 1987. Simple perturbation method for convex molecule fluids. J. Chem. Phys., 87: 1751-1756. Douslin, D.R., Harrison, R.H., Moore, R.T. and McCullough, J.P., 1964. Pu T relations for methane. J. Chem. Eng. Data, 9: 358-363. Fischer, J., 1980. Perturbation theory for the free energy of two-center Lennard-Jones liquids. J. Chem. Phys., 72: 5371-5377. Fischer, J., Lustig, R., Breitenfelder-Manske, H. and Lemming, W., 1984. Influence of intermolecular potential parameters on orthobaric properties of fluids consisting of spherical and linear molecules. Mol. Phys., 52: 485-497. Fischer, J., Saager, B., Bohn, M. and Haile, J.M., 1987. Specific heat of simple liquids. Mol. Phys., 62: 1175-1185. Friend, D.G., Ely, J.F. and Ingham, H., 1989. ‘I’hermophysical properties of methane. J. Phys. Chem. Ref. Data, 18: 583-638. Gear, C.W., 1971. Numerical initial value problems in ordinary differential equations, Prentice Hall, Englewood Cliffs, NJ. Goodwin, R.D. and Prydz, R., 1972. Densities of compressed liquid methane, and the equation of state. J. Res. NBS, 76A: 81-101. Gupta, S., Sediawan, W.B. and McLaughlin, E., 1988. Computer simulation of benzene using the modified Gaussian overlap and six-site potentials. Mol. Phys., 65: 961-975. Haile, J.M., 1980. A primer on the computer simulation of atomic fluids by molecular dynamics, private communication. Hanley, H.J.M. and Watts, R.O., 1975a. The self-diffusion coefficient of liquid methane. Mol. Phys., 29: 1907-1917. Hanley, H.J.M. and Watts, R.O., 1975b. Molecular dynamics calculation of the thermodynamic properties of methane. Aust. J. Phys., 28: 315-324. Kleinrahm, R., Duschek, W. and Wagner, W., 1986. Pressure, density, temperature measurements in the critical region of methane. J. Chem. Thermodyn., 18: 1103-1114. Kortbeek, P.J., Biswas, S.N. and Trappeniers, N.J., 1986. poT and sound velocity measurements for CH, up to 10 kbar. Physica B, 139/140: 109-112. Kortbeek, P.J. and Schouten, J.A., 1989. Compressibility and sound velocity measurements in CH, up to 1 GPa, revisited. Private communication. Luckas, M. and Lucas, K., 1989. Thermodynamic properties of fluid carbon dioxide from the SSR-MPA potential. Fluid Phase Equilibria, 45: 7-23. Lustig, R., 1986. A thermodynamic perturbation theory for non-linear multicentre LennardJones molecules with an anisotropic reference system. Mol. Phys., 59: 173-194. Lustig, R., 1987. Application of thermodynamic perturbation theory to multicentre LennardJones molecules. Results for CF,, Ccl,, neo-C,H,, and SF, as tetrahedral and octahedral models. Fluid Phase Equilibria, 32: 117-137. Lustig, R. and Steele, W.A., 1988. On the thermodynamics of liquid propane. A molecular dynamics study. Mol. Phys., 65: 475-486. Lustig, R., Toro-Labbe, A. and Steele, W.A., 1989. A molecular dynamics study of the thermodynamics of liquid ethane. Fluid Phase Equilibria, 48: I-10.

46 MacRury, T.B., Steele, W.A. and Beme, B.J., 1976. Intermolecular potential models for anisotropic molecules with applications to Nz, CO, and benzene, J. Chem. Phys., 64: 1288-1299. McDonald, I.R. and Singer, K., 1972. An equation of state for simple liquids. Mol. Phys., 23: 29-40. Murad, S. and Gubbins, K.E., 1978. ACS Symp. Ser., 86: 62. Murad, S., Evans, D.J., G&bins, K.E., Street, W.B. and Tildesley, D.J., 1979. Molecular dynamics simulation of dense fluid methane. Mol. Phys., 37: 725-736. Murthy, C.S., Singer, K. and McDonald, I.R., 1983. Modeling of simple nonpolar molecules for condensed phase simulations. In: J.M. Haile and G.A. Mansoori (Eds.), Molecular based study of fluids, Advances in Chemistry Series 204, American Chemical Society, Washington, DC, p. 189. Saager, B., Lotfi, A., Bohn, M., Nguyen Van Nhu and Fischer, J., 1990. Prediction of gas PvT data with effective intermolecular potentials using the Haar-Shenker-Kohler equation and computer simulations. Fluid Phase Equilibria, 54: 237-246. Setzmann, U., 1989. Dr.-Ing. thesis supervised by W. Wagner, Ruhr-Universitat Bochum, F.R.G. Singer, K., Taylor, A. and Singer, J.V.L., 1977. Thermodynamic and structural properties of liquids modelled by ‘2-Lennard-Jones centres’ pair-potentials. Mol. Phys., 33: 1757-1795. Trappeniers, N.J., Wassenaar, T. and Abels, J.C., 1979. Isotherms and thermodynamic properties of methane at temperatures between 0 ’ C and 150 ‘C and at densities up to 570 amagat. Physica A, 98: 289-297. van Waveren, G.M., Michels, J.P.J. and Trappeniers, N.J., 1986. Molecular dynamics simulation of CH, in the dense fluid phase. Physica B, 139/140: 144-147. Williams, D.E., 1967. Nonbonded potential parameters derived from crystalline hydrocarbons. J. Chem. Phys., 47: 4680-4684. Younglove, B.A., 1974. The specific heats cp and c, of compressed and liquified methane. J. Res. NBS, 78A: 401-410.