Routes to cubic ice through heterogeneous nucleation.

The freezing of water into ice is one of the most important processes in the physical sciences. However, it is still not understood at the molecular level. In particular, the crystallization of cubic ice ([Formula: see text])-rather than the traditional hexagonal polytype ([Formula: see text])-has become an increasingly debated topic. Although evidence for [Formula: see text] is thought to date back almost 400 y, it is only in the last year that pure [Formula: see text] has been made in the laboratory, and these processes involved high-pressure ice phases. Since this demonstrates that pure [Formula: see text] can form, the question naturally arises if [Formula: see text] can be made from liquid water. With this in mind, we have performed a high-throughput computational screening study involving molecular dynamics simulations of nucleation on over 1,100 model substrates. From these simulations, we find that 1) many different substrates can promote the formation of pristine [Formula: see text]; 2) [Formula: see text] can be selectively nucleated for even the mildest supercooling; 3) the water contact layer's resemblance to a face of ice is the key factor determining the polytype selectivity and nucleation temperature, independent of which polytype is promoted; and 4) substrate lattice match to ice is not indicative of the polytype obtained. Through this study, we have deepened understanding of the interplay of heterogeneous nucleation and ice I polytypism and suggest routes to [Formula: see text] More broadly, the substrate design methodology presented here combined with the insight gained can be used to understand and control polymorphism and stacking disorder in materials in general.

Ice I has two distinct polytypes: stable hexagonal ice () and metastable cubic ice (). Under conditions ordinary on Earth, the free energy difference between these polytypes has been experimentally determined to be just 30 to 50 J in favor of (1–3). Pristine has not directly and unambiguously been observed in nature, yet there is much indirect evidence of its existence. For example, observations of a rare halo of the sun, called Scheiner’s halo, have provided evidence of in the Earth’s atmosphere (4). Similarly, studies of the topology of snowflakes provide evidence of growth from an initial nucleus of (5, 6), and crystals sampled from the polar stratosphere appear to have characteristics (7). In addition, cirrus clouds have been observed to have a relatively high water vapor supersaturation, possibly due to the presence of (8).

These observations have led to a large body of experimental work related to . Metastable forms of ice with cubic diffraction patterns have been made in various ways, such as water vapor deposition on cold substrates (9–11), freezing confined water in mesopores (12), or heating low-density amorphous ice (8). However, it has only become clear in the last few years that what was widely referred to as is actually stacking disordered ice (), an ice I structure composed of interlaced cubic and hexagonal sequences (13). An example of an stacking sequence is shown in Fig. 1, as well as images of the and polytypes.

Fig. 1.

(A) The three polytypes of ice I: cubic ice (), hexagonal ice (), and an example of stacking disordered (). consists of interlaced layers of and , which stack along the (111) and basal plane of the respective polytypes. Different stacking sequences are possible; the image shown is just a single example. For simplicity, is referred to as its own polytype in this article. (B) Example of a simulation before and after nucleation. Water molecules (blue) and substrate (gray) are shown. (C) Workflow for creating novel systems that nucleate a desired polytype, illustrated on . Different Miller indices for faces unique to the polytype are taken from the unit cell to make slab substrates (atoms are colored by height). Water is placed on the substrates, and cooling ramps are performed. From this, substrates that nucleate are identified. Successful substrates can then be pruned, or the ice contact layer that forms on them can be used to make additional substrate designs (as discussed in the text).

The realization that pristine had not been directly observed sparked renewed interest in and attempts to make it without stacking faults (i.e., with 100% “cubicity”). Significantly, in 2020, del Rosso et al. (14) and Komatsu et al. (15) reported the formation of pure in the laboratory (16). del Rosso et al. (14) made it by heating a powder of O ice XVII to 160 K under vacuum, and Komatsu et al. (15) made it by decompression of hydrogen hydrate at 100 K. There can be little doubt that these studies represent important milestones in the field. However, the approaches taken to obtain pristine are rather complex, involving very well-defined experimental conditions that will not readily be encountered in nature. With this in mind, it is interesting to ask: “Are there other routes to pristine ?” In particular, given that in nature, or invariably forms through a heterogeneous nucleation process at the interface with foreign materials, the question arises if (completely lacking in stacking faults) can form through a heterogeneous nucleation process. Previous studies lend support to this possibility. For instance, was formed in simulations of supercooled water confined within graphene-based nanostructures (17) and between two slabs of “ice 0” (18), which shows that desired polytypes of ice I can be achieved through heterogeneous nucleation. has also been observed in simulation studies upon the application of external static electric fields (19). More broadly, in pharmaceuticals foreign substrates are used as an experimental means of controlling polymorphism and crystal structure (20–24). Another example is elemental silicon where conversely, the cubic polytype is stable over the hexagonal; hexagonal silicon was first achieved experimentally via epitaxial growth upon a foreign substrate in 2015 (25).

