Homogeneous ice nucleation in an ab initio machine-learning model of water

Ice crystallization from supercooled liquid water is one of the most emblematic phase transformations to be found in nature. It is of key importance in the regulation of our planet’s climate (1) and in many applications, such as artificial cloud seeding, cryopreservation, and food processing. Molecular simulations have proven an invaluable tool to obtain insight into molecular-level details of this process and to make predictions at conditions not readily accessible to experiments. For instance, Lupi et al. (2) considered the effect of stacking disorder (i.e., the presence of alternate layers of hexagonal and cubic ice) on the nucleation rates, and Sanz et al. (3) used systems of more than 100,000 molecules in order to compute nucleation rates at low supercoolings.

However, simulations of ice nucleation carried out so far have employed relatively simple empirical models, such as the coarse-grained monoatomic model of water mW (4) or the four-site rigid TIP4P water models (5). A different route to study this phenomenon is using ab initio molecular dynamics (MD) (6). In this technique, the forces acting on the atomic nuclei are derived from electronic structure calculations. At variance with empirical models, the ab initio potential energy surface does not rely on empirical information, captures complex bonding behavior between atoms, and describes the formation and breaking of chemical bonds. The solution of the many-body electronic Schrödinger equation is, in general, not tractable, and a widely used approximation in this context is Kohn–Sham density-functional theory (7) (DFT). The application of ab initio MD, however, has been limited for several decades to the simulation of relatively small systems (∼1,000 atoms) and short times (∼100 ps) due to its high computational cost. This limitation has precluded the study of ice nucleation from first principles.

A solution to this conundrum has been the use of machine-learning algorithms that are able to learn the energies and forces derived from DFT data (8). The machine-learning interatomic models constructed in this fashion reproduce the ab initio potential energy surface with high fidelity, are several orders of magnitude faster than DFT, and also show linear scaling with the number of nuclei. Such models have recently been applied to the study of crystal nucleation in silicon (9) and gallium (10). Previous simulations using first-principles models, however, explored only relatively large supercoolings, for which systems of a few thousand atoms are able to contain the required crystalline cluster.

Here, we compute ice-nucleation rates using an ab initio machine-learning model of water. We employ the seeding technique (3) and systems of up to 300,000 atoms in order to obtain nucleation rates in a broad range of supercoolings. Our results allow us to compare predictions from a model derived from first principles with direct experimental measurements of nucleation rates. Although we only simulate explicitly clusters of hexagonal ice, we take into account the effect of stacking disorder using a model for the chemical potential of ice with stacking disorder.

During homogeneous ice nucleation, an ice cluster is formed within bulk liquid water. Typically, this phenomenon takes place below the melting temperature, and, thus, there is a driving force for the formation of ice. However, the formation of an ice cluster in the liquid creates a liquid–solid interface with an associated energetic penalty. The competition between a favorable bulk term and an unfavorable surface term leads to a free-energy barrier that the system must surmount in order to proceed with the transformation. The existence of a free-energy barrier makes nucleation a rare event and severely hinders the ability to study the phenomenon directly by using molecular simulations. Although there have been attempts to study ice nucleation by using straightforward molecular simulations (11), in general, the problem must be tackled by using more sophisticated techniques.

A possible route to study ice nucleation on the computer is rare-event techniques, such as path sampling (2, 12), forward-flux sampling (13, 14), or metadynamics (15, 16). These approaches can provide valuable insights into the nucleation mechanism, albeit at a high computational cost. A simpler alternative is the seeding technique (3), which is based on performing a series of relatively short simulations at different temperatures, starting from a configuration that contains an ice cluster embedded in liquid water. The aim of these simulations is to find the temperature $T^{*}$, for which the chosen cluster is critical—that is, to say, at $T^{*}$, the cluster has equal probabilities of growing and thawing. This information is then used in combination with the equations of classical nucleation theory (CNT) (17) to calculate the nucleation rate. This approach has several potential pitfalls that can affect the calculated rates, such as the appropriate choice of an order parameter to calculate the cluster size, and the applicability of CNT to the nucleation process under study. These limitations have been carefully considered in the literature (18), and the seeding technique has been shown to provide nucleation rates in good agreement with other methods (19, 20).

