Crystalline clusters in mW water: stability, growth, and grain boundaries

TL;DR

This study uses numerical simulations to analyze the stability, growth, and grain boundaries of crystalline clusters in mW water, revealing size-dependent stability and the influence of grain boundaries on polymorph selection.

Contribution

It introduces a detailed analysis of size-dependent stability and grain boundary effects in crystalline water clusters, highlighting mechanisms behind polymorph competition.

Findings

Small clusters favor Ice 0 phase stability.

Large clusters favor Ice I phase stability.

Grain boundaries hinder growth of metastable phases.

Abstract

With numerical simulations of the mW model of water, we investigate the energetic stability of crystalline clusters for both Ice I (cubic and hexagonal ice) and for the metastable Ice 0 phase as a function of the cluster size. Under a large variety of forming conditions, we find that the most stable cluster changes as a function of size: at small sizes the Ice 0 phase produces the most stable clusters, while at large sizes there is a crossover to Ice I clusters. We further investigate the growth of crystalline clusters with the seeding technique and study the growth patterns of different crystalline clusters. While energetically stable at small sizes, the growth of metastable phases (cubic and Ice 0) is hindered by the formation of coherent grain boundaries. A five-fold symmetric twin boundary for cubic ice, and a newly discovered coherent grain boundary in Ice 0, that promotes cross nucleation of cubic ice. Our work reveals that different local structures can compete with the stable phase in mW water, and that the low energy cost of particular grain boundaries might play an important role in polymorph selection.

Equations

$${\mathcal{H}}^{\prime }=\mathcal{H}+|\Delta \mu |n^{2/3}(n^{1/3}-3n_c^{1/3}/2)$$

$$t_{\text{fp}}(n)=\frac{1}{2kV}\left\{1+\text{erf}\left[c(n-n_c)\right]\right\}$$

$$\Delta E_{\text{c0}}=e_I-\frac{e_0+e_c}{2}$$

Full text

I Introduction

The pioneering work of Sir Charles Frank Frank (1952) revealed that the lowest energy local structure can differ from the arrangement found in the stable crystalline phases. In the Lennard-Jones fluid, for example, Frank found that the 13-atom icosahedral cluster has an energy which is $8.4\%$ lower than 13-atom clusters corresponding to nearest-neighbour environments in close-packed crystals (both the fcc and hcp crystalline phases). The idea of competing local structures has had a profound influence on our understanding of supercooled liquid states and of their non-equilibrium extension, i.e. glasses Tanaka et al. (2019); Tanaka (2012); Shintani and Tanaka (2006); Leocmach and Tanaka (2012); Royall and Williams (2015); Royall et al. (2018); Russo et al. (2018a).

The study of locally favoured structures has also played an important role in the understanding of the behaviour of water Russo and Tanaka (2016). Differently from simple liquids with isotropic interactions, water is characterized by strong directional interactions that favor a local tetrahedral arrangement of particles. It is possible to classify these tetrahedral environments in two different groups: amorphous and crystalline local structures.

The first group is composed of amorphous local structures. Amorphous structures are characterized by a lack of full translational and rotational symmetry, meaning that there is a high degree of degeneracy with respect to the relative position and/or orientation of the tetrahedral units. One of the key properties of water, that sets it apart from non-anomalous liquids, is that its amorphous structures can be further grouped into two states Tanaka (2000). The first state, which we call the $S$-state, is a low density/energy configuration in which a central water molecule is surrounded by four tetrahedra that are translationally ordered (i.e. they form a well-defined second shell of nearest neighbours), but orientationally disordered (meaning that the relative orientation of the tetrahedra is disordered) Russo and Tanaka (2014). The second state, which we call the $\rho$-state, is a high-density/energy configuration with a collapsed second shell of nearest neighbours. The definition of an appropriate order parameter has revealed bimodality in liquid structures up to ambient conditions, and has allowed the definition of a microscopic two-state model that accounts for all water thermodynamic anomalies Russo and Tanaka (2014). Locally favoured structures also impose strong constraints on the dynamics of water molecules: more specifically, a water molecule has a slower diffusive dynamics the more $S$-states surround it. A hierarchical two-state model was built with this idea in mind, and was found to be able to accounts for water’s dynamic anomalies Shi et al. (2018a), while also predicting a strong-to-strong crossover at supercooled conditions Shi et al. (2018b).

The idea of different states playing a pivotal role in determining water’s behaviour is further supported by the behaviour of water in supercooled and out-of-equilibrium conditions. In particular, several simulation works have shown that water can exist in two distinct liquid forms: the LDL (low density liquid) and HDL (high density liquid) phases, which become indistinguishable at a liquid-liquid critical point Poole et al. (1992); Mishima and Stanley (1998); Gallo and Sciortino (2012); Gallo et al. (2016). A non-equilibrium counter part of this transition is also found in water glasses, where low-density (LDA) and high density (HDA) amorphous glasses have been found (with the proposal of an additional very-high density VHDA amorphous glass Handle et al. (2017))Mishima and Stanley (1998); Martelli (2019).

The second group of local structures is composed of crystalline local structures. While the anomalous behaviour of the liquid phase appears to be strongly linked with the competition of two different amorphous local structures, the process of ice nucleation is associated with a different type of local structure, i.e. crystalline local structures. These structures are defined as clusters of molecules whose nearest-neighbour environment is the same as those found in bulk crystals. They differ from the amorphous local structures for their high degree of local translational and rotational order. They can thus extend to arbitrarily large sizes, and their growth becomes irreversible above the critical nucleus size. The relationship of these crystalline local structures to amorphous local structures formed in a supercooled state is one of the important factors controlling the ease of crystallization Tanaka (2012); Russo et al. (2014); Russo and Tanaka (2016); Russo et al. (2018a). These clusters originates from fluctuations that are very localised in space and time. For this reason, computer simulations are giving a fundamental contribution to our understanding of this important process Brukhno et al. (2008); Moore and Molinero (2011a); Gallo et al. (2016); Russo and Tanaka (2016); Russo et al. (2018a).

