Homogeneous ice nucleation rates for mW and TIP4P/ICE models through Lattice Mold calculations

When liquid water is supercooled below its melting temperature, freezing occurs through a two-step mechanism: the first step being the emergence of an ice embryo and the second its subsequent growth driving crystallization of the whole system. Generally, impurities in water assist the formation of a solid nucleus that, once it reaches a certain critical size, promotes irreversible crystal growth.^1–3 However, in the absence of impurities, the liquid phase can remain metastable for a certain amount of time (depending on the supercooling conditions and the sample size) before undergoing the liquid-to-crystal transition. This phenomenon is termed homogeneous nucleation, and it is an activated process in which the system must overcome a free energy barrier to form the critical cluster.^4–6 

Understanding ice formation is extremely relevant in different fields such as atmospheric science,^7–12 food industry,^13,14 physics,^5,15–17 and geology.^18,19 One magnitude which is important to understand and characterize is the nucleation rate, which accounts for the typical time that it will take the system at some given conditions to form the critical nucleus. More technically, the nucleation rate (J) is defined as the number of critical clusters formed per unit of volume and time. Experiments on homogeneous ice nucleation are typically performed with micrometer-sized water droplets at temperatures ranging between 230 and 242 K.^17,20–26 Below this range, water freezes instantaneously (within timescales of the order of milliseconds or less), while at higher temperatures, freezing requires inaccessibly long timescales.²⁷ However, some experiments with nanosized water droplets have been able to measure nucleation rates in the so-called “no-man’s land,” which represents an extremely low temperature range from 200 to 230 K.^28–31 

Computer simulations have proved to be useful in evaluating nucleation rates for many different systems, ranging from Hard Spheres,^32–37 Lennard-Jones,^34,38–45 NaCl,^46–49 and different ionic salts,^50,51 among many other systems.^52–55 Specifically, and due to its aforementioned importance, ice nucleation has been extensively studied by means of computer simulations.^{27,34,56–61} Previous efforts have made use of different computational methods, such as Forward Flux Sampling (FFS),⁶² Transition Path Sampling (TPS),⁶³ Metadynamics,^64,65 Umbrella Sampling (US),⁶⁶ and Seeding.²⁷ Some of them are exact (FFS and TPS) or almost exact (US and Metadynamics) and provide rigorous values of J. Although these techniques use an order parameter as reaction coordinate to follow or drive the liquid-to-solid transition,^{33,58,64,65,67} their values of J do not depend, in practice, on the choice of the order parameter. However, they are computationally very expensive.^49,68–70 In contrast, we have recently proposed another technique denoted as Seeding, which is computationally very efficient and relies on the Classical Nucleation Theory (CNT).^{34,35,44,52,71–75} This technique provides reliable results of J provided that the order parameter is chosen with care. However, admittedly, the results of Seeding depend on the choice of the order parameter,⁷⁶ which is in contrast with the rigorous techniques described above.

In this work, we employ the Lattice Mold (LM) technique⁷⁷ to compute homogeneous ice nucleation rates. By evaluating the required free energy of creating a precritical cluster, and the rate of such cluster to induce crystallization, J can be estimated in a computationally efficient manner that does not rely on any theory (i.e., CNT) or choice of the local order parameter (in contrast with the Seeding technique). Nevertheless, although in LM, no traditional order parameters are used, the reaction coordinate to drive crystallization consists of a mold of potential wells (where the crystal density and lattice mold positions assume a given structure) that promotes the formation of a precritical cluster. Admittedly, LM is more expensive than Seeding, but it is still computationally efficient. Previously, LM has been successfully validated for the calculation of nucleation rates of Hard Spheres and the Tosi–Fumi NaCl model,⁷⁷ and here, it is adopted to measure ice nucleation rates within a temperature range that includes both micrometer- and nanometer-sized water droplet experiments.^{17,20–22,24–26,28,29,31} We first apply the LM technique to the case of ice nucleation in the monoatomic water (mW) model,⁷⁸ a coarse-grained description of water that mimics the hydrogen-bonding structure of water by means of a nonbonded angular dependent term that encourages tetrahedral configurations. Once we compare and validate the rates for the mW potential, we extend this technique to the TIP4P/ICE⁷⁹ model, which describes water atomistically as a nonpolarizable rigid molecule. The TIP4P/ICE model provides outstanding predictions of the solid and supercooled liquid properties of water, including the phase diagram, equations of state, ice growth rate, and interfacial free energy for the ice Ih–water interface,^61,79,80 which is remarkable considering the anomalous behavior of supercooled water.^81–83 We apply the LM technique to shed light on the nucleation rate of the TIP4P/ICE model at 230 K and normal pressure, which is currently under debate.^60,68,84 Overall, with our calculations, we aim to find a consensus on the ice nucleation rates of these two popular water models by means of a novel and computationally efficient rare event method.

