Anharmonic nuclear motion and the relative stability of hexagonal and cubic ice

We use extensive first-principles quantum mechanical calculations to show that, although the static lattice and harmonic vibrational energies are almost identical, the anharmonic vibrational energy of hexagonal ice is significantly lower than that of cubic ice. This difference in anharmonicity is crucial, stabilising hexagonal ice compared with cubic ice by at least 1.4 meV/H2O, in agreement with experimental estimates. The difference in anharmonicity arises predominantly from molecular O-H bond stretching vibrational modes and is related to the different stacking of atomic layers.


I. INTRODUCTION
Six-fold symmetric snow crystals are formed from hexagonal ice (Ih), which covers about 10% of the Earth's surface and plays a prominent role in determining its climate [1][2][3]. A cubic form of ice also occurs in nature, but is very rare [4,5]. The structures of pure hexagonal and cubic ice differ only in the stacking of layers of tetrahedrally coordinated water molecules (see Fig. 1). Yet, hexagonal ice is thermodynamically more stable than cubic ice, with experiments indicating the difference in stability to lie in the meV/H 2 O range [6][7][8][9][10][11][12][13][14][15]. There is a growing realisation that real "cubic ice" typically contains many stacking faults and is not pure cubic ice (Ic) as originally suggested by König [16]. Stacking faulted ice is a highly complex material, whose nature and properties depend heavily on the free energy difference between pure Ih and Ic.
The very similar free energies of Ih and Ic have so far prevented state-of-the-art first-principles quantum mechanical calculations from explaining the stability of Ih. Both density functional theory (DFT) and diffusion quantum Monte Carlo studies have found that Ih and Ic are almost degenerate in energy when nuclear motion is neglected [17]. Our calculations show that the harmonic zero-point vibrational energies of Ih and Ic are large, at roughly 700 meV/H 2 O, but they are almost identical. Consequently, when averaged over different protonorderings, the two phases are almost degenerate when harmonic vibrations are included (see Fig. 2). However, the small mass of hydrogen gives rise to large amplitude vibrations and large anharmonic effects.
A substantial body of theoretical work exists on water and ice, based on force-field path-integral molecular dynamics (PIMD), first-principles classical molecular dynamics (MD), and first-principles vibrational calculations in the quasi-harmonic approximation. This work has led to significant successes in understanding the important role of quantum vibrations and anharmonicity for various phenomena observed experimentally. Examples include, but are not limited to, (a) the isotope effects, e.g., on the melting temperature, of Ih upon going from protiated to deuterated ice [18][19][20], (b) accurate O-H bond lengths and infrared O-H stretching frequencies in water [21][22][23], (c) reproduction of the anomalous thermal expansion, and the isotope effect on the volume in Ih [24] and (d) the heat capacity of water [25]. Attempts to calculate the relative stability of Ih and Ic have either relied on empirical force-fields such as TIP4P [26,27] or have lacked an accurate description of anharmonicity [17,28]. TIP4P has since been shown to produce incorrect protonordering energetics and an incorrect static lattice energy difference between Ih and Ic compared to highly accurate diffusion Monte Carlo methods [29]. Moreover, no successful attempts to explain the origin of the greater stability of Ih have been made. With our fully anharmonic, first-principles DFT study we show that the inclusion of accurate anharmonic quantum nuclear motion is decisive in stabilising Ih with respect to Ic, and relate the difference in stability to the different stacking of the atomic layers.
The water molecules in both Ih and Ic are tetrahedrally coordinated, each donating and accepting two hydrogen bonds and thus satisfying the "Bernal-Fowler ice rules" [31]. The oxygen sublattices of Ih and Ic arise from an ABAB and an ABC stacking of puckered layers of oxygen atoms, respectively. Correspondingly, the basic building blocks of Ih are chair-and boat-form hexamers, while Ic is built exclusively from chair-form hexamers (Fig.  1). Pure Ic has so far proven elusive. Experimentally synthesised "Ic" typically contains many stacking faults which strongly affect its physical and chemical properties [6,13,32]. Ice containing cubic sequences interlaced with hexagonal sequences is commonly known as stacking-disordered ice (Isd). Unlike Ih and Ic, which refer to a unique stacking arrangement of puckered layers, Isd refers to the infinite set of possible stacking sequences. This set smoothly connects Ih as one end member to Ic as the other. Isd has trigonal P 3m1 symmetry [33][34][35]. As of yet, it is unclear whether stacking disorder is kinetically or thermodynamically driven. A full understanding of Isd will require understanding the properties of Ih and Ic, including their relative stability. At the most basic level, the free energy difference between Ih and Ic is required to understand why Isd is found to anneal to Ih rather than Ic.
The fraction of cubic stackings of layers, or "cubicity", typically does not exceed around 60% [32] in experiments and both the fraction itself as well as the nature FIG. 1. The oxygen sublattices in (a) hexagonal ice, Ih, and (b) cubic ice, Ic, consist of ABAB and ABC stacked puckered layers (red, blue and green colours), respectively. Projections of the periodic structures along the direction orthogonal to the puckered layers are shown on the lower faces of the diagrams. Ih is built of the chair-and boat-form hexamers shown in a-I and a-II, respectively, while Ic is built exclusively from chairform hexamers as in b-I. The structures were visualised using VESTA [30].
of the stacking arrangements depend heavily on the synthesis pathway [13,32,36]. Crucially, ice synthesised via both homogeneous and heterogeneous freezing of (supercooled) water has a random stacking of cubic and hexagonal layers [32,36], which is consistent with a layer-bylayer growth mechanism. Heterogeneous freezing in particular is central to atmospheric and climate physics and, due to random stacking-disorder, depends vitally on the free energy difference between Ih and Ic. It is clear that Isd is an extremely important and highly complex material.
Molecular dynamics and Monte Carlo simulations using empirical ice potentials to model ice nucleation processes have successfully reproduced stacking-disorder (see [27,36,37] and references therein). However, they struggle, amongst other issues, with the accuracy of the empirical potentials in describing the melting temperature and the relative stability of Ih and Ic [38].
In the following we limit ourselves to the study of pure Ih and Ic.