Homogeneous ice nucleation occurs by the growth of the stable Ice I phase from a few nucleation sites within the supercooled liquid state. Ice I has two different polytypes: hexagonal ice (space group $P6_{3}/mmc$), and cubic ice (space group $Fd\bar{3}m$). Hexagonal ice is the thermodynamic stable form, but nucleation at supercooled conditions often proceeds through metastable phases, and deviations from six-fold rotational symmetry have been found in diffraction patterns for a vast array of crystallization scenarios (see Ref. Malkin et al. (2015) for a review). These include re-crystallization from either amorphous ices Shilling et al. (2006); Dowell and Rinfret (1960); Kohl et al. (2000); Loerting et al. (2006); Geiger et al. (2014), or high-pressure ices Arnold et al. (1968); Kuhs et al. (1987); Hansen et al. (2008a, b). This metastable phase was long identified with the cubic phase, but it has recently being pointed out that due to the high amount of stacking disorder Murray and Bertram (2007); Carignano (2007); Pirzadeh and Kusalik (2010); Rozmanov and Kusalik (2012); Seo et al. (2012); Malkin et al. (2012); Choi et al. (2014); Quigley (2014) this phase should be more appropriately called stacking disordered ice, ice $I_{\text{sd}}$, which retains a three-fold rotational symmetry and can be identified with the space group $P3m1$ Kuhs et al. (2012); Hansen et al. (2008a).

The formation of ice $I_{\text{sd}}$ can be understood in terms of cross-nucleation Desgranges and Delhommelle (2006); Nguyen and Molinero (2014), which occurs when two polymorphs share a common crystallographic plane, i.e. the (111) plane for cubic ice, and the basal plane (0001) of hexagonal ice. In this case the two crystals match perfectly at the interface plane (both in symmetry and lattice constant), without discontinuities across the two phases. This causes a great reduction of the interface free energy. In the language of interphase interfaces, the two Ice I polytypes form a coherent interface. Very recently, computer simulations of the mW model have shown that ice crystals with stacking disorder are about $14$ KJmol*-1* more stable than hexagonal crystals for sizes up to 100000 molecules at $T=230$ K Lupi et al. (2017), nicely showing that the driving force for crystallization is size dependent.

Cubic ice and ice $I_{\text{sd}}$ are the only metastable phases found in laboratory experiments at ambient pressure. But there are different candidate crystalline structures, which retain the local tetrahedral arrangement of water molecules, and have additional properties not found in the known forms of Ice Engel et al. (2018). We focus here on the so-called Ice [math] form, which is currently the only known phase (of the mW potential) that is metastable in the regime of ambient pressure and low temperatures that is relevant for homogeneous ice nucleation Russo et al. (2014). The unit cell of Ice [math] contains 12 Oxygen atoms, and is tetragonal where the ratio between the long and the short sides of the cell is $c/a=1.81$.

Several simulation studies have addressed the problem of homogeneous nucleation of water, both in molecular models, such as TIP4P Matsumoto et al. (2002); Radhakrishnan and Trout (2003); Malkin et al. (2012) and the coarse grained mW model Molinero and Moore (2008); Moore and Molinero (2010, 2011a, 2011b), which sacrifices orientational degrees of freedom in exchange of a considerable faster crystallization dynamics. To enhance the sampling of crystalline cluster, several techniques like Umbrella Sampling, Metadynamics, and Forward Flux Sampling have been employed Quigley and Rodger (2008); Li et al. (2011); Reinhardt et al. (2012); Reinhardt and Doye (2013). An understanding of ice nucleation requires accurate studies of interfacial properties Benet et al. (2014a, b); Cheng et al. (2015); Espinosa et al. (2016a, b); Ambler et al. (2017) which are commonly used within Classical Nucleation Theory parametrizations Koop and Murray (2016); Cheng and Ceriotti (2017); Cheng et al. (2017). The behaviour of water at interfaces is central in understanding heterogeneous nucleation Lupi et al. (2014); Reinhardt and Doye (2014); Cox et al. (2015); Sosso et al. (2016); Lupi et al. (2016), nucleation with free surfaces Johnston and Molinero (2012); Li et al. (2013); Hudait and Molinero (2014); Haji-Akbari et al. (2014); Malek et al. (2015); Haji-Akbari and Debenedetti (2017), in nano-confinement Moore et al. (2010); Gonzalez Solveyra et al. (2011), and crystallization in solutions Bullock and Molinero (2013); Soria et al. (2017); Conde et al. (2017).

In the present work we focus on crystalline clusters for the mW model of water Molinero and Moore (2008). We consider both polytypes of Ice $I$, i.e. the cubic and hexagonal ice phases, and the metastable Ice [math] phase. We first examine their energetic stability under a wide variety of conditions (temperature and density). We will show that for small clusters, the Ice $I$ crystalline clusters are not the most energetically stable. We instead show that small clusters with the nearest neighbour environment of Ice 0 have a lower total energy due to the abundance of $5$-membered rings. While Ice I clusters (both hexagonal and cubic forms) grow isotropoically in all directions, Ice [math] grows along layers composed of $5$-membered rings. We then perform seeding simulations at supercooled conditions, and examine the growth behaviour of the different clusters. In particular we highlight the importance of grain boundaries in the kinetic pathway of growth of the different crystalline clusters. We focus here on three coherent grain boundaries, whose low free energy penalties make them particularly abundant. We will estimate these free energy penalties, and try to rationalize their impact on ice nucleation.

II Methods

Here we perform Monte Carlo simulations of the mW model of water Molinero and Moore (2008). We employ isobaric simulations at ambient pressure, and focus on highly supercooled conditions, $T=218$ K, which is low enough to have accessible nucleation barriers, but high enough that nucleation is a rare process (not spontaneously occurring in unbiased simulations) and that dynamics can still be observed over $1000$ s of Brownian times. Ice crystals are then grown both via biased simulations, with the CNT-US method, and unbiased simulations, with the seeding technique.

The CNT-US method is a variant of Umbrella Sampling introduced in Ref. Russo et al. (2014), where the expression for the free-energy barrier of Classical Nucleation Theory is used as a biasing potential in simulations. The new Hamiltonian is then written as

[TABLE]