Understanding the formation of ice is important not just from a fundamental desire to understand one of the two polytypes of ice I but because the existence of different polytypes and stacking fractions affects important physiochemical properties of ice including crystal shape (26), light scattering (26), vapor pressure (8), and surface chemistry (27). In addition, unlike , is believed to be benign to biological tissue, yielding potential in cryopreservation (28, 29). Considering this, the aim of the current study is to explore routes to pristine through heterogeneous nucleation. This requires developing the understanding of the interplay of ice polytypism and heterogeneous ice nucleation and answering a more general question: “To what extent can templating achieve heterogeneous nucleation of desired polytypes?” To tackle this, we use molecular dynamics on model slab substrate plus water systems using the coarse-grained monoatomic water (mW) water model (30). Since effective ice nucleating agents (INAs) are many and varied in nature [e.g., AgI (31), feldspar (32–34), kaolinite (35–37), organic crystals (38), and biological materials (39–44)], we simulate a broad range of model substrates (over 1,100) to minimize overreliance upon a single case study and to boost the chances of extracting general insights. From these simulations, we find that 1) many different substrates can promote the formation of pristine ice I polytypes, including the elusive ; 2) pristine can be selectively nucleated for even the most mild supercooling; 3) the water contact layer’s resemblance to a face of ice is the key factor determining the polytype selectivity and nucleation temperature (Tn), independent of which polytype is promoted; and 4) substrate lattice match to ice is not indicative of the polytype obtained. In addition, the substrate design methodology presented here and the insights obtained can be used to control polymorphism and stacking disorder in materials beyond water–ice.

Designing Substrates to Control Ice I Polytypism

Substrate designs with polytype selectivity were discovered for , , and using the general workflow illustrated in Fig. 1. Starting with the bulk unit cell of the desired polytype, slabs of substrate are generated by taking different Miller index cuts of the faces unique to the polytype (templating shared faces will lead to stacking disorder—which can also be done if this is the desired outcome). The substrate is modeled as either a Lennard–Jones (LJ) system (over a range of interaction parameters) or as an arrangement of OH groups. Water is then placed upon the substrate, and cooling ramps are performed, as illustrated in Fig. 1. Successful systems, where the desired polytype forms, are then identified and used to generate new substrate designs through two routes: 1) substrate pruning via manually deleting atoms and 2) extraction of the ice contact layer with the substrate, which is subsequently used as a new substrate. The former was done with the aim of making the substrate simpler and/or flatter since if successful, candidates would be easier to engineer experimentally. The latter was done to utilise the deformations present in ice contact-layers to widen the variety of substrate designs considered. Through these approaches, substrates with different surface densities, structures, and symmetries from the original ice lattice can be identified. Cooling ramps are then performed on these new designs, and the process is repeated for systems that are successful. In total, 219 systems were generated this way, which combined with the 900 taken from prior work (), gives a total of 1,119 systems in the database.

Ice I Polytypism: Prevalent but Controllable

Fig. 2 summarizes the results of the cooling ramps performed on the database of over 1,100 slab substrate and water systems. Tn provides a quantitative measure of the nucleation ability of a substrate, with higher Tn for more potent INAs. A substrate’s polymorph/type selectivity is given by its ability to consistently nucleate a polymorph/type. The full range of nucleating abilities is present, from systems where nucleation is not observed at all to potent ice nucleation ability requiring only very mild supercooling. The nucleation of all three polytypes is also seen for many different systems, with the relative proportions of the different polytypes found shown in Fig. 2, Inset.

Fig. 2.

Histogram of nucleation temperature (Tn) for all systems in the database where nucleation was observed. For each system, five cooling ramps are performed, and Tn is the average temperature at which nucleation is observed. Systems where nucleation is not observed are omitted from the histogram. Inset is a donut plot showing the proportion of systems in the database for which a polytype (or no nucleation) was observed. A system’s polytype is labeled “varied” when no single polytype is obtained 100% of the time.

In , we explore the polytype distribution with Tn and find the nucleation of , , and for even very mild supercooling (temperatures up to 272.5 K). The presence of not only , which contains individual layers of the metastable , but also, pristine crystals at relatively high temperatures is a surprising result. Due to its metastability, the existence of cubic layers of ice has been believed to be only prominent at strong supercooling. In heterogeneously formed ice, it has only been observed up to 257 K in experiment (13), and samples of will eventually transform to pristine upon heating above 240 K (45–47). However, our results suggest ice polytypism could be prevalent in nucleation not only for highly supercooled conditions, but for all levels of supercooling. Before explaining this observation and the general insights obtained into polytype selectivity, we now focus on the nucleation of given its previously elusive nature.

