Further Insight and Future Perspectives.

By direct simulation of heterogeneous ice nucleation, we have revealed the detailed interplay between ice I polytypism and heterogeneous ice nucleation. Our results explicitly show, for flat smooth substrates, that both the temperature at which nucleation is induced and the polytype promoted by a substrate is determined by the water contact layer’s resemblance to a face of ice. The two metrics introduced here, similarity and exclusivity, capture the geometric bias underlying this. In addition, the structure within the contact layer is shown to be more informative than the substrate lattice match. A liquid structure mechanism for polytype/morph selectivity has also been suggested for homogeneous nucleation of other systems such as hard spheres (60, 61), soft spheres (62), carbon (63), and molybdenum (64). Thus, the insight gained here for heterogeneous nucleation could extend beyond ice.

The results of this study reinforce the point that future nucleation studies should take into account the possibility of $$ {\mathrm{I}}_{\text{h}}$$, $$ {\mathrm{I}}_{\text{c}}$$, and $$ {\mathrm{I}}_{\text{sd}}$$, instead of just the thermodynamically stable $$ {\mathrm{I}}_{\text{h}}$$. Evidence of this is provided by the prevalence of all of the ice I polytypes seen here, combined with both the fact that the stable $$ {\mathrm{I}}_{\text{h}}$$ and metastable $$ {\mathrm{I}}_{\text{c}}$$ are achieved in the same manner and the surprising ease in which $$ {\mathrm{I}}_{\text{c}}$$ is promoted. We therefore highlight the suggestion of prior work (65), which formulated heterogeneous classical nucleation theory taking into account the polytype selectivity of the substrate. We will now place our results in a broader context and discuss their implications and links to experiment.

Connections to Experiment.

In this study, we propose that polytypism is prevalent at even the most mild supercooling (in initial ice nuclei). Our evidence is that 1) we see the formation of not only $$ {\mathrm{I}}_{\text{sd}}$$ but also, pristine $$ {\mathrm{I}}_{\text{c}}$$ for temperatures up to $$ \sim $$272.5 K on many systems; 2) promotion of $$ {\mathrm{I}}_{\text{c}}$$ at this temperature is shown to be extremely reliable via many simulations on one of our most promising designs; and 3) analysis of the water contact layer shows the nature of the INA’s templating effect and not the temperature of the system is the determining factor for polytype selection.

However, this is in conflict with the view that ice polytypism is important only at strong supercooling. del Rosso et al. (14) and Komatsu et al. (15) found $$ {\mathrm{I}}_{\text{c}}$$ to $$ {\mathrm{I}}_{\text{h}}$$ transition temperatures of 215 and 250 K, respectively. The difference between these results leaves the relative stability of $$ {\mathrm{I}}_{\text{c}}$$ unclear, but it appears to be more stable at strong supercooling. With regard to $$ {\mathrm{I}}_{\text{sd}}$$, experimental studies have shown it is energetically unfavorable to have stacking disorder in ice, and above 240 K, it will eventually anneal to form pristine $$ {\mathrm{I}}_{\text{h}}$$ (10, 4546–47). Thus, $$ {\mathrm{I}}_{\text{sd}}$$ could be expected to be more prevalent at lower temperatures. Yet, as noted by Murray et al. (26), while this is true for, for example, crystals at $$ -80°$$C sampled from the tropical tropopause where 50% of crystals had the trigonal structure indicative of $$ {\mathrm{I}}_{\text{sd}}$$, the trend is less clear for warmer temperatures. For instance, studies by Yamashita (66) showed for crystals grown at around $$ -15°$$C, no trigonal crystals formed, whereas at $$ -7°$$C, 69% were trigonal. Other experimental studies have also seen evidence of $$ {\mathrm{I}}_{\text{sd}}$$ for very mild supercooling (67, 68), yielding the conclusion that initial presence of cubic layers in ice could be relevant over a very wide range of temperatures (26).

Our results support this conclusion and help resolve this debate by showing that in fact the promotion of any of the ice I polytypes is possible for even the most mild supercooling. Moreover, we see that unless a unique face of $$ {\mathrm{I}}_{\text{h}}$$ is strongly templated in the nucleation event, then cubic layers of ice will be present in the initial ice nuclei. It therefore appears that both hexagonal and cubic layers of ice are ubiquitous in heterogeneous ice nucleation. Even as cubic layers anneal at such temperatures, their presence in the initial nuclei can affect the properties of the resulting macroscopic crystals (5, 6). Ice polytypism could therefore be more important than previously thought for many fields, a prime example being atmospheric science. Thus, weight is added to the claims of $$ {\mathrm{I}}_{\text{c}}$$’s potential importance in phenomena such as the topology of snowflakes (5, 6), cryopreservation (28, 29), and the high water vapor supersaturation of cirrus clouds (8) as well as its general existence in the atmosphere (7).

Routes to Polymorphs/Types.

