Simulations of Vapor–Liquid Equilibria: Routine versus Thoroughness

Review pubs.acs.org/jced

Simulations of Vapor−Liquid Equilibria: Routine versus Thoroughness Ivo Nezbeda* Faculty of Science, J. E. Purkinje University, 400 96 Ú stí nad Labem, Czech Republic ABSTRACT: Selected vapor−liquid equilibrium literature data have been analyzed revealing in many cases very large fluctuations or very low precision, and inconsistencies, with some of them being even incorrect. Consequently, use of such data may result in a questionable estimate of other properties as, for example, the critical point location. It is shown that plotting the vapor compressibility factor as a function of temperature along its equilibrium curve may provide a simple but strong test of precision and correctness of such data. Examples are given for several models of water and also for three selected complex organic fluids for which simulation data have been reported in literature.

often so easy to prove it. For example, in 1979 we proved,4 using theoretical arguments, that the equation of state of a hard dumbbell fluid published by Aviram et al.5 obtained from simulations was incorrect but yet their equation has been subsequently used by many researchers. A similar story applies to extensive data for the EXP-6 potential published by Belonoshko and Saxena.3,6,7 As a recent example of potential problems with simulations we mention recently reported severe discrepancies involving results of solubilities of aqueous NaCl solution.8 On the basis of this experience we published some time ago certain recommendations for anyone carrying out simulation data about how to keep running simulations under control and to check her/his results.3 To summarize, when running simulations there are several pitfalls of which any simulator should be aware, and these were summarized in a recent excellent comprehensive review by Frenkel.9 It is also obvious that in addition to general rules, the simulation of specific properties also allows the possibility to use other criteria suited particularly for the particular property considered. Vapor−liquid equilibria (VLE) are likely the most frequently reported data of simulations on fluids. There are three commonly used methods: Monte Carlo (MC) or molecular dynamics (MD) simulation with an explicit interface, Gibbs ensemble (GE) MC simulation, and Grand Canonical MC simulation.1,2 All these methods provide, for pure fluids at a given temperature, T, the equilibrium liquid and vapor densities, ρL and ρV, respectively, and the orthobaric pressure, P, of varying precision in their dependence on the details of the

1. INTRODUCTION In view of rapidly advancing computer technology, molecular simulations1,2 have become a common tool to study a large number of problems at the molecular level and are thereby also becoming a common tool in applied science and technology. Moreover, with respect to the availability of both commercial and freeware software, simulation carrying out computer experiments is as easy as never before and practically anyone can “study” nowadays a system of her/his own interest by means of such software, even without proper understanding of the underlying methods. However, this general availability is also accompanied by a potential danger of producing low quality or even erroneous data. In addition to specifying parameters defining the system or problem at hand, for example, the thermodynamic conditions and force fields, simulations also require specification of a number of technical parameters such as the number of molecules in the simulation box, cutoff parameters, corrections for finite size effects, and length of equilibration run, etc., all of which affect the results. Properly setting all of the parameters may not be always straightforward and may require some experience. There are also certain auxiliary quantities3 which should be monitored during simulations to keep their development under control; however, the use of such quantities has been only rarely reported. Finally, papers reporting results should provide all necessary information on the simulation details and the estimation of the precision in the results, however, this is not always the case. Consequently, unless such simulations are not repeated independently in another lab, the inaccuracy or even gross errors of the reported data is difficult to detect and may be subsequently used by other trusting researchers. Over the years of our own involvement in developing and implementing various simulation methodologies, we have found errors in a number of published results but it is not © XXXX American Chemical Society

Special Issue: Proceedings of PPEPPD 2016 Received: June 29, 2016 Accepted: September 5, 2016

A

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Review

