The performance of several popular water models (TIP3P, TIP4P, TIP5P and TIP4P/2005) is analyzed. For that purpose the predictions for ten different properties of water are investigated, namely: 1. vapour-liquid equilibria (VLE) and critical temperature; 2. surface tension; 3. densities of the different solid structures of water (ices); 4. phase diagram; 5. melting-point properties; 6. maximum in the density of water at room pressure and thermal coefficients alpha and KT; 7. structure of liquid water and ice; 8. equation of state at high pressures; 9. self-diffusion coefficient; 10. dielectric constant. For each property, the performance of each model is analyzed in detail with a critical discussion of the possible reason of the success or failure of the model. A final judgement on the quality of these models is provided. TIP4P/2005 provides the best description of almost all properties of the list, the only exception being the dielectric constant. In second position, TIP5P and TIP4P yield a similar performance overall, and the last place with the poorest description of the water properties is provided by TIP3P. The ideas leading to the proposal and design of the TIP4P/2005 are also discussed in detail. TIP4P/2005 is probably close to the best description of water that can be achieved with a non-polarizable model described by a single Lennard-Jones (LJ) site and three charges.

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 Celsius for the basal plane, at about -80 Celsius for the primary prismatic plane and at about -70 Celsius for the secondary prismatic plane.

Molecular dynamics simulations have been performed to estimate the three-phase (solid hydrate-liquid water-gaseous methane) coexistence line for the water-methane binary mixture. The temperature at which the three phases are in equilibrium was determined for three different pressures, namely, 40, 100, and 400 bar by using direct coexistence simulations. In the simulations water was described by using either TIP4P, TIP4P/2005, or TIP4P/Ice models and methane was described as simple Lennard-Jones interaction site. Lorentz-Berthelot combining rules were used to obtain the parameters of the cross interactions. For the TIP4P/2005 model positive deviations from the energetic Lorentz-Berthelot rule were also considered to indirectly account for the polarization of methane when introduced in liquid water. To locate the three-phase coexistence point, two different global compositions were used, which yielded (to within statistical uncertainty) the same predictions for the three-phase coexistence temperatures, although with a somewhat different time evolution. The three-phase coexistence temperatures obtained at different pressures when using the TIP4P/Ice model of water were in agreement with the experimental results. The main reason for this is that the TIP4P/Ice model reproduces the melting point of ice I(h).

In this work the high pressure region of the phase diagram of water has been studied by computer simulation by using the TIP4P/2005 model of water. Free energy calculations were performed for ices VII and VIII and for the fluid phase to determine the melting curve of these ices. In addition, molecular dynamics simulations were performed at high temperatures (440 K) observing the spontaneous freezing of the liquid into a solid phase at pressures of about 80,000 bar. The analysis of the structure obtained lead to the conclusion that a plastic crystal phase was formed. In the plastic crystal phase the oxygen atoms were arranged forming a body center cubic structure, as in ice VII, but the water molecules were able to rotate almost freely. Free energy calculations were performed for this new phase, and it was found that for TIP4P/2005 this plastic crystal phase is thermodynamically stable with respect to ices VII and VIII for temperatures higher than about 400 K, although the precise value depends on the pressure. By using Gibbs-Duhem simulations, all coexistence lines were determined, and the phase diagram of the TIP4P/2005 model was obtained, including ices VIII and VII and the new plastic crystal phase. The TIP4P/2005 model is able to describe qualitatively the phase diagram of water. It would be of interest to study if such a plastic crystal phase does indeed exist for real water. The nearly spherical shape of water makes possible the formation of a plastic crystal phase at high temperatures. The formation of a plastic crystal phase at high temperatures (with a bcc arrangements of oxygen atoms) is fast from a kinetic point of view occurring in about 2 ns. This is in contrast to the nucleation of ice Ih which requires simulations of the order of hundreds of ns.