where $\mathcal{H}$ is the given by the mW potential, $\Delta\mu$ is the chemical potential difference between Ice I and supercooled water, $n$ is the size of the largest nucleus in the system, and $n_{c}$ is the critical nucleus size. The advantage of this method over standard Umbrella Sampling is that the whole free energy landscape can be explored in equilibrium in one independent simulation, avoiding kinetic traps that are often incurred when restricting simulations to small sampling windows.

To distinguish the nucleation barriers between different phases, e.g. Ice $I_{h}$, $I_{c}$, and nuclei with stacking disorder we perform standard Umbrella Sampling simulations in which we initialize the simulations with a nuclei of the desired phase. These simulations only preserve the biasing potential of the CNT-US scheme, but, as in traditional US schemes, the range of the order parameter (the size of the largest nucleus) is divided in many small windows, each one allowing fluctuations of the order parameter between $N_{i}$ and $N_{i}+\Delta N$. Since the initialization of each window occurs with a seed of the correct type, fluctuations to a seed with a different composition occur very rarely. We thus monitor the trajectory in each window and perform a check of the nucleus composition every 1000 Monte Carlo (MC) cycles to ensure that the nucleus has preserved the correct composition. If not, we perform a rejection move of the trajectory of the last 1000 MC time steps. We typically set $\Delta N=5$, but we have checked that the results are insensitive to this choice for $\Delta N=1$ (where the technique is formally equivalent to the Successive Umbrella Sampling scheme Rovigatti et al. (2018)) and $\Delta N=10$.

We then monitor the different simulations to make sure that the original crystal structure is preserved during the sampling, and ensure that no additional stacking faults have originated.

For unbiased simulations, we use the seeding technique Knott et al. (2012); Zimmermann et al. (2015); Espinosa et al. (2016c, d), which consists in the insertion of perfect crystalline nuclei in the melt in order to determine at which conditions the nucleus is critical, and then use Classical Nucleation Theory to estimate parameters like the surface tension and the crystallization rate. While being successful in the estimation of the crystallization rate of water at moderate supercooling, with critical nuclei composed of several thousand molecules and where the approximations of Classical Nucleation Theory are most likely to hold, potential problems with the extension of this technique to deeply supercooled conditions have been recently highlighted Lifanov et al. (2016). On the contrary, for systems seeded with stable Ice I crystals at deeply supercooled conditions, we will find that the critical nucleus sizes are in excellent agreement with the ones we obtain independently from Umbrella Sampling simulations.

To distinguish the different crystalline phases, and to bias the Umbrella Sampling simulations we use the standard bond-orientational order parameters. In particular we employ the $Q_{12}$ bond orientational order parameter, where a molecule is identified as having a crystal-like environment if it has more than 12 connected neighbours among its first 16 nearest neighbour. A connected neighbour of particle $i$ is defined as a neighbour $j$ for which the normalized product satisfies $Q_{12}(i)\cdot Q_{12}(j)>0.75$. See Ref. Tanaka et al. (2019) for a recent review of this methodology.

Where not differently specified, the system size of seeding simulations is $N=10000$ water molecules, while Umbrella Sampling simulations are run with $N=4000$ molecules. Since we only focus on small nuclei (up to $N=100$ molecules), no strong system size dependence of our results is expected.

III Results

III.1 Energy of crystalline clusters

In order to understand crystal clusters competition in supercooled water, we compare the energy of clusters of perfect crystals in three different ways to account for different scenarios of cluster formation: (i) at the equilibrium density $\rho$, corresponding to a specific temperature $T$, which depends on the structure under consideration (see Fig.1) ; (ii) at fixed density corresponding to that of the liquid phase at equilibrium; (iii) at fixed bond distance, in analogy with the work of Frank Frank (1952). Afterwards, (iv) we extend the approach (i) to account for clusters equilibration in the liquid phase, and (v) we consider also the energetic contribution coming from the hydration layer.

In Fig.1 we show the $\rho$ versus $T$ relation at constant pressure $P=0$ atm for the liquid phase (L), cubic ice (Ic), hexagonal ice (Ih) and ice 0 (Ice 0) as obtained from Monte Carlo simulations of the mW model in the NPT ensemble for $N=5376$ particles with periodic boundary conditions. We can notice the well known temperature of maximum density (TMD) anomaly of liquid water found at around T=250 K in the mW model at normal conditions Molinero and Moore (2008); Russo et al. (2018b). We also note crystal densities being lower than the liquid density above 207 K (see Fig.1).

Clusters of mW particles are built by adding particles to a crystalline lattice (or from a configuration of the liquid at equilibrium in the case of disordered clusters) considering the following energetic approach: starting from any position in the crystalline lattice, we add one particle at each step choosing among the lattice positions that result in the minimum energy cluster. This involves computing the bond energy of every particle on the surface of the cluster with its potential neighbours, and choosing to add the particle that generates the cluster with minimum energy. In case of degenerate lattice positions, i.e. different clusters with the same size and same energy, the choice of which lattice position to fill is done randomly. We repeat this procedure until we get a cluster of the desired size. Since clusters will have different morphologies (i.e. different numbers of inter-particle bonds), we compute the average energy of clusters over 50 independent growth processes.

In Fig.2(a),(b) we show the energy per particle of clusters, $E/n$, as a function of the cluster size $n$ for the crystals and liquid phases at the conditions described above (ii), (iii) and (i), respectively.

From Fig.2(a),(b) we can notice that clusters of Ic and Ih have always the same energy for any size $n$, as expected. We also notice that Ice 0 has a lower energy respect to Ice I for small clusters size up to the critical value $n^{*}\sim 40$ for the three conditions (i), (ii) and (iii) (we computed the value of $n=n^{*}$ at which curves associated to Ice 0 and Ic/Ih cross for other values of $\rho$ and $T$ following methods (ii) and (iii) and found always $n^{*}\sim 40$).

To extend the previous results, we account for clusters equilibration in the liquid phase (method (iv)). In order to do so we perform Monte Carlo simulations sampling the dynamics according to a fixed topology structure, where the energetic neighbors of each particle of a cluster are constrained to be the same during the dynamics and their bond elongation bounded by a cutoff set to $\simeq~{}1.3~{}{\AA}$. In Fig.3 we show the energy per particle of clusters, $E/n$, as a function of the cluster size $n$ for the crystals obtained following method (iv). The curves are obtained averaging for each size $n$ over 10 different uncorrelated simulation times and 50 independent growth processes. From Fig.3 we can see that clusters of Ice 0 are more stable respect to clusters of ice I up to the size $n_{eq}^{*}\simeq 30$.

