The melting temperature of the most common models of water

The melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E,TIP3P,TIP4P, TIP4P/Ew, and TIP5P) is obtained from computer simulations at p = 1 bar. Since the melting temperature of ice I(h) for the TIP4P model is now known [E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phys. Rev. Lett. 92, 255701 (2004)], it is possible to use the Gibbs-Duhem methodology [D. Kofke, J. Chem. Phys. 98, 4149 (1993)] to evaluate the melting temperature of ice I(h) for other potential models of water. We have found that the melting temperatures of ice I(h) for SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P models are T = 190 K, 215 K, 146 K, 232 K, 245 K, and 274 K, respectively. The relative stability of ice I(h) with respect to ice II for these models has also been considered. It turns out that for SPC, SPC/E, TIP3P, and TIP5P the stable phase at the normal melting point is ice II (so that ice I(h) is not a thermodynamically stable phase for these models). For TIP4P and TIP4P/Ew, ice I(h) is the stable solid phase at the standard melting point. The location of the negative charge along the H-O-H bisector appears as a critical factor in the determination of the relative stability between the I(h) and II ice forms. The methodology proposed in this paper can be used to investigate the effect upon a coexistence line due to a change in the potential parameters.     

Relation between the melting temperature and the temperature of maximum density for the most common models of water

Water exhibits a maximum in density at normal pressure at 4 degrees above its melting point. The reproduction of this maximum is a stringent test for potential models used commonly in simulations of water. The relation between the melting temperature and the temperature of maximum density for these potential models is unknown mainly due to our ignorance about the melting temperature of these models. Recently we have determined the melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P). In this work we locate the temperature of maximum density for these models. In this way the relative location of the temperature of maximum density with respect to the melting temperature is established. For SPC, SPC/E, TIP3P, TIP4P, and TIP4P/Ew the maximum in density occurs at about 21-37 K above the melting temperature. In all these models the negative charge is located either on the oxygen itself or on a point along the H-O-H bisector. For the TIP5P and TIP5P-E models the maximum in density occurs at about 11 K above the melting temperature. The location of the negative charge appears as a geometrical crucial factor to the relative position of the temperature of maximum density with respect to the melting temperature.     

The thickness of a liquid layer on the free surface of ice as obtained from computer simulation

Molecular dynamic simulations were performed for ice I(h) with a free surface by using four water models, SPC/E, TIP4P, TIP4P/Ice, and TIP4P/2005. The behavior of the basal plane, the primary prismatic plane, and of the secondary prismatic plane when exposed to vacuum was analyzed. We observe the formation of a thin liquid layer at the ice surface at temperatures below the melting point for all models and the three planes considered. For a given plane it was found that the thickness of a liquid layer was similar for different water models, when the comparison is made at the same undercooling with respect to the melting point of the model. The liquid layer thickness is found to increase with temperature. For a fixed temperature it was found that the thickness of the liquid layer decreases in the following order: the basal plane, the primary prismatic plane, and the secondary prismatic plane. For the TIP4P/Ice model, a model reproducing the experimental value of the melting temperature of ice, the first clear indication of the formation of a liquid layer, appears at about -100 degrees C for the basal plane, at about -80 degrees C for the primary prismatic plane, and at about -70 degrees C for the secondary prismatic plane.     

The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface

In this work we present an implementation for the calculation of the melting point of ice I(h) from direct coexistence of the solid-liquid interface. We use molecular dynamics simulations of boxes containing liquid water and ice in contact. The implementation is based on the analysis of the evolution of the total energy along NpT simulations at different temperatures. We report the calculation of the melting point of ice I(h) at 1 bar for seven water models: SPC/E, TIP4P, TIP4P-Ew, TIP4P/ice, TIP4P/2005, TIP5P, and TIP5P-E. The results for the melting temperature from the direct coexistence simulations of this work are in agreement (within the statistical uncertainty) with those obtained previously by us from free energy calculations. By taking into account the results of this work and those of our free energy calculations, recommended values of the melting point of ice I(h) at 1 bar for the above mentioned water models are provided.     

