ResearchGate See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/233597508 Inadequacy of the Lorentz-Berthelot Combining Rules for Accurate Predictions of Equilibrium Properties by Molecular Simulation Article in Molecular Physics ■ April 2001 DO 1:10.1080/00268970010020041 CITATIONS READS 70 192 2 authors, including: /*S v Jerome Delhommelle University of North Dakota 96 PUBLICATIONS 1,746 CITATIONS SEE PROFILE All in-text references underlined in blue are linked to publications on ResearchGate, lettingyou access and read them immediately. Available from: Jerome Delhommelle Retrieved on: 04 August 2016 Molecular Physics, 2001, Vol. 99, No 8, 619-625 Inadequacy of the Lorentz-Berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation JEROME DELHOMMELLE1-2* and PHILIPPE MILLIE2 1 Research School of Chemistry, Australian National University, Canberra ACT 0200, Australia 2 DSM/DRECAM/SPAM (Laboratoire Francis Perrin), CEA-CE Saclay, F-91191 Gif-sur-Yvette, France (Received 28 September 2000; accepted 8 November 2000) Though molecular beam experiments have revealed deficiencies in the Lorentz-Berthelot combining rules, these rules are still used widely to parametrize effective pair potential models or to calculate the thermodynamic properties of mixtures. Gibbs ensemble Monte Carlo and isothermal isobaric Monte Carlo simulations were used to compute the liquid-vapour phase equilibria and the liquid properties of binary mixtures of rare gases modelled by effective pair potentials. Three sets of simple combining rules were tested in this work: the commonly used Lorentz-Berthelot rules, the Kong rules (Kong, J., 1973, /. chem. Phys., 59, 2464) and the Waldman-Hagler rules (Waldman, M., and Hagler, A. T, 1993, /. comput. Chem., 14, 1077). These three sets of rules do not require any additional parameter. It is shown that: (1) the choice of a set of combining rules has a significant effect on the thermodynamic properties, (2) using the Lorentz-Berthelot rules yields significant deviations from experiment and (3) the Kong rules provide a much better description of the mixture properties both for coexistence curves and liquid properties We therefore recommend the use of the Kong rules instead of the Lorentz-Berthelot when parametrizing potential models. 1. Introduction The application of molecular simulation methods to an understanding of and to the prediction of properties of practical interest requires the development and the parametrization of accurate potential models. Several force fields have been developed in the past few years [1-7]. Within most of these force fields, a molecule is modelled as a chain of Lennard-Jones (or exponential-6 [4]) sites and possibly of electrostatic sites. The determination of the Lennard-Jones (or exponential-6) effective parameters is achieved from fits to properties of pure compounds. During this parametrization, interactions between unlike sites are estimated using a set of combining rules. Generally a geometric mean rule (the Berthelot rule) is applied for the well depth parameter while an arithmetic mean rule (the Lorentz rule) is applied for the exclusion diameter (in [1], a geometric mean rule is also applied for the diameter). Although the Berthelot rule is known to overestimate the well depth parameter significantly [8, 9], the effect of the Lorentz-Berthelot rules on the parametrization of these models and on the prediction of the properties of mixtures has been overlooked. Author for correspondence e-mail: delhom@rsc.anu.edu.au Two main reasons can account for this. First, more accurate combining rules often require the prior determination of additional parameters such as polarizabil-ities or ionization potentials [10-12]. Second, it is believed that the choice of the combining rules has little influence on the predicted properties. When all the Lennard-Jones sites are of similar sizes, the various sets of combining rules will yield similar unlike-pair potential parameters. This was observed in a previous work on alkane-H2S mixtures [13], where the CH2, CH3 or H2S groups were modelled by 'united atom' sites of similar sizes. However, in order to increase the accuracy of the predictions of molecular simulation, 'all-atoms' force fields (where the considerably smaller hydrogen atoms are treated explicitly) have been developed recently [1, 2]. In these models, each atom is represented by one interaction site. The atomic Lennard-Jones parameters associated with H, C or S take very different values. The choice of the combining rule will then be of great significance: according to the chosen set of combining rules, the value taken by the Q dispersion coefficient can vary by 9% for the interaction between a carbon atom and a hydrogen atom within the OPLS-all atoms model [1]. United atoms models parameters are also sometimes dissimilar enough to be sensitive to the choice of the set of combining rules. Potoff et al. [14] Molecular Physics ISSN 0026-8976 print/ISSN 1362-3028 online © 2001 Taylor & Francis Ltd http://www.tandf.co.uk/journals DOI: 10.1080/0026897001002004 1 620 J. Delhommelle and P. Millie showed recently that the Lorentz-Berthelot rules perform better for the CH3OH-hexane and H20-ethane mixtures, whereas the Kong rules [15] yield better results for the C02-alkane systems. However, two points prevented Potoff et al. from drawing definite conclusions. First, the potential parameters were determined by fitting the properties of pure compounds using the Lorentz-Berthelot rules to estimate the unlike-pair potential parameters between groups of like molecules. All those calculations were therefore partly based on the Lorentz-Berthelot rules. Second, polarization effects, which play an important role in those mixtures [13], were not taken into account. The aim of this work is to provide a clear picture of the performance of several combining rules when used together with common force fields. Simple fluids (rare gases), involving no electrostatic or polarization interactions, were studied in this work. We insist on the fact that our goal is not to give the most accurate description of the phase equilibria of mixtures of rare gases. Our goal is to study the effects of combining rules on the interactions between rare gases atoms in order to draw conclusions on such effects when atomic or united atoms sites of complex molecules interact with each other. It was shown recently that using accurate 'true' (two-body) pair potentials [16] in conjunction with three-body terms yielded a good description of the liquid-vapour phase equilibria of pure rate gases [17, 18]. However, since molecules are commonly described by 'effective' pair potential models, we chose to model the rare gases studied in this work by effective Lennard-Jones pair potentials. Three sets of combining rules, which do not require any additional parameters, were then chosen: the commonly used Lorentz-Berthelot rules (LB), the Kong rules (KG) [15] and the Waldman-Hagler rules (WH) [19]. They were applied to the calculation of liquid-vapour phase equilibria and compressed-liquid properties of mixtures of rare gases. The calculation of coexistence curves allowed us to test the effect of the sets of combining rules across a wide range of composition and pressure. A comparison of simulation results with experimental data enabled us then to assess the relative performance of each set of combining rules and even- tually to recommend the use of a particular set regarding this type of calculation. 2. Potential models and combining rules Each atom was modelled by a Lennard-Jones site. In such a model, two sites i and j interact with each other according to the following expression: r / \ 12 / \ 6l where r,y, e,y and cr,y are the site-site separation, the Lennard-Jones well depth and the Lennard-Jones diameter for the sites i and j. The Lennard-Jones parameters for each atom were determined by using the law of corresponding state, i.e. by matching the experimental critical point. We have previously estimated a reduced critical temperature of 77* = Tc/s =1.31 and a reduced critical density of p* = pca3 = 0.311 from Gibbs ensemble Monte Carlo (GEMC) simulations of a Lennard-Jones fluid [20]. These values are in good agreement with other simulation results [21-23]. By replacing Tc and pc by the experimental values of the critical temperature and critical density in these two expressions, it is straightforward to obtain the values of the Lennard-Jones parameters for each atom. This procedure ensures that the simulated GEMC coexistence curve matches the experimental coexistence curve for each pure compound. We have applied this procedure to Ne, Ar and Kr. The parametrization yielded unsatisfactory results for He, presumably because of quantum effects. Mixtures containing He were therefore discarded in the remainder of this work. The effective parameters so obtained are given in table 1 (the experimental critical points were taken from [24]). We point out that these parameters are within the range of values used in most force fields (for instance, cH = 2.5 A and cc = 3.5 A in the OPLS-AA model [1]). Therefore the conclusions drawn on mixtures containing these atoms will be valid for most force fields. 'True' (two-body) pair potential parameters obtained from molecular beam experiments [25] are also given in table 1. The effective e's are smaller than the experimental values, whereas the effective c's are Table 1. Lennard-Jones parameters and atomic masses (in [25] several molecular beam experimental results are given for atom. We report here the range of variations of these data). Atom Mass/gmol 1 This work Ref. [25] (e/*B)/K a/A (e/*B)/K a/A Ne 20.183 33.89 2.79 41-42 2.75-2.76 Ar 39.944 115.17 3.38 142-143 3.35-3.36 Kr 83.800 159.85 3.62 199-200 3.58-3.59 Inadequacy of Lorentz-Berthelot combining rules 621 slightly larger than the experimental er's. The effective pair potentials determined in this work therefore are less attractive than the true pair potential. This is consistent with the findings in [17, 18]. Both works show that the repulsive three-body contribution must be added to the two-body interactions when using a true pair potential to obtain a good agreement with the experimental coexistence curve. If i and j refer to unlke atoms, the Lennard-Jones parameters for these interactions are determined using a set of combining rules. It is quite common to use a geometric mean for e (the Berthelot rule) and an arithmetic mean for a (the Lorentz rule): 1/2 = U<7Ü + <7jj). (2) The Berthelot rule is known to overestimate the well depth parameter [8]. On the basis of the London formula, a geometric mean combining rule is more likely to be successful for the C6(= 4ecr6) dispersion coefficient I21i (3) 6 / 6 6x1/2 The KG and WH sets of combining rules both use equation (3) to determine the well depth. Using the Berthelot rule (equation (2)) and equation (3), the ratio between the C6 coefficients from the two combining rules for the dispersion part can be expressed in the following way: ^KG or WH C6_ (4) Equation (4) shows that the choice of a set of combining rules will have a small influence provided that the parameters a of unlike groups are close to each other [13]. Kong has then developed a combining rule for a according to the atomic distortion model [26]. In this model, the repulsive interaction between unlike atoms can be expressed as the sum of the deformation energies of each atom: C/7(r, + r/)=I[C/-P(2r,)+^p(2r/)], (5) where 2r, and 2rj are the separations at which the slopes of the potential energy curves of like species are equal, i.e. at the distortion plane. This condition can be expressed as dR dU^ (R) R=2r, dR (6) R=2r, By introducing equation (6) in (5), Kong obtained the following expression for the 12-6 Lennard-Jones potential: Table 2. Influence of the set of combining rules on the unlike-pair potential Lennard-Jones parameters. KG LB WH Ne- -Ar (e/kB)IK a/A 55.46 3.13 62.47 3.09 53.39 3.15 Ne- -Kr (e/kB)IK a/A 59.71 3.29 73.60 3.21 55.71 3.33 12 1 1/13- (7) The combining rule derived by Waldman and Hagler was based on a careful graphical analysis of rare gas data. They showed that an arithmetic mean combining rule for a6 rather than for a yielded accurate results for the unlike-pair diameter. The second combining rule is therefore in this case 1/6 (8) The unlike-pair potential parameters obtained with each set of combining rules are given in table 2. The KG and WH approaches yield parameters which are very close to each other. Using the KG rules, we obtain higher values of unlike-pair well depths and smaller diameters. The WH rules will therefore result in less attractive unlike-pair interactions. The Berthelot rules yield much higher values of the well depth (by approximately 15%). We note that the Lorentz rule yields the lowest values for the diameter. The Lorentz rule therefore tends to enhance the effect of the Berthelot rule: the application of the LB rules results in strongly attractive unlike-pair interactions. 3. Simulation details The Gibbs ensemble Monte Carlo (GEMC) method [27, 28] was used to compute the phase equilibria of the Ne-Ar and Ne-Kr binary mixtures. These calculations were performed for constant temperature, pressure and total number of molecules (GEMC-NPT simulations). The ratio of the different moves used in calculations was: 60% particle displacement, 1% volume fluctuation and 39% particle transfer between the boxes. Simulations were run on systems containing a total number of 400 atoms. An initial equilibration period of 107 steps was followed by 2 x 107 steps where thermodynamic properties were averaged. Isothermal isobaric (NPT) Monte Carlo simulations were also performed for compressed liquid mixtures of the Ne-Ar binary mixtures. Simulations were run on systems of a total number of 400 atoms. In this case, 622 J. Delhommelle and P. Millie Table 3. Coexisting Ne mole fraction in Ne-Ar mixture at T = 110.78 K (statistical uncertainties are of 0.7% for the Ne mole fraction in the liquid phase and of 2% in the vapour phase). Exp. [31] KG LB WH P = 27.34 bar % Ne (vap) 67.01 63.6 63.3 66.2 % Ne (liq) 2.46 2.6 4.5 2.3 p = 41.64 bar % Ne (vap) 75.08 73.7 72.0 72.2 % Ne (liq) 4.22 4.8 7.6 3.7 P = 55.30 bar % Ne (vap) 78.71 75.8 76.5 76.7 % Ne (liq) 5.93 6.4 11.1 5.2 P = 69.22 bar % Ne (vap) 81.03 79.2 79.4 79.3 % Ne (liq) 7.64 8.6 15.0 7.3 an initial equilibration period of 106 steps was followed by 107 steps where thermodynamic properties were averaged. In all the simulations, a spherical cutoff radius of 12 A was applied for the Lennard-Jones interactions. Average long range corrections were applied for distances beyond this cutoff radius [29]. 4. Application to the prediction of equilibrium properties of binary mixtures GEMC simulations were performed to determine the coexistence curve of a binary mixture of Ne-Ar at T = 110.78 K. The simulation results and experimental data are given in table 3 and plotted in figure 1. The coexisting Ne mole fractions in the vapour phase are correctly described by the three sets of combining rules. Simulated Ne mole fractions agree with experimental results with an accuracy of 1.9% for the KG rules, 3.7% for the LB rules and 1.9% for the WH rules (all rrms, relative root mean-square, errors). The LB rules are found to overestimate the Ne mole fraction significantly in the liquid phase (rrms of 91%). Unlike the LB rules, the KG and WH rules both agree with experimental results within statistical uncertainty (rrms at 11% and 9%, respectively). GEMC simulations were also performed to determine the coexistence curve of a binary mixture of Ne-Kr at T = 150 K. The simulation results as well as experimental data are given in table 4 and plotted in figure 2. The three sets of combining rules yield a good description of the coexisting vapour phase. Simulated Ne mole fractions agree with experimental results with an accuracy of 1.4% for the KG rules, 2.5% for the LB rules and 1.5% for the WH rules (all rrms errors). However, only the KG rules yield accurate results for the liquid phase (deviation with experimental Ne mole fractions of 3.5%). The WH rules underestimate the Ne mole fraction in the liquid phase by Figure 1. Coexisting Ne mole fraction in an Ne-Ar binary mixture at T = 110.78 K. Filled circles are experimental data [31], diamonds are simulation results obtained with the KG rules, triangles represent results obtained with the LB rules and squares those obtained with the WH rules. /u 1 CD1 A 1 1 1 1 1 ' ' i I ' 1 1 c» 1 i—i— 60 - ■ 50 - Q- 40 - 30 - AO - 20 ■ ■ ■ 1 ■ ■ ■ i......... , * 0.2 0.4 0.6 Ne mole fraction 0.8 Inadequacy of Lorentz-Berthelot combining rules 623 Figure 2. Coexisting Ne mole fraction in an Ne-Kr binary mixture at T = 150 K (details as for figure 1; the experimental data are taken from [32]). 110 100 90 80 CO CL 30 1 1 1 1 1 1 1 1 . 1 _ £9 A o J ■ . ........ . . i ... i ... : 0.2 0.4 0.6 Ne mole fraction 0.8 Table 4. Coexisting Ne mole fraction in Ne-Kr mixture at T = 150K (statistical uncertainties are of 0.7% for the Ne mole fraction in the liquid phase and of 2% in the vapour phase). Exp. [32] KG LB WH P = 27.34 bar % Ne (vap) 78.1 77.0 76.5 78.6 % Ne (liq) 2.21 2.2 4.7 1.9 P = 61.79 bar % Ne (vap) 83.8 82.8 81.9 82.4 % Ne (liq) 3.26 3.4 7.5 2.5 P = 81.14 bar % Ne (vap) 86.2 84.9 83.8 84.7 % Ne (liq) 4.55 4.3 10.2 3.3 P= 103.12 bar % Ne (vap) 87.6 86.4 85.2 86.2 % Ne (liq) 5.64 5.6 12.7 3.8 28.7% whereas the LB rules overestimate significantly the Ne mole fraction in the liquid phase (rrms of more than 100%). We have then extended these results by studying liquid binary mixtures. The atoms parameters were determined by matching the coexistence curve (i.e. by fitting the experimental critical point). We have first checked that these parameters allowed us to obtain accurate results for compressed liquid densities. Therefore we performed NPT Monte Carlo simulations of pure Ne and pure Ar with the parameters given in table 1. The results showed that the parameters for Ne and Ar could be used to study compressed liquids. The simulated densities were of 0.399 ± 0.01 gem-3 (1.374±0.01 gem-3) whereas the experimental density is of 0.405gem"3 [21] (1.366gem"3 [30]) for pure Ne (Ar). These values are for state points (r = 110K, P = 200 bar for Ne and r = 110.78 K, P = 344.7 bar for Ar) which are representative of the state conditions studied thereafter for the mixtures. NPT Monte Carlo simulations of compressed liquid mixtures of Ne-Ar Table 5. Density of compressed liquid mixtures of Ne-Ar for an Ne mole fraction of 50% at T = 110.78 K (densities are given in gem-3; statistical uncertainties are of 0.010 gcm~3). Pressure/bar Exp. [31] KG LB WH 296.5 0.979 0.985 1.031 0.966 344.7 1.021 1.019 1.066 1.005 413.7 1.068 1.066 1.109 1.049 482.6 1.107 1.103 1.146 1.090 were then performed for two different compositions at T = 110.78 K. The results obtained for Ne mole fractions of 50% and 22.5% are given in tables 5 and 6, respectively. All these results are plotted in figure 3, and they confirm what we have observed in the calculation of coexistence curves. If we consider a mixture of Ne-Ar with an Ne mole fraction of 50%, we find that the predicted densities obtained using the KG rules are in excellent agreement with the experimental data (rrms of 0.37%). The WH rules underestimate the liquid den- 624 J. Delhommelle and P. Millie 1.35 r-r Figure 3. Density of compressed liquid mixtures of Ne-Ar at 7 = 110.78 K. Results are presented from bottom to top for Ne mole fractions of 50% and 22.5% (details as for figure 1; the experimental data are taken from [31]). o 3 IS) c Q 200 250 300 350 400 P (bar) 450 500 550 Table 6. Density of compressed liquid mixtures of Ne-Ar for an Ne mole fraction of 22.5% at T = 110.78 K (densities are given in gem-3; statistical uncertainties are of 0.01 gcm~3). Pres sure/bar Exp. [31] KG LB WH 180.6 1.147 1.152 1.185 1.143 275.8 1.202 1.208 1.234 1.196 413.7 1.265 1.265 1.288 1.255 551.6 1.316 1.311 1.331 1.303 sities by 1.6%. The LB rules are found again to significantly overestimate the experimental data (rrms of 4.3%). Again, we observe that the KG rules yield accurate predictions for these compressed liquid mixtures. In view of the unlike-pair parameters listed in table 2 and of these last simulation results, it seems that the LB rules yield excessively attractive unlike-pair potentials whereas the WH rules give unlike-pair potentials that are slightly too weak. The same trends are observed also for the other Ne mole fraction, although the effects are of course weaker. It should be noted that even when one component is in large excess (Ne mole fraction of 22.5%), the LB rules significantly overestimate (by 2.3%) the liquid densities. combining rules, which do not require any additional parameter, were tested: the conventional Lorent-Berthelot and two sets of rules due to Kong and Waldman-Hagler. The choice of a set of combining rules is shown to have a significant effect on the prediction of thermodynamic properties of simple mixtures. The Lorentz-Berthelot rules are shown to be unable to predict accurately the coexistence curve of these simple binary mixtures. The Kong and Waldman-Hagler rules give results which are in good agreement with experimental data. The Kong rules are found to yield more accurate predictions. Simulations of compressed liquids confirm these findings. The Lorentz-Berthelot rule results in overly attractive unlike-pair potential parameters and overestimates significantly the liquid densities. On the other hand, the Waldman-Hagler rules yield slightly too weakly attractive unlike-pair parameters and underestimates the liquid densities. The Kong rules yield densities in excellent agreement with experiment for the two compositions studied. We believe that these results provide clear evidence of the inadequacy of the Lorentz-Berthelot combining rules. From the results obtained in this work, it seems that the Kong rules should be used instead. We hope this result will be considered in developments of future force fields. 5. Conclusion We have carried out GEMC simulations of the vapour-liquid phase equilibria of mixtures of Ne-Ar and Ne-Kr and NPT Monte Carlo simulations of compressed liquid mixtures of Ne-Ar for several Ne mole fractions. The sensitivity of the results to the choice of a set of combining rules was studied. Three sets of simple References [1] Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J., 1996, /. Amer. chem. Soc, 108, 6638. [2] Cornell, W. D., Cieplak, P., Bayly, C, Gould, I. R., Merz, K. M., Ferguson, D. M., Spellmeyer, D. C, Fox, T., Caldwell, J. W., and Kollman, P. A., 1995, /. Amer. chem. Soc, 117, 5179. Inadequacy of Lorentz-Berthelot combining rules 625 [3] Martin, M. G., and Siepmann, J. I., 1997, /. Amer. chem. Soc, 119, 8921. [4] Errington, J. R., Boulougouris, G. C, Economou, I. G., Panagiotopoulos, A. Z., and Theodorou, D. N., 1998, /. phys. Chem. B, 102, 8865. [5] Nath, S. K., Escobedo, F. A., and De Pablo, J. J., 1998, /. chem. Phys., 108, 9905. [6] Ungerer, P., Beauvais, C, Delhommelle, J., Boutin, A., Rousseau, B., and Fuchs, A. H., 2000, /. chem. Phys., 112, 5499. [7] Delhommelle, J., Tschirwitz, C, Ungerer, P., Granuccl G., Millie, P., Pattou, D., and Fuchs, A. H., 2000, /. phys. Chem. B., 104, 4745. [8] Maitland, G. C, Rigby, M., Smith, E. B., and Wakeham, W. A., 1981, Intermolecular forces (Oxford: Clarendon Press). [9] Stone, A. J., 1996, The Theory of Intermolecular Forces (Oxford: Clarendon Press). [10] Pena, M. D., Pendo, C, and Renuncio, J. A. R., 1981, /. chem. Phys., 76, 333, [11] Tang, K. T., and Toennies, J. P., 1986, Z. Phys. D, 1, 91. [12] Ihm, G., Cole, M. W., Toigo, F., and Klein, J. R., 1990, Phys. Rev. A, 42, 5244. [13] Delhommelle, J., Millie, P., and Fuchs, A. H., 2000, Molec. Phys., 98, 1895. [14] Potoff, J. J., Errington, J. R., and Panagiotopoulos, A. Z., 1999, Molec. Phys., 97, 1073. [15] Kong, C. L., 1973, /. chem. Phys., 59, 2464. [16] Aziz, R. A., 1993, /. chem. Phys., 99, 4518. [17] Anta, J. A., Lomba, E., and Lombardero, M., 1997, Phys. Rev. E, 55, 2707. [18] Marcelli, G., and Sadus, R., 1999, /. chem. Phys., 99, 4518. [19] Waldman, M., and Hagler, A. T., 1993, /. comput. Chem., 14, 1077. [20] Delhommelle, J., Boutin, A., Tavitian, B., Mackie, A. D.,and Fuchs, A. H., 1999, Molec. Phys., 96, 1517. [21] Smit, B., 1992, /. chem. Phys., 96, 8639. [22] Panagiotopoulos, A. Z., 1994, Int. J. Thermophus., 15, 6. [23] Martin, M. G., and Siepmann, J. I., 1998, /. phys. Chem. B, 102, 2569. [24] Vargaftik, N. B., 1975, Handbook of Physical Properties of Liquids and Gases (Washington, DC: Hemisphere). [25] Aziz, R. A., 1984, Inert Gases, edied by M. L. Klein (Berlin: Springer-Verlag) Chap. 2. [26] Smith, F. T., 1972, Proc. R. Soc. Lond. A, 5, 1708. [27] Panagiotopolous, A. Z., 1987, Molec. Phys., 61, 813. [28] Panagiotopolous, A. Z., Quirke, N., Stapleton, M., and Tildesley, D. J., 1988, Molec. Phys., 63, 527. [29] Allen, M. P., and Tildesley, D. J., 1987, Computer Simulation of Liquids (Oxford University Press). [30] Street, W. B., 1967, /. chem. Phys., 46, 3282. [31] Street, W. B., 1965, /. chem. Phys., 42, 500. [32] Miller, R. C, Kidnay, A. J., and Hiza, M. J., 1972, /. chem. Thermodynam., 4, 807.