Another crucial ingredient in the simulation of ice nucleation is an accurate description of the interatomic interactions. Here, we derive the forces between nuclei from first-principles calculations. In particular, we use DFT adopting the Strongly Constrained and Appropriately Normed (SCAN) (21) exchange and correlation functional. SCAN is arguably one of the best semilocal functionals available, and many properties of ice and water have been studied using this functional—e.g., in refs. 22 and 23. Driving the dynamics directly using DFT forces would be unduly costly, and, instead, we use a machine-learning model trained on DFT data. The model is based on deep neural networks and was constructed by using the deep potential methodology developed by Zhang et al. (24). Below, we refer to this model as SCAN-ML (i.e., SCAN-trained, machine-learning-based model). The SCAN-ML model was carefully trained to reproduce data over a vast region of the phase diagram of water (25). SCAN-ML has been used to provide evidence of the existence of a liquid–liquid transition at deeply supercooled conditions (26) and to study the ice I${}_{\text{h}}$–ice XI transition (27). The thermodynamic properties of this model relevant to ice nucleation were thoroughly characterized in ref. 23. The model has a melting temperature of 312 K, around 40 K larger than the experimental value. The density change upon melting is 6% in the model, somewhat smaller than the 9% found in experiments. Another important property is the relative stability between ice I${}_{\text{h}}$ and ice I${}_{\text{c}}$, which are the two competing polymorphs during ice nucleation at ambient pressure. The SCAN-ML model correctly predicts that ice I${}_{\text{h}}$ is more stable than ice I${}_{\text{c}}$, in agreement with experiments. Ref. 23 also analyzed the ability of the SCAN-ML model to reproduce SCAN energies in configurations that contain atomic environments compatible with both liquid water and ice and found that the model is a faithful representation of SCAN with deviations of less than 1.3 meV per H₂O molecule. We provide in Table 1 a summary of the properties of the SCAN-ML model, and we compare them with experimental data and results using the empirical water models TIP4P/Ice and mW.

Before describing the results of our simulations, we briefly discuss the advantages of SCAN-ML over empirical models. SCAN-ML is an all-atom, fully flexible model at variance with empirical potentials such as mW, which is a coarse-grained model, and TIP4P/Ice, which is an all-atom rigid model. Since SCAN-ML reproduces the DFT potential energy surface, the flexibility of the OH bonds depends on the environments, while in flexible empirical models, such as TIP4P/2005f (30), the flexibility of the bonds is modeled by using simple functional forms and a few parameters that do not depend on the environment. Another property that depends on the environment is the dipole moment of the water molecule. For instance, the dipole moment is different in liquid water and ice (22), but can also exhibit more subtle changes with the environment (31, 32). SCAN-ML is polarizable and able to capture the effects connected to changes in dipole moment (27, 32). SCAN-ML is also fully reactive and can describe the proton-transfer process in water. This model captures many-body interactions beyond two- and three-body, while mW is limited to three-body interactions, and TIP4P/Ice is based only on two-body interactions. SCAN-ML and TIP4P/Ice can both describe an important feature of ice I${}_{\text{h}}$—namely, proton disorder, which is absent in the coarse-grained mW due to the lack of protons.

Our simulations based on SCAN-ML also have several limitations. While the electronic degrees of freedom are treated quantum-mechanically, the dynamics of the nuclei are based on the equations of motion of classical mechanics. Therefore, we ignore nuclear quantum effects (NQEs) that could be modeled by using path-integral MD. Another disadvantage is that SCAN-ML is around 1 to 2 orders of magnitude more computationally expensive than empirical models. Also, the properties of SCAN-ML differ somewhat from experimental properties, and this shows the limitations of the SCAN functional in the description of water and ice. Lastly, the model is short-ranged, with an interaction cutoff of 6 Å. It thus cannot capture the long-range electrostatic interactions (present, for instance, in TIP4P models) or van der Waals forces beyond this range. Long-range electrostatic interactions could be modeled by using the recently introduced deep potential long-range scheme (33).