Lattice Mold is a computational approach that can provide reliable estimations of the nucleation rate as previously shown for hard spheres and molten NaCl.⁷⁷ Within the LM framework, nucleation is induced by introducing a mold of potential energy wells with the crystal lattice of a given solid phase within the metastable liquid, as illustrated in Fig. 1(a)—similar in methodology to the Mold Integration technique for the calculation of coexistence liquid–solid interfacial free energies.^35,85,86 The interactions between the particles of the system (water) and the “virtual particles” of the mold is described by a square-well (SW) potential (a continuous version of it to obtain forces analytically) that induces the ordering of water molecules into a solid structure (ice in this case). The technique consists of two distinct steps: First, the reversible work to form a precritical cluster is computed by gradually switching on the interaction between the mold and the liquid molecules. Second, independent simulations are conducted to measure the average required time for the induced precritical cluster to complete crystallization (i.e., overcome the remaining nucleation free energy barrier). Initially, the free energy difference between the liquid and the liquid with such precritical cluster is

where ϵ_m is the maximum well depth to evaluate this integral, which ensures that all the sites in the mold will be occupied, $⟨N(ϵSW)⟩$ is the average number of mold sites occupied by liquid molecules at a given ϵ_SW value, and N_w is the total number of wells. This integral is evaluated numerically by performing several Molecular Dynamics (MD) simulations at different well depths (ϵ_SW). We compute the average number of occupied sites by dividing the potential energy of the water–well interaction by the well depth (ϵ_SW). It is crucial to make sure that the integral is reversible so that the cluster induced by the mold cannot be sustained in the absence of potential wells and quickly dissolves once the interaction between the mold and the liquid molecules is switched off. Since in LM simulations the position and orientation of the mold are imposed, this constraint must also be reflected in the calculation of the free energy such that

where k_B is the Boltzmann constant, T is the temperature, ρ_f is the fluid number density, and V_w is the volume of a well. The second and third terms on the right-hand side of Eq. (2) account for the translational and orientational corrections to the free energy for imposing a fixed position and orientation to the induced cluster, respectively. The probability per unit volume of finding a crystal cluster of the size of the one induced by the mold (P) is given by

For the second step, we carry out independent simulations with the precritical cluster sustained by the mold and calculate the induction time of nucleation. This second step of the technique gives an average induction time $⟨t⟩$⁠, which in combination with the first step provides an estimation of the nucleation rate (J),

It is crucial that the cluster formed by the wells is precritical so that the integral of Eq. (1) can be evaluated reversibly (i.e., not inducing crystallization of the system along the integration pathway). This is initially checked by performing simulations with different mold sizes (i.e., different number of wells) at ϵ_m. Then, we choose the smallest mold that results in complete crystallization after an induction time that is computationally affordable. By applying this protocol, we gurantee: (i) reversibility in Eq. (1); and (ii) minimization of the imposed constraint across the transition state. Nevertheless, as long as the chosen mold presents induction time before promoting system crystallization, the choice of the mold size does not alter the obtained nucleation rate, as discussed in Ref. 77. Taken together, the two steps of this technique provide the free energy of forming a precritical cluster and the time required for such a cluster to overcome the nucleation free energy barrier.

In this work, we employ the mW⁷⁸ and TIP4P/ICE⁷⁹ water models to study homogeneous ice nucleation at normal pressure. Intermolecular potential details can be found in the respective model references. Simulations were performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Molecular Dynamics package⁸⁷ in the NVT ensemble (with N equal to the number of particles, V the volume, and T the temperature), keeping T constant with the Nosé–Hoover thermostat.^88–90 Additionally, we compute equations of state in order to obtain mold structures with a relaxed density by performing NpT simulations, fixing the pressure (p) with the Nosé–Hoover barostat.^88–90 We integrate the equations of motion using the velocity-Verlet integrator. The simulation time step chosen is 3 fs for the mW model and 2 fs for the TIP4P/ICE model, and the thermostat relaxation time is 0.1 ps for the mW potential and 2 ps for the TIP4P/ICE model. In NpT simulations, the barostat relaxation time is 0.5 and 2 ps for mW and TIP4P/ICE models, respectively.