We aimed to understand if heterogeneous nucleation of pristine $$ {\mathrm{I}}_{\text{c}}$$ could be a valid route to its formation and found much evidence indicating it is. A methodology to design substrates with desired polytype selectivity leads to the discovery of many promising designs that promote $$ {\mathrm{I}}_{\text{c}}$$. No dependence in a substrate’s nucleation capacity upon the polytype being promoted is observed, and the nucleation of pristine $$ {\mathrm{I}}_{\text{c}}$$ is shown to have a high degree of adaptability and robustness to both substrate deformations and system size. This shows the small free energy difference between $$ {\mathrm{I}}_{\text{h}}$$, $$ {\mathrm{I}}_{\text{c}}$$, and that of $$ {\mathrm{I}}_{\text{sd}}$$ can be overcome via substrate design. Thus, we suggest that heterogeneous nucleation could be a valid route to the formation of pristine $$ {\mathrm{I}}_{\text{c}}$$.

Synthesizing in the laboratory one of the substrate designs proposed here, or discovering them in nature, yields the possibility of achieving pristine $$ {\mathrm{I}}_{\text{c}}$$ from supercooled water. This possibility is further enhanced by the observation of $$ {\mathrm{I}}_{\text{c}}$$ forming at its “most metastable temperature” (most mild supercooling), and the fact that many hydroxyl patterns were discovered (as this gives a great number of candidate materials). We note that experiments with the aim of achieving $$ {\mathrm{I}}_{\text{c}}$$ should be carried out at strong supercooling where $$ {\mathrm{I}}_{\text{c}}$$ is more stable. This will help avoid annealing to $$ {\mathrm{I}}_{\text{h}}$$ after the initial nucleation. If a suitable material is discovered, utilizing, for example, the techniques of water vapor deposition could be used to control the growth of crystals to larger sizes, analogous to the widespread use of chemical vapor deposition to make, for example, laboratory-grown diamonds. Experiments should consider the front propagation speed of different ice faces, as anisotropy in ice crystal growth has been seen (69, 70), and a faster front propagation speed could encompass defects, resulting in stacking disorder. The effects of confinement on water can also be utilized. First, in, for example, the case of 2D confinement, two material surfaces could impose a stronger orientation upon the resulting crystal. Second, confinement can affect the relative stabilities of polytypes; for example, this was reported to result in $$ {\mathrm{I}}_{\text{sd}}$$ with high cubicity being observed in alumina pores (71). With regard to candidate materials, the approximate structural-based guideline presented earlier can help their discovery. Interesting avenues would be groups of materials with hydroxylated surfaces, where examples of effective INAs are already known. These include the myriad of clays and minerals (3233343536–37), organic crystals [such as alcohols (38), amino acids, and steroids (39, 40)], and biological matter (414243–44). The latter is particularly interesting as there exists the possibility of utilizing the technique of sequence-controlled polymerization to make biomaterials with desired hydroxyl group structures (72, 73). Furthermore, synthesis of the 2D designs of Fig. 3A could be achieved by designing self-assembled monolayers. A design with a hexamer chair-like configuration is shown in Fig. 3A. This templating could be mimicked via adsorption of molecules with such structures upon flat substrates. In addition, materials with a close lattice match to ice that possess the cubic structure could be promising, such as the cubic forms of AgI and CuI (74).

Finally, the substrate design methodology presented here can easily be extended to the plethora of materials that show polymorphism and stacking disorder in materials science, ranging from pharmaceuticals to close-packed metals and alloys. As shown in Fig. 1, the key steps of the approach are to take unique faces from the desired polymorph/type, simulate nucleation of the liquid/gas form upon this, and generate new designs by pruning and/or the extraction of contact layers. This can be utilized to discover substrate designs for nucleating pristine polymorphs/types or even for desired levels of stacking faults to give higher-order polytypes with desired stacking sequences. This would enable the fine tuning of physiochemical properties (e.g., solubilities, mechanical, band gaps, dielectric properties) by giving access to the full range of possible structures (7576–77). This is highly desirable as stacking disorder can affect things such as the solubility of pharmaceuticals and even the band gap and dielectric constant of diamond (75). Ritonavir is a well-known example of the former, a drug that has found use in treating HIV and was part of the World Health Organization’s “SOLIDARITY” global trial to treat COVID-19 (7879–80). For carbon, the metastable hexagonal polymorph, Lonsdaleite, has been calculated to have superior mechanical properties to diamond but is yet to be synthesized without stacking faults (76, 77). Systems that can nucleate desired crystal structures can be used to grow to macroscopic scales via, for example, vapor deposition. Moreover, the connection between polytype selectivity and nucleation ability seen in this study means that the substrate design methodology could instead be used to discover potent nucleating agents, even in situations where controlling the polymorph/type obtained is not required.

SI Appendix.

We provide additional material in SI Appendix about 1) the polytype distribution with T_n for the substrate database, 2) a selection of promising substrate designs to nucleate pristine $$ {\mathrm{I}}_{\text{c}}$$, 3) the large-scale simulation of the nucleation of $$ {\mathrm{I}}_{\text{c}}$$, and 4) the 2D similarity calculation methodology.