Melting and superheating of sI methane hydrate: molecular dynamics study

Melting and decay of the superheated sI methane structure are studied using molecular dynamics simulation. The melting curve is calculated by the direct coexistence simulations in a wide range of pressures up to 5000 bar for the SPC/E, TIP4P/2005 and TIP4P/Ice water models and the united-atom model for methane. We locate the kinetic stability boundary of the superheated metastable sI structure that is found to be surprisingly high comparing with the predictions based on the classical nucleation theory.     

Simulated surface tensions of common water models

Initial simulated values of the surface tension for the SPC/E water model have indicated excellent agreement with experiment. More recently, differing values have been obtained which are significantly lower than previous estimates. Here, we attempt to explain the differences between the previous studies and show that a variety of simulation conditions can affect the final surface tension values. Consistent values for the surface tensions of six common fixed charge water models (TIP3P, SPC, SPC/E, TIP4P, TIP5P, and TIP6P) are then determined for four temperatures between 275 and 350 K. The SPC/E and TIP6P models provide the best agreement with experiment.     

Prediction of self-diffusion coefficient and shear viscosity of water and its binary mixtures with methanol and ethanol by molecular simulation

Density, self-diffusion coefficient, and shear viscosity of pure liquid water are predicted for temperatures between 280 and 373 K by molecular dynamics simulation and the Green-Kubo method. Four different rigid nonpolarizable water models are assessed: SPC, SPC/E, TIP4P, and TIP4P/2005. The pressure dependence of the self-diffusion coefficient and the shear viscosity for pure liquid water is also calculated and the anomalous behavior of these properties is qualitatively well predicted. Furthermore, transport properties as well as excess volume and excess enthalpy of aqueous binary mixtures containing methanol or ethanol, based on the SPC/E and TIP4P/2005 water models, are calculated. Under the tested conditions, the TIP4P/2005 model gives the best quantitative and qualitative agreement with experiments for the regarded transport properties. The deviations from experimental data are of 5% to 15% for pure liquid water and 5% to 20% for the water + alcohol mixtures. Moreover, the center of mass power spectrum of water as well as the investigated mixtures are analyzed and the hydrogen-bonding structure is discussed for different states.     

A comparison of the value of viscosity for several water models using Poiseuille flow in a nano-channel

The viscosity-temperature relation is determined for the water models SPC/E, TIP4P, TIP4P/Ew, and TIP4P/2005 by considering Poiseuille flow inside a nano-channel using molecular dynamics. The viscosity is determined by fitting the resulting velocity profile (away from the walls) to the continuum solution for a Newtonian fluid and then compared to experimental values. The results show that the TIP4P/2005 model gives the best prediction of the viscosity for the complete range of temperatures for liquid water, and thus it is the preferred water model of these considered here for simulations where the magnitude of viscosity is crucial. On the other hand, with the TIP4P model, the viscosity is severely underpredicted, and overall the model performed worst, whereas the SPC/E and TIP4P/Ew models perform moderately.     

Computer simulation study of metastable ice VII and amorphous phases obtained by its melting

Molecular dynamics simulations of metastable ice VII and cubic ice Ic are carried out in order to examine (1) the ability of commonly used water interaction potentials to reproduce the properties of ices, and (2) the possibility of generating low-density amorphous (LDA) structures by heating ice VII, which is known to transform to LDA at approximately 135 K at normal pressure [S. Klotz, J. M. Besson, G. Hamel, R. J. Nelmes, J. S. Loveday, and W. G. Marshall, Nature (London) 398, 681 (1999)]. We test four simple empirical interaction potentials of water: TIP4P [W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983)], SPC/E [H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. B 91, 6269 (1987)], TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)], and ST2 [F. H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 (1974)]. We have found that TIP5P ice VII melts at 210 K, TIP4P at 90 K, and SPC/E at 70 K. Only TIP5P water after transition has a structure similar to that of LDA. TIP4P and SPC/E have almost identical structures, dissimilar to any known water or amorphous phases, but upon heating both slowly evolve towards LDA-like structure. ST2 ice VII is remarkably stable up to 430 K. TIP4P and SPC/E predict correctly the cubic ice collapse into a high-density amorphous ice (HDA) at approximately 1 GPa whereas TIP5P remains stable up to approximately 5 GPa. The densities of the simulated ice phases differ significantly, depending on the potential used, and are generally higher than experimental values. The importance of proper treatment of long-range electrostatic interactions is also discussed.     