For consistency with previous nucleation rate calculations of both models,^{15,34,56,57,59–61,68,91} we set the potential cutoff at 8.5 and 4.3 Å for TIP4P/ICE and mW, respectively. Moreover, for the TIP4P/ICE model, long-range Coulombic interactions are treated with the pppm/tip4p particle–particle particle-mesh (PPPM) solver in LAMMPS⁹² with a relative error in forces of 10⁻⁵. Long-range corrections to the Lennard-Jones part of the potential for energy and pressure are also applied. While in the mW model the potential cutoff (4.3 Å) is an intrinsic model parameter that cannot be modified, we check that for the TIP4P/ICE model, the neglected potential energy contribution of employing a short cutoff (8.5 Å instead of a long cutoff such as 15 Å) is lower than 0.05% when evaluating the potential energy of a bulk liquid phase applying PPPM and long-range corrections (for further insights on potential truncation effects on ice nucleation, please see Ref. 93). Finally, for the TIP4P/ICE model, we keep the O–H bond length and H–O–H angle values constant with the Shake algorithm implemented in LAMMPS.⁹⁴ 

To implement the LM technique, we simulate the potential wells’ interaction with water molecules using a square-well like (i.e., continuous) potential (SW)⁸⁶ of the following form:

where ϵ_SW is the depth of the potential energy well, r_w is the radius of the attractive well, α controls the steepness of the well, and r is the distance between water molecules and the center of the wells. We choose α = 0.017 Å as discussed in Ref. 86. In practice, for the TIP4P/ICE model, this is the effective interaction between oxygen atoms and the wells, while hydrogen atoms have no interaction with the wells and spontaneously arrange into the correct orientation of the lattice. It is interesting to point out that the wells’ positions are fixed (i.e., they act as an external field) and we only integrate the positions of water. LM simulations are performed in the NVT canonical ensemble (to keep fixed the crystal lattice positions of the mold along the calculation), and the density of the system is chosen as that of a bulk liquid at the temperature and pressure (1 bar in this work) of interest. We make the simulation box large enough so that the emergence of the precritical nucleus barely alters the desired pressure. We also note that although in this work we compute nucleation rates at different temperatures along a given isobar (1 bar), the LM technique can be indistinctly applied to evaluate pressure-dependent rates (i.e., across an isotherm) as long as the metastable liquid and the crystal mold display their characteristic equilibrium density of the state of interest (p, T).

We first focus on the crystal nucleation of mW at 220 K and 1 bar by means of the LM technique. As discussed in Ref. 77, the use of molds consisting of potential wells of a limited width constrains the solid structure, hence, leading to a moderate overestimation of the free energy barrier (ΔG). Such an effect can be minimized by calculating ΔG with wells that are in the limit of fitting two water molecules at the same time, to avoid constraining the crystal lattice and providing correct estimates of the nucleation rate. To find the optimal well width that does not constrain the induced nuclei, we perform simulations of the liquid phase with one well of ϵ_SW = 8k_BT (sufficient well depth to completely fill the well) and we measure the average occupation of such a well with different potential widths. Reaching the maximum well radius before the occurrence of double-particle occupation ensures negligible translational constraint of confining a water molecule within a given position in the mold crystal lattice. In Fig. 1(b), we plot the average occupation of a well inserted into an equilibrated liquid at 220 K and 1 bar as a function of r_w. We find that the optimal width is 1.32 Å, the limit at which the average occupation of a single well is roughly 100%.

