Phase diagram of the TIP4P/Ice water model by enhanced sampling simulations

The phase diagram of water is diverse and includes one liquid phase, one vapor phase, and 20 known ice polymorphs^1–3 as well as hypothesized low- and high-density phases of supercooled water.^4–6 Molecular simulations are an important tool to obtain insight into the phase behavior of water on scales and conditions inaccessible to experiments.⁶ Furthermore, the ability to calculate the phase diagram using molecular simulations and compare it with its experimental counterpart provides a stringent test for the capability of a given molecular model to represent water accurately across different thermodynamic conditions. For these reasons, the development of methods to accurately determine the phase diagram of water has received considerable attention over the years.

Many methods have been employed to study phase equilibria involving liquid water and ice. One of the most common approaches is thermodynamic integration, which calculates free energies by performing a series of simulations that link a phase of known free energy to the phase of interest. For crystalline solids, two well-known flavors of thermodynamic integration are the Einstein crystal method developed by Frenkel and Ladd⁷ and the Einstein molecule variant proposed by Vega and Noya.^8,9 Another simple, yet powerful approach is the direct coexistence technique, where melting points are determined either by a series of NPT simulations or by the equilibrium temperature of a single coexistence simulation within the NPH ensemble.¹⁰ The methods described above are often complemented by the Gibbs–Duhem integration technique developed by Kofke.¹¹ In this method, the Clausius–Clapeyron equation is integrated starting from a known coexistence point, and a continuous and smooth coexistence curve can be traced from it. In a seminal work, Sanz et al. used the Einstein molecule method and Gibbs–Duhem integration to calculate a large region of the phase diagram of different empirical water models that provided an unprecedented insight into the ability of water models to reproduce the experimental phase diagram.¹² Very recently, a similar approach has also been used to calculate the phase diagram of machine learning models trained on energies and forces derived from ab initio calculations. For instance, Zhang et al. studied the phase diagram of a water model derived from density-functional theory (DFT) calculations based on the SCAN functional,¹³ and Reinhardt and Cheng followed a similar procedure to calculate the phase diagram of a water model based on a hybrid density functional with van der Waals corrections.¹⁴ However, the application of the Einstein crystal and Einstein molecule methods is not devoid of challenges. An important feature of many ice phases is proton disorder, i.e., the existence of many proton configurations compatible with oxygen atoms in their regular crystallographic positions. A delicate approximation is typically made in the Einstein crystal and Einstein molecule methods in order to account for the entropic contribution of proton disorder. This approximation is critical in the cases of partially disordered ice phases, such as ice III and V.¹⁵ Indeed, for the TIP4P/2005 model, a dramatic increase in the region of stability of ice III and an increase of about 25 K in the melting temperature were observed when the melting line was computed using the direct coexistence technique rather than the Einstein molecule method.¹⁶ 

Enhanced sampling methods provide an alternative route to calculate the free energies of solids. These methods are based on introducing a bias potential along with a suitably defined collective variable. For many years, the lack of structure-specific collective variables precluded the application of enhanced sampling to the accurate calculation of free energies of solids. Recently, two different approaches have been proposed to determine structure-specific collective variables. One is based on using particular peaks of the structure factor,^17,18 and the other relies on a similarity kernel (environment similarity) in order to compare environments in the simulation to reference environments of the target structure.^19,20 In the interface pinning method¹⁸ introduced by Pedersen et al., one considers a liquid–solid coexistence configuration and uses an umbrella bias potential, a function of a structure factor peak, to avoid growth or melting. From such a simulation one can extract liquid–solid chemical potential differences. This method was applied by Cheng et al. to calculate the chemical potential of ice I_h.²¹ Another approach is based on performing a simulation in which the bulk liquid and the bulk solid are reversibly interconverted. This alternative, together with the “environment similarity” collective variable, was used by Piaggi and Car to determine the chemical potentials of ice I_h.²² As in the case of the direct-coexistence method, the entropic contribution of proton disorder is automatically included in this case due to the explicit simulation of the crystallization process,²² thus providing a key advantage over thermodynamic integration. The environment similarity collective variable was also recently used to calculate the ice Ih–liquid interfacial free energy.²³ Until now, however, enhanced sampling methods have not been applied to the calculation of complex phase diagrams that involve multiple crystal polymorphs.