The importance of polarizability in the modeling of solubility: quantifying the effect of solute polarizability on the solubility of small nonpolar solutes in popular models of water

In recent work by Paschek [J. Chem. Phys. 120, 6674 (2004)] and others [see H. Docherty et al., J. Chem. Phys. 125, 074510 (2006) for a review] it has been suggested that, when coupled to a simple Lennard-Jones model for various small nonpolar solute molecules, the most common models of water (e.g., SPC/E and TIP4P) fail to reproduce quantitatively the solubility of small nonpolar solute molecules in water due in part to failing to account for polarization of the solute molecule. Given the importance of such systems as test-case prototype models of the solubility of proteins and biomolecules, in this work, we investigate the impact of using a polarizable solute model with the SPC/E, TIP3P, TIP4P, TIP4P-Ew, and TIP4P/2005 rigid water models. Specifically we consider Ne, Ar, Kr, Xe, and methane as solutes. In all cases we observe that the use of a polarizable solute improves agreement between experiment and simulations, with the best agreement seen for the largest solutes, Kr, CH(4), and Xe and the modern reparametrizations of the TIP4P model, i.e., the TIP4P-Ew and TIP4P/2005 models.     

Thermodynamic properties of methane/water interface predicted by molecular dynamics simulations

Molecular dynamics simulations have been performed to examine the thermodynamic properties of methane/water interface using two different water models, the TIP4P/2005 and SPC/E, and two sets of combining rules. The density profiles, interfacial tensions, surface excesses, surface pressures, and coexisting densities are calculated over a wide range of pressure conditions. The TIP4P/2005 water model was used, with an optimized combining rule between water and methane fit to the solubility, to provide good predictions of interfacial properties. The use of the infinite dilution approximation to calculate the surface excesses from the interfacial tensions is examined comparing the surface pressures obtained by different approaches. It is shown that both the change of methane solubilities in pressure and position of maximum methane density profile at the interface are independent of pressure up to about 2 MPa. We have also calculated the adsorption enthalpies and entropies to describe the temperature dependency of the adsorption.     

Monte Carlo calculation of the thermodynamic properties of water

The Monte Carlo method and parallel computing are used to calculate the thermodynamic properties of water (density, heat capacity, compressibility, thermal expansion coefficient, and static dielectric constant) in a wide range of temperatures (from 70 K to 530 K) at constant (atmospheric) pressure. Four groups of computational experiments are carried out, each for its own model of the water molecule: TIP3P (Jorgensen et al., 1983), SPC/E (Berendsen et al., 1987), TIP4P/2005 (Abascal&Vega, 2005), and TIP5P-E (Rick, 2004). An additional calculation based on the replica exchange method is conducted for the TIP4P/2005 model. A comparison of the calculated properties of water with experimental data suggests that the TIP4P/2005 model can provide highly realistic computer simulation results for water and aqueous solutions.     

Structural properties and thermodynamics of water clusters: a Wang-Landau study

The temperature dependence of structural properties and thermodynamic behavior of water clusters has been studied using Wang-Landau sampling. Four potential models, simple point charge/extended (SPC/E), transferable intermolecular potential 3 point (TIP3P), transferable intermolecular potential 4 point (TIP4P), and Gaussian charge polarizable (GCP), are compared for ground states and properties at finite temperatures. Although the hydrogen bond energy and the distance of the nearest-neighbor oxygen pair are significantly different for TIP4P and GCP models, they approach to similar ground state structures and melting transition temperatures in cluster sizes we considered. Comparing with TIP3P, SPC/E model provides properties closer to that of TIP4P and GCP.     