We now turn to discuss the results of the seeding simulations. We studied ice I${}_{\text{h}}$ clusters of around 200; 700; and 4,500 molecules embedded in liquid water, and the corresponding total number of water molecules in the simulation boxes were around 4,000; 12,000; and 100,000, respectively. The choice of system size is discussed in detail in SI Appendix. The initial, equilibrated configurations of such ice clusters are shown in Fig. 1 A–C. We refer the reader to Materials and Methods for information about the equilibration procedure. The clusters are nearly spherical, an observation that will be important when CNT is used to calculate several physical properties (see below). Some faceting of the clusters can be observed, and the hexagonal shape of the clusters is compatible with the sixfold symmetry of the basal plane of ice I${}_{\text{h}}$.

MD simulations were performed at different temperatures, starting from the equilibrated configurations. The change in cluster size as a function of time is shown in SI Appendix, Fig. S1 for the three cluster sizes at different temperatures. From these simulations, we identified the temperatures $T^{*}$ at which these clusters have equal probabilities of growing and thawing. In Fig. 1D, we show $T^{*}$ for the three cluster sizes studied here with the SCAN-ML model. In order to determine the cluster size $N^{*}$, we must choose a local order parameter. The results depend somewhat on this choice (18), and, here, we employ a criterion similar to the one used by Espinosa et al. (19) in order to compare our results with the data reported therein (see Materials and Methods for details about our criterion to identify ice-like molecules). We also include in Fig. 1D the comparison with the results of Espinosa et al. (36) for two widely used empirical models—namely, mW and TIP4P/Ice. The results show that the critical cluster sizes are fairly independent of the model. At the highest supercooling studied here, around 50 K, the dynamics of liquid water are very slow, and thermal equilibration might not have been reached (see SI Appendix for a detailed discussion on the relaxation times of liquid water in the SCAN-ML model).

In order to illustrate the ability of SCAN-ML to capture subtle quantum-mechanical polarization effects, we calculated the average dipole moment of the water molecule as a function of the distance from the center of the ice I${}_{\text{h}}$ clusters. The quantum-mechanical molecular dipole moment was computed according to the modern theory of polarization (37, 38), adopting the formulation in terms of Wannier centers (39). The dependence of the Wannier centers on the coordinates of the atoms in the system was described by a deep neural network, as described in refs. 31 and 32 (see SI Appendix for further details of the calculation).The results are shown in Fig. 1 E–G. The average dipole moment changes from around 3.25 D in the ice I${}_{\text{h}}$ cluster to around 3.1 D in the liquid water surrounding the cluster. We also show in Fig. 1 E–G the reference values for bulk ice I${}_{\text{h}}$ and liquid water (dashed lines), and the agreement with the dipoles in the cluster configurations is very good. There is also good agreement between the bulk dipole moments calculated here and reference values obtained with SCAN DFT (22). The experimental dipole moment of liquid water at 298 K is 2.9 ± 0.6 D and is reproduced relatively well by SCAN (22). We note that the average dipole moment of the water molecule is a function of the temperature. Since each cluster has been equilibrated at a different temperature, the reference bulk values differ in Fig. 1 E–G. Furthermore, in the configurations with the ice cluster embedded in liquid water, the average dipole moment transitions smoothly from the bulk ice I${}_{\text{h}}$ value to the bulk liquid value, and the water molecules at the interface have, on average, intermediate values of the dipole moment.

We now turn to assess the performance of the SCAN-ML model to describe ice nucleation. We calculate nucleation rates by combining the information obtained from the seeding simulations with CNT. The predictions of CNT rest on various assumptions (17)—for instance, CNT assumes that clusters are spherical and show bulk ice I${}_{\text{h}}$ properties. Within CNT, the nucleation rate (nuclei per unit time per unit volume) is,where ρ_l is the density of the liquid; Z is the Zeldovich factor, which represents the probability of a critical cluster to cross the energy barrier; f is the attachment rate; $\beta =1/(k_{B}T)$; T is the temperature; and k_B is the Boltzmann constant. The nucleation free-energy barrier, $\mathrm{\Delta }{G}^{*}$, can be calculated by using the CNT formula,