Here, we study the phase diagram of the TIP4P/Ice water model using enhanced sampling molecular dynamics (MD) simulations. We employ a method based on simulating the coexistence between an ice polymorph and liquid water through the application of a bias potential that allows for reversibly sampling the melting and growth of ice layers. We target five different ice polymorphs, namely, ice I_h, II, III, V, and VI, with appropriately chosen “environment similarity” collective variables. This procedure allows us to calculate ice–liquid water chemical potential differences and, thus, determine the corresponding coexistence points. We subsequently employ Gibbs–Duhem integration to trace continuous coexistence lines starting from these points.

Our goal is to calculate the melting lines of ice I_h, II, III, V, and VI. To this end, we follow a protocol based on calculating a small number of points along the melting line of each polymorph which are then used within the Gibbs–Duhem integration method to trace the corresponding melting lines. In Sec. II A, we describe the biased coexistence approach that we use to calculate melting points, and in Sec. II B, we provide specific details about our implementation of the Gibbs–Duhem integration method.

We consider a system in which liquid water is in direct coexistence with an ice polymorph. In a conventional coexistence simulation, if the temperature is below the melting temperature, the crystal grows irreversibly at the expense of the liquid, and if the temperature is above the melting temperature, the crystal melts irreversibly. For our biased coexistence simulations, we take inspiration from the interface pinning technique¹⁸ and introduce a bias potential such that the number of ice-like molecules in the system reversibly increases and decreases by the number of molecules in one ice layer. This procedure allows us to compute chemical potentials and thus coexistence points. As discussed in Sec. II A 1, a key element in these simulations is the use of order parameters that are able to drive the crystallization of each ice polymorph.

To construct our order parameters, we consider the atomic environments within a prescribed cutoff r_c of each water molecule in a given ice polymorph. The number m of distinct environments is equal to the number of molecules in the basis of the crystal structure. The environments X = χ₁, …, χ_m are the starting point for the construction of the order parameter of each polymorph and are shown in Fig. 1 for ice I_h, II, III, V, and VI. The environments are centered at the oxygen atoms, and only the oxygen atoms are considered as neighbors. Protons are not included in the environments in order for the proton configuration to form spontaneously during the simulation without any bias toward a particular configuration.²⁴ The number m of distinct environments and the cutoff r_c for each ice polymorph is summarized in Table I.

The “environment similarity” kernel^19,20 is then used to compare the atomic environments χ_l ∈ X with a generic environment χ,

where $ρχl(r)$ and ρ_χ(r) are the atomic densities corresponding to the environments χ_l and χ, respectively. If the densities are represented by the sums of Gaussians centered at the neighbors’ positions with spread σ, the kernel becomes

where n is the number of neighbors in the environment χ_l and $ril$ and r_j are the positions of the neighbors in environments χ_l and χ, respectively. In Eq. (2), we have added a normalization such that $kχl(χl)=1$⁠. The m similarity kernels derived from Eq. (2) allow us to identify if a given environment is compatible with one of the environments of an ice polymorph. However, since it is preferred to have a single similarity measure between a given environment and any of the m reference environments of a polymorph, we define the best-match kernel

that compares the environments χ with all the environments X of a given polymorph. k_X(χ) as defined in Eq. (3) is not continuous and not differentiable and cannot be used directly in enhanced sampling MD simulations. For this reason, we use a smooth approximation to the maximum function,

where λ has to be set to a large value. In our simulations, we employed λ = 100.

Since there is one value of k_X(χ) per oxygen atom, any given bulk configuration has a distribution of this quantity. Figure 2 shows the distributions of k_X(χ) in liquid water and each ice polymorph. Since the distributions have very small overlap, k_X(χ) can be used to distinguish liquid water from any ice polymorph. We have chosen the spread of the Gaussians σ such that the overlap between the two distributions is minimized (the overlap between distributions as a function of σ is shown in Fig. S1 of the supplementary material).

The similarity kernel defined in Eq. (3) provides a way to characterize the environments in a given configuration as being compatible with the environments in a given ice polymorph. For a system of N water molecules, there are N oxygen–oxygen environments χ¹, χ², …, χ^N. We, thus, define a global order parameter n_ice that represents the number of environments consistent with an ice polymorph,