Considering also the contribution to the energy of clusters coming from the hydration layer (method (v)), we observe (Inset of Fig.3) that the maximum size up to which clusters of Ice 0 are more stable respect to ice I is reduced further respect to the previous cases up to the value $n_{eq-hy}^{*}\simeq 24$.

This behavior is similar to what has been observed for the stability of 13-particles icosahedron respect to 13-particles fcc and hcp clusters interacting through the Lennard-Jones potential when embedding clusters in liquidlike environment Mossa and Tarjus (2003).

III.1.1 Clusters morphology

In Fig.4 we show typical clusters of Ices Ic, Ih and Ice 0 projected in the xy (left panels) and xz (right panels) planes for clusters of size $n=40$ that is the typical value at which the curves of energy per particle associated to Ice 0 and Ic/Ih cross. Particles are colored according to the growth sequence such that the seeding particle is the darkest one and the last particle attached to the cluster is the lightest one. For Ic the xy plane corresponds to the direction [111], while xz to [1̄10]. For Ih xy is the basal face corresponding to the Miller indices {0001}, while xz is the secondary prism face corresponding to the Miller indices {11-20}. For Ice 0 the long side of the unit cell $c$ is aligned along the z axis.

The first important difference between Ice 0 and Ic/Ih and L is that Ice 0 grows preferentially along the (100) plane, while other clusters grow with spherical symmetry. We understand the planar symmetry of growing clusters of Ice 0 in terms of its neighbouring bond lengths. While Ic and Ih are characterized by a unique length, which is the distance between two neighbouring particles, Ice 0 is characterized by three different neighbouring distances ($d_{s}$, $d_{i}$, $d_{l}$), very close to each other. Two of these three lengths ($d_{s}$ and $d_{l}$) are involved in the formation of planar structures. In the snapshots of Ice 0 in Fig.4, white bonds are associated to the short distance ($d_{s}\simeq 2.67$ Åfor $\rho=0.976$ g/cm3), while blue bonds to the long distance ($d_{l}\simeq 2.73$ Åfor $\rho=0.976$ g/cm3, which is about $2\%$ longer than $d_{s}$). The intermediate distance ($d_{i}\simeq 2.68$ Åfor $\rho=0.976$ g/cm3) connects particles along the z direction forming distorted hexagonal rings and then it doesn’t play any role in the formation of clusters built following the minimum energy rule, as discussed below. From Fig.4 we can see that clusters of Ice 0 grow forming rings of pentagons composed of $d_{s}$ (white) and $d_{l}$ (blue) bonds. For this reason the energy of Ice 0 is lower than Ic/Ih already for clusters of $n=5$ particles which form pentagons in Ice 0, that is a ring, while Ic/Ih can form the first ring only for $n=6$. Each ring forming in the cluster gives an extra term in the calculation of its total energy. Indeed, the number of bonds $N_{b}$ forming in a cluster is $N_{b}=n+N_{r}-1$, where $N_{r}$ is the number of rings. Counting the number of bonds for the different growing clusters, we find that on average Ic and Ih can form more bonds than Ice 0 only for $n\gtrsim 40$, which corresponds to the typical value at which the curves of energy per particle associated to Ice 0 and Ic/Ih cross (highlighted by the orange vertical band in Fig.2).

Instead of adding the strongest bonded particle to a cluster, we can add particles with increasing radial distance from a central particle. In this case we get for Ice 0 clusters that are similar to those obtained following the minimum energy rule up to a size of $n=13$, while for $n>13$ we get a spherical symmetric structure whose energy becomes in average larger than the energy of clusters of Ic/Ih (computed with this same rule) for $n\gtrsim 22$ (see Inset of Fig.2(a)). This value of $n$ is smaller than the critical value of $n$ ($\gtrsim 40$) found for clusters of Ice 0 growing according to the minimum energy rule because for equal number of particles, a spherical symmetric cluster forms a smaller number of rings than a planar symmetric cluster for the Ice 0, while for Ic/Ih the growth is always spherical symmetric following both growing rules.

The reason why the number of bonds $N_{b}$ can be used to compute the difference in the energy content between clusters of Ic/Ih and Ice 0 is related to the two characteristic lengths in Ice 0, $d_{s}$ and $d_{l}$, being very close to each other (i.e. $(d_{l}-d_{s})/d_{s}\ll 1$). If the characteristic length of neighbouring particles forming a cluster is only one, as in the case of Ic and Ih, this length will be the minimum of the potential energy (note that for perfect crystals the three-body term of the mW potential is zero). If, on the other hand, the characteristic lengths of neighbouring particles forming a cluster are more than one, in general they will be different from the length corresponding to the minimum of the potential energy and then the eventual increase in the number of bonds will not balance the increase in energy of each individual bond. In the case of Ice 0, since the two characteristic lengths $d_{s}$ and $d_{l}$ forming clusters are close to each other and close to the distance corresponding to the minimum of the potential energy (which is 2.71 Åfor $\rho=0.976$ g/cm3), the energy of each bond is very close to that of bonds in Ic/Ih clusters such that the value of $N_{b}$ is a clear indication of the energy content of a cluster.

The mechanism responsible for the lower energy content of small clusters of Ice 0 respect to clusters of Ic/Ih now described is the same causing the icosahedron to have a lower energy content than fcc and hcp clusters in simple one-component liquids with particles interacting through spherically symmetric pair-potentials, like the Lennard-Jones. Indeed, fcc and hcp structures are characterized by a single length, which is the minimum of the potential. The icosahedron is characterized by two different lengths: the distance between the central particle and its 12 neighbours, and the edge of the 20 equilateral triangles forming the surface of the icosahedron where the vertices are the 12 neighbouring particles, which is about $5\%$ longer. The total number of bonds forming the icosahedron is 42. In the case of fcc or hcp structures, the 13 particles (central particle plus its 12 first near neighbours) composing them can form only 36 bonds. Although the 42 bonds formed in the icosahedron have a slightly higher energy respect to the 36 bonds formed in fcc and hcp, the total energy of a 13-particles cluster of icosahedron is smaller (by 8.4$\%$ in the case of a Lennard-Jones liquid Frank (1952)) than that of a 13-particles cluster of fcc or hcp. Another analogy between Ice 0 and the icosahedron is their tendency to form more clusters as the temperature is lowered Steinhardt et al. (1983); Mossa and Tarjus (2003). On the other hand, an important difference between Ice 0 and the icosahedron is that the former can form a crystal, even if it is metastable, while the latter, due to its five-fold rotational symmetry which is not compatible with translational periodicity, cannot tile the space and form a crystal so that it leads to frustration in simple one-component liquids Tarjus et al. (2005).

