7412
J. Phys. Chem. B 2004, 108, 7412-7414
Vapor-Liquid Equilibria in Five-Site (TIP5P) Models of Water Martin Lı´sal,†,‡ Ivo Nezbeda,*,†,‡ and William R. Smith§ E. Ha´ la Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences of the Czech Republic, 165 02 Prague 6, Czech Republic, Department of Physics, J. E. Purkyneˇ UniVersity, 400 96 U Ä stı´ nad Labem, Czech Republic, and Faculty of Science, UniVersity of Ontario Institute of Technology, 2000 Simcoe Street N, Oshawa, Ontario L1H 7K4, Canada ReceiVed: February 3, 2004; In Final Form: March 22, 2004
Using the Gibbs ensemble Monte Carlo simulations, we have determined vapor-liquid equilibria in a new, reparametrized five-site model (TIP5P-E) of water and compared them with the original TIP5P model of Mahoney and Jorgensen and TIP4P water. It is shown that for vapor-liquid equilibria properties, the new model provides only a marginal improvement over the original model and both models are considerably inferior to the TIP4P model.
1. Introduction The need to understand and predict properties of water has stimulated work toward its quantitative molecular modeling. Because of computational expense of quantum mechanical methods, the majority of models employ van der Waals forces with point-charge electrostatics embedded in a rigid body. Typical examples of such commonly used models are the simple point-charge models (SPC1 or SPC/E2) or transferable intermolecular potentials models (TIPn3); for an exhaustive review on water potential models, see ref 4. As an attempt toward improving fixed-charge models, Mahoney and Jorgensen5 have recently extended the family of TIP models by a five-site model with the ambition to develop a superior effective pair potential. The parameters of the resulting potential, called TIP5P, were optimized primarily with respect to the temperature of density maximum at 1 atm and the structural and thermodynamic properties at ambient conditions. To simplify and speed up the optimization process, they used in simulations a simple spherical cutoff at 9 Å. Not only does this makes the model dependent on technical details of simulations,6 but it may also be a considerable obstacle for its application to mixtures (and heterogeneous systems in general) when other potentials are not necessarily developed using cutoffs. These facts have prompted Rick to reparametrize the model by accounting properly for the long-range electrostatic interactions using the Ewald summation.7 This new model, referred to as TIP5P-E model, reproduces the same set of experimental data as the original TIP5P model. A general idea behind all the intermolecular potential models is the possibility (necessity) to predict various properties of the considered fluid at various thermodynamic conditions, not necessarily at those used in the parametrization. Success or failure to do so may also serve, indirectly, as an indicator of the usefulness of the developed model. In connection with the new TIP5P-E model, two immediate questions arise: (i) how does the new model perform in comparison with the original TIP5P model outside the range of conditions used for param* To whom correspondence should be addressed. E-mail: IvoNez@ icpf.cas.cz. † Academy of Sciences of the Czech Republic. ‡ J. E. Purkyne ˇ University. § University of Ontario Institute of Technology.
TABLE 1: The Lennard-Jones Parameters, E and σ, of the Oxygen Site and Partial Charges qi of the TIP5P and TIP4P Potential Models; kB Is Boltzmann’s Constant TIP5P [5] TIP5P-E [7] TIP4P [3]
O/kB [K]
σO [Å]
qO
qH
80.57 89.63 78.02
3.120 3.097 3.154
0 0 0
0.241 0.241 0.52
etrization and (ii) how does it perform for properties not considered in parametrization. Two of the most important properties required by the engineering community in practical applications are vaporliquid equilibrium densities and the vapor pressure, which are not usually considered in the process of optimization of parameters for models of water. The SPC and TIP4P3 models have been successfully used for this purpose,8 and it is therefore appealing to test the five-site models also against equilibrium data. This is the goal of this article. 2. Models and Computational Details The five-site potential functions are descendants of the ST2 model9 and thus involve a rigid water monomer with positive charges on two hydrogen H-sites, and negative charges on two M-sites mimicking the lone electron pair. Denoting by O the neutral oxygen site, the potential function assumes the form
u(1, 2) ) uLJ(rOO) + uCoul(1, 2) ) 4O
[( ) ( ) ] σO
rOO
12
-
σO
6
+
rOO qiqj
∑∑
i∈{1}j∈{2}
rij
(1)
where rij denotes the intermolecular i-j charge-charge separation, O and σO are the oxygen Lennard-Jones (LJ) interaction parameters, and qi is the partial charge. The monomers of both TIP5P and TIP5P-E models have the same form of a distorted tetrahedral geometry, with H-O-H and M-O-M angles being 104.52° and 109.47°, respectively, and oxygen-charge separations O-H ) 0.9572 Å and O-M ) 0.70 Å. Furthermore, the models also have the same charges on the sites and differ only in the LJ parameters (Table 1). We compare the results for the TIP5P models with the TIP4P model, which has also two H-sites
10.1021/jp0495242 CCC: $27.50 © 2004 American Chemical Society Published on Web 05/04/2004
Vapor-Liquid Equilibria in Five-Site Models of Water
J. Phys. Chem. B, Vol. 108, No. 22, 2004 7413
TABLE 2: Vapor-Liquid Equilibrium Data (Excess Internal Energy, u, Density, G, and Vapor Pressure, Psat) of TIP5P-E Water Obtained from Gibbs Ensemble Monte Carlo Simulationsa
a
T [K]
ug [kJ/mol]
ul [kJ/mol]
Fg [kg/m3]
Fl [kg/m3]
Psat [bar]
298.15 325 350 375 400 425 450 475 500 525
-0.18 ( 0.12 -0.42 ( 0.18 -0.61 ( 0.11 -1.14 ( 0.24 -1.76 ( 0.22 -2.58 ( 0.19 -3.47 ( 0.19 -5.24 ( 0.28 -7.05 ( 0.36 -9.04 ( 0.75
-41.13 ( 0.10 -38.73 ( 0.10 -36.76 ( 0.07 -35.00 ( 0.10 -33.19 ( 0.06 -31.46 ( 0.07 -29.64 ( 0.09 -27.89 ( 0.09 -25.93 ( 0.10 -23.54 ( 0.25
0.052 ( 0.008 0.21 ( 0.02 0.57 ( 0.03 1.42 ( 0.06 3.03 ( 0.91 6.13 ( 0.20 11.13 ( 0.42 21.14 ( 1.01 37.69 ( 2.42 71.38 ( 8.16
999.74 ( 6.00 984.98 ( 3.95 958.26 ( 2.32 927.19 ( 3.33 888.77 ( 2.42 848.98 ( 3.03 796.78 ( 5.04 743.51 ( 5.02 674.48 ( 8.05 578.53 ( 16.24
0.070 ( 0.008 0.327 ( 0.037 0.882 ( 0.037 2.274 ( 0.124 5.028 ( 0.310 10.01 ( 0.27 17.90 ( 0.54 31.26 ( 0.63 50.04 ( 2.25 78.87 ( 5.81
Superscripts g and l denote the vapor and liquid phases, respectively.
TABLE 3: Comparison of the Critical Properties of TIP4P and Five-Site Models of Water (as Estimated from the Wegner Expansion, Eq 2, and the Antoine Equation, Eq 3) with Experimental Data16 and Constants A, B and C of the Antoine Equation TIP5P TIP5P-E TIP4P exptl
Tc [K]
Pc [bar]
Fc [kg/m3]
A
B
C
521.3 540.9 588.2 647.096
85.6 102.2 149.3 220.640
337 328 315 322.000
13.310679 13.642988 11.774845
-4642.0526 -5073.9387 -3563.6212
2.5861483 21.888973 -61.71318
but only one negatively charged M-site on the H-O-H bisector 0.15 Å away from the oxygen site; the parameters of this model are also given in Table 1. As a consequence of the different number of sites and their different arrangement in the TIP5P and TIP4P models, the respective monomers bear slightly different dipole moments: 2.29 D for the TIP5P model and 2.18 D for the TIP4P model. To study vapor-liquid equilibria (VLE), we used the standard Monte Carlo simulations with N ) 756 molecules in an NVT Gibbs ensemble (GE).10 The simulations were organized in cycles as follows. Each cycle consisted of three steps: nD translational and rotational moves, nV volume moves, and nT particle transfers. The types of moves were selected at random with fixed probabilities, chosen so that the appropriate ratio of moves was obtained; nD:nV:nT was set at N:1:7600. The acceptance ratios for translational and rotational moves and for volume changes were adjusted to approximately 30%. After an initial equilibration period of 2 × 104 cycles, we generated 5 × 104 cycles to accumulate averages of the desired quantities. The precision of the simulated data was obtained using block averages,11 with 5000 cycles per block. The long-range interactions were treated by means of the reaction field method12 with the dielectric constant RF ) 1 for vapor box and RF ) ∞ for liquid box.
the minus sign to the gas phase, and subscript e at β indicates that the exponent is considered as an effective quantity. For
3. Results and Conclusions Results of the GE simulations are given in Table 2 and in Figure 1, where the equilibrium densities and vapor pressures of the five-site models are compared with those of the TIP4P model.13 As we see, the TIP5P-E model provides only a marginal improvement over the original TIP5P model. Liquid density agrees with experimental data only at lowest temperatures, and concerning the vapor pressure, it is considerably off at all temperatures. We have also attempted to locate the critical point from the obtained data using the standard procedure based on an expansion in powers of ∆T ) | T - Tc|:14
1 F( ) Fc + C2∆T ( A0∆Tβe 2
(2)
where F is the density, T is the temperature, subscript c refers to the critical point, the plus sign refers to the liquid phase and
Figure 1. Equilibrium vapor-liquid densities (top) and vapor pressures (bottom) of different models of water (symbols) and their comparison with experimental data (solid line).
estimating the critical pressure, Pc, we fitted the vapor pressure data to the Antoine equation:15
ln Psat ) A +
B T+C
(3)
7414 J. Phys. Chem. B, Vol. 108, No. 22, 2004 The estimated critical properties along with the parameters of the Antoine equation are given in Table 3. As one could expect from the performance for the VLE, both five-site models give the critical properties grossly off both the correct values and the prediction of the TIP4P model. The success of any potential model may be assessed by its ability (i) to maintain its accuracy outside the range of conditions used in optimization and (ii) to predict correctly the properties which were not considered in optimization. There seems to be no doubt that at low temperatures the TIP5P-E model is the best choice for estimating a number of thermodynamic and structural properties. However, when considering vapor-liquid equilibria, the vapor pressure is considerably off even at these conditions. Equilibrium liquid densities agree well at the lowest temperature, but with increasing temperature they rapidly deviate from experimental data. Consequently, for these properties the five-site models perform considerably worse in comparison with the four-site TIP4P model. To summarize, another attempt to describe the properties of water by means of an effective pairwise potential has not been completely successful. The fundamental question, to what extent this is possible by enlarging the set of properties and range of thermodynamic conditions considered, or impossible in principle, thus remains open. Acknowledgment. This research was supported by the Grant Agency of the Academy of Sciences of the Czech Republic (Grant No. IAA4072309), the Grant Agency of the Czech
Lı´sal et al. Republic (Grant No. 203/02/0764), and the National Research Council of Canada (Grant OGP 1041). We would like to thank Dr. S. W. Rick, University of New Orleans, for providing us with the manuscript of ref 7 prior to publication. Calculations were made on the SGI ALTIX 3300 in the Supercomputing Centre of the Czech Technical University, Prague. References and Notes (1) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (2) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (3) (a) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (b) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381. (4) Guillot, B. J. Mol. Liquids 2002, 101, 219. (5) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (6) Lı´sal, M.; Kolafa, J.; Nezbeda, I. J. Chem. Phys. 2002, 117, 8892. (7) Rick, S. W. J. Chem. Phys. 2004, 120, 6085. (8) Chialvo, A. A.; Cummings, P. T. AdV. Chem. Phys. 1999, 109, 115. (9) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (10) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, CA, 1996. (11) Flyvbjerg, H.; Petersen, H. G. J. Chem. Phys. 1989, 91, 461. (12) Neumann, M. Mol. Phys. 1983, 50, 841. (13) Lı´sal, M.; Smith, W. R.; Nezbeda, I. Fluid Phase Equilib. 2001, 181, 127. (14) Wegner, F. J. Phys. ReV. B 1972, 5, 4529. (15) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids; McGraw-Hill: New York, 2000. (16) NIST Webbook. http://webbook.nist.gov/chemistry/fluid.