Nucleation of : Highly Robust and Promoted by Numerous Diverse Substrates

The workflow illustrated in Fig. 1 produced many substrates that nucleated . In total, 70 substrate setups, 56 of which are unique structures (where the difference is not just a varied water–substrate interaction), were found to nucleate 100% of the time. Most effective substrates are composed of OH groups, and several active substrates lack any strong similarity to the lattice. Some examples of effective and structurally simple substrates are shown in Fig. 3. Structure files for a selection of the most promising substrates are provided in .

Fig. 3.

Hydroxyl group substrate designs to nucleate . (A) Top and side views of some promising designs discovered that nucleate , labeled with percentage chance nucleates and Tn. OH groups are colored by height, with distances for one unit cell given. (B) Heat map exploring the adaptability of a substrate’s nucleation characteristics to atomic displacements within the substrate. Specifically, Tn is plotted against mean absolute displacement of substrate atoms in the xy plane and z axis for the substrate indicated by the asterisk in A. Atoms are randomly displaced by vectors derived from a Gaussian distribution centered at 0.5 and with SD of 0.2. The axes show the mean displacement, and the maximum displacements present will be approximately twice this. Each marker is a pie chart, whose filled area indicates the percentage chance of obtaining —calculated from 6 cooling ramps over five different substrates, giving 30 cooling ramps per pie chart.

The example substrates shown in Fig. 3 have specific arrangements of OH groups and surface unit cells, as is indicated for one of the substrates. If precisely this specific arrangement of OH groups is required to nucleate , then the chances of finding such a structure in nature or creating it in a laboratory would be extremely slim, and even thermal vibrations might render it ineffective in exclusively nucleating . To explore how robust to change a substrate is, we took the two-dimensional (2D) hydroxyl group design marked with an asterisk in Fig. 3, made random displacements of the substrate atoms, and examined ice nucleation upon it. The results of this analysis are shown in Fig. 3, from which it can be seen that the nucleation ability of this substrate is indeed highly adaptable and robust to random displacement. Pristine is achieved for 100% of simulations with mean absolute displacements of up to 0.45 Å in the xy plane and for at least 0.75 Å in z (maximum displacements present will be approximately twice this). Therefore, there is a high degree of resilience in this system. is also obtained for a large range of Tn, from 272 to 210 K.

Considering the possible introduction of stacking faults leading to , it is natural to think of system size effects. Larger systems will have a greater chance for the introduction of stacking faults and will involve ice growth farther away from the substrate. It has been suggested in prior work that the resulting ice will always be despite any initial nucleation of a particular polytype (13). To address this, we simulated nucleation of a computationally very large system of 600,000-mW molecules (which is at the limit of affordability) upon a perturbed 2D substrate from Fig. 3 (i.e., a more realistic substrate than a pristine one), at constant temperature (above the system’s Tn to ensure the nucleation event is not spinodal decomposition). More details of this simulation are given in . The system fully crystallized to pristine . The present-day mismatch between available scales in computational and experimental work means it is not possible to fully understand the scale limit here, but this result is very encouraging and shows this 2D design has a remarkable ability to impose a specific crystal orientation over significant length scales.

To further guide the search for appropriate materials to nucleate pristine , we suggest the following approximate structural-based guidelines.

OH–OH separation: The symmetry of the substrate can differ from that of ice; however, generally, the substrate’s radial distribution function (rdf) contains a peak in the region of 2.5 to 2.7 and/or 4.3 to 4.6 . This corresponds to the first and second peaks of the rdf of bulk water, respectively.

Water–substrate interaction: Most substrates that are effective in promoting have a minimum water–substrate bond length of 2.6 to 3.0 . None had a minimum bond beyond 3.5 . The first coordination shell of the water–substrate rdf also has a well-defined peak, indicating a well-defined local structure.

Water contact layer: The structure of the water contact layer is not required to match that of the substrate.

In addition to these structure-based guidelines, we also show later that a high Tn is a good indicator that a strong crystal orientation is being imposed upon the water.

Structure of Water Contact Layer Determines Polytype Selectivity

