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 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.     

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.     

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 points and thermal expansivities of proton-disordered hexagonal ice with several model potentials

A method of free energy calculation is proposed, which enables to cover a wide range of pressure and temperature. The free energies of proton-disordered hexagonal ice (ice Ih) and liquid water are calculated for the TIP4P [J. Chem. Phys. 79, 926 (1983)] model and the TIP5P [J. Chem. Phys. 112, 8910 (2000)] model. From the calculated free energy curves, we determine the melting point of the proton-disordered hexagonal ice at 0.1 MPa (atmospheric pressure), 50 MPa, 100 MPa, and 200 MPa. The melting temperatures at atmospheric pressure for the TIP4P ice and the TIP5P ice are found to be about T(m)=229 K and T(m)=268 K, respectively. The melting temperatures decrease as the pressure is increased, a feature consistent with the pressure dependence of the melting point for realistic proton-disordered hexagonal ice. We also calculate the thermal expansivity of the model ices. Negative thermal expansivity is observed at the low temperature region for the TIP4P ice, but not for the TIP5P ice at the ambient pressure.     

Kinetic aspects of the thermostatted growth of ice from supercooled water in simulations

In experiments, the growth rate of ice from supercooled water is seen to increase with the degree of supercooling, that is, the lower the temperature, the faster the crystallization takes place. In molecular dynamics simulations of the freezing process, however, the temperature is usually kept constant by means of a thermostat that artificially removes the heat released during the crystallization by scaling the velocities of the particles. This direct removal of energy from the system replaces a more realistic heat-conduction mechanism and is believed to be responsible for the curious observation that the thermostatted ice growth proceeds fastest near the melting point and more slowly at lower temperatures, which is exactly opposite to the experimental findings [M. A. Carignano, P. B. Shepson, and I. Szleifer, Mol. Phys. 103, 2957 (2005)]. This trend is explained by the diffusion and the reorientation of molecules in the liquid becoming the rate-determining steps for the crystal growth, both of which are slower at low temperatures. Yet, for a different set of simulations, a kinetic behavior analogous to the experimental finding has been reported [H. Nada and Y. Furukawa, J. Crystal Growth 283, 242 (2005)]. To clarify this apparent contradiction, we perform relatively long simulations of the TIP4P/Ice model in an extended range of temperatures. The temperature dependence of the thermostatted ice growth is seen to be more complex than was previously reported: The crystallization process is very slow close to the melting point at 270 K, where the thermodynamic driving force for the phase transition is weak. On lowering the temperature, the growth rate initially increases, but displays a maximum near 260 K. At even lower temperatures, the freezing process slows down again due to the reduced diffusivity in the liquid. The velocity of the thermostatted melting process, in contrast, shows a monotonic increase upon raising the temperature beyond the normal melting point. In this case, the effects of the increasing thermodynamic driving force and the faster diffusion at higher temperatures reinforce each other. In the context of this study, we also report data for the diffusion coefficient as a function of temperature for the water models TIP4P/Ice and TIP4P/2005.     

Computer simulation of two new solid phases of water: Ice XIII and ice XIV

NpT Monte Carlo simulations have been performed for two recently discovered solid phases of water which have been denoted as ice XIII and ice XIV C. G. Salzmann et al. [Science311, 1758 (2006)]. Several potential models of water were considered, namely, the traditional SPC/E, TIP4P, and TIP5P and the more recent TIP5P-E, TIP4P-Ew, TIP4P/Ice, and TIP4P/2005 models. Significant differences in density and oxygen-oxygen radial distribution functions are found between the predictions of the SPC/E, TIP5P, and the models of the TIP4P family. The models TIP4P/Ice and TIP4P/2005 provide the best estimates of the density.     

Melting temperature of ice Ih calculated from coexisting solid-liquid phases

We carried out molecular-dynamics simulations by using the two-phase coexistence method with the constant pressure, particle number, and enthalpy ensemble to compute the melting temperature of proton-disordered hexagonal ice I(h) at 1-bar pressure. Four models of water were considered, including the widely used TIP4P [W. L. Jorgensen, J. Chandrasekha, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys.79, 926 (1983)] and TIP5P [M. W. Mahoney and W. L. Jorgensen J. Chem. Phys.112, 8910 (2000)] models, as well as recently improved TIP4P and TIP5P models for use with Ewald techniques-the TIP4P-Ew [W. Horn, W. C. Swope, J. W. Pitera, J. C. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys.120, 9665 (2004)] and TIP5P-Ew [S. W. Rick, J. Chem. Phys.120, 6085 (2004)] models. The calculated melting temperature at 1 bar is T(m) = 229 +/- 1 K for the TIP4P and T(m) = 272.0 +/- 0.6 K for the TIP5P ice I(h), both are consistent with previous simulations based on free-energy methods. For the TIP4P-Ew and TIP5P-Ew models, the calculated melting temperature is T(m) = 257.0 +/- 1.1 K and T(m) = 253.9 +/- 1.1 K, respectively.     

Quantum effects in liquid water and ice: model dependence