Since it would be incorrect to evaluate the nucleation rate by Eq. (4) with a mold of potential wells in the limit of double-particle occupancy, we measure the nucleation rate at different well widths lower than the optimal one, and we later proceed to extrapolate to 1.32 Å. We choose to apply the LM technique at r_w values of 0.85, 0.99, and 1.12 Å and solve Eq. (1) by performing several simulations at different ϵ_SW. To achieve this, we run simulations of a liquid with a mold of potential wells forming a spherical cluster with the equilibrium positions of the stable crystal phase at such conditions (ice Ih). The integrand of Eq. (1) is the average number of occupied wells (⁠$〈$N(ϵ_SW)$〉$⁠) at different values of ϵ_SW. In Fig. 1(c), we show the time-evolution of the mold average occupancy for three different simulations exhibiting low, middle, and high mold occupancy, at values of ϵ_SW of 0.1, 1.25, and 8k_BT, using a potential well width of r_w = 1.12 Å. By evaluating the average number of wells as a function of ϵ_SW, we construct the curves shown in Fig. 1(d) and find the area below the curves, which provides ΔG* through Eq. (1). We have checked that changing the grid intervals of ϵ_SW from 0.1 to 0.25k_BT along the steepest part of the integrand [i.e., from ϵ_SW = 1 to 2k_BT in Fig. 1(d)] does not result in significant differences in ΔG* (e.g., ∼0.3k_BT). Nevertheless, it is always recommendable to employ a thin grid in the transition region between mostly empty molds (ϵ_SW < 1k_BT) and practically full ones (ϵ_SW > 2k_BT) to ensure that such region does not add uncertainty to the computed free energy barriers. Next, we employ Eq. (2) to obtain the free energy of forming a precritical cluster with freedom to translate and rotate. As it can be seen in Fig. 1(c), intermediate ϵ_SW values along the integrand [Fig. 1(d)] can result in strong fluctuations of the mold occupancy N(ϵ_SW) over time, clearly denoting reversibility along the integration pathway.

Finally, to obtain the nucleation rate we must also compute the rate of the formed precritical cluster to overcome the nucleation free energy barrier and induce crystallization across the whole system. We conduct simulations with the wells enabled at maximum depth, ϵ_m = 8k_BT, and find the average time $⟨t⟩$ it takes for the cluster to irreversibly grow. In Fig. 1(e), we show the potential energy (just water–water interactions) evolution of ten different trajectories (with different initial velocities) for r_w = 1.12 Å. The dramatic drop in potential energy marks the transformation of the liquid into a solid, therefore, indicating the time required for irreversible growth of a given trajectory.⁹⁶ In Table I, we provide the obtained values of ΔG and $⟨t⟩$ for the studied well widths (r_w). We acknowledge that cluster size and irreversible growth can be tracked by other means, taking for example a local order parameter such as the $q̄4−q̄6$ proposed by Lechner and Dellago.⁹⁷ However, it is not strictly necessary to monitor the system evolution through a local order parameter as we demonstrate in Fig. 1(e) (using a global indicator as the system potential energy). Please note that when the potential energy reaches values that indicate irreversible growth [i.e., −45 kJ mol⁻¹ in Fig. 1(e)], the pressure has only increased 30 bars with respect to the initial pressure (p = 1 bar), negligibly affecting our nucleation rate calculations,⁹⁸ hence making the NVT ensemble suitable for the LM technique as long as the total system size is reasonably large (see Table II). Once that average time has been obtained, we can compute the nucleation rate by means of Eq. (4). Finally, in Fig. 1(f), we show the nucleation rate for the three different well widths studied and the extrapolation to the optimal width (see also Table I).

We now apply the LM technique at different temperatures to obtain a broad picture of the mW nucleation scenario at normal pressure (1 bar). We achieve this by inserting spherical molds with an ice Ih lattice that can induce precritical clusters from a temperature range between 215 and 240 K. The inserted molds are spherical and display a perfect ice lattice with the equilibrium density of the bulk crystal structure at each temperature, which was obtained by performing independent NpT simulations of bulk ice. Additionally, we also prepare molds of ice Ic and a stacking mixture of ice Ic/Ih as shown in Fig. 2(a) and Table II. Since ice Ic and Ih share the basal plane [(111) in ice Ic and (0001) in ice Ih], the Ic/Ih mixture was generated by alternating AB layers from ice Ih and ABC from ice Ic stacked along the basal plane, generating a hybrid crystal structure. All mold sizes, number of water molecules in each system, solid and liquid densities, and nucleation rates for the different states are reported in Table II. In Fig. 2(a), we show snapshots of a typical simulation box for the LM technique, along with closer views of three different molds employed at T = 240 K (pure Ih, pure Ic, and Ih/Ic stacking mixture) surrounded by water molecules. In Fig. 2(b), we plot the nucleation rate obtained by using LM for the three different structures compared with previous results obtained through different rare event techniques and brute force simulations.^{15,34,56,58,59,91} We find good agreement with the results obtained by Li et al.⁵⁶ by means of Forward Flux Sampling and by both Russo et al.⁵⁸ and Cheng et al.⁵⁹ through Umbrella Sampling simulations. In contrast, some disagreement is found with the nucleation rates obtained by Haji-Akbari et al. using Forward Flux Sampling calculations at 230 and 235 K.⁵⁷ Nevertheless, those calculations were recently revisited in Ref. 91 by the same author providing nucleation rates in excellent agreement with Refs. 56, 58, and 59 as well as with our present calculations (Fig. 2). The extrapolated trend of our results toward higher supercooling is also consistent with the nucleation rate of Moore and Molinero¹⁵ obtained by means of brute force simulations at low temperatures (200 K). Seeding results obtained by Espinosa et al.³⁴ are 4 orders of magnitude below the ones obtained through LM; however, this discrepancy lies within the uncertainty of Seeding calculations, which is of about 5 orders of magnitude due to the reliance of Seeding simulations on local order parameters to identify the number of particles in the critical cluster.⁷⁶ The excellent agreement between the predicted nucleation rates obtained through LM and previous calculations^{15,34,56,58,59,91} indirectly evidences that the choice of the mold size barely alters the obtained nucleation rates by this technique as discussed in Ref. 77.