and 30 million configurations were generated in the production runs, with the actual number depending on the thermodynamic state. A low acceptance rate for molecule transfer to the liquid state was experienced, being about only 1% for temperatures below 423 K. During the simulation, the vapor box was alternatively empty or populated by only a few molecules, which caused large fluctuations in properties of the gas phase. The reported values were evaluated from 5 to 15 blocks of equal size. For the purpose of a fair comparison with other results, we use the standard error of the mean (SEM), considering averaging over 15 blocks at low temperatures and 5 blocks at higher temperatures. This corresponds to a roughly 90% confidence interval. Boulougouris et al. provide only scarce information on the details of their GE simulations. They used between 200 and 250 molecules and for the majority of the state points they ran two independent simulations at one state. Although in the table presenting the results there are listed uncertainties in both pressure and densities, no details on their meaning or how they were obtained are provided; with respect to the common practice it may also be assumed that the uncertainties mean the SEM. Finally, NIST data were obtained by grand-canonical Wang−Landau/Transition-matrix Monte Carlo and histogram reweighting.13−15 The uncertainties of the results were obtained from three independent simulation runs and represent 95% confidence limits. Concerning the likely most accurate currently available nonpolarizable model of water, the TIP4P/2005 model,16 we are aware of only one set of data by Vega et al.17 These authors used a number of techniques to perform cross checks to ensure that the produced numbers were correct; the final reported values were then obtained by the Gibbs−Duhem methodology18 starting from equilibrium values at T = 450 K and using 360 particles. The use of the Gibbs−Duhem methodology is particularly useful at low temperatures, because it bypasses the problems associated with the low number of particles in the gas phase and low probability of transfer of molecules to the dense liquid phase in the GE. Once the orthobaric pressures were determined, the coexistence densities were then obtained from long run simulations in an NPT ensemble. The equilibrium densities were obtained by averaging over 10 blocks, and the declared uncertainties represent the SEM. For completeness we also consider another member of the 4-site family of water models, TIP4P.16 Results for this model were reported by Vorholz et al. along with the data for SPC/E vapor and there were only minor differences in the technical details between the simulations on these two different models. From a large number of various polarizable models of water, those that use the Gaussian charge rather than point charge model seem to be currently preferred because of their more realistic description of physical reality. Typical examples of these models are the Gaussian Charge Polarizable Model (GCPM) of Paricaud et al.19 and the BK3 model of Kiss and Baranyai.20 Paricaud et al. and Moucka et al.21 used the GE to obtain VLE properties of GCPM and BK3 water, respectively, whereas Kiss and Baranyai used one simulation cell with two phases, that is, with an explicit interface, to study BK3 water. Paricaud et al. provide a table of the orthobaric densities and pressure but without any estimate of their precision and other details, and the same applies also to Kiss and Baranyai’s results. One likely reason for omission of this information is that their papers dealt with the development of the respective force fields and the VLE data were only peripheral results for demonstrating the accuracy of the developed models. Moucka

employed method. In the overwhelming majority of such simulations, results for the liquid phase and of interfacial properties are the main interest, with little attention paid to thermodynamics of the vapor phase. This is, first, because properties of the vapor phase can be estimated by the virial expansion, and also by the frequently unspoken shared conviction that the saturated vapor density is the upper limit of conditions of interest which is obviously not fully true. For example, in the turbine industry accurate data on supersaturated (dense) steam are of very high importance, particularly due to the lack of any experimental data in this region, with molecular simulations being thus the only available “experimental” tool; reliable and accurate data on at least the stable water vapor phase are thus highly demanded. Moreover, an equally important fact is that vapor densities obtained from simulations may be used for the evaluation of other thermodynamic properties, for example, for an estimation of the critical point, another property of great importance in industrial applications; if the vapor densities are inaccurate or wrong, the location of the critical point must also be inaccurate or wrong. In connection with our recent interest in supersaturated steam, we used literature data for the most common models of water to calculate the compressibility factor along the vapor coexistence line and to our surprise we discovered that some of the data were, literally, useless. This finding has motivated us to further investigate this problem and we found that a similar problem applies also to literature data found in the literature for various other compounds. In this short communication we thus propose and demonstrate a very simple test to check the consistency/precision of the VLE data resulting from any simulation method employed. This is illustrated using only data available from the literature, not only for water but also for several other compounds. In all such cases we found data which, for a given force field, cannot be completely correct. In the next section we present an analysis of the simulation details of the selected data and point out their differences, particularly with respect to estimates of their precision, and in the subsequent section we illustrate the problem and suggest how the consistency of simulation VLE data can be examined by computing the vapor compressibility factor and plotting it along the vapor saturation curve.