Higher-order virial coefficients of water models

We use the Mayer sampling method, with both direct and overlap sampling, to calculate and compare classical virial coefficients up to B6 for various water models (SPC, SPC/E, MSPC/E, TIP3P, and TIP4P). The precision of the computed values ranges from 0.1% for B2 to an average of 25% for B6. When expressed in a form scaled by the critical properties, the values of the coefficients for SPC water are observed to greatly exceed the magnitude of corresponding coefficients for the simple Lennard-Jones model. We examine the coefficients in the context of the equation of state and the Joule-Thomson coefficient. Comparisons of these properties are made both to established molecular simulation data for each respective model and to real water. For all models, the virial series up to B5 describes the equation of state along the saturated vapor line better than the series that includes B6. At supercritical temperatures, however, the sixth-order series often describes pressure-volume-temperature behavior better than the fifth-order series. For example, the sixth-order virial equation of state for SPC/E water predicts the 673 K isotherm within 8% of published molecular simulation values up to a density of 9 mol/L (roughly half the critical density of SPC/E water).     

Molecular dynamics simulations of ice nucleation by electric fields

Molecular dynamics simulations are used to investigate heterogeneous ice nucleation in model systems where an electric field acts on water molecules within 10-20 ? of a surface. Two different water models (the six-site and TIP4P/Ice models) are considered, and in both cases, it is shown that a surface field can serve as a very effective ice nucleation catalyst in supercooled water. Ice with a ferroelectric cubic structure nucleates near the surface, and dipole disordered cubic ice grows outward from the surface layer. We examine the influences of temperature and two important field parameters, the field strength and distance from the surface over which it acts, on the ice nucleation process. For the six-site model, the highest temperature where we observe field-induced ice nucleation is 280 K, and for TIP4P/Ice 270 K (note that the estimated normal freezing points of the six-site and TIP4P/Ice models are ～289 and ～270 K, respectively). The minimum electric field strength required to nucleate ice depends a little on how far the field extends from the surface. If it extends 20 ?, then a field strength of 1.5 × 10(9) V/m is effective for both models. If the field extent is 10 ?, then stronger fields are required (2.5 × 10(9) V/m for TIP4P/Ice and 3.5 × 10(9) V/m for the six-site model). Our results demonstrate that fields of realistic strength, that act only over a narrow surface region, can effectively nucleate ice at temperatures not far below the freezing point. This further supports the possibility that local electric fields can be a significant factor influencing heterogeneous ice nucleation in physical situations. We would expect this to be especially relevant for ice nuclei with very rough surfaces where one would expect local fields of varying strength and direction.     

Sources of the deficiencies in the popular SPC/E and TIP3P models of water

Motivated by the results of Vega et al. [J. Phys. Condens. Matter 20, 153101 (2008)] about the phase diagram of water, and by the results of Kiss and Baranyai [J. Chem. Phys. 131, 204310 (2009)] about the properties of gas-phase clusters, we carried out a comparative study of the structure modeled by SPC∕E and TIP3P interactions in ambient liquid water. The gas-phase clusters of SPC∕E and TIP3P models show erroneous structures, while TIP4P-type models, either polarizable or not, provide qualitatively correct results. The trimers of SPC∕E and TIP3P are planar in gas phase, contrary to experimental and TIP4P-type models. The aim of this study was to see whether traces of these false geometries characteristic to SPC∕E and TIP3P in gas phase can also be found in the liquid phase. For this purpose we selected trimers formed by adjacent neighbors of water molecules in the liquid and calculated their geometrical features. We determined angles formed by the HO bonds of the molecules with OO vectors and with the normal vector of the OOO plane in the selected trimers. Our results showed that, despite high temperature, the SPC∕E and TIP3P water contains larger number of planar arrangements than other TIP4P-type models. Although structural differences presented in this study are small, they are accurately detectable. These results weaken the reliability of studies obtained by the SPC∕E or TIP3P models even in the liquid phase.