Importantly, when comparing the different ice structures Ic, Ih, and the Ic/Ih stacking mixture, we roughly find the same nucleation rate for all structures within the uncertainty at the distinct studied temperatures: between pure ice Ic and Ih nuclei at 235 K and among pure ice Ic and Ih as well as with an ice Ic/Ih stacking mixture mold at 240 K [see Table II and Fig. 2(b)]. Please note that we only generate a mold of a stacking mixture of Ic/Ih at 240 K so that the mold, and consequently the induced precritical cluster, is large enough to display ice stacking disorder across the basal plane of the two competing crystal structures while maintaining the least possible surface/volume ratio (i.e., spherical shape). A higher fraction of interfacial particles in smaller clusters (i.e., molds of ∼150 wells at 235 K) would lead to minimal stacking differences between pure ice Ic or Ih molds and stacking mixture molds (with negligible differences at even lower temperatures). Results reported by Cheng et al.⁵⁹ already suggested that the nucleation rate of ice Ic and Ih are very similar, while that of the Ih/Ic stacking mixture is moderately higher (i.e., 2 orders of magnitude). In our case, we find that at 240 K and 1 bar, the nucleation rate of a stacking mixture nucleus is within the uncertainty of those of pure ice Ic and Ih. For the rest of previous calculations,^56–58,91 the employed rare event methods were not biasing a priori the nucleating ice polymorph; hence, their calculations are ascribed to mixtures of both polymorphs, although ice Ih is usually the dominant phase. Overall, our results in Fig. 2(b), where we compare LM results with other previous calculations for the mW model, reaffirm the validity of this technique, which had previously been proven for Hard Spheres and NaCl Tosi–Fumi models.⁷⁷ 

Once the results of J for the mW potential are computed, we apply the same methodology to evaluate the nucleation rate of the TIP4P/ICE model,⁷⁹ an accurate potential for describing supercooled water and the different ice phases.^83,99–102 We carry out our calculations at 230 K and normal pressure, conditions at which the homogeneous nucleation rate has been previously reported, making use of Metadynamics,⁶⁰ FFS,^68,91 and Seeding⁸⁴ calculations. The reported values of J obtained by means of these different techniques significantly disagree in several orders of magnitude; with the results obtained by Niu et al.⁶⁰ being 9 orders of magnitude greater than those of Haji-Akbari and Debenedetti.⁶⁸ The revisited work of this calculation by Haji-Akbari⁹¹ did not calculate the nucleation rate directly, but it estimates a value of J 4 orders of magnitude higher than the original published value in Ref. 68, thus reducing the discrepancy to 5 orders of magnitude. Moreover, Seeding results obtained by Espinosa et al.⁸⁴ were extrapolated towards lower temperatures using CNT-like fits³⁴ and provided a nucleation rate 5 orders of magnitude greater than those reported by Niu et al.⁶⁰ Although such differences lie within the uncertainty of both calculations (Metadynamics has ∼2 orders of magnitude of uncertainty and Seeding ∼5 orders of magnitude when the mislabeling criterion is used to identify molecules as liquid-like or solid-like³⁴), there is a noticeable difference of 10 orders of magnitude between the Seeding results and those of Haji-Akbari⁹¹ (far beyond the uncertainty between both methodologies, with the reported uncertainty of FFS being ∼1 order of magnitude⁹¹).