2. SIMULATION DETAILS AND UNCERTAINTY ESTIMATES OF THE SELECTED LITERATURE DATA Although the goal of this paper is neither the presentation of new simulation results nor an analysis of available data, it is worth providing, both for comparison and general assessment, some technical details of simulations which produced the data used in this paper. When presenting simulation data, in all cases authors describe in detail how they computed the energy of the studied system, particularly how they handled the long-range interactions and how they accounted for finite size effects but little attention is usually paid to processing the obtained results. So it is not unusual that the “experimental” (i.e., simulation) data are given either without any uncertainty estimate at all or if it is provided then without its detailed description. Because of its wide use in various applications, we focus first on the SPC/E nonpolarizable model of water for illustrative purposes. The VLE data for this model we consider are those by Vorholz et al.,10 Boulougouris et al.,11 and NIST.12 Vorholz et al. obtained the VLE data by GE simulations with the total number of particles varying from 140 to 350, and between 5 B

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Review

et al. implemented the Multi-Particle-Move MC22 into the GE and started their simulations with 250 particles in the liquid phase, and the number of particles in the gas phase varied from 5 at T = 298 K to 110 at 590 K. They provide an estimate of the SEM based on the block method, wherein the number of blocks varied in the dependence on the thermodynamic state. To examine the proposed consistency test also for other fluids, we have chosen from the available literature sources three organic compounds: benzene, n-decane, and 1,2ethanediol (ethylene glycol). The VLE data of benzene were determined by Janecek et al.23 using MC simulations with 343 molecules and an explicit interface. The density profile was discretized into slabs of thickness of 0.25 Å, and the densities of both phases were determined as averages of the respective local densities. Two different methods were used to determine pressure: (1) as the average value of the normal pressure profile within the vapor phase (method P1) and (2) from the zzcomponent of the virial tensor (method P2). Method P1 provides data of higher precision and also agrees well with the experimental values. The reported error bars were estimated as three times the standard deviation of a set of five blocks. n-Decane was studied by Muller and Meija.24 They performed MD simulations with an explicit interface to obtain both VLE data and interfacial properties. They considered three force fields with the main goal of their study being an assessment of these force fields. Density profiles were calculated by dividing the system in 250 slabs along the z direction, and the pressure was determined as the z component of the pressure tensor. The results were reported along with an estimate of the SEM. VLE data of ethylene glycol were determined by Huang et al.25 and also by Stubb et al.26 Huang et al. used the Grand Equilibrium method,27 a variant of the NPT+Test Particle method.28 In this method the liquid and gas phases are simulated separately. Using an NPT ensemble with 500 molecules, the chemical potential is determined first for a series of pressures and only one simulation of the gas phase then follows. The declared uncertainty means the SEM. Stubb et al. used the GE methodology with the coupled-dicoupled configurational bias. They used 250 molecules and about 50000 cycles in production runs. The declared errors are assumed to be the SEM. The quantity used in this paper to assess quality of the simulation data is the compressibility factor of vapor, z = PV/ NkBT ≡ P/ρkBT, where V is the volume, N is the number of molecules, and kB is the Boltzmann’s constant. Simulations provide pressure and density, which may be correlated and this must be taken into account when the uncertainty of the compressibility factor is estimated, ⎛ uρ ⎞2 ⎛ u ⎞2 uρuP uz = z ⎜ ⎟ + ⎜ P ⎟ − 2 corr(ρ , P) ⎝ ⎠ ρP P ⎝ρ⎠

we assume that they represent SEM unless explicitly specified otherwise.