To get further insight into the connection between substrate structure and polytype selectivity, we look again at the full database of systems. Analyzing the water contact-layer structure on each substrate shows a strong connection between its similarity to an ice face and the polytype obtained. The similarity measurement is defined in Eq. and returns a value between zero and one for the match between 2D patterns. The contact layers were computed using the same simulation setup as described in but run at a constant temperature of 285 K for 2 ns. The water contact-layer structure is extracted every 0.5 ps, giving independent snapshots that are used to extract the average structure. Fig. 4 shows the maximum similarity to and and the “exclusivity” of the structure to that polytype (equal to maximum similarity to polytype minus maximum similarity to a different polytype). Positive values of exclusivity correspond to templates that look most similar to the polytype, and negative values indicate templates that look more similar to a different polytype.

Fig. 4.

Resemblance of the water contact layer of systems to polytypes of ice. The similarity in the xy plane to each polytype is calculated (details are in ). (Left) Maximum result obtained for similarity to . (Right) Maximum result obtained for similarity to . Points plotted against exclusivity, which is equal to the maximum similarity to the polytype minus the maximum similarity to a different polytype. (A) Points labeled as a polytype if it was obtained for 100% of simulations of that system. (B) The same data as in A but with the points colored by Tn.

From Fig. 4, it is clear that both the water contact layer’s similarity to an ice face and its exclusivity to a polytype are related to polytype selectivity. A strong connection between this and nucleation ability can be seen in Fig. 4. Physically, strong similarity and exclusivity can be understood as an effective geometric bias of the system toward forming a particular ice face. Stronger similarity of the water contact layer to an ice face increases the chances of ice clusters forming, which in turn, increases the chance a cluster reaches the size of the critical nucleus. This determines the systems nucleation ability/temperature. Polytype selectivity can be seen as a level of detail greater than nucleation ability alone, where the geometric bias is strong enough to ensure that only one polytype forms.

Similar trends are seen for and (Fig. 4, Left and Right, respectively). As the water contact layer template becomes less exclusively similar to an ice face of the respective polytype, is obtained. This should not be confused for some global preference for the faces [basal and (111) of and , respectively]. is obtained for systems that still have a reasonably strong resemblance to unique or faces due to the introduction of stacking faults. It is known that ice forming homogeneously lacks any imposed crystal orientation, resulting in with 50:50 randomly arranged hexagonal/cubic stacking (48–50). It has also been shown that the homogeneous critical nucleus (at 230 K) with stacking faults has a lower free energy than that of pure (51). Thus, the promotion of ice without stacking faults requires a strong crystal orientation to be imposed. Bi et al. (17) achieved this by using a concave wedge to give a multidimensional structural match to the water that reduced the ability of stacking faults to enter. Here, we see that individual substrates can promote ice without stacking faults by imposing a strong orientation upon the water contact layer. Specifically, the water contact layer needs to 1) be similar to a face of the desired polytype (similarity measurement) and 2) not be similar to a face of a different polytype (exclusivity measurement). The requirement of exclusivity is in agreement with the idea that a substrate that does not promote the basal face can lead to ice without stacking disorder. (35, 52–54) As the quality of the template weakens (similarity and exclusivity drop), the water contact layer imposes a weaker orientation upon the crystal, resulting in the greater likelihood of stacking faults entering and thus, yielding the formation of . This yields the intermixing of polytypes promoted and poorer nucleation ability that can be seen in Fig. 4. Furthermore, both the metastable and the stable show a very similar behavior in Fig. 4. This shows that for heterogeneous ice nucleation, it may not be the very small free energy difference between these two polytypes, and that of , but instead the geometric biasing of the water contact-layer that determines both Tn and the ability of a polytype to form. Consequently, no significant difference in nucleation ability nor polytype selectivity is seen between substrates that promote and those that promote .

We note that when the water contact layer strongly resembles [basal and (111) face of and , respectively], nucleation upon the faces is also seen in the exact same manner. The systems have high Tn and result in ice with that crystal orientation. Although we do not focus on in this study, we note that it could be possible to use the techniques presented here to design potent INAs for or achieve “higher-order polytypes” with desired stacking sequences [e.g., the polytype, which consists of perfectly alternating layers of cubic and hexagonal ice].

Weak Correlation between the Substrate and Water Contact Layer Template

Predicting the nucleation characteristics of substrates a priori is highly desirable and still not currently possible. Given our finding that the water contact layer’s similarity to ice is key, one might think the requirement is simply a high substrate lattice match with ice. Lattice match of substrates has been discussed widely since the discovery of the nucleation ability of AgI (55). However, much recent work has shown lattice match alone is not enough to predict a substrates nucleation ability (56–59). Less work has investigated the connection to polytype selectivity, but a similar (if not stronger) requirement would be natural to infer. This apparent contradiction in the understanding of the water contact layer leads to an interesting question: What is the connection between the similarity of the substrate to ice and the similarity of the resulting water contact layer to ice?