II. COMPUTATIONAL MODEL
Atomistic simulations of proton-disordered systems such as Ih and Ic require sets of explicit atomic positions and a calculation must be performed for each protonordering studied. The number of proton-ordered, energetically quasi-degenerate structures allowed by the ice rules increases exponentially with the size of the simulation cell. This leads to Pauling's residual configurational entropy [39], S config , which has been confirmed experimentally [40,41]. For large systems with negligible surface effects the associated configurational free energies of Ih and Ic, ∆G config = −T S config , are almost identical [42,43]. We therefore neglect ∆G config in the following.
To gain an understanding of the effects of protonordering on the vibrational properties of ice, we consider 16 distinct proton-ordered eight-molecule Ih configurations as constructed by Hirsch and Ojamäe [44], and 11 distinct proton-ordered eight-molecule Ic configurations [17]. We also consider the "conventional" hexagonal, 12molecule P 6 3 cm Ih and quasi-cubic, eight-molecule P 4 3 Ic structures (numbers 13 and 1 in Fig. 2, respectively). Details of the proton-ordered structures and their numbering are provided in Supplementary Table 1.
We performed electronic structure calculations using plane-wave pseudopotential DFT as implemented in the castep code [45] (version 7.02). We employed the Perdew-Burke-Ernzerhof (PBE) [46] generalised gradient approximation to the exchange-correlation functional, and on-the-fly generated ultrasoft pseudopotentials [47] with core radii of 0.7Å and 0.8Å for the hydrogen and oxygen atoms, respectively. Supplementary Section V describes results obtained with other density functionals. We used a plane-wave energy cut-off of 1600 eV and Monkhorst-Pack reciprocal space grids of spacing less than 2π × 0.04Å −1 for all total energy calculations and geometry optimisations. The resulting energy differences between frozen-phonon configurations are converged to below 10 −4 eV/H 2 O, the atomic positions are converged to within 10 −5Å , and the residual forces to within 10 −4 eV/Å.
We obtained harmonic vibrational free energies from the harmonic frequencies calculated using the k-space Fourier interpolated dynamical matrix. The latter was obtained by Fourier transforming the real-space matrix of force constants constructed using a finite displacement method. Anharmonic vibrational free energies were calculated using the method described in [48], which has so far been successfully applied to high-pressure systems [49,50]. As in [48], we describe the 3N -dimensional BO energy surface (where N is the number of atoms in the simulation cell) by mapping 1D subspaces along the harmonic normal mode axes up to large amplitudes of four times the harmonic root-mean-square (RMS) displacements, where anharmonicity is important. We then reconstruct the 3N -dimensional BO surface from the 1D subspaces. The resultant representation of the BO energy surface is an approximation to the true 3N -dimensional BO energy surface. This approximation only weakly affects the free energy difference between Ih and Ic (see Supplementary Section VI). The 1D energy surfaces were fitted using cubic splines. The anharmonic vibrational Schrödinger equation was solved within a vibrational selfconsistent field (VSCF) framework. The vibrational wave function was expanded in a basis of simple harmonic oscillator eigenstates, and the inclusion of 25 states for each vibrational degree of freedom was found sufficient to obtain converged results.