3. CONSISTENCY TEST AND DISCUSSION As previously mentioned, the reported VLE values for pure fluids are pressure, temperature, and phase densities. When these results are individually compared with either experimental or literature data, they are typically in agreement. However, when they are combined, for example, to calculate the compressibility factor, inaccuracies or imprecision may be compounded and become discernible. In Figure 1 we show the

Figure 1. Compressibility factor of water vapor along its saturation curve. For the sake of clarity, no error bars are shown for the NIST and experimental data.

compressibility factor of water vapor along its saturation curve reported for the SPC/E model of water. At low temperatures the compressibility factor must start from a value around unity (ideal gas limit) and then be monotonically decreasing as temperature rises, up to the critical point value. As it is seen, only the NIST simulation data12 and those of Boulougouris et al.11 satisfy this requirement, with the latter exhibiting a large uncertainty. On the other hand, the data of Vorholz et al.10 are incorrect at low temperature and at higher temperatures they display large uncertainties. It is obvious that when these data are used to estimate the critical temperature, the result must be subject, at minimum, to a very large error or even to be erroneous. To further demonstrate that it is possible to generate fully consistent results over a very large temperature range we show in Figure 2 results for TIP4P/2005 water17 along with data by Vorholz et al.10 for the TIP4P water model, and compare both results also with experiment. Again, the former data form a smooth curve that also compares quite well with the experimental results, whereas the data of Vorholz et al. are subject to considerable fluctuations. Similar conclusions can be drawn also from an examination of data for polarizable models shown in Figure 3. The data of Paricaud et al.19 for the GCPM model demonstrate that even for such a complex model it is possible to obtain a smooth curve. On the other hand, the data of Kiss and Baranyai20 for the BK3 model are very chaotic and the data by Moucka et al.21 for the same model display an incorrect trend with decreasing temperature. These observa-

(1)

where u denotes uncertainty and “corr” is the correlation coefficient. Following Dinpajooh et al.,29 in the vapor phase this coefficient may be set to unity to a good approximation; if pressure and densities are uncorrelated then corr = 0. Pressure and densities obtained from the GE or explicit interface simulations are always correlated, in other cases they may be usually considered as uncorrelated. The above analysis points to diversity in the presentation of simulation results. Published data are certain average values, and when their errors are also given without any specification, C

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Review

Figure 4. Compressibility factor of benzene vapor along its saturation curve. Pi refers to different methods used to determine the orthobaric pressures.

Figure 2. Compressibility factor of water vapor along its saturation curve. Error bars for the data of ref 17 are comparable with the size of the symbols.

Figure 3. Compressibility factor of water vapor along its saturation curve. Error estimates of data from refs 19 and 20 are not available.

tions may also be one of the reasons for the reported difference in the critical point location for this model.21 The main source of error in the vapor compressibility factor is a very low vapor density and its large fluctuations. However, inaccuracies in the estimation of the orthobaric pressure may also contribute to large discrepancies as illustrated in Figure 4. Here we show results for benzene modeled by the TraPPE force-field simulated with an explicit interface (inhomogeneous simulations).23 In this case two different methods were used to determine pressure: (1) as the average value of the normal pressure profile within the vapor phase (method P1) and (2) from the zz-component of the virial tensor (method P2). As is seen, method P1 provides data of high precision and are also quite accurate when compared with experiment. On the other hand, the results of method P2 are accurate only at temperatures T > 450 K. Moreover, at low temperatures they are subject to very large errors. We have also checked the reported data for a number of other complex liquids and in Figure 5 we show results for n-

Figure 5. Compressibility factor of n-decane vapor (upper graph) and of ethylene glycol vapor (lower graph) along their saturation curves. For higher temperatures the declared errors of data for ethylene glycol from ref 25 are of the same size as the symbols.

decane24 and ethylene glycol vapors.25,26 As is obvious, the data of Muller and Mejia for n-decane, and those of Huang et al. for ethylene glycol seem practically useless. For the former liquid, the compressibility factor is monotonically decreasing for T < D

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Review