The procedure used within the LM framework for the TIP4P/ICE is identical to that followed for the mW model. By inserting a mold that creates a precritical cluster in the liquid [214 wells, see snapshot in Fig. 3(a)], we are able to perform multiple simulations at different well depths (ϵ_SW) ranging from 0 to 10k_BT, and consequently solve Eq. (1). Once we have obtained the free energy of formation of such precritical cluster induced by the mold, we calculate the average time it takes to overcome the nucleation free energy barrier. Gathered this information, we apply Eq. (4) to compute J. In the same way as we do for the mW model, we compute the nucleation rate at r_w values of 0.95, 1.02, 1.09, 1.1, and 1.23 Å and extrapolate J to the optimal well width of 1.36 Å in this case.

A snapshot of the employed mold in TIP4P/ICE LM calculations is shown in Fig. 3(a). A summary of these calculations, providing the number of wells in the mold, nucleation rate, and density of the stable (Ih) and metastable (liquid) phases, can be found in Table II. Moreover, further details on the free energy needed to form a precritical cluster with the employed mold as well as the induction time to complete nucleation are reported in Table III. Additionally, in Fig. 3(b), we plot the number of water molecules in the crystal nucleus as a function of time for eight different trajectories at a well radius of r_w = 1.09 Å, through which we obtain $⟨t⟩$ [required for Eq. (4)]. The result of the nucleation rate against temperature is shown in Fig. 3(c) along with the results obtained by Niu et al. through Metadynamics,⁶⁰ by Espinosa et al. through Seeding,⁸⁴ and by Haji-Akbari and Debenedetti⁶⁸ (also including the revisited correction from Ref. 91) via FFS. We find excellent agreement between our results and those of Niu et al.⁶⁰ We predict a nucleation rate of log₁₀ (J/m⁻³ s⁻¹) = 14.1(3.0), which matches within the uncertainty of the Metadynamics logarithm of the nucleation rate of 14.8(2.7) reported by Niu et al.⁶⁰ Seeding results are 5 orders of magnitude apart, which is a reasonable agreement given the nature of the Seeding technique and the fact that 230 K is an extrapolated state from direct Seeding measurements at lower supercooling.^34,84 Additionally, we find a moderate disagreement with the revisited FFS calculation by Haji-Akbari,⁹¹ which lies 5 orders of magnitude lower than our LM calculations and the recently reported nucleation rates using Metadynamics simulations.⁶⁰ 

Remarkably, we note that much larger discrepancies between nucleation rates from different authors can be found for the TIP4P/ICE model [Fig. 3(c)] than for the mW potential [Fig. 2(b)]. We hypothesize that the reason behind such fact is the much better sampling (due to a faster relaxation time) that can be attained with the mW potential compared to the atomistic TIP4P/ICE model.⁶¹ Generally, nucleation rate discrepancies among different rare event techniques have emerged at low supercooling states where extremely large sampling is required.⁹¹ For instance, in FFS calculations, inconsistencies have emerged when a notorious number of milestones needs to be employed to connect the metastable state with the critical cluster. That usually happens at low metastabilities, and it is especially accentuated in computationally expensive models with long relaxation times such as the TIP4P/ICE model here compared to the mW potential.⁶⁸ That FFS is more prone to underestimate J at low supercooling is supported by the nucleation rates computed via FFS at 230 and 235 K from Ref. 57 [shown in Fig. 2(b)] where inconsistent values compared to the rest of techniques were reported only at moderate supercooling (i.e., 230 and 235 K at 1 bar), while at deep supercooling (i.e., from 215 to 225 K at 1 bar), good agreement with different independent calculations was found. Moreover, it has been recently noted in oppositely charged colloids that such lack of extreme sampling under conditions of low metastability may lead to underestimated values of J via FFS calculations.¹⁰³ Hence, the better agreement between independent techniques in mW calculations compared to TIP4P/ICE nucleation rates might be due to the less demanding computational effort required to sample the pathway connecting the metastable fluid and the critical cluster in the mW model.⁹¹ 