The stability range in size of clusters of Ice 0 respect to ice Ic/Ih will depend on the particular water model adopted. For example, the mW model respect to molecular models underestimates the lattice energy of Ice 0 relative to ice Ic/Ih Quigley et al. (2014). However, since in our study the stability of Ice 0 respect to ice Ic/Ih originates from topology (that is the number of bonds formed in a cluster), we expect that a different water model should give a similar result for the energy per particle difference between Ice 0 and ice Ic/Ih (that is Ice 0 more stable than ice Ic/Ih up to a specific cluster size), while the energy per particle of each ice could change respect to the mW model.

Another aspect to take into account when considering fully atomistic realizations of clusters is the configuration of the hydrogen bond network. In the case of small clusters, the HB network rearrangements take place in the hydration layer, as verified for small aggregated domains in the case of an atomistic model of water in Ref. Matsumoto et al. (2007).

In conclusion, we find that even a tetrahedral liquid like water can be characterized by competing crystalline clusters. In our case clusters of Ice 0 are even more stable than isolated clusters of Ic/Ih up to a critical size of $n^{*}\sim 40$ for clusters growing following the minimum energy rule (see Fig.2), while up to a critical size of $n^{*}\sim 22$ for clusters growing following the closest distance rule (see Fig.2).

III.2 Cluster growth and grain boundaries

In this section we address the growth of different crystalline seeds, both from critical nuclei (with the seeding technique and Umbrella Sampling simulations), and as planar interfaces. We will show how the growth of each polymorph is associated with a specific coherent grain boundary. Ice Ih is hindered by stacking disorder formed via a twin boundary between its basal plane and the (111) surface of Ice Ic, giving rise to stacking disordered Ice I${}_{\text{sd}}$. Ice Ic growth is instead hindered by a five-fold coherent grain boundary. Ice 0 has a coherent grain boundary between its (100) face and the (310) plane of cubic ice.

III.2.1 Growth of hexagonal ice

Nucleation of the stable Ice I phase in the mW model has been extensively addressed in previous studies Moore and Molinero (2011b); Li et al. (2011); Reinhardt and Doye (2012); Quigley (2014); Espinosa et al. (2016c); Hudait et al. (2016); Lupi et al. (2017); Cheng et al. (2018). The free energy difference between the two polytypes, cubic and hexagonal ice, was found to be very small, $\Delta G_{\text{ch}}=2\pm 1.5\,\text{J\,mol}^{-1}$ at $T=200$ K Quigley (2014). No free energy difference was indeed found in similar studies Romano et al. (2014); Russo et al. (2014). The difference in nucleation barriers between an hexagonal and cubic nucleus at $T=240$ K was also found to be negligibly small, $2\pm 2\,\text{J\,mol}^{-1}$ Cheng et al. (2018), which is equivalent to a difference of about $1\,k_{B}T$.

Figure 5 shows the result for our calculations of the nucleation barriers at $T=218$ K with the Umbrella Sampling method detailed in the Methods section. The height of the nucleation barrier is the same ($\beta\Delta G=27.2\pm 0.25$, where $\beta=1/k_{B}T$) for both a nucleus of Ice Ih and Ic. We see that the nucleation barrier difference between the two ices is indeed negligible at lower temperatures, and within the Umbrella Sampling method’s precision.

Both Ice Ih and Ic are made of layers composed of six-membered rings, and differ only in the way these layers are stacked. In Ice Ih the oxygen positions are the same between all layers (parallel to the basal plane), while in Ice Ic each layer is shifted by a distance equal to half the diameter of a ring. Given the low free energy difference between the cubic and hexagonal polytypes, the growth of ice Ih along the direction perpendicular to the basal plane is controlled by the free energy cost of introducing a twin grain boundary that changes the order of the stacking. It was recently pointed out that the presence of stacking disorder gives rise to a new crystal, ice I${}_{\text{sd}}$ which retains only a three-fold rotational symmetry instead of the full six-fold symmetry of ice Ih Malkin et al. (2012). The signature of stacking disorder was found in many X-ray and neutron diffraction patterns which were previously attributed to cubic ice Malkin et al. (2012), leading to the idea that at low temperatures ice always nucleates as Ice I${}_{\text{sd}}$ rather than cubic ice Malkin et al. (2012); Kuhs et al. (2012). We will challenge this idea in the next Section.

The free energy cost of forming a stacking fault is offset by an increase in configurational entropy coming from the different ways to arrange the layers. Previous studies have addressed both these terms Quigley (2014); Lupi et al. (2017); Cheng et al. (2018). The most recent estimates Cheng et al. (2018) put the stacking fault free energy per unit area between $\gamma_{\text{sf}}=0\,\text{mJ m}^{-2}$ at $T=180$ K and $\gamma_{\text{sf}}=0.11\,\text{mJ m}^{-2}$ at the melting point.