where $\mathrm{\Delta }\mu$ is the difference in chemical potential between liquid water and ice I${}_{\text{h}}$, and $N^{*}$ is the number of water molecules in the critical cluster. Eq. 2 provides a convenient way to calculate rates from $N^{*}$ and $T^{*}$ obtained in the seeding simulations. ρ_l and $\mathrm{\Delta }\mu$ of the SCAN-ML model were calculated in ref. 23, and further details about the determination of $N^{*}$, Z, f, ρ_l, and $\mathrm{\Delta }\mu$ can be found in Materials and Methods and SI Appendix.

The nucleation rates thus calculated are shown in Fig. 2, together with results of experiments (40–46) and simulations using empirical models mW (36) and TIP4P/Ice (14, 16, 36). Most experiments are performed on micrometer-sized droplets and yield nucleation rates in the supercooling range 35 to 40 K (41, 42). However, experiments in the last 10 y have also used nano-sized droplets to reach much deeper supercoolings (40, 43). We have not included in our plot the experimental results of Laksmono et al. (49) since there are discrepancies between their rates and most measurements. Furthermore, it has been argued that more than one nucleus could have formed in those experiments (50). We have included in Fig. 2 a horizontal line representing the experimental homogeneous nucleation limit—i.e., the rate at which a micrometer-sized droplet freezes in 1 s (19). We note that there are significant differences between the nucleation rates calculated from simulations using different methods. For instance, in the case of TIP4P/Ice, estimates from forward-flux sampling (14), metadynamics (16), and seeding (36) span about 10 orders of magnitude, which is, however, the typical error bar of the seeding technique. Nucleation rates of the SCAN-ML model are in good agreement with experimental measurements within the uncertainty of our calculation. Furthermore, the rates of SCAN-ML are intermediate between those of mW and TIP4P/Ice. Therefore, the performance of SCAN-ML is similar to that of the best available semiempirical models. As we shall see later, the inclusion of stacking disorder makes rates faster and reduces to some extent the discrepancy with experiment.

Another quantity that can be easily obtained from the seeding simulations is the interfacial free energy averaged over all orientations $\overline{\gamma }$. For this purpose, we employ the CNT expression,where the symbols have the same meaning as in Eq. 1 and Eq. 2, and ${\rho }_{\text{ice}}$ is the density of ice I${}_{\text{h}}$. The results of this calculation are shown in Fig. 3A. Data for the mW and TIP4P/Ice models obtained from seeding simulations (36) are also shown. The dependence of $\overline{\gamma }$ on supercooling is a consequence of two different factors. The first is that the interfacial free energy of a flat interface depends on the temperature, and the second is that each cluster has a different size, and this will affect $\overline{\gamma }$, as shown by the Tolman equation (17).

In order to validate the results obtained using seeding simulations, we also calculated the interfacial free energy γ at coexistence for flat interfaces using advanced sampling simulations. For this purpose, we computed γ for the most relevant interfaces in ice I${}_{\text{h}}$—namely, the prismatic ($1\overline{1}00$), secondary prismatic ($11\overline{2}0$), and basal (0001) planes. The method to compute γ for flat interfaces at coexistence is based on the reversible interconversion of the liquid and the respective liquid–ice I${}_{\text{h}}$ interface. This is achieved by a suitably designed bias potential that increases the probability of observing the high free-energy interfacial configuration. A schematic of the interface sampled during the simulation of the secondary prismatic plane is shown in Fig. 3B. Further details of this approach and its validation can be found in Materials and Methods. We also computed the interfacial free energy averaged over all orientations ($\overline{\gamma }$) as the mean of the three studied interfaces (19). The results of the free-energy calculations are summarized in Table 2, and $\overline{\gamma }$ is shown in Fig. 3A. As seen in this figure, the agreement between $\overline{\gamma }$ obtained from advanced sampling calculations and seeding is very good. For reference, we show in Table 2 results for the models mW and TIP4P/Ice, as reported in ref. 52.