Here, κ is a watershed between values of k_X(χⁱ) consistent with the liquid and those consistent with the solid. Appropriate values of σ and κ for each ice polymorph are shown in Table I. The order parameter defined in Eq. (5) can be made continuous and differentiable using

where f is a switching function that is $∼0$ and $∼1$ for k_X(χⁱ) values consistent with the liquid and the solid, respectively. In particular, we choose

where y = (k_X(χⁱ) − k₁)/(k₂ − k₁), and k₁ and k₂ can be chosen as the peaks of the liquid and ice distributions of k_X(χⁱ), respectively.

In order to build a bias potential, we use the on-the-fly probability enhanced sampling (OPES) algorithm^25,26 that shares many features with the well-known metadynamics technique.²⁷ The bias potential is a function of n_ice and targets a multiumbrella ensemble²⁵ that produces an approximately uniform distribution in an interval from $nicelow$ to $nicehigh$⁠. The limits of the interval are chosen as integer multiples of the number of molecules in a crystal layer, Δn_layer, i.e., $nicelow=NlΔnlayer$ and $nicehigh=(Nl+1)Δnlayer$ with $Nl∈N$⁠. The functional form of the bias potential is

where β is the inverse temperature and $niceλ$ are N_λ points uniformly spaced in the interval $nicelow$ to $nicehigh$⁠. In the simulations, the spacing between $niceλ$ is one molecule. ΔF_λ is the difference in free energy between the unbiased system and the system with an umbrella potential at $niceλ$ and can be determined iteratively as described in Ref. 25. A schematic plot of the target probability, bias potential, and free energy in the multiumbrella OPES method is shown in Fig. 3. We note that other methods, such as umbrella sampling or metadynamics, can be used instead of OPES to construct the bias potential.

In Fig. 4(a), we report an example of the time evolution of the number of ice-like molecules for ice I_h at 270 K in one of our biased simulations. After an initial transient period during which a suitable bias potential is determined, the simulation reversibly samples the full targeted range of ice-like molecules. This is in stark contrast with the behavior of an unbiased, direct coexistence simulation, which will have a single transition that either melts or freezes the sample. As such, the biased simulation can be continued as long as needed, providing us with the statistics on n_ice required to estimate the difference in free energy between the ice polymorph and the liquid.

Once the bias potential is converged, the equilibrium properties of the system can be calculated through the reweighting technique. In particular, the unbiased probability of observing an ice slab composed of n_ice molecules is

where δ is the Dirac delta function and ⟨·⟩₀ denotes an average over the unbiased ensemble (e.g., NPT) and n_ice(R) is the number of ice-like molecules in a configuration characterized by atomic coordinates R. This probability can be computed from the biased simulations using the formula

where ⟨·⟩_V is an average over the biased ensemble and V is the bias potential. In practice, the application of Eq. (10) corresponds to constructing a histogram with weights e^βV. The free-energy surface for n_ice is obtained simply by using the relation

In Fig. 4(b), we show an example of free-energy profiles for ice I_h–liquid coexistence at atmospheric pressure. All free-energy profiles show a linear behavior, and, as expected, the overall trend is for the slope to increase with temperature due to the higher stability of the liquid over ice I_h.

The free energy of ice in direct coexistence with the liquid can be written as

where μ are the chemical potentials, n is the number of molecules, S(n_ice) is the contribution of the surfaces to the free energy, and the subscripts _ice and _liq denote ice and liquid. Using the fact that n_ice + n_liq = N, with N being the total number of water molecules in the system, we can express the free energy as

The behavior of S(n_ice) depends on the characteristics of the specific surface considered and may introduce non-linear behavior in ΔG as a function of n_ice. To extract the chemical potential, we evaluate the slope of ΔG over an entire layer Δn_layer,

In general, the surface term will be periodic over a layer. Using the property S(n_ice + Δn_layer) = S(n_ice) for the surface term, we find

which implies that the slope of ΔG(n_ice) with respect to n_ice over a full layer is equal to the difference in chemical potential between the two phases.

In practice, we determine μ_liq→ice = μ_ice − μ_liq by extracting the slope from ΔG(n_ice) calculated using Eq. (11). The slope is obtained by performing an error-weighted linear least-squares fit. Once the chemical potential differences are known, the melting points T_m are calculated using the condition μ_liq→ice(T_m) = 0 at constant pressure. T_m can be easily obtained by an error-weighted linear least-squares fit of μ_liq→ice vs temperature.