The phase diagram of water at negative pressures as obtained from computer simulations for two models of water, TIP4P/2005 and TIP5P is presented. Several solid structures with lower densities than ice Ih, so-called virtual ices, were considered as possible candidates to occupy the negative pressure region of the phase diagram of water. In particular the empty hydrate structures sI, sII, and sH and another, recently proposed, low-density ice structure. The relative stabilities of these structures at 0 K was determined using empirical water potentials and density functional theory calculations. By performing free energy calculations and Gibbs-Duhem integration the phase diagram of TIP4P/2005 was determined at negative pressures. The empty hydrates sII and sH appear to be the stable solid phases of water at negative pressures. The phase boundary between ice Ih and sII clathrate occurs at moderate negative pressures, while at large negative pressures sH becomes the most stable phase. This behavior is in reasonable agreement with what is observed in density functional theory calculations.

The recently proposed Einstein molecule approach is extended to compute the free energy of molecular solids. This method is a variant of the Einstein crystal method of Frenkel and Ladd[J.Chem. Phys. 81, 3188 (1984)]. In order to show its applicability, we have computed the free energy of a hard-dumbbells solid, of two recently discovered solid phases of water, namely, ice XIII and ice XIV, where the interactions between water molecules are described by the rigid non-polarizable TIP4P/2005 model potential, and of several solid phases that are thermodynamically stable for an anisotropic patchy model with octahedral symmetry which mimics proteins. Our calculations show that both the Einstein crystal method and the Einstein molecule approach yield the same results within statistical uncertainty. In addition, we have studied in detail some subtle issues concerning the calculation of the free energy of molecular solids. First, for solids with non-cubic symmetry, we have studied the effect of the shape of the simulation box on the free energy. Our results show that the equilibrium shape of the simulation box must be used to compute the free energy in order to avoid the appearance of artificial stress in the system that will result in an increase of the free energy. In complex solids, such as the solid phases of water, another difficulty is related to the choice of the reference structure. As in some cases there is not an obvious orientation of the molecules, it is not clear how to generate the reference structure. Our results will show that, as long as the structure is not too far from the equilibrium structure, the calculated free energy is invariant to the reference structure used in the free energy calculations. Finally, the strong size dependence of the free energy of solids is also studied.

In this note we present results for the heat capacity at constant pressure for the TIP4PQ/2005 model, as obtained from path integral simulations. The model does a rather good job of describing both the heat capacity of ice I h and of liquid water. Classical simulations using the TIP4P/2005, TIP3P, TIP4P, TIP4P-Ew, SPC/E and TIP5P models are unable to reproduce the heat capacity of water. Given that classical simulations do not satisfy the third law of thermodynamics, one would expect such a failure at low temperatures. However, it seems that for water, nuclear quantum effects influence the heat capacities all the way up to room temperature. The failure of classical simulations to reproduce C p points to the the necessity of incorporating nuclear quantum effects to describe this property accurately.

The three phase equilibrium line (hydrate-liquid water-liquid carbon dioxide) has been estimated for the water + carbon dioxide binary mixture using molecular dynamics simulation and the direct coexistence technique. Both molecules have been represented using rigid nonpolarizable models. TIP4P/2005 and TIP4P/Ice were used for the case of water, while carbon dioxide was considered as a three center linear molecule with the parameterizations of MSM, EPM2, TraPPE, and ZD. The influence of the initial guest occupancy fraction on the hydrate stability has been analyzed first in order to determine the optimal starting configuration for the simulations, paying attention to the influence of the two different cells existing in the sI hydrate structure. The three phase coexistence temperature was then determined for a pressure range from 2 to 500 MPa. The qualitative shape of the equilibrium curve estimated is correct, including the high pressure temperature maximum that determines the hydrate re-entrant behaviour. However, in order to obtain quantitative agreement with experimental results, a positive deviation from the classical Lorentz-Berthelot combining rules must be considered.