We have also included in Fig. 3A and in Table 2 experimental results for $\overline{\gamma }$ at the melting temperature (53). There is no direct experimental measurement of $\overline{\gamma }$ at other temperatures, and estimates based on CNT differ significantly (53). For this reason, we have not included them in our analysis. The spread of the experimental results at the melting temperature is relatively large (∼20 mJ/m²) and has a mean value ∼31.5 mJ/m² after removing outliers. It has also been argued (19, 53) that the experiments of Hardy (51) based on the shape of the grain-boundary groove provide the most reliable estimate, with a value of 29.1 ±0.8 mJ/m². $\overline{\gamma }$ for SCAN-ML is well within the region of uncertainty of the experimental measurements. However, the interfacial free energy of SCAN-ML is higher than the average experimental estimate. This behavior can be rationalized by taking into account that SCAN-ML has a melting temperature and enthalpy of fusion higher than both the corresponding numbers for real water (experiments) and TIP4P/Ice. Turnbull (54) observed that there is a strong correlation between the interfacial free energy and the enthalpy of fusion, and Laird (55) has made a similar observation for the correlation between the interfacial free energy and the melting temperature. It is thus expected that the interfacial free energy of SCAN-ML should be higher than in the experiment and in TIP4P/Ice [melting temperature ∼270 K (5)]. In SI Appendix, Fig. S4, we show that, indeed, the interfacial free energy correlates very well with the melting temperature in the TIP4P family and SCAN-ML. An option to account for the different melting temperatures of the models is to compare $\overline{\gamma }$ using units of $k_{B}T$ for the energy. A plot of $\overline{\gamma }$ in units of $k_{B}T/{\text{m}}^2$ vs. supercooling is shown in SI Appendix, Fig. S11. One could also estimate the value of $\overline{\gamma }$ in mJ/m² that SCAN-ML would have if its melting temperature were the experimental one. An appropriate rescaling of $\overline{\gamma }$ is $\overline{\gamma }\prime =\overline{\gamma }\hspace{0.17em}{T}_m^{\text{exp}}/T_m^{\text{SCAN}\text{-}\text{ML}}$, where $T_m^{\text{exp}}$ and $T_m^{\text{SCAN}\text{-}\text{ML}}$ are the melting temperatures in the experiment and SCAN-ML, respectively. In this way, we obtain a interfacial free energy for SCAN-ML of $\overline{\gamma }\prime =31.5$ mJ/m² at the melting temperature.

We now turn to analyze the thermodynamic properties of the models that affect nucleation rates by assuming the validity of CNT. The nucleation barrier $\mathrm{\Delta }{G}^{*}$ controls nucleation rates at low and moderate supercoolings since it is exponentiated in Eq. 1. The CNT expression for $\mathrm{\Delta }{G}^{*}$ is,

Therefore, the central physical quantities that govern nucleation rates at low and intermediate supercoolings are 1) the difference in chemical potential between liquid water and ice I${}_{\text{h}}$ ($\mathrm{\Delta }\mu$), 2) the interfacial free energy of ice I${}_{\text{h}}$ with liquid water ($\overline{\gamma }$), and 3) the density of ice (${\rho }_{\text{ice}}$). In the next paragraphs, we analyze these quantities for SCAN-ML, TIP4P/Ice, and mW.

In Fig. 4A, we show the difference between $\left|\mathrm{\Delta }\mu \right|$ in different models and in the experiment $\left|\mathrm{\Delta }{\mu }^{\text{exp}}\right|$ as a function of supercooling. $\left|\mathrm{\Delta }{\mu }^{\text{exp}}\right|$ cannot be measured directly, and its calculation from experimentally measured heat capacities of liquid water and ice I${}_{\text{h}}$ (56) is described in SI Appendix. $\left|\mathrm{\Delta }\mu \right|-\left|\mathrm{\Delta }{\mu }^{\text{exp}}\right|$ is reported in units of $k_{B}T$ in Fig. 4A in order to compare models with different melting temperatures. At 35 K of supercooling, $\left|\mathrm{\Delta }\mu \right|$ is underestimated by 9% in SCAN-ML. The performance of SCAN-ML in describing this property is somewhat better than that of TIP4P/Ice, which underestimates $\mathrm{\Delta }\mu$ by 17% at the same supercooling. The mW model is the most accurate among the models considered here, with $\left|\mathrm{\Delta }\mu \right|$ at 35 K within 1% of the experimental value. However, mW changes from a underestimation of $\left|\mathrm{\Delta }\mu \right|$ at low supercoolings to an overestimation at large supercoolings. This is a consequence of a much weaker deviation of $\mathrm{\Delta }\mu$ from a linear dependence with temperature than the other models (SI Appendix, Fig. S8).