(5) Aviram, I.; Tildesley, D. J.; Streett, W. B. The virial pressure in a fluid of hard polyatomic molecules. Mol. Phys. 1976, 34, 881−885. (6) Belonoshko, A.; Saxena, S. K. A molecular dynamics study of the pressure-volume-temperature properties of super-critical fluids: 1. H2O. Geochim. Cosmochim. Acta 1991, 55, 381−387. (7) Belonoshko, A.; Saxena, S. K. A molecular dynamics study of the pressure-volume-temperature properties of super-critical fluids. 2. CO2, CH4, CO, and H-2. Geochim. Cosmochim. Acta 1991, 55, 3191− 3208. (8) Nezbeda, I.; Moučka, F.; Smith, W. R. Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility. Mol. Phys. 2016, 114, 1665−1690. (9) Frenkel, D. Simulations: The dark-side. Eur. Phys. J. Plus 2013, 128, 10−21. (10) Vorholz, J.; Harismiadis, V. I.; Rumpf, B.; Panagiotopoulos, A. Z.; Maurer, G. Vapor+liquid equilibrium of water, carbon dioxide, and the binary system, water+carbon dioxide, from molecular simulation. Fluid Phase Equilib. 2000, 170, 203−234. (11) Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N. Engineering a molecular model for water phase equlibrium over a wide temperature range. J. Phys. Chem. B 1998, 102, 1029−1035. (12) SAT-TMMC: Liquid-Vapor coexistence properties - SPC/E Water. http://www.nist.gov/mml/csd/informatics_research/spce_ water_sat_lrc.cfm (accessed August, 16, 2016). (13) Errington, J. R. Vapour liquid equilibrium of a pure fluid from test particle method in combination with NpT molecular dynamics simulations. J. Chem. Phys. 2003, 118, 9915−9925. (14) Shell, M. S.; Debenedetti, P. G.; Panagitopoulos, A. Z. An improved Monte Carlo method for direct calculation of the density of states. J. Chem. Phys. 2003, 119, 9406−9411. (15) Rane, K. S.; Murali, S.; Errington, J. R. Monte Carlo Simulation Methods for Computing Liquid-Vapor Saturation Properties of Model Systems. J. Chem. Theory Comput. 2013, 9, 2552−2566. (16) Vega, C.; Abascal, J. L. F. Simulating water with rigid nonpolarizable models: a general perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663−19688. (17) Vega, C.; Abascal, L. F.; Nezbeda, I. Vapor-liquid equilibria from the triple point up to the critical point for the new generation of TIP4P-like models: TIP4P/Ew, TIP4P/2005, and TIP4P/ice. J. Chem. Phys. 2006, 125, 034503−034511. (18) Kofke, D. A. Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line. J. Chem. Phys. 1993, 98, 4149−4162. (19) Paricaud, P.; Předota, M.; Chialvo, A. A.; Cummings, P. T. From dimer to condensed phases at extreme conditions: Accurate predictions of the properties of water by a Gaussian charge polarizable model. J. Chem. Phys. 2005, 122, 244511−244524. (20) Kiss, P.; Baranyai, A. A systematic development of a polarizable potential of water. J. Chem. Phys. 2013, 138, 204507−204523. (21) Moučka, F.; Nezbeda, I. Gibbs ensemble simulation on polarizable models: Vapor-liquid equilibrium in Baranyai-Kiss models of water. Fluid Phase Equilib. 2013, 360, 472−476. (22) Moučka, F.; Nezbeda, I.; Smith, W. R. Computationally efficient Monte Carlo simulations for polarisable models: multi-particle move method for water and aqueous electrolytes. Mol. Simul. 2013, 39, 1125−1134. (23) Janeček, J.; Krienke, H.; Schmeer, G. Inhomogeneous Monte Carlo simulation of the vapor-liquid equilibrium of benzene between 300 and 530 K. Condens. Matter Phys. 2007, 10, 415−423. (24) Muller, E. A.; Mejia, A. Comparison of united-atom potentials for the simulation of vapor-liquid equilibria and interfacial properties of long-chain n-alkanes up to n-C100. J. Phys. Chem. B 2011, 115, 12822−12834. (25) Huang, Y. L.; Merker, T.; Heilig, M.; Hasse, H.; Vrabec, J. Molecular modeling and simulation of vapor-liquid equlibria of ethylene oxide, ethylene glycol and water as well as their binary mixtures. Ind. Eng. Chem. Res. 2012, 51, 7428−7440.