From the estimated coexistence points, we calculate coexistence lines using the Gibbs–Duhem integration.¹¹ This is done by numerically integrating the Clausius–Clapeyron equation,

where Δh and Δv are per-molecule differences in enthalpy and volume between the two phases under consideration, respectively, and T_m is the coexistence temperature at pressure P. We integrate Eq. (16) using the fourth-order Runge–Kutta method,²⁸ evaluating the integrand for each integration step by performing one MD simulation for each phase. While we only calculate liquid–ice coexistence points, we can locate triple points from the intersection between liquid–ice coexistence lines. From each liquid–ice–ice triple point, we then conduct a new Gibbs–Duhem simulation with the two ice phases. We find new triple points from these coexistence lines and repeat the procedure until the whole phase diagram is computed. We note that while the ice II coexistence line is only metastable with liquid water, it is essential as its intersection with the coexistence line of ice I_h provides the last triple point necessary to construct the full phase diagram of all five phases of ice considered here.

The MD simulations were carried out using LAMMPS²⁹ patched with the PLUMED 2 enhanced sampling plugin.³⁰ The input files to reproduce all simulations are available on PLUMED-NEST³¹ as part of a collective effort to improve the transparency and reproducibility of enhanced molecular simulations. We employed a time step of 2 fs for the integration of the equations of motion. The temperature was controlled using the stochastic velocity rescaling algorithm with a 0.1 ps relaxation time.³² The pressure was maintained by the Parrinello–Rahman barostat with a 1 ps relaxation time.³³ The bond lengths and angles were kept fixed using the SHAKE algorithm. A cutoff of 0.85 nm was used for the Lennard-Jones and Coulomb interactions. Long-range Coulomb interactions beyond this cutoff were computed with the particle–particle particle-mesh (PPPM) solver. Tail corrections to the pressure and energy were included to take into account long-range effects neglected due to the Lennard-Jones potential truncation.

The initial configurations for the different ice phases were prepared using the GENICE software.^34,35 This ensures the correct proton order/disorder and zero total dipole moment, while simultaneously respecting the Bernard–Fowler ice rules.³⁶ Starting from these ice structures, the initial coexistence configurations were created as follows: First, the ice box was equilibrated with an anisotropic barostat for 1 ns. This was followed by 300 ps of simulation during which the motion of half of the water molecules on one side of the box was kept frozen, while the water molecules on the other half of the box were melted at 450 K. The whole system was then equilibrated to 270 K and 1.0133 × 10⁻⁴ GPa, 250 K and 0.303 atm, 230 K and 0.303 GPa, 250 K and 0.5066 GPa, and 300 K and 1.0133 GPa for ice I_h, II, III, V, and VI, respectively. Finally, following best practices for ordinary direct coexistence simulations,³⁷ we performed enhanced coexistence simulations in the Np_xT ensemble (fluctuating box dimensions in the x-direction and constant area in the yz-directions), with the box dimensions along the yz-directions set to the average equilibrium values of the ice phase at the thermodynamic condition of the enhanced coexistence simulation.

Following this procedure, we prepared the enhanced coexistence simulations for the systems specified in Table II. In Fig. 5, we show the interfaces employed in our simulations with the oxygen atoms colored using the environment similarity order parameter. We conducted enhanced coexistence simulations using four walkers (four simultaneous simulations with different random seeds that share the same bias potential). We first carried out 0.5 ns of sampling in the Np_xT ensemble with harmonic wall-bias potentials at the outer ice ranges (see Table II). This step allows the coexistence configuration to relax toward the thermodynamically favored outer ice range, which tends to speed up the convergence of the bias potential. After this step, the enhanced sampling within the multiumbrella ensemble was initiated. The environments for the construction of the collective variables were obtained using the Environment Finder software.³⁸ We refer the reader to our digital repository, which includes all the input files to reproduce our simulations, for further details.