The connection between the similarity of the substrate to ice and the similarity of the resulting water contact layer to ice is reported in Fig. 5. The maximum similarity of water to any polytype is extracted, and the polytype is noted, which we call SW. The similarity of the substrate to this polytype is then extracted to give Ssub. We see that SW and Ssub do not correlate well, with an R2 value of only 0.35. As discussed earlier, the structure of the water contact layer determines much of a substrate’s nucleation ability and polytype selectivity. The substrate lattice match is thought to imply the structure of the water contact layer; however, we show here that it does not. It is therefore necessary to separate these two ideas in the understanding of heterogeneous ice nucleation. Substrates with a poor lattice match to ice can give water contact layers very similar to ice, yielding potent ice nucleation ability. Conversely, substrates with high lattice match can give water contact layers very dissimilar to ice, yielding poor ice nucleation ability. Therefore, for the systems studied here, the water contact layer structure is a better descriptor of ice nucleating ability than the substrate structure. Observing Fig. 5, the descriptive improvement this gives can be understood by comparing the regions (dashed outlines in Fig. 5) of large SW but low Ssub, where Tn is large, and of large Ssub but low SW, where Tn is small.

Fig. 5.

Substrate’s and water contact layer’s similarity to ice, colored by Tn. Maximum similarity of water to any polytype is extracted, and the polytype is noted, giving SW. The similarity of the substrate to this polytype is then extracted to give Ssub. The correlation between SW and Ssub is weak, with an R2 value of 0.35. Comparing the dashed outlined regions of large SW but low Ssub, where Tn is large, and large Ssub but low SW, where Tn is small, illustrates the difference in these metrics’ ability to describe Tn.

Discussion

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 , , and , instead of just the thermodynamically stable . 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 and metastable are achieved in the same manner and the surprising ease in which 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 but also, pristine for temperatures up to 272.5 K on many systems; 2) promotion of 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 to transition temperatures of 215 and 250 K, respectively. The difference between these results leaves the relative stability of unclear, but it appears to be more stable at strong supercooling. With regard to , 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 (10, 45–47). Thus, 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 C sampled from the tropical tropopause where 50% of crystals had the trigonal structure indicative of , the trend is less clear for warmer temperatures. For instance, studies by Yamashita (66) showed for crystals grown at around C, no trigonal crystals formed, whereas at C, 69% were trigonal. Other experimental studies have also seen evidence of 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 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 ’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 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 . No dependence in a substrate’s nucleation capacity upon the polytype being promoted is observed, and the nucleation of pristine 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 , , and that of can be overcome via substrate design. Thus, we suggest that heterogeneous nucleation could be a valid route to the formation of pristine .

Synthesizing in the laboratory one of the substrate designs proposed here, or discovering them in nature, yields the possibility of achieving pristine from supercooled water. This possibility is further enhanced by the observation of 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 should be carried out at strong supercooling where is more stable. This will help avoid annealing to 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 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 (32–37), organic crystals [such as alcohols (38), amino acids, and steroids (39, 40)], and biological matter (41–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. 3 could be achieved by designing self-assembled monolayers. A design with a hexamer chair-like configuration is shown in Fig. 3. 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 (75–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 (78–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.

.

We provide additional material in about 1) the polytype distribution with Tn for the substrate database, 2) a selection of promising substrate designs to nucleate pristine , 3) the large-scale simulation of the nucleation of , 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 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).

Materials and Methods

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. 1. 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. 1. The simulation boxes are 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 . For each system, five cooling ramps are performed, and Tn 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 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 or 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, 90–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 compared with (30, 100) as well as the thermodynamics of stacking faults. Both the calculated free energy cost of a growth fault in 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 to and from to , 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 .

The similarity between an extracted, f, and reference, g, pattern was calculated asThe first term on the right-hand side determines whether the pattern g is present in f regardless of additional noise and whether g is a weaker signal than f (giving 0.5 if so), whereas the second term on the right-hand side determines whether f is present in g, regardless of additional noise and whether f is a weaker signal than g (giving 0.5 if so). By combining the two, one can determine if f is the same as g, as a value of one will be obtained only if f is exactly g, both in pattern and intensity [i.e., f(x, y) = g(x, y)]. Effectively, this is a calculation of the similarity between two images, a topic of much work with numerous methods (102). Here, this method is chosen for its simplicity, which maintains a strong physical connection and avoids obscuring interpretation.