450 K instead of increasing toward the limiting ideal gas value at low temperatures as demonstrated by the REPPROP correlation.30 The data for ethylene glycol of Huang et al. resemble randomly scattered numbers without any general trend. This is in a sharp contrast to the results of Stubb et al. which form a smooth curve going to the correct low temperature limit. It is also quite stunning, in the light of the graph, that the Huang et al. declare a very high precision of their data.

4. CONCLUSION In this paper we have attempted to draw the attention of researchers who produce VLE data by molecular simulations to a simple test to check the precision/consistency of their data. Plotting the vapor compressibility factor along its coexistence curve may reveal considerable imprecision (inconsistency) in the data and also warn the simulator not to use them for the evaluation of other thermodynamic properties. It also suggests that the compressibility factor be evaluated directly in the course of simulations (if the chosen methodology makes it possible). We have presented only several selected examples, but many more could be presented. It follows that users of published VLE data should be cautious when using them and should test them before further use. In any case, since the low temperature data are, in general, subject to larger uncertainty, it is recommended that their use for the estimation of the critical point location be avoided. At a more general level, the presented analysis of available data should be taken as a message to the simulation community to pay more attention to presentation of simulation results. Molecular simulations, referred to sometimes also as pseudoexperiments, are a modern counterpart of real laboratory experiments. Whereas it is unthinkable to publish experimental data without a detailed description of the experimental procedure, data estimation and their error uncertainty, this culture does not seem to be fully rooted in simulation community. The lack of such information makes it then very difficult to independently repeat a simulation and verify its result, for example, for the purpose of checking simulation codes. Furthermore, the comparison of different data sets and their assessment may be problematic.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +420 220390296. Funding

Support for this work was provided by the Czech Science Foundation (Grant No. 16-02647S). Notes

The author declares no competing financial interest.

■

REFERENCES

(1) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press; Oxford, 1987. (2) Frenkel, D.; Smit, B. Understanding Molecular Simulations; Academic Press; San Diego CA, 2002. (3) Nezbeda, I.; Kolafa, J. The use of control quantities in computer simulation experiments: application to the exp-6 potential fluid. Mol. Simul. 1995, 14, 153−163. (4) Nezbeda, I.; Smith, W. R.; Boublik, T. Conjectures on fluids of hard spherocylinders, dumbells, and spheres. Mol. Phys. 1979, 37, 985−989. E

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Review

(26) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-atom Description for Ethers, GLycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596−17605. (27) Vrabec, J.; Hasse, H. Grand Equilibrium: vapor-liquid equilibria by a new molecular simulation method. Mol. Phys. 2002, 100, 3375− 3383. (28) Moller, D.; Fischer, J. Vapour liquid equilibrium of a pure fluid from test particle method in combination with NpT molecular dynamics simulations. Mol. Phys. 1990, 69, 463−473. (29) Dinpajooh, M.; Bai, P.; Allan, D. A.; Siepmann, J. I. Accurate and precise determination of critical properties from Gibbs ensemble Monte Carlo simulations. J. Chem. Phys. 2015, 143, 114113−114125. (30) Lemmon, E. W.; Huber, M. L.; McLinden, M. O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, version 9.1; National Institute of Standards and Technology: Gaithersburg, MD, 2013.

F

DOI: 10.1021/acs.jced.6b00539 J. Chem. Eng. Data XXXX, XXX, XXX−XXX