We estimate the value of $\gamma_{\text{sf}}$ from Umbrella Sampling simulations for an Ice I${}_{\text{sd}}$ nucleus in which the cubic and hexagonal layer alternate between each other (see, e.g., Fig. 9). The barrier is shown with green diamond symbols in Fig. 5. We emphasise that this is not the nucleation barrier for Ice I${}_{\text{sd}}$, as it does not take into account the entropy contribution coming from different realizations of the stacking disorder. The curve in Fig. 5 is instead the nucleation barrier for a single realization of stacking disorder, from which we can estimate the cost of forming the twin boundary defect. The free energy difference between the stacking disordered nucleus and the hexagonal (or cubic) nucleus is shown in the inset of Fig. 5 and corresponds to a difference of $0.3\pm 0.1\,k_{B}\,T$. The estimation of the difference in free energies between stacking disordered and cubic and hexagonal ice was done by fitting the nucleation barriers with the classical nucleation theory expression in the range $60<N<120$ around the top of the barrier. The expression is $\Delta G=a\,x+b\,x^{2/3}$. We obtain maxima at $\Delta G=27.6$ for the stacking disordered seed, and $\Delta G=27.3$ for both cubic and hexagonal seeds. The propagated error from fitting the classical nucleation theory expression is low ($~{}1\%$), but to account for systematic effects coming from the Umbrella Sampling method, we have increased this error to be of the same order of the deviation between barriers computed from independent trajectories, which is of the order of $0.1\,k_{B}\,T$.

To the value $0.3\pm 0.1\,k_{B}\,T$, one would need to subtract the bulk free energy difference between the cubic and hexagonal crystals, but as one can see from the inset of Fig. 5, this difference falls below our resolution at $T=218$ K, so we simply set it to zero. We then measure the surface tension between the hexagonal and cubic layers during the simulation and estimate it as $\gamma_{\text{sf}}=0.16\pm 0.05\,\text{mJ m}^{-2}$. Ignoring the bulk free energy difference makes our value an upper bound to the true surface penalty, so overall our estimate is in agreement with the results of Ref. Cheng et al. (2018).

III.2.2 Growth of cubic ice

We first consider unbiased simulations from cubic ice seeds following the seeding technique. We first determine the equilibrium densities of both the crystal and supercooled liquid phase, and prepare crystalline nuclei at different system size. We then run isobaric simulations of $N=10000$ water molecules at $T=218$ K and find the size that is critical, where the nucleus has the same probability to either grow or shrink. The inset of Fig. 6 shows the time evolution of different runs for seeds of size $n=110$ molecules, where half the trajectories (continuous blue lines) see the nucleus melting, and half the trajectories (lines with symbols) see the nucleus growing. We thus determine the size $n\sim 110$ to be critical, in excellent agreement with the Umbrella Sampling simulations shown in Fig. 5, and study the growth process of the initial ice I${}_{\text{c}}$ seed. Via bond orientational order parameters we determine the fraction of cubic ice present in the nucleus, and compute the cubicity parameter, defined as $f_{\text{cubic}}=n_{\text{cubic}}/(n_{\text{cubic}}+n_{\text{hexagonal}})$. The main panel in Fig. 6 plots the cubicity parameter as a function of the growth steps for the three trajectories that experience crystal growth, showing that cubicity steadily decreases during growth in all cases. These results confirm that at these sizes, crystal growth involves the formation of stacking faults.

In order to study how the process of stacking disorder proceeds, in Fig. 7 we plot the $\left<110\right>$ projection of the crystal growth process from an initial cubic seed (panel (a)). In the growth stages (panels (b)-(d)) we notice the development of stacking disorder. A crucial difference between Ice I${}_{\text{sd}}$ and Ice I${}_{\text{c}}$ is that the former has stacking disorder only in one direction, while the latter, due to its cubic symmetry, has a family of equivalent planes $\{1,1,1\}=(111),(\bar{1}11),(1\bar{1}1),(11\bar{1})$. The growth along two of such planes is highlighted in Fig. 7(d), where the order of the stacking is indicated with “c” for cubic arrangement, and “h” for hexagonal arrangement.

The appearance of random stacking along different crystallographic directions poses a problem with cubic seeds trying to grow. As seen in Fig. 7, the interface between the two stacking directions is incompatible with crystalline order, i.e. a grain boundary forms. We observe the formation of such grain boundaries in all 10 trajectories that are started with a cubic seed.

While Ice I${}_{\text{sd}}$ only pays a free energy penalty for stacking, ice I${}_{\text{c}}$ must additional incur in the free energy penalty associated with the formation of grain boundaries. For incoherent boundaries, this penalty would be too high, and cubic ice would then play no role in crystal nucleation. Instead we always observe the formation of a coherent grain boundary with five-fold symmetry. This particular structure was first noted in Ref. Li et al. (2011). We can already get a hint of this process by considering the angle $\theta$ at which the $\{1,1,1\}$ planes meet in Fig. 7d, which is $\theta=\cos^{-1}(1/3)$. We notice that this angle, $\theta\sim 70.5^{\circ}$, is very close to the value $72^{\circ}$ that corresponds to five-fold symmetry. The five-fold grain boundary is depicted in Fig. 8. It is a defect in which five cubic ice grains with hexagonal stacking on co-planar {111} planes merge together in a coherent way.

It is important to notice that this particular grain only appears as a consequence of stacking disorder of ice Ic, and cannot originate when nucleation of Ice Ih or I${}_{\text{sd}}$ occurs instead (since there is only one stacking direction in these forms of ice). The formation of this grain boundary is thus an indirect piece of evidence towards nucleation of ice Ic compared to Ice Ih or I${}_{\text{sd}}$. In the following we will show that nucleation simulations indicate the presence of this structure in at least $50\%$ of trajectories.

To study nucleation events without seeding, we use CNT-US simulations at $T=218$ K. According to Eq. 1, the two unknown parameters that need to be determined are the driving force for crystallization, $\Delta\mu$, and the critical nucleus size, $n_{c}$. The first one is obtained via thermodynamic integration from the melting point of mW ice I, $T=275.6$ K. The parameter $n_{c}$ is instead determined as follows: several simulations at different values of $n_{c}$ are run independently, and the fluctuations of the largest nucleus size is monitored. When $n_{c}$ is close to the true value of the critical size, one observes large fluctuations of the size $n$, with the system sampling from small to large crystalline sizes in the course of the same simulation. From these simulations one can reconstruct the free energy barrier, and thus obtain the true value of $n_{c}$. We run $10$ independent trajectories at the pre-determined values of $\beta\Delta\mu$ and $n_{c}$.