III. RESULTS
Our calculations show that the static lattice energies of Ih and Ic, E static , vary by up to 5 meV/H 2 O with protonordering. This agrees with Refs. [17,44,51] and, more importantly, with Ref. [52], which evaluates DFT static lattice energies for 16 8-molecule orthorhombic, 14 12molecule hexagonal and 63 48-molecule orthorhombic Ih proton-orderings. This strongly suggests that our sets of proton-orderings provide a good representation of the distribution of energies in disordered ice.
We also find that the harmonic vibrational contributions to the free energies of different proton-orderings, ∆G har , vary by up to 2 meV/H 2 O. We have employed the VSCF method described above to calculate the anharmonic contribution to the vibrational free energy, ∆G anh , finding a variation between proton-orderings of up to 5 meV/H 2 O The free energies of Ih and Ic at the harmonic vibrational level, G har = E static + ∆G har , averaged over the proton-orderings, are virtually indistinguishable: also positive, are systematically lower than those of the Ic configurations, so that the total free energy, G anh = E static + ∆G har + ∆G anh , averaged over the different proton-orderings, is significantly lower for Ih than for Ic: The values obtained for ∆ Ic→Ih av G har and ∆ Ic→Ih av G anh depend significantly on the method of averaging. For example, using a Boltzmann distribution for the free energy of the proton-orderings leads to values of ∆ Ic→Ih Boltzmann G har ≈ 5.5 meV/H 2 O and 6.1 meV/H 2 O at 10 and 100 K, respectively.
It is noteworthy that, given cell volumes that are reasonably close to experiment, the differences in E static , ∆G har and ∆G anh between Ih and Ic depend only weakly on the details of the DFT calculations and in particular on the choice of exchange-correlation functional (see Fig.  3 and Supplementary Section V).
The most stable proton-ordered configurations of Ih and Ic, referred to as XIh and XIc, display free energy differences very similar to ∆ Ic→Ih av G har and ∆ Ic→Ih av G anh . However, the inclusion of anharmonic vibrational energies changes the relative stability of the different protonorderings. XIh is experimentally and theoretically known to have space group Cmc2 1 [44,53] (structure 15 in Fig.  2), a result confirmed by our calculations. A structure of space group I4 1 md has been proposed [17] (structure 2) for XIc on theoretical grounds. Our calculations support this proposal at the harmonic vibrational level, but the inclusion of anharmonic contributions suggests that Ic P c, P ca2 1 and P na2 1 (structures 7, 9 and 12) may also be strong candidates for XIc. Ref. [54] reports experimental evidence for Ic I4 1 md and P na2 1 in partially proton-ordered Ic via Fourier transform infrared spectroscopy. This lends significant support to our result that anharmonicity provides the decisive contribution to the energy differences between proton-orderings, since P na2 1 has a high free energy at the static lattice and harmonic vibrational levels and only becomes a lowfree-energy structure when anharmonicity is taken into account.
At typical experimental temperatures of below 100 K, the proton-ordering is largely frozen in. Consequently one cannot expect to measure a change in free energy corresponding to a transition from Ic to Ih in thermal equilibrium, but rather (assuming the Ic sample is annealed at low temperatures and consists mostly of XIc) a change in free energy corresponding to transitions from XIc to a proton-ordering of Ih, which is likely to be smaller than ∆ Ic→Ih av G anh . As indicated in Fig. 2, we evaluate a lower bound on the free energy difference as This lower bound is consistent with, but on the high side of, experimentally measured free energy differences of 0.3 − 1.6 meV/H 2 O [7][8][9][10][11][12][13][14][15]. Notably, the experimental value is rather uncertain, mainly because the free energy difference is very small, and because the Ic samples are typically not fully characterised in terms of stacking faults or proton ordering. As shown in Fig. 4 We have investigated the convergence of ∆G har and ∆G anh with respect to the size of the simulation cell using cells containing up to 192 molecules for Ih P 6 3 cm and 128 molecules for Ic P 4 3 as shown in Fig. 6.
The vibrational energies of the different protonorderings were calculated using 64-molecule cells for the harmonic vibrational energy and eight-molecule cells for the anharmonic energies. The latter is justified by the short-range nature of anharmonicity, evidence for which is described in section IV. Calculations using cells with up to 192 molecules for the P 6 3 cm Ih structure and 128 molecules for the P 4 3 Ic structure indicate that the difference between ∆G har for Ih and Ic is converged to within 0.1 meV/H 2 O using 64-molecule simulations cells. The anharmonic corrections converge analogously, and the difference between ∆G anh for the two structures is converged to within 0.2 meV/H 2 O using 64-molecule simulations cells. More details may be found in Supplementary Section II.