Results for the interfacial free energy are presented in Fig. 3A. There is limited experimental information to ascertain the deviation of the interfacial free energy with respect to experiments. However, the values of $\overline{\gamma }$ for TIP4P/Ice and SCAN-ML (adjusted for the different melting temperature) are in relatively good agreement with most experimental results and most likely within a 5% error. Instead, $\overline{\gamma }$ in the mW model is around 35 mJ/m², which is higher than the most reliable experimental estimates of $\overline{\gamma }$ and is most likely overestimated by around 10%. It is also possible to characterize the temperature dependence of $\overline{\gamma }$ using the interfacial entropy, $S_{\gamma }=-\partial \overline{\gamma }/\partial T$, that can be estimated from the slope of $\overline{\gamma }$, with respect to temperature in Fig. 3A. We observe that mW has a lower slope $\partial \overline{\gamma }/\partial T$ than SCAN-ML and TIP4P/Ice, and that the latter two models have a similar slope. This indicates that the interfacial entropy of the coarse-grained mW model is higher than in the TIP4P/Ice and SCAN-ML all-atom models that include protons explicitly. Following ref. 28, the interfacial entropy can also be calculated by using,where $\mathrm{\Delta }{H}_f$ is the enthalpy of fusion, and $\mathrm{\Delta }{C}_p(T_m)$ is the difference in heat capacity between liquid water and ice I${}_{\text{h}}$ at the melting temperature T_m. Using the values of $\mathrm{\Delta }{H}_f$ and $\overline{\gamma }$, reported in Tables 1 and 2, and $\mathrm{\Delta }{C}_p(T_m)=49$ J/(mol$\hspace{0.17em}·\hspace{0.17em}$K) (23), we obtain $S_{\gamma }\approx -232 \text{μ}$J$\hspace{0.17em}·\hspace{0.17em}$m${}^{-2}\hspace{0.17em}·\hspace{0.17em}$K⁻¹ for SCAN-ML. This result is in good agreement with $S_{\gamma }$ calculated from experimental data ($-215 \hspace{0.17em}\text{μ}$J$\hspace{0.17em}·\hspace{0.17em}$m${}^{-2}\hspace{0.17em}·\hspace{0.17em}$K⁻¹) (28). A similar analysis for TIP4P/Ice gives $S_{\gamma }\approx -226 \hspace{0.17em}\text{μ}$J$\hspace{0.17em}·\hspace{0.17em}$m${}^{-2}\hspace{0.17em}·\hspace{0.17em}$K⁻¹, also in good agreement with the experiment. The mW model has a $S_{\gamma }$ of $-44 \hspace{0.17em}\text{μ}$J$\hspace{0.17em}·\hspace{0.17em}$m${}^{-2}\hspace{0.17em}·\hspace{0.17em}$K⁻¹ (28) that is around a fifth of the experimental value. These results are in agreement with the discussion above based on the slopes of the lines in Fig. 3A. From this discussion, we deduce that an all-atom description seems essential to capture $\overline{\gamma }$ and its temperature dependence.

The density of ice I${}_{\text{h}}$ in the different models considered here is shown in Table 1 and in SI Appendix, Fig. S8 (data from refs. 29, 36, and 58). The density of SCAN-ML ice I${}_{\text{h}}$ is around 3% higher than in the experiment, and, according to Eq. 4, this would partially compensate the somewhat low $\left|\mathrm{\Delta }\mu \right|$ in this model. For TIP4P/Ice, the density of ice I${}_{\text{h}}$ is around 1% lower than in the experiment, and we expect it to have a negligible effect compared to other errors. Finally, the mW model overestimates the density of ice I${}_{\text{h}}$ with respect to experiment by ∼7%, and this might compensate in part for a large $\overline{\gamma }$.