Figure 9 shows representative configurations in which only nuclei of size $n\sim 300$ are plotted, and colored according to their phase (orange for the hexagonal crystal, and blue for the cubic crystal). The results reveal three possible outcomes of nucleation. In the top row we see the formation of a five-fold symmetric grain boundary of Fig. 8. In the bottom row (left panel) we see nuclei with stacking disorder in at least two directions, which also originates from the cubic phase. Finally, the bottom row (right panel) shows nuclei with stacking in only one direction, which corresponds to Ice I${}_{\text{sd}}$ nucleation.

The abundance of stacking disorder in two or more directions, and in particular the presence of the five-fold symmetric grain, provides indirect evidence that nucleation indeed involves also the cubic phase. The cubic phase is able to avoid paying the full free-energy cost of an incoherent grain boundary by forming coherent five-fold symmetric grain boundaries, in which the hydrogen bond network is only distorted and not disrupted. The five-fold symmetric grain boundary forms a coherent interface at the expense of some elastic energy. We will now attempt a first free energy estimate of the grain boundary cost.

We tried obtaining a free energy barrier for crystals containing the five-fold grain boundary, but the nucleus fluctuations always favored growth on one of the cubic sides, until eventually the grain boundary disappeared. This could be remedied by introducing additional biasing fields that stabilize the grain boundary, but here instead we use a simpler approach, and leave a more detailed estimate for future studies.

We decrease the temperature to $T=204$ K, in a regime where nucleation occurs spontaneously in the system. We then estimate the free energy barrier of nucleation via mean first passage theory Wedekind and Reguera (2008); Russo et al. (2013). The mean first passage time $t_{\text{fp}}(n)$ is defined as the average time elapsed until the appearance of a nucleus of size $n$ in the system. For homogeneous nucleation in the steady-state, the mean first passage time is given by

[TABLE]

where $k$ is the nucleation rate, $n_{c}$ is the critical nucleus size, erf is the error function, and $c$ is related to the curvature of the free energy barrier at its maximum (i.e. the Zeldovich factor), $c=\sqrt{\Delta G^{\prime\prime}(n_{c})/k_{B}\,T}$.

We run $100$ homogeneous nucleation trajectories at $T=204$ K until the appearance of a super-critical nucleus. We then divide the trajectories in two sets depending on whether the nucleus has formed a five-fold symmetric grain boundary or not. To identify the grain boundary we use bond orientational order parameters to detect the atoms in the inner five-membered ring of the defect (depicted in cyan in Fig. 8). These atoms are identified when the following two conditions are satisfied: $Q_{4}<0.1$ and $Q_{12}>0.2$. Both $Q_{4}$ and $Q_{12}$ refer to coarse-grained bond orientational parameters Lechner and Dellago (2008) (see Ref. Tanaka et al. (2019) for a review of these methods).

The computed mean first passage times for nuclei of size up to $n=100$ are plotted in Fig. 10, distinguishing between nuclei with the five-fold grain (red squares) and without it (black circles). We notice a very clear distinction between the two different sets, with nuclei with the five-fold grain having shorter mean first passage times. This already reveals that the nucleation barrier for these states is lower compared to the nuclei without the grain. By using Eq. 2 we fit both curves, and display the results as continuous lines in Fig. 10. The results for the fit are (the superscript $5f$ refers to structures with the five-fold grain boundary): $n_{c}^{5f}=41$, $n_{c}=44$, $k^{5f}/k=1.44$, $c^{5f}=0.040$, and $c=0.038$. We notice that both the critical size and the free energy curvature are about the same in both cases, while the main difference is in the nucleation rate. We can thus estimate the free energy difference between the nuclei with and without a five-fold grain boundary $\Delta G^{5f}$ via the following relation: $\beta\Delta G^{5f}=-\log(k^{5f}/k)=-0.37$, or $\Delta G^{5f}=-0.6KJ/mol$. As discussed above, for small nuclei (here $n_{c}\sim 40$) the energetic cost of forming the five-fold grain boundary is very small, so the free energy advantage comes from the entropic gain in the coherent stacking in multiple directions. For much bigger nuclei, we expect the strain on the bonds far from the grain core to make $\Delta G^{5f}$ big and positive.

Overall, the abundance of five-fold grain boundaries in both direct nucleation and umbrella sampling simulations, and the absence of a significant free-energy penalty for small nuclei, shows that the nucleation of the cubic phase occurs distinctly and at least with the same probability of the nucleation of Ice Ih and I${}_{\text{sd}}$. The free-energy cost of the five-fold grain boundary is expected to grow considerably with increasing nucleus size, due to the strain in the bonds far from the grain core, and this is likely to stop the stacking of cubic ice in multiple directions. It should thus be more advantageous, for cubic seeds above a certain size, to have stacking disorder in only one direction, thus converting them into Ice I${}_{\text{sd}}$. This suggests another mechanism, besides just the difference in bulk free energy between the cubic and hexagonal phases, for the conversion of cubic ice in Ice I${}_{\text{sd}}$, that is governed by the size dependence of the free energy cost of the five-fold grain boundary.

III.2.3 Growth of Ice 0

In this section we consider the growth of Ice 0 nuclei with the seeding technique at $T=208$ K. We run simulations with an Ice 0 seed and determine the critical size to be around $n_{c}=400$ molecules. We then run 100 independent simulations and focus on the trajectories where the size of the nucleus grows. We observe two different scenarios, occurring with similar probability (number of trajectories that follow one growth mechanism or the other): i) growth of a defect free ice 0 crystal; ii) cross-nucleation of ice Ic on the ice [math] seed. A snapshot from a configuration following scenario ii) is shown in Fig. 11. While the first scenario is the expected growth mode for a metastable nucleus, the observation of cross nucleation of cubic ice poses some interesting questions. Firstly, cross nucleation at such small sizes implies that the surface penalty between ice [math] and cubic ice should be considerably low, which is unexpected given the differences between the two different crystalline structures. Secondly, the question remains why only the cubic form is cross-nucleated, and not the hexagonal phase.

To answers these questions, we run growth simulations with the direct coexistence method, where a slab of Ice 0 is prepared in contact with the supercooled liquid phase. Also in this case we observe the cross-nucleation of defect free cubic ice from the (100) plane of Ice 0. The (100) plane is the same that minimizes the cluster energy, as we described in Fig. 4.