For the Gibbs–Duhem integration, we carried out a series of MD simulations, with an equilibration time of 0.125 ns followed by 0.375 ns of production simulation. To reduce the required equilibration time, we used the final configuration of the previous integration step as the starting configuration. For steep coexistence lines in the temperature and pressure plane (negative pressure part of L-I_h; L-III, L-V, L-VI), the integration was carried out over the pressure using an integration step of 0.0253 GPa and 0.0126 GPa for ice II and ice III, respectively. For near-flat coexistence lines in the temperature and pressure plane (positive pressure part of I_h-L, I_h-III, II-I_h, III–V, V–VI), the integration was instead carried out over the temperature with an integration step of 5 K.

We carried out enhanced coexistence simulations at different temperatures along the isobars. The time evolution of the number of ice-like molecules for all simulations is reported in Sec. S2 of the supplementary material but shares the essential features of the example provided in Fig. 4 for ice I_h. Each simulation required different simulation times to reach adequate convergence. In detail, we performed four parallel simulations lasting 5 ns for ice V and 10 ns for ice I_h, II, III, and VI.

The free-energy profiles for all the studied ice polymorphs and thermodynamic conditions are reported in Sec. S3 of the supplementary material. Except for ice VI, we obtain approximately linear free-energy profiles as in the example shown in Fig. 4(b). These free-energy profiles are characteristic of the continuous or rough growth mechanism.³⁹ For the (100) interface of ice VI (see Fig. 6), there is a clear deviation from linearity and we observe a free-energy barrier to form the new layer of ice. This free-energy profile is characteristic of the birth and spread growth mechanism.³⁹ It is possible that this mechanism is observed only for ice VI because the simulation cell is large enough to manifest it. The nearly linear profiles for ice I_h, ice II, ice III, and ice V suggest that one may sample these interfaces on smaller ice ranges than one layer, as done in the interface pinning method.¹⁸ This may reduce the sampling time required to obtain the free-energy profiles and allows for using a single umbrella potential. Nonetheless, by probing the collective variable on a full layer, we obtain information on the overall slope.

From the free-energy profile and error associated with each point, we determined the chemical potential using an error-weighted linear least-squares fit. The estimated chemical potentials for each ice polymorph as a function of temperature follow a linear behavior [see the example in Fig. 4(c) and Sec. S4 of the supplementary material] within a range of ∼25 K.

The melting points of the ice polymorphs at various pressures were calculated from linear fits of the chemical potentials as a function of temperature as described in Sec. II and are listed in Table S1 of the supplementary material. The accuracy of our estimates can be assessed through comparisons with analogous results reported in the literature for TIP4P/Ice (see Table III). Our estimate of the melting point of ice I_h at ambient pressure is 269.1 K, which is in close agreement with previous estimates, deviating by 0.7 K from the value obtained from direct coexistence simulations using very large system sizes.¹⁶ Although additional sampling could improve our estimate, we considered a deviation smaller than 1 K to be sufficient for the purposes of this study. Furthermore, finite-size effects can result in systematic errors of the same order of magnitude.^40,41

Using the liquid–solid coexistence points, we traced the coexistence lines using the Gibbs–Duhem integration method. While a single liquid–solid coexistence point per polymorph provides sufficient information for the Gibbs–Duhem integration, we determined a total of 19 points to check the internal consistency of our estimates. Figure 7(a) shows the coexistence curves constructed starting from the points reported in Table S1 of the supplementary material. Also shown in Fig. 7(a) are the corresponding coexistence lines reported in Refs. 42 and 46. With the exception of ice III, our results are in agreement with those reported in Ref. 42. The discrepancy found for ice III was already observed in Ref. 46, which also reported new results for the liquid–ice III coexistence line that are in agreement with our predictions. A similar discrepancy between thermodynamic integration and direct coexistence estimates of the liquid–ice III coexistence line was also observed for TIP4P/2005.¹⁶ We suggest that the origin of the discrepancies seen in Fig. 7 for the liquid–ice III coexistence line lies in the approximation adopted for the residual entropy of ice III in the Einstein molecule approach used in Ref. 42. We also note that there is excellent agreement between the coexistence points calculated here using biased coexistence simulations and the lines traced using the Gibbs–Duhem integration, demonstrating the internal consistency of our calculations.

We calculated the full phase diagram of TIP4P/Ice using a series of Gibbs–Duhem integration simulations starting from the triple points listed in Table IV. We note that, while the ice II–liquid coexistence line is not present in the equilibrium phase diagram of water, its calculation is essential to obtain the liquid–ice V–ice II triple point. The full phase diagram for TIP4P/Ice is shown in Fig. 7(b). Compared to the phase diagram reported in Ref. 42, we observe a significant increase in the size of the ice III stability region. A similar increase of the ice III stability region was reported for TIP4P/2005 by Conde et al. when the phase diagram was computed using direct-coexistence simulations instead of the Einstein molecule method.¹⁶ 