We then computed the nucleation free-energy barriers using Eq. 4, the values for $\mathrm{\Delta }\mu$ and ${\rho }_{\text{ice}}$ reported in SI Appendix, Fig. S8, and the linear fits to $\overline{\gamma }$ shown in Fig. 3A. The results are shown in Fig. 4B. We have included barriers from refs. 16, 29, and 57 obtained by using umbrella sampling and metadynamics free-energy calculations. For SCAN-ML, we expect that the barrier should be overestimated since $\left|\mathrm{\Delta }\mu \right|$ is underestimated. This is compatible with the nucleation rates being somewhat slower than the experiment (Fig. 2). In the TIP4P/Ice water model, $\left|\mathrm{\Delta }\mu \right|$ is underestimated more than in SCAN-ML, and, therefore, we would expect an overestimation of the nucleation barrier and nucleation rates slower than in the experiment. At variance with this prediction, the seeding nucleation rates (36) of TIP4P/Ice seem to agree relatively well with the experiments. The rates calculated by Niu et al. (16) and Haji-Akbari and Debenedetti (14) for TIP4P/Ice, however, are slower than the experimental measurements. In the case of the mW water model, $\mathrm{\Delta }\mu$ is in very good agreement with the experiment. For this reason, we surmise that the slow nucleation rates in this model can be traced back to an overestimation of $\overline{\gamma }$ not fully compensated by the overestimation of ${\rho }_{\text{ice}}$.

Another important aspect of ice nucleation is stacking disorder. There is significant experimental (59) and computational (2) evidence that nucleating ice clusters contain stacking faults—i.e., alternating layers of ice I${}_{\text{h}}$ and ice I${}_{\text{c}}$—and the solid polymorph that exhibits this feature is called ice I${}_{\text{sd}}$ (60). The prevalence of stacking faults in ice at equilibrium depends on two thermodynamic properties—namely, the difference in chemical potential between ice I${}_{\text{h}}$ and ice I${}_{\text{c}},\hspace{0.17em}\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$, and the interfacial free energy between these two polymorphs, ${\gamma }_{I_{\text{h}}\to I_{\text{c}}}$. The experimental evidence on the value of $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ is limited due to the difficulty in obtaining pure ice I${}_{\text{c}}$, although very recently, it has become possible to prepare samples with high structural purity (61). The available experimental data put $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ in the range from 0 to ∼200 J/mol (see ref. 62 for a review). An alternative point of view is provided by Lupi et al. (2), who argue that the experimental results reported in ref. 63 put an upper limit to $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ at 16.5 ±1.7 J/mol. On the computational side, the TIP4P/Ice and mW models have very small values of $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ of ∼0 (64) and ∼5 J/mol (29, 65), respectively. In ref. 23, we have found a $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ for the SCAN-ML model of 65 ±37 J/mol. As we shall see, the precise value of $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ has an influence on rates, and further experimental and computational efforts are needed to shed light on its value.

We described the effect of stacking disorder using a model for the chemical potential of ice I${}_{\text{sd}}$ that rests on the following assumptions: 1) The entropy of mixing of ice I${}_{\text{c}}$ and ice I${}_{\text{h}}$ layers is ideal; 2) the interfacial free energy is negligible: and 3) stacking is only relevant in one direction—namely, the direction perpendicular to the basal plane of ice I${}_{\text{h}}$. It can be shown that the first two assumptions give a lower bound for the chemical potential of ice I${}_{\text{sd}}$. Since the effect of stacking disorder is more relevant when the chemical potential of ice I${}_{\text{sd}}$ is lower, then our model gives an upper bound for the possible effects of stacking disorder. A more sophisticated two-dimensional model has been used by Lupi et al. (2), and it was found that the simplified one-dimensional model underestimates the entropic stabilization due to stacking disorder. We also note that a similar model has been used by Pronk and Frenkel (66). Further details can be found in Materials and Methods. In Fig. 5A, we show the difference in chemical potential between ice I${}_{\text{sd}}$ and ice I${}_{\text{h}},\hspace{0.17em}\mathrm{\Delta }{\mu }_{I_{\text{sd}}\to I_{\text{h}}}$, as a function of supercooling, as obtained from our model. The model takes as input the difference in chemical potential between ice I${}_{\text{c}}$ and ice I${}_{\text{h}},\hspace{0.17em}\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$. We used two different values for this quantity, one compatible with the free energy of the mW model, $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}=5$ J/mol, and another one compatible with the SCAN-ML model, $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}=$ 65 J/mol. In both cases, $\mathrm{\Delta }{\mu }_{I_{\text{sd}}\to I_{\text{h}}}$ becomes negligible as the supercooling goes to zero and the critical cluster size goes to infinity. This reflects the fact that ice I${}_{\text{h}}$ is the most stable phase in the thermodynamic limit. At larger supercoolings, the finite size effects start to be important, and ice I${}_{\text{sd}}$ becomes progressively more stable against ice I${}_{\text{h}}$. The magnitude of the stabilization, around 50 J/mol, is small compared with $\left|\mathrm{\Delta }\mu \right|$, even at relatively large supercoolings.