IV. DISCUSSION
The origin of the difference in the anharmonicities of Ih and Ic can be traced back to the libration (80 − 140 meV) and, predominantly, the molecular symmetric and antisymmetric O-H bond stretching modes (365 − 380 meV and 395−410 meV, respectively) indicated in Fig. 7. The vibrational density of states (DoS) and the distribution of anharmonic corrections over the vibrational frequencies show (Fig. 8) that the dominant anharmonic contributions arise from the O-H bond stretching modes, which correspond to large amplitude displacements of the hydrogen atoms relative to their neighbouring (essentially stationary) oxygen atoms. The comparatively small role of the oxygen atoms is confirmed by studying Ih and Ic analogues with fixed oxygen positions, which recover the value of ∆G anh observed for real Ih and Ic to within 5% (see Supplementary Table III). The O-H bond stretching modes contribute > 2/3 of the difference in anharmonicity between Ih and Ic, when averaged over protonorderings. Note that the energy difference between Ih and Ic shown in Fig. 8 increases at high frequencies.
Variations in vibrational frequencies with cubicity have recently been identified in infrared absorption experiments [6]. Carr et al. observed that the O-H stretching modes were shifted to higher vibrational frequencies with increasing cubicity of their Isd samples. They also observed an increasing broadening of the absorption peak. According to Carr et al., both trends are thought to be associated with the stacking disorder, which peaks at a cubicity of 50%. For samples with cubicities of 50% Carr et al. observed shifts of 28 ± 2 cm −1 and 13 ± 2 cm −1 for protiated and deuterated Isd, respectively. Our calculations reproduce the widening of the O-H stretching peak  scales roughly linearly with the degree of cubicity rather than the amount of stacking disorder, and is largest in pure Ic. In accurately reproducing the ratio between the blueshifts of the molecular stretching frequencies for the protiated and deuterated phases, our results also provide a good account of the effects of isotopic substitution.
The role of the high energy modes can be further illuminated by considering the H-H radial distribution functions ( Fig. 9 (b)) of Ih and Ic and RMS displacements of the protons (see Supplementary Table IV). The latter show that the harmonic vibrational amplitudes in Ic are about 1% smaller than in Ih. Yet, at the harmonic level the H-H radial distribution functions of Ih and Ic are essentially identical for distances less than about 3Å. Beyond 3Å the static structures of the two phases differ. Unlike the RMS displacements of the protons, the RDFs measure a two-particle quantity that describes correlated motions of pairs of protons. At the harmonic level the protons in Ic move less with respect to their equilibrium positions than in Ih, but they move by just as much with respect to each other. We note that the O-H and H-H RDFs for Ih shown in Figs. 9 (a) and (b) agree well with, e.g., the experimental RDFs in [55].
While the protons in both phases feel the same local environment, differences occur starting with the fourthnearest-neighbour protons (see Fig. 10). Also, for systems as small as 8 to 12 molecules ∆G anh is already about 3/4 of the converged value ( Supplementary Fig.  1). This system size limits the wavelength of the vibrational modes responsible for the difference in anharmonic energies to roughly the same distance as the separation of fourth-nearest-neighbour protons. Together these observations indicate that the influence of more distant nuclei   The RDFs for Ih P 63cm (blue) and Ic P 43 (red) are sampled from the anharmonic nuclear wavefunctions. For radii < 3Å the difference between Ih and Ic at the harmonic level (green), ∆g har ≡ g Ih har −g Ic har , are minimal. The difference at the anharmonic level (black), ∆g anh ≡ g Ih anh −g Ic anh , is small, but non-negligible. The features in ∆g har and ∆g anh for radii beyond around 3Å originate predominantly from differences in the static structures of Ih and Ic.
is small. Allowing for both chair-and boat-form hexamers, there are 12 distinct arrangements of fourth-nearestneighbour pairs of protons. Out of these 12 arrangements only 8 are realised in the proton-orderings we have considered. These are shown in Fig. 11. Three of these arrangements (numbers VI -VIII) are associated with boat-form hexamers of H 2 O molecules, which only exist in Ih, and 5 (numbers I -V) are associated with the chair-form hexamers of H 2 O molecules found in Ih and Ic. The RDFs in Fig. 11 show that, on going from Ic to Ih, arrangements II and IV are depopulated. For arrangements II and IV the displacement of the first proton of the pair from its equilibrium position along its hydrogen bridge bond leads to a large displacement relative to the second proton. Conversely, arrangements VI -VIII, for which the same displacement of the first proton leads to a far smaller displacement relative to the second, are populated. This explains why on average the protons in Ic move less with respect to their equilibrium positions than in Ih, while moving by just as much with respect to each other, resulting in the same G har as in Ih.
Going beyond the harmonic approximation, anharmonicity reduces the RMS vibrational amplitudes in Ih and Ic by around 1% and 2.5%, respectively, localising the nuclear wavefunction more in ice Ic than Ih. On the level of the collective vibrational modes, the localisation is typically driven by a strong quartic contribution to the respective BO energy surfaces (see Supplementary  Fig. 6). Moreover, at the anharmonic level the RDFs show that the protons in ice Ic move less with respect to each other than in Ih, instead of just moving less with respect to their equilibrium positions, as they do at the harmonic level. This difference in the relative motion of pairs of protons is the origin of the difference in ∆G anh between Ih and Ic. The larger effect of anharmonicity in Ic is again due to the stronger geometric "coupling" between pairs of protons.

V. PERSPECTIVE
Ih and Ic are important in various branches of science. Examples include climate modelling and the simulation of ice nucleation and formation, where cubic ice plays an important role and for which ∆G anh is an essential input parameter. As an example of the relevance in biological sciences, the benign shape of cubic ice crystals is of potential interest for cryopreservation [56]. Here we have demonstrated that accounting for anharmonic nuclear vibrations is central to understanding and correctly predicting the free energy difference between Ih and Ic. However, the importance of anharmonic vibrations in hydrogen-bonded systems reaches far beyond ice. An accurate treatment of anharmonicity is likely to be crucial in correctly describing the energy differences between very similar polymorphs of hydrogen-bonded molecular crystals which are important in, e.g., pharmaceutical materials science.
Calculating anharmonic vibrational energies in solids is a challenging computational task, which has only recently been successfully achieved using first-principles quantum mechanical methods. Anharmonic effects are particularly important for light elements, such as the hydrogen atoms in H 2 O. Anharmonic vibrations are also expected to be important at the surfaces of ice, and when impurities or other defects are present.