Dataset S1.

A selection of 15 substrates is provided that are promising in both their ability to nucleate $$ {\mathrm{I}}_{\text{c}}$$ and their potential for discovery/synthesis. The coordinates are in a single “txt” file, which gives the OH-group (x, y, z) coordinates for each substrate (units are angstroms).

Molecular Dynamics Setup.

To model heterogeneous nucleation, water molecules are placed above a frozen slab of crystalline surface, periodic in (x, y). An example of this before and after nucleation can be seen in Fig. 1B. In total, 1,119 different substrate–water systems were used in this study. This database consists of the LJ systems from ref. 53; graphitic and graphite oxide surfaces modeled in a manner similar to those in refs. 81 and 82, respectively; and the OH-group patterns from ref. 57 plus over 200 systems generated for this work using the methodology presented in Fig. 1C. The simulation boxes are $$ \sim $$45 to 60 Å in (x, y); 4,000 to 6,000 water molecules are used giving layers of 35- to 60-Å thickness, thick enough to yield a region of water within the film in which the bulk water density is recovered. For the water–substrate interaction, various Stillinger–Weber or LJ (12/6) parameters are used, depending on the nature of the substrate.

All simulations are performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) code (83); and sampled the constant number of particles, constant volume, and constant temperature (NVT) canonical ensemble, using chains of 10 Nosé–Hoover thermostats with a relaxation times of 0.5 ps; and integrated the equations of motion with a time step of 10 fs. Simulations are initially equilibrated for at least 50 ps with a time step of 5 fs, prior to the production run. To observe nucleation, cooling ramps are used following the procedures in, for example, refs. 57, 82, and 84, with the temperature reduced from 275 K at 1 K $$ {\mathrm{n}\mathrm{s}}^{-1}$$. For each system, five cooling ramps are performed, and T_n is extracted by the drop in the systems potential energy. Nucleation events were checked by hand, and there is no indication that any were homogeneous. The polytype formed upon nucleation is identified using the local $$ {\text{q}}_3$$ order parameter of Li et al. (85), in combination with the CHILL+ algorithm of Nguyen and Molinero (86). These are computed using the plugin for metadynamics 2 (PLUMED2) (87, 88) and the open visualization tool (OVITO) (89) software, respectively. Ice structures are labeled as $$ {\mathrm{I}}_{\text{h}}$$ or $$ {\mathrm{I}}_{\text{c}}$$ only if they are completely lacking the other polytype.

This study employed the coarse-grained mW model. This model has inherent limitations; for example, no distinction is made between acceptors and donors of hydrogen bonds, atomic-scale resolution is lost, and water dissociation is not accounted for. It also does not have charge, meaning the effects of surface charge on water in contact with substrates are not included. However, the mW model is an appropriate choice for our study for several reasons. First, it has been widely and successfully used in studying water, especially the nucleation of ice (30, 51, 53, 56, 81, 82, 84, 85, 90919293949495969798–99). Second, it accurately describes important properties of water such as the density, structure, and melting temperature (30). Third, it captures the small metastability of $$ {\mathrm{I}}_{\text{c}}$$ compared with $$ {\mathrm{I}}_{\text{h}}$$ (30, 100) as well as the thermodynamics of stacking faults. Both the calculated free energy cost of a growth fault in $$ {\mathrm{I}}_{\text{h}}$$ and the cost to grow a pair of cubic layers have been shown to be in very close agreement with experiment, with comparative calculation and experimental values from $$ 15.3\pm 2.3$$ to $$ 16.5\pm 1.7{\text{J mol}}^{-1}$$ and from $$ 9.7\pm 1.9$$ to $$ 8.0{\text{J mol}}^{-1}$$, respectively (51, 91, 101). Last, but not least, the low computational cost of the model—compared with atomistic models—allows for nucleation simulations to be performed on an extensive database of substrates, as is done here. Extending this study of polytype selection to atomistic models would be interesting, as it is foreseeable that effects such as the hydrogen ordering induced by a substrate could be important and potentially even increase the polytype selectivity. The effect of the faster kinetics of coarse-grained models could also then be investigated (we note that this could be small, as the faster kinetics of mW should be present for all polytypes; thus, the relative trends observed here should hold). However, capturing the energetics between polytypes is a fundamental requirement for the model, and we are not aware of any atomistic model that does this and is affordable for nucleation studies.

Similarity Calculations.

In this study, similarities between 2D images of the water contact layers (and also the substrates) and faces of ice (from the three ice I polytypes) are calculated. These structures are represented as 2D histograms on a fine grid with spacing 0.1 Å, upon which a rotationally symmetric normalized 2D Gaussian is placed at the position of each atom. Here, the general concepts of the similarity calculations are explained. More details are given in SI Appendix.

The similarity between an extracted, f, and reference, g, pattern was calculated as