A comparison between computational and experimental estimates of the nucleation rate is critical to elucidate how good our understanding of homogeneous ice nucleation is. In Fig. 4, we show the nucleation rate with respect to supercooling conditions (ΔT = T_melting − T, where T_melting is 273 K for the mW model,¹⁰⁴ 270 K for the TIP4P/ICE model,⁶¹ and 273.15 K for real water) for the two studied models against experimental nucleation rates. Experiments performed with microdroplets^{17,20–22,24–26} and nano-droplets^28,29,31 are colored in dark and light green, respectively, while results for the mW water model^{15,34,56,58,59,91} are colored in blue and those obtained from the TIP4P/ICE model^60,68,84 in red. As depicted in Fig. 4, our understanding of ice nucleation at normal pressure spans over a wide range of temperatures. However, within such range, we can still find discrepancies between different experimental and computational approaches. At supercooling conditions ranging between 30 and 40 K, there is good agreement between all experimental results reported.^{17,20–22,24,25} Moreover, estimations of the nucleation rate in the same range for the TIP4P/ICE model also match experimental results, while the mW potential severely underestimates the nucleation rate within ∼10–15 orders of magnitude. This can be explained by the higher liquid–crystal interfacial free energy of the mW model at that supercooling conditions when compared to the TIP4P/ICE model.⁶¹ However, at greater supercooling conditions (above ΔT = 45 K), the mW model appears to recover the experimental J trend.

It is worth noting that, as discussed in Ref. 84, in nano-droplet experiments (depicted by light green symbols in Fig. 4), the pressure inside such drops is larger than the atmospheric one by virtue of the Laplace pressure, and over ΔT = 50 K, metastable liquid nano-droplets can reach pressures up to 500 bars. The Seeding curve in Fig. 4 accounts for this effect; however, direct comparison between calculations performed at normal pressure and those of nanoscale experiments can be affected by the higher pressure in the latter, which has a decelerating effect on ice nucleation.⁹⁸ Under high supercooling conditions (ΔT $>$ 40 K), TIP4P/ICE predictions reported by Niu et al.⁶⁰ moderately underestimate the experimental nucleation rates^17,28,29,31 by 4–6 orders of magnitude. Moreover, in the range 35 K < $Δ$T $<$ 50 K, the mW potential underestimates the experimental nucleation rates by ∼5–10 orders of magnitude. We hypothesize that the difference in melting enthalpy of both models (and the high value of the ice–liquid interfacial free energy in the case of the mW model) with respect to that of real water might be behind such moderate nucleation rate underestimation (Fig. 4). Whereas the experimental melting enthalpy (ΔH_m) at normal pressure is 1.44 kcal/mol, the same magnitude for both potentials is 1.29 and 1.26 kcal/mol for the TIP4P/ICE and mW models, respectively. Given that the ice–water chemical potential difference (Δμ) along an isobar can be approximated until moderate supercooling conditions through Δμ = ΔH_m(1 − T/T_m), where T indicates temperature and T_m the melting temperature, we hypothesize that the lower melting enthalpies of both potentials can translate into lower chemical potential differences and, thus, lower driving force to nucleate. Furthermore, it has been suggested that ice nucleation might be favored at the subsurface of water droplets,¹⁰⁶ which would enhance the experimental nucleation rate. Even though there is intense debate in that regard (and still it is not clear whether nucleation may be favored or not; see Refs. 57 and 107), the fact that in computer simulations ice nucleation rates are evaluated in bulk liquid conditions (i.e., absence of interfaces) might also explain the lower nucleation rate predictions of both TIP4P/ICE and mW models.

Importantly, the most noticeable disagreement between experiments also occurs within this temperature range where Laksmono et al.²⁶ using microdroplets to measure J disagree with nanodroplet nucleation rates reported by Manka et al.,²⁹ Hagen et al.,³¹ and Amaya and Wyslouzil²⁸ by 8 orders of magnitude. In contrast, TIP4P/ICE Seeding estimations⁸⁴ in this regime match the experimental rates reported by Manka et al.,²⁹ Hagen et al.,³¹ and Amaya and Wyslouzil²⁸ (although probably due to an overestimation in J by virtue of the employed local order parameter). The interference of ice growth with ice nucleation in the experiments of Laksmono et al.²⁶ has been hypothesized to be the cause of such discrepancy.⁶¹ Interestingly, we note that under very high supercooling conditions (ΔT $>$ 55 K), mW nucleation rates become higher than those of TIP4P/ICE (considering that the J values reported by Niu et al. are more accurate than those obtained from Seeding simulations). In fact, spontaneous crystallization of the mW model has been reported by Moore and Molinero¹⁵ near 200 K. The moderately higher J values of the mW model compared to those of the TIP4P/ICE model at deep supercooling conditions (Fig. 4) on top of the significantly higher crystal growth rate of the mW model (∼2 orders of magnitude higher than that of TIP4P/ICE⁶¹) enables the direct observation of water freezing in mW simulations at large supercooling conditions without the need of rare event techniques. Finally, in Fig. 4, we show the nucleation rate obtained by means of Seeding simulations using SCAN-ML, a recently developed machine learning trained model,¹⁰⁵ which greatly resembles TIP4P/ICE results obtained by Niu et al.⁶⁰ Remarkably, both TIP4P/ICE and SCAN-ML recover most of the experimental rates within 5–7 orders of magnitude for a wide range of supercooling conditions (Fig. 4).