It is also interesting to evaluate how $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$ affects the cubicity of the nucleating clusters. We evaluated the cubicity as a function of the size of the cluster, and the results are reported in Fig. 5B. We found that for $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}=5$ J/mol, the cubicity drops below 1% for clusters of around 100,000 molecules, in excellent agreement with findings of Lupi et al. (2) for the mW water model. Instead, for $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}=65$ J/mol, the cubicity drops below 1% at around 2,000 molecules. Therefore, the extent to which stacking disorder is relevant in small clusters depends largely on $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$.

$\mathrm{\Delta }{\mu }_{I_{\text{sd}}\to I_{\text{h}}}$ can be used within CNT to estimate the nucleation rates of ice I${}_{\text{sd}}$. In Fig. 5C, we show the nucleation rates of ice I${}_{\text{sd}}$ calculated in this fashion. For clarity, we only show the results using a low value of $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}$—namely, 5 J/mol, since this gives the greatest effect for the rates and is easier to visualize. For this reason, it should be considered an upper bound for the effect, rather than the most reliable estimate. Despite the systematic choices we have made to obtain the maximum possible influence of stacking disorder on the rates, the effect of stacking disorder is around 2 or 3 orders of magnitude at deep supercoolings. This relatively small change, however, improves somewhat the agreement of SCAN-ML with the experiment. We also calculated nucleation rates using $\mathrm{\Delta }{\mu }_{I_{\text{h}}\to I_{\text{c}}}=65$ J/mol, and the results are shown in SI Appendix, Fig. S9.

This work shows that the latest advances in ab initio MD allow studies of complex phenomena such as ice nucleation from first principles. Our findings indicate that nucleation rates predicted based on SCAN DFT are in reasonably good agreement with experiment. The rates are similar to those estimated with the TIP4P/Ice model and somewhat faster than the rates of the mW model. The nucleation rate is a complex quantity that depends on many different properties of an atomistic model, such as the density of liquid water and ice, the water–ice interfacial free energy, and the difference in chemical potential between water and ice. We have performed a careful analysis of these properties, and we have also compared them to the results of empirical models. SCAN-ML gives a balanced description of these properties that results in good agreement between the calculated rates and the experimental measurements.

We have also highlighted limitations of the SCAN functional in describing some of the properties of liquid water and ice. More accurate functional approximations and/or higher-level quantum chemical data are expected to improve the description of the properties of water. It has also been shown that the MB-pol model (67) based on Coupled Cluster CCSD(T) calculations reproduces experimental properties with high accuracy and is thus an interesting model to study in the future. Furthermore, in this work, we have neglected NQEs that may have an important impact on some properties. For instance, the difference in chemical potential between liquid water and ice is influenced by heat capacities, and the latter are affected significantly by NQEs (68). Despite their possible relevance, modeling NQEs through path-integral MD is still computationally impractical for nucleation simulations that employ large system sizes, such as the ones considered here. Understanding the impact of NQEs on ice nucleation is an interesting direction for future work. Finally, ab initio machine-learning models of water can be extended to describe a substrate in order to simulate heterogeneous ice nucleation, a process of direct relevance to atmospheric science and climate modeling. This would allow one to include hitherto-neglected phenomena, such as the effect of pH and the spontaneous hydroxylation of surfaces. For these reasons, we foresee continued progress in the simulation of ice nucleation from first principles and the prediction of rates that are in progressively improved agreement with experiments, as a result of an accurate description of the thermodynamic and kinetic properties of water and ice.