In Fig. 12 we plot a snapshot from such configurations. Panel (d) shows that a coherent grain boundary is formed between the face (100) of Ice 0 (panel (b)) and the (310) face of cubic ice (panel (a)). The two faces join coherently (without breaking of hydrogen bonds), as shown for example in panel (c). From these equilibrium simulations we measure a very small strain (around $2\%$) of the cubic unit cell compared to its equilibrium value close to the interface. This signals a very low surface free energy cost of the defect. The cost of coherent grain boundaries between two solid surfaces of the same material (mW) originates almost entirely from the strain energy. We can estimate this cost as

[TABLE]

where $e_{I}$ is the energy of an interface atom, while $e_{0}$ and $e_{c}$ are respectively the energy of atoms in the bulk Ice 0 and Ic phases. We use our simulations to estimate these energies as: $e_{I}=94.71$ kJ/mol, $e_{0}=96.23$ kJ/mol, $e_{c}=97.18$ kJ/mol, for $\Delta E_{\text{c0}}=1.99$ kJ/mol. $e_{I}$ is obtained by measuring the average energy of the molecules at the interface between ice 0 and cubic ice from the direct-coexistence simulations . Since the interface is coherent, it is very sharp (one particle width), and there is no ambiguity in the selection of the correct particles. $\gamma_{\text{c0}}$ is obtained by assuming (as usually done for coherent interfaces) that the dominant contribution to the free energy is energetic (given by $\Delta E_{\text{c0}}$), and then dividing the total surface energy for the area of the surface. This leads to a surface energy of stacking between the (100) plane of Ice 0 and the (310) plane of cubic ice of $\gamma_{\text{c0}}=29.67$ mJ/m2. This energy is much larger than the stacking surface energy of cubic and hexagonal layers in Ice I (that we estimated earlier as $\gamma_{\text{sf}}\sim 0.16$ mJ/m2 at similar thermodynamic conditions), but is well below the limit of $200$ mJ/m2 that is conventionally used to separate coherent from semi-coherent grains.

We observe that the low cost of the grain boundary also makes it difficult to obtain free energy barriers of Ice 0 (as shown for the other phases in Fig. 5), as we always observe the appearance of cubic ice in the Umbrella Sampling simulations.

One important consequence of this grain boundary is the fact that the (100) face of Ice 0 appears to be an ideal template for the cross nucleation of Ice Ic. Our simulations always show the absence of hexagonal stacking, and the formation of a crystal with cubicity of $c=1$. The reason for this is that the grain boundary involves a high-order face of cubic ice, the (310) face, and so the new cubic phase grows along a direction incompatible with stacking faults.

IV Conclusions

In this article we have focused on crystalline clusters in the mW model of water. In the spirit of Frank’s original argument, we have shown that there is a competition between different crystalline clusters, and that the most energetic clusters at small sizes include five-membered rings. Here we have demonstrated this with the Ice 0 structure, which is a metastable crystal for the mW potential. Five-membered rings allow small clusters to increase the ratio between number of bonds and number of particles in the cluster, without significant penalties to the bond energy. As the clusters grow in size, six-membered rings (as found in the stable Ice I forms) become favourable due to their optimal bond energy. Frustration effects between local structures can suppress the nucleation of crystals, and are thus linked to glass-forming ability Ronceray and Harrowell (2017); Russo et al. (2018a); Wei et al. (2018); Tanaka et al. (2019). The competition between different local structures also plays an important role in polymorph selection, i.e. which metastable phase is nucleated.

To study the growth of crystalline clusters in supercooled conditions we have used both direct simulations, seeding, and umbrella sampling biasing techniques. In particular we have highlighted the role of grain boundaries in the growth process of the clusters. The lowest energy grain boundary is the stacking fault between hexagonal and cubic layers that can appear both on the basal plane of Ice Ih and the (111) plane of Ice Ic, which we estimate as $\gamma_{\text{sf}}=0.16\pm 0.05$mJm*-2* at $T=218$ K. Given the negligible bulk free energy difference between Ic and Ih, this predicts, at low temperatures, the formation of Ice I${}_{\text{sd}}$ as the first stage of nucleation, favored by the configurational entropy associated with its disordered stacking sequence. But our results show also that cubic ice still plays an important role in ice nucleation, thanks to another coherent grain boundary. Differently from Ice Ih and I${}_{\text{sd}}$, cubic ices can generate stacking faults in four different directions. The associated cost of intersecting these stacking faults, which would otherwise suppress the nucleation of small cubic crystallites, is instead accommodated by the formation of a coherent five-fold symmetric grain boundary. At deep supercooled conditions ($T=208$ K) we estimate the free energy difference between nuclei with and without the five-fold grain boundary to be of the order of $0.37\,k_{B}T$ in favor of the first ones. The abundance of configurations with the five-fold grain boundary is an indirect measure of the importance of the nucleation of cubic ice. Due to its five-fold symmetric nature, the grain boundary cost will increase rapidly with nucleus size, due to the bond strain far from the core of the boundary. We can thus expect that the grain boundary will also be a contributing factor in selecting one growth direction for cubic ice, which would then gradually transform in Ice I${}_{\text{sd}}$, in accordance with Ostwald’s step rule.

We have also considered seeding simulations of Ice 0 nuclei, and found that its growth is also controlled by a coherent grain boundary. The new grain boundary is formed between the (100) face of Ice 0 and the (310) face of cubic ice. We estimated the surface free energy of the grain boundary to be $\gamma_{\text{c0}}=29.67$ mJ/m2. The fact that the grain boundary involves a high order face of cubic ice, makes the (100) face of Ice 0 an ideal template for the cross-nucleation of pure cubic ice. Indeed our simulations show growth of cubic ice without stacking faults.

AcknowledgementsJ.R. and F.L. acknowledge support from the European Research Council Grant DLV-759187. J.R. acknowledges support from the Royal Society University Research Fellowship and a travel support from Institute of Industrial Science (IIS), University of Tokyo for a short stay at IIS. This study was partially supported by Grants-in-Aid for Specially Promoted Research (Grant No. JP25000002) and Scientific Research (A) (18H03675) from the Japan Society for the Promotion of Science (JSPS).