In summary, here we have computed ice nucleation rates for two popular water models through the Lattice Mold technique, an efficient and simple method that does not rely on any theoretical assumption (i.e., CNT). Our results show good agreement with some previously reported data for both models obtained through FFS,⁵⁶ Seeding,^34,84 US,^58,59 and Metadynamics.⁶⁰ The computational performance achieved through our method is also of high relevance: More than 21 × 10⁶ central processing unit (CPU) hours were required to compute the nucleation rate of the TIP4P/ICE at 230 K via FFS calculations,⁶⁸ while through the LM technique, less than 1 × 10⁶ of CPU hours were required. Despite the fact that LM calculations are still far more expensive than Seeding simulations—where evaluating the TIP4P/ICE nucleation rate at 230 K costs ∼150 000 CPU hours⁶¹—the technique still has a very good cost balance, especially considering the advantages with respect to Seeding, such as non-reliance on the Classical Nucleation Theory or in regard to the usage of complex local order parameters to determine the critical cluster size.

Overall, multiple publications can provide relatively disparate estimations of the nucleation rate between experiments and simulations for a molecule as essential as water, hence, demonstrating the complexity of these calculations. Despite this, our understanding on the time required to freeze water has considerably improved in the last few decades (Fig. 4), and here we have presented the Lattice Mold technique, aiming to reach consensus on the value of the homogeneous ice nucleation rate for two distinct water models that describe molecular interactions at different levels of resolution. Establishing a solid and computationally efficient methodology to compute nucleation rates is fundamental, since in the near future, the development and validation of more accurate water models will require feasible and reliable measurements of this critical magnitude controlling atmospheric phenomena and cryopreservation processes.

This project received funding from the Oppenheimer Research Fellowship of the University of Cambridge. A.R.T. was funded by the Universidad Politécnica de Madrid (Grant No. ESTANCIAS-PIF 20-TYOSR8-13-4N0WPQ) and the Oppenheimer Fellowship. I.S.-B. acknowledges funding from the Derek Brewer scholarship of Emmanuel College and EPSRC Doctoral Training Programme studentship, Grant No. EP/T517847/1. M.M.C. and J.R. acknowledge funding from the Spanish Ministry of Economy and Competitivity (Grant No. PID2019-105898GA-C22) and the Madrid Government (Comunidad de Madrid-Spain) under the Multiannual Agreement with Universidad Politécnica de Madrid in the line Excellence Programme for University Professors, in the context of the V PRICIT (Regional Programme of Research and Technological Innovation). J.R.E. also acknowledges funding from the Roger Ekins Research Fellowship of Emmanuel College. This work was performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk) funded by EPSRC Tier-2 capital Grant No. EP/P020259/1. The authors gratefully acknowledge the Universidad Politécnica de Madrid (www.upm.es) for also providing computing resources on Magerit Supercomputer. We thank Pablo Piaggi for sending us the preprint of Ref. 105 prior to publication. E.S. and C.V. acknowledge funding received under Grant No. PID2019-105898GB-C21.

The authors have no conflicts to disclose.

Ignacio Sanchez-Burgos: Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Andres R. Tejedor: Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – review & editing (equal). Carlos Vega: Conceptualization (equal); Methodology (lead); Writing – review & editing (supporting). Maria M. Conde: Methodology (equal); Writing – review & editing (supporting). Eduardo Sanz: Conceptualization (equal); Methodology (lead); Writing – review & editing (equal). Jorge Ramirez: Conceptualization (equal); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – review & editing (supporting). Jorge R. Espinosa: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.