This paper explores the influence of choice of potential model on the quantum effects observed in liquid water and ice. This study utilizes standard rigid models and a more formal context for the rigid-body centroid molecular dynamics methodology used to perform the quantum simulations is provided. Quantum and classical molecular dynamics simulations are carried out for liquid water and ice Ih at 298 and 220 K, respectively, with the simple point charge/extended and TIP4P-Ew water models. The results obtained for equilibrium and dynamical properties are compared with those recently reported on TIP4P [L. Hernandez de la Pena and P. G. Kusalik, J. Chem. Phys. 121, 5992 (2004); L. Hernandez de la Pena et al., J. Chem. Phys 123, 144506 (2005)]. For the liquid, an energy shift of about 8% and an average molecular uncertainty of about 11 degrees were found independently of the water model. The self-diffusion coefficient consistently increases by more than 50% when going from the classical to the quantum system and quantum dynamics are found to reproduce the experimental isotopic shifts with the models examined. The ice results compare remarkably well with those previously reported for the TIP4P water model; they confirm that quantum effects are considerable and that the quantum mechanical uncertainty and the energy shifts due to quantization are smaller in ice than in liquid water. The relevance of these findings in the context of the construction of water models is briefly discussed.     

Predicting the melting temperature of ice-Ih with only electronic structure information as input

The melting temperature of ice-Ih was calculated with only electronic structure information as input by creating a problem-specific force field. The force field, Water model by AFM for Ice and Liquid (WAIL), was developed with the adaptive force matching (AFM) method by fitting to post-Hartree-Fock quality forces obtained in quantum mechanics∕molecular mechanics calculations. WAIL predicts the ice-Ih melting temperature to be 270 K. The model also predicts the densities of ice and water, the temperature of maximum density of water, the heat of vaporizations, and the radial distribution functions for both ice and water in good agreement with experimental measurements. The non-dissociative WAIL model is very similar to a flexible version of the popular TIP4P potential and has comparable computational cost. By customizing to problem-specific configurations with the AFM approach, the resulting model is remarkably more accurate than any variants of TIP4P for simulating ice-Ih and water in the temperature range from 253 K and 293 K under ambient pressure.     

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.     

The ice/water interface

The TIP4P model of water is used to simulate the ice/water interface. A stable interface is constructed with the basal plane of ice 1h exposed to liquid water. Density and orientation profiles, and diffusion constants, are calculated as a function of the distance perpendicular to the interface.     

A general purpose model for the condensed phases of water: TIP4P/2005

A potential model intended to be a general purpose model for the condensed phases of water is presented. TIP4P/2005 is a rigid four site model which consists of three fixed point charges and one Lennard-Jones center. The parametrization has been based on a fit of the temperature of maximum density (indirectly estimated from the melting point of hexagonal ice), the stability of several ice polymorphs and other commonly used target quantities. The calculated properties include a variety of thermodynamic properties of the liquid and solid phases, the phase diagram involving condensed phases, properties at melting and vaporization, dielectric constant, pair distribution function, and self-diffusion coefficient. These properties cover a temperature range from 123 to 573 K and pressures up to 40,000 bar. The model gives an impressive performance for this variety of properties and thermodynamic conditions. For example, it gives excellent predictions for the densities at 1 bar with a maximum density at 278 K and an averaged difference with experiment of 7 x 10(-4) g/cm3.     

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.     

A potential model for the study of ices and amorphous water: TIP4P/Ice

The ability of several water models to predict the properties of ices is discussed. The emphasis is put on the results for the densities and the coexistence curves between the different ice forms. It is concluded that none of the most commonly used rigid models is satisfactory. A new model specifically designed to cope with solid-phase properties is proposed. The parameters have been obtained by fitting the equation of state and selected points of the melting lines and of the coexistence lines involving different ice forms. The phase diagram is then calculated for the new potential. The predicted melting temperature of hexagonal ice (Ih) at 1 bar is 272.2 K. This excellent value does not imply a deterioration of the rest of the properties. In fact, the predictions for both the densities and the coexistence curves are better than for TIP4P, which previously yielded the best estimations of the ice properties.     

Anisotropy in the crystal growth of hexagonal ice, I(h)

Growth of ice crystals has attracted attention because ice and water are ubiquitous in the environment and play critical roles in natural processes. Hexagonal ice, I(h), is the most common form of ice among 15 known crystalline phases of ice. In this work we report the results of an extensive and systematic molecular dynamics study of the temperature dependence of the crystal growth on the three primary crystal faces of hexagonal ice, the basal {0001} face, the prism {1010} face, and the secondary prism {1120} face, utilizing the TIP4P-2005 water model. New insights into the nature of its anisotropic growth are uncovered. It is demonstrated that the ice growth is indeed anisotropic; the growth and melting of the basal face are the slowest of the three faces, its maximum growth rates being 31% and 43% slower, respectively, than those of the prism and the secondary prism faces. It is also shown that application of periodic boundary conditions can lead to varying size effect for different orientations of an ice crystal caused by the anisotropic physical properties of the crystal, and results in measurably different thermodynamic melting temperatures in three systems of similar, yet moderate, size. Evidence obtained here provides the grounds on which to clarify the current understanding of ice growth on the secondary prism face of ice. We also revisit the effect of the integration time step on the crystal growth of ice in a more thorough and systematic way. Careful evaluation demonstrates that increasing the integration time step size measurably affects the free energy of the bulk phases and shifts the temperature dependence of the growth rate curve to lower temperatures by approximately 1 K when the step is changed from 1 fs to 2 fs, and by 3 K when 3 fs steps are used. A thorough investigation of the numerical aspects of the simulations exposes important consequences of the simulation parameter choices upon the delicate dynamic balance that is involved in ice crystal growth.