Journal of Molecular Liquids 129 (2006) 120 – 124 www.elsevier.com/locate/molliq
On the re-engineered TIP4P water models for the prediction of vapor–liquid equilibrium Ariel A. Chialvo a,⁎, Albert Bartók b , András Baranyai b a b
Chemical Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831-61 10, USA Department of Theoretical Chemistry, Eötvös University, 1518 Budapest 112, PO Box 32, Hungary Available online 7 September 2006
Abstract We perform extensive Gibbs Ensemble Monte Carlo simulations to study the capability of some recently re-parameterizations of the original TIP4P model intended to predict accurately the vapor–liquid coexistence envelope of water, its critical point, and its temperature dependence for the vapor pressure and second virial coefficient, and complement this analysis with the characterization of some specific crystalline faces of ice. We also disclose some trends between the resulting dipole moment of the models and the Lennard–Jones parameters, the location of the negative charge, as well as the estimated critical temperature. Finally, we discuss the inability of these models to predict accurately and simultaneously the melting temperature and the temperature of maximum density. © 2006 Elsevier B.V. All rights reserved. Keywords: Gibbs Ensemble Monte Carlo; Water models; Ices; Melting point; Density maximum
1. Introduction After more than three decades of effort behind the development of intermolecular potential models for water, since the pioneering work of Ben-Naim and Stillinger [1], it becomes clear that despite their simplicity, these models have been able to capture a great deal of the features behind the complex behavior of water [2]. This is obviously a remarkable feat considering that these models involve frequently 2–3 forcefield parameters, their parameterization is usually done at specific state conditions, and yet they are able to predict many water properties away from those conditions [3,4]. Among these simple models, the SPCE [5] and the TIP4P [6] models have become the most successful and frequently used in the molecular simulation studies of water at extreme conditions, i.e., either high temperature fluid—[4,7] or low temperature solid-state conditions [8]. A renewed effort in the modeling of water has recently emerged, one centered around the original version of the TIP4P ⁎ Corresponding author. Fax: +1 865 574 4961. E-mail address:
[email protected] (A.A. Chialvo). 0167-7322/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2006.08.018
model [6] and thrust by the refitting of the force-field, namely, (a) to refine the original model by an elaborate fitting procedure and to account properly for the long-range electrostatic interactions via an Ewald sum, i.e., the TIP4P-EW [9], as well as (b) to predict more accurately the melting temperature, i.e., TIP4P-ice [10], and (c) to predict more accurately the temperature of the density maximum as well as other properties, i.e., TIP4P-2005 [11]. These activities have also provided a larger dataset of predicted properties on a wider range of state conditions, which in turn, are giving us the unprecedented opportunity to unravel some additional consequences regarding the limitations of these re-engineered TIP4P models.
Table 1 Parameters of the TIP4P model and its variants Model
qH (e)
ℓOH (Å)
∡HOH (deg)
ℓOM (Å)
ε/k (K)
σ(Å)
μ (d)
TIP4P TIP4p-ice TIP4p-2005 TIP4p-I
0.52 0.5897 0.5564 0.5673
0.9572 0.9572 0.9572 0.9572
104.52 104.52 104.52 104.52
0.15 0.1577 0.1546 0.1556
78.12 106.1 93.20 93.90
3.153 3.1668 3.3589 3.1630
2.177 2.425 2.305 2.345
A.A. Chialvo et al. / Journal of Molecular Liquids 129 (2006) 120–124
121
Table 2 Critical conditions, melting and maximum density temperatures as predicted by the TIP4P models Model
Tc (K)
ρc (g/ Pc cm3) (bar)
Ref.
Tm#
Tmd
Ref.
TIP4p TIP4p-ice TIP4p-2005 TIP4P-I Exp.
558 705 617 639 647.1
0.313 0.305 0.309 0.315 0.322
[14] [13] This work This work [17]
232.0 272.2 252.1 – 273.15
253, 258 295 278 – 277
[15,16] [16] [11] – [17]
#
112.0 184.0 121.0 118.0 220.6
For the Ih ice polymorph.
In this research note we first determine the vapor–liquid phase envelope for the TIP4P-2005 water model and its corresponding critical point, its temperature dependence of the vapor pressure, and its second virial coefficient to complete the characterization of the model. Then, we calculate the configurational internal energy and density of the hexagonal ice (Ih) and ice VII and ice VIII as examples of low- and high-pressure phases, respectively. Finally, we rationalize the findings in terms of the trend in the corresponding force-fields, and highlight the limitations of this type of parameterization aimed at the more accurate simultaneous description of larger number of water properties.
Fig. 2. Temperature dependence of the second virial coefficient of water as predicted by the three models in comparison with the experimental data of Ref. [18].
Our analysis comprises the original TIP4P water [6] and its two recent variations, i.e., TIP4P-ice [10] and TIP4P-2005 [11]. For the sake of clarity we display the corresponding model parameters in Table 1, and the resulting critical conditions based on the Gibbs Ensemble Monte Carlo simulated orthobaric curves in conjunction with the Wegner expansion in Table 2 (see Ref. [13] for details). Details of the simulation methodologies as well as for the prediction of the second virial coefficients have
already been given in our previous publications [12,13] and therefore they will not be repeated here. According to Fig. 1, the TIP4P-2005 model describes the vapor–liquid phase envelope with the same level of accuracy as for the TIP4P-EW (see Fig. 2 and Table 2 of Ref. [13]), a not too surprising result given the fact that the two models predict practically the same temperature dependence for both the second virial coefficient and the vapor pressure (see Figs. 2 and 3). In fact, the TIP4P-2005 model under-predicts the temperature dependence of the orthobaric densities along the entire vapor branch, and above 400 K for the liquid branch. Consequently, the temperature dependence of the vapor pressure becomes also under-predicted, following approximately the same trend as that for the TIP4P-EW model, i.e., between the slight over-prediction
Fig. 1. Vapor–liquid equilibrium envelope of the TIP4P-2005 water model in comparison with the experimental data of Ref. [17].
Fig. 3. Coexistence vapor pressure of the TIP4P, TIP4P-ice, TIP4P-2005, and TIP4P-I water models in comparison with the experimental data of Ref. [17].
2. Models, methodology, and results
122
A.A. Chialvo et al. / Journal of Molecular Liquids 129 (2006) 120–124
Fig. 6. Second virial coefficient of the TIP4P-I in comparison with the other three model variations.
of the original TIP4P model and the large under-prediction of the TIP4P-ice model. Notice that the magnitude of the dipole moment of these models is represented by μ = 2qH(ℓOHcosϑ − ℓOM) with ϑ ≡ 0.5 ∡ HOH. Therefore, at fixed geometrical conditions for ℓOH, and ϑ, the dipole moment μ goes linear with ℓOM as clearly depicted in Fig. 4. 3. Discussion and final remarks
Fig. 4. Dipole moment dependence of the adjusted force-fields for the TIP4P model variations.
Fig. 5. Vapor–liquid equilibrium envelope of the TIP4P-I water model in comparison with the experimental data of Ref. [17].
The simulation analysis provides the strongest evidence yet about the inability of these TIP4P model variations to predict accurately and simultaneously the melting temperature and the temperature of the density maximum of water at ambient pressure. In other words, because of the simplicity, and consequently
Fig. 7. Dependence of the configurational internal energies of Ih ice on the dipole moment of the water model (■), ice VII (●) and ice VIII (▴). Full lines at 5 K. Dashed lines at the temperature of measurements, 250 K, 295 K, and 10 K, respectively (See Table 1 1.2 in Ref. [19]).
A.A. Chialvo et al. / Journal of Molecular Liquids 129 (2006) 120–124
Fig. 8. Dependence of the resulting ice densities on the dipole moment of the water model using the same notation as in the previous figure. Due to the small difference in temperature for ice VIII the calculated values are practically indistinguishable.
small number of adjustable parameters in these models, we are limited to only a few choices for the relevant properties to be fitted. The underlying consequence is that the adjusted forcefields, the model geometry and the resulting properties (such as electrical multipoles and critical conditions) become strongly correlated as clearly illustrated in Fig. 4. In fact, based on the data of Tables 1 and 2, and the quasilinear behavior depicted in Fig. 4, we could have predicted the sets of force-fields for which the critical temperature of the resulting model were that of the real water. As an illustration, a simple interpolation based on the dipole moment as the independent variable, allowed us to define the TIP4P-I (interpolated) model (see Table 1), to determine the complete vapor–liquid equilibrium envelope, and finally to estimate the corresponding critical conditions (see Table 2, and Fig. 5). As we would expect it, the TIP4P-I model suffers from the same defects as the other three members of the TIP4P family, since their parameters scale linearly with the strength of the dipole moment of the corresponding model. This is also clearly captured by the second virial coefficient, since the dipole–dipole interactions are the leading contributions to the total potential energy (Fig. 6). The quasi-linear variation of properties in terms of the dipole moment of the molecule is also manifested in the energies and the densities of the ice phases studied as shown in Figs. 7 and 8. Moreover according to Table 2 and Fig. 4, it becomes evident that these models are not able to capture some essential physics behind the anomalous density maximum of water, in that, most models investigated in the literature (see also Vega and Abascal [16].) exhibit an approximately constant temperature gap (≈ 20 K) between the melting temperature of the Ih ice and the temperature of the maximum density. According to the results of Abascal and Vega [11] the densities of the ice phases, as described by the TIP4P-2005 model, are predicted with reasonable accuracy and the variation of the liquid density with temperature is a remarkably good agreement with experiments. (Although, our experience indicated that to achieve a four-digits
123
accuracy in the resulting density, as they reported, for an isothermal-isobaric simulation requires extremely long simulations). However, our results in Fig. 7 for the dipole dependence of the configurational energies of the ice phases indicate that the energy differences between ices VII and Ih at 5 K are too large [12], ∼ 16 kJ/mol, for all models. In addition, the best estimate for the configurational energy for the ice Ih is given by the TIP4P-2005 model. Note that this value is much smaller than the corresponding value predicted by the TIP4P-Ice model, even though the latter was parameterized to get the correct melting temperature of this ice. This finding would suggest that, in this respect, TIP4P-I might show a better performance. In summary, our findings indicate that any attempt to produce a good match between the relevant water properties and the corresponding predictions from planar nonpolarizable models is unattainable, even when constraining the study to ice polymorphs, as opposed to the entire water phase diagram. The inability of these models to decouple these predictions of the normal melting temperature of ice Ih, and the temperature of maximum density might be providing a hint on what might be missing in their formulation. Acknowledgements This research was sponsored by the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences under contract number DE-ACO5-OOOR22725 with Oak Ridge National Laboratory, managed and operated by UT-Battelle, LLC. AB gratefully acknowledges the financial support from OTKA grant TO43542. References [1] A. Ben-Naim, F.H. Stillinger, in: R.A. Horne (Ed.), Structure and Transport of Processes in Water and Aqueous Solutions, WileyInterscience, New York, 1972. [2] B. Guillot, Journal of Molecular Liquids 101 (1–3) (2002) 219. [3] Y. Guissani, B.J. Guillot, Journal of Chemical Physics 98 (1993) 8221. [4] A.A. Chialvo, P.T. Cummings, in: S.A. Rice (Ed.), Advances in Chemical Physics, vol. 109, Wiley & Sons, New York, 1999, p. 115. [5] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, Journal of Physical Chemistry 91 (1987) 6269. [6] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Journal of Chemical Physics 79 (1983) 926. [7] A.G. Kalinichev, Reviews in Mineralogy and Geochemistry 42 (2001) 83. [8] C. McBride, C. Vega, E. Sam, L.G. MacDowell, J.L.F. Abascal, Molecular Physics 103 (1) (2005) 1; S.C. Gay, E.J. Smith, A.D.J. Haymet, Journal of Chemical Physics 116 (20) (2002) 8876; L.A. Baez, P. Clancy, Journal of Chemical Physics 103 (22) (1995) 9744; I.M. Svishchev, P.G. Kusalik, Chemical Physics Letters 239 (4–6) (1995) 349; P.H. Poole, U. Essmann, F. Sciortino, H.E. Stanley, Physical Review E 48 (6) (1993) 4605. [9] H.W. Horn, W.C. Swope, J.W. Pitera, J.D. Madura, T.J. Dick, G.L. Hura, T. Head-Gordon, Journal of Chemical Physics 120 (20) (2004) 9665. [10] J.L.F. Abascal, E. Sam, R.G. Fernandez, C. Vega, Journal of Chemical Physics 122 (23) (2005) 234511. [11] J.L.F. Abascal, C. Vega, Journal of Chemical Physics 123 (23) (2005) 234505. [12] A. Baranyai, A. Bartok, A.A. Chialvo, Journal of Chemical Physics 123 (5) (2005) 054502.
124
A.A. Chialvo et al. / Journal of Molecular Liquids 129 (2006) 120–124
[13] A. Baranyai, A. Bartok, A.A. Chialvo, Journal of Chemical Physics 124 (7) (2006) 074507. [14] A. Baranyai, A. Bartok, A.A. Chialvo, Journal of Molecular Liquids. (submitted for publication). Special issue on the EMLG-2005 Meeting (Prague, Sep. 2005). [15] W.L. Jorgensen, C. Jenson, Journal of Computational Chemistry 19 (10) (1998) 1179. [16] C. Vega, J.L.F. Abascal, Journal of Chemical Physics 123 (14) (2005) 144504.
[17] L. Haar, J.S. Gallagher, G.S. Kell Steam, Steam Tables, Hemisphere Publishing Corporation, New York, 1984. [18] P.G. Hill, R.D.C. Macmillan, Industrial and Engineering Chemistry Research 27 (5) (1988) 874. [19] V.F. Petrenko, R.W. Whitworth, Physics of Ice, Oxford University Press, 2002.