The topology of the TIP4P/Ice phase diagram predicted by our simulations is in qualitative agreement with the experimental data, with the different phase regions being correctly located with respect to each other. However, the large stability region occupied by ice III significantly reduces the ability of TIP4P/Ice to reproduce the actual phase diagram of water. We suggest that the presence of the large region of stability of ice III in the TIP4P/Ice model is likely due to the explicit fit to the range of existence of ice III during the parameterization of the model.⁴² 

We have calculated the phase diagram of the TIP4P/Ice water model using enhanced coexistence simulations. Our phase diagram differs from the previously published phase diagram by having a substantially larger stability region for ice III. While the phase diagram that we have calculated correctly positions the stability regions of each polymorph with respect to each other, the qualitative agreement with experiments deteriorated with respect to the phase diagram reported in Ref. 42. Despite this shortcoming, the TIP4P/Ice model can still be considered one of the best empirical potentials for water.

In computing the phase diagram of TIP4P/Ice, we have provided a state-of-the-art methodology for determining complex phase diagrams. An important ingredient of our method is an order parameter that is able to drive the growth and melting of all the ice polymorphs studied here. While the Einstein crystal method, requiring about 8 ns for estimating the free energy per thermodynamic state point,^47,48 is somewhat faster than our methodology and direct coexistence, we highlight several features that might provide advantages for our method over existing approaches. First, compared to direct coexistence, our approach allows us not only to find coexistence points but also to calculate differences in chemical potential that are crucial for assessing the stability of phases. Second, the explicit simulation of the growth and melting process is useful to obtain an insight into these processes. For example, from the shape of the free-energy profiles, we detected a birth and spread growth mechanism for the (100) plane of ice VI. Third, our approach allows us to systematically improve our estimates by extending the sampling time and, therefore, increasing the number of instances of creation and destruction of a layer of ice until statistical convergence is reached. This can be an advantage in small systems in which the accuracy of the direct coexistence methodology is limited by the stochasticity of the growth and melting events. However, for small systems, both the direct coexistence technique and our method would suffer from the same systematic finite-size effects in the melting points. Finally, compared to the Einstein crystal method, our approach does not rely on approximations for the residual entropy associated with proton disorder and can thus be used to study both proton-disordered and partially proton-ordered ice polymorphs. As such, our methodology fills a gap between the Einstein crystal method and direct coexistence, by providing both free energy estimates without approximating proton disorder.

We have shown that we are able to study the phase diagram of water with individual simulations lasting no longer than 5–10 ns. Therefore, the methodology is well suited for machine learning models of water trained on DFT calculations and data-driven many-body models such as MB-pol, which has shown a remarkable ability to reproduce the properties of liquid water and ice over a wide range of thermodynamic conditions.^49–56 

See the supplementary material for details about the environment similarity kernel for water and ice polymorphs, additional analyses of the coexistence between water and ice polymorphs, free energy differences, chemical potentials, list of melting points for the ice polymorphs, and Gibbs–Duhem coexistence lines.

We are grateful to Carlos Vega for a discussion on the effect of partial proton ordering on the stability of ice polymorphs. We thank Christoph Salzmann for providing us with the data for drawing the experimental phase diagram. S.L.B. and F.P. were supported by the Air Force Office of Scientific Research under Award No. FA9550-20-1-0351. P.M.P. and R.C. acknowledge support from the center Chemistry in Solution and at Interfaces funded by the U.S. Department of Energy under Award No. DE-SC0019394. Computational resources were provided by the Department of Defense High Performance Computing Modernization Program (HPCMP) and the Triton Shared Computing Cluster (TSCC) at the San Diego Supercomputer Center (SDSC).

The authors have no conflicts to disclose.

MD input, analysis scripts, and key results are publicly available on GitHub. MD trajectories are available from the corresponding author upon request. Input files for PLUMED are openly available on PLUMED-NEST (https://www.plumed-nest.org/), the public repository of the PLUMED consortium, as plumID:22.021.