Violations of conventional ice rules lead to a quasi-one-dimensional hydrogen bonded ice phaseIn bulk ice, the Bernal-Fowler ice rules dictate that each water molecule in a stable ice phase should participate in four hydrogen bonds with its neighboring water molecules: 2 donated and 2 accepted hydrogen bonds25. Previous work on nanoconfined water at the DFT level has suggested that the thermodynamically stable monolayer phase is a square structure that still exhibits this 4-fold coordination9,12. However, more recent first-principles-level simulations have suggested a wide variety of stable phases of nanoconfined water spanning a broad range of different temperature and lateral pressure conditions14,15,16,17. Kapil et al.16 report a hexagonal phase stable below 0.1 GPa, a pentagonal phase stable from 0.1 to 0.5 GPa, and a flat-rhombic phase stable above 0.5 GPa for a confinement width of 5 Å. Reference16 also confirms the greater thermodynamic stabilities of these phases with respect to the square phase using the chemically accurate and robust Quantum Monte Carlo (QMC) calculation setup35,36. References15,17 explored even higher lateral pressures for a confinement width of 6 Å. These studies found that the stable phase beyond 0.5 GPa is a zigzag monolayer ice (ZZMI) phase, which is topologically equivalent to the flat-rhombic phase, except that the larger confinement width leads to a small buckling of the oxygen atoms. In the 15–20 GPa regime, Lin et al. report a zigzag quasi bilayer ice (ZZ-qBI) structure as the stable phase17. This phase also resembles the flat-rhombic phase, although the buckling in this phase results in two nearly distinct layers of water molecules.As shown in Fig. 1a and Table 1, none of the above thermodynamically stable nanoconfined ice phases satisfy conventional ice rules: all of these phases display fewer than four hydrogen bonds per water molecule. However, the water molecules in each of these phases violate the ice rules in different ways. In the low-pressure hexagonal phase, all water molecules exhibit a 3-fold hydrogen bonding coordination: half of the water molecules donate two and accept one hydrogen bond, and the other half donate one and accept two. Water molecules in the intermediate-pressure pentagonal phase exhibit diverse hydrogen bonding motifs, with coordination numbers ranging from 2 to 4. This results in an average hydrogen bond coordination of 3.2 in the pentagonal crystal structure. The high-pressure flat-rhombic (or the ZZMI) phase exhibits a single hydrogen-bonding motif in which each water molecule donates one and accepts one hydrogen bond. The resulting hydrogen bonding network contains just two hydrogen bonds per water molecule, deviating significantly from the conventional ice rules. The ZZ-qBI phase exhibits a similar network, with just two hydrogen bonds per water molecule. Since this phase is stable even at lateral pressures much higher than the flat-rhombic phase, this suggests that such extreme violations of the bulk ice rules are likely to occur in high lateral pressure conditions.Fig. 1: Nanoconfined ice phases and the role of vdW stacking of hydrogen bonded chains.a The four nanoconfined ice phases considered in this work from refs. 16,17, along with the lateral pressures used for simulating each phase. Solid lines between water molecules indicate hydrogen bonds. We highlight chains of hydrogen-bonded water molecules in each phase. For estimating Echain, we use the unique highlighted chains in the case of the flat-rhombic and ZZ-qBI phases. For the hexagonal and pentagonal phases, we report the average value and the standard deviation as the error bar, where the standard deviation is computed across different choices of chains, as described in the text. b For each phase, we compute the stabilization energy of a single molecule in the 0 K crystal structure, Elattice; the (average) stabilization energy of a water molecule within the hydrogen bonded chain(s), Echain; and the stabilization energy between chains in the crystal structure, Estack. The hatched bars show the contribution of vdW interactions to the corresponding stabilization energy. Source data are provided as a Source Data file.Table 1 First-principles calculations for four nanoconfined ice phasesThe flat-rhombic phase is thermodynamically stable in a sizeable temperature-pressure region of the phase diagram15,16,17 despite having just two hydrogen bonds per molecule. Considering that ice in general is primarily stabilized by a dense network of hydrogen bonds24,37, this stability over a broad range of conditions is highly unusual. To understand this unexpected behavior, we examine the structure of the flat-rhombic phase and look for similarities with other types of materials. As depicted in Fig. 1a, the flat-rhombic phase is characterized by a zigzagging hydrogen bond network that runs along quasi-one-dimensional chains15,17. These chains are stacked alongside each other without any hydrogen bonding between chains. This structure is similar to quasi-one-dimensional vdW materials32,33, which exhibit unusual electronic, optical, and vibrational properties. However, while these materials consist of covalently bonded chains held together by weak inter-chain vdW interactions, the chains in the flat-rhombic phase are stabilized by hydrogen bonding.To determine whether the flat-rhombic phase can be classified as a quasi-one-dimensional vdW material, we calculate the binding energy of the one-dimensional chains and study the role of vdW interactions in stabilizing these chains. In Fig. 1b and Table 1, we report three different types of binding energies of a water molecule in each phase: the stabilization energy of a water molecule in the corresponding crystal structure, Elattice; the stabilization energy of a water molecule within a continuous chain of hydrogen-bonded water molecules, Echain (see highlighted chains of water molecules in Fig. 1a); and the remaining stabilization energy from interactions between chains of water molecules in each phase, Estack. This means that the stacking energy is the stabilization energy of the lattice that is not captured by the chains alone: Estack = Elattice − Echain. In the hexagonal phase, we consider two distinct chain types in the ‘zigzag’ and ‘armchair’ directions and report the average binding energies. As can be seen in the pentagonal phase snapshot in Fig. 1a, there is no clear choice for a unique hydrogen bond chain in the pentagonal phase. This shows qualitatively that the pentagonal phase cannot clearly be separated into distinct hydrogen bond chains. Nonetheless, to perform the quantitative comparison shown in Fig. 1b, we averaged our energies over three different hydrogen bond chains in the pentagonal phase. The energies for these different chains were all extremely close to each other, suggesting that this quantitative comparison is robust with respect to the exact choice of hydrogen bond chain in the pentagonal phase. For the flat-rhombic and ZZ-qBI phases, we considered a single unique chain.In Fig. 1b, we observe that Echain contributes to approximately 50% or less of the crystal stabilization energy in both hexagonal and pentagonal phases. Therefore, these phases have high stacking energies between chains, stabilized primarily by hydrogen bonds, with vdW interactions (indicated by hatched regions) playing only a secondary role in stabilizing the chains. Conversely, in the flat-rhombic phase, the primary crystal stabilization source is from the chains. Here, the stacking energy is low, with roughly 90% of the stacking energy resulting from vdW interactions. This difference indicates that the interactions stabilizing water molecules in the flat-rhombic phase contrast starkly with those in the hexagonal and pentagonal phases. The flat-rhombic phase appears to form hydrogen bonds in one dimension while being stabilized by weak vdW forces in the other, classifying it as a quasi-one-dimensional vdW crystal. Similarly, the stacking energy of the ZZ-qBI phase is almost entirely driven by vdW interactions, meaning that vdW forces also stabilize the interactions between this phase’s hydrogen-bonded chains. In fact, in the absence of vdW stacking, the lattice energy is positive, meaning that the ZZ-qBI crystal structure would be thermodynamically unstable. This suggests that these vdW interactions between hydrogen-bonded chains play an essential role in stabilizing nanoconfined crystal structures at high lateral pressures. As our MLP isn’t trained for pressures beyond 8 GPa, we are unable to comment on the dynamical stability or metastability of the ZZ-qBI phase. Similarly, in order to comment on the true dynamical stability or metastability of any of these phases in the absence of vdW interactions, we would need to retrain an MLP at the revPBE0 level (in the absence of the D3 dispersion correction).Anomalous finite-temperature dependence of hydrogen bondingThe characteristic interactions in the flat-rhombic phase lead to anomalous finite-temperature hydrogen bonding structure and dynamics. In bulk or conventionally hydrogen-bonded ice phases, such as the hexagonal and pentagonal monolayer phases, the most probable molecular orientations are the ones that maximize the number of hydrogen bonds25, and the protons maximize their residence time in orientations with the maximum possible hydrogen bonds. However, in the flat-rhombic phase, even at finite temperatures (extending past 300 K), the most probable molecular orientations are not those that maximize hydrogen bonds but instead, those that maintain a distinct network of quasi-one-dimensional chains.To characterize the rotational motion of individual water molecules, we define Ï• and θ angles for molecular dipoles. These are depicted in Fig. 2a and b, respectively. These figures also indicate the coordinate axes used to define these angles. The \(\phi=\frac{\pi }{2}-\arccos \left(\hat{{{{\boldsymbol{\mu }}}}}\cdot \hat{{{{\bf{z}}}}}\right)\) angle captures the alignment between the molecular dipole vector μ and the z-axis, with hats above vectors indicating normalized unit vectors. A water molecule with ϕ = ± π/2 will have its dipole vector aligned with the  ± z-axis. We also wrap around the Ï• angles computed such that Ï• always lies in the range \([-\frac{\pi }{2},+ \frac{\pi }{2}]\). The \(\theta=\arctan \left({{{\bf{c}}}}\cdot \hat{{{{\bf{z}}}}}/{{{\bf{c}}}}\cdot \hat{{{{\bf{x}}}}}\right)\) angle captures a molecule’s rotation about its center of mass around the y axis with \({{{\bf{c}}}}={{{{\bf{r}}}}}_{{{{{\rm{OH}}}}}_{1}}\times {{{{\bf{r}}}}}_{{{{{\rm{OH}}}}}_{2}}\) and \({{{{\bf{r}}}}}_{{{{{\rm{OH}}}}}_{1}}\) and \({{{{\bf{r}}}}}_{{{{{\rm{OH}}}}}_{2}}\) being the OH bond vectors of an individual water molecule. θ measures the rotation of c around the y-axis. Since the choice of H1 and H2 is arbitrary, we wrap around the θ angles so that they always lie in the range [0, π].Fig. 2: Depiction of the orientations of water molecules.a The Ï• angle captures in-plane and out-of-plane fluctuations of water molecules’ dipole moments. b The θ angle captures in-plane rotations of water molecules. The mathematical definitions of the Ï• and θ angles are provided in the text. c The left and right images depict the molecular orientations at the probability maxima in Fig. 3a, while the middle image depicts the molecular orientation at the maximum for the number of putative hydrogen bonds in Fig. 3b.To illustrate the flat-rhombic phases’ apparent hesitancy to form additional hydrogen bonds, we computed the probabilities of molecular orientations defined by these θ and Ï• angles. These two-dimensional log probabilities are juxtaposed with the mean number of hydrogen bonds for each two-dimensional bin determined by the θ and Ï• angles, as depicted in Fig. 3. As shown in Supplementary Note V, the nanoconfined hexagonal and pentagonal phases exhibit the expected trend: molecular orientations with more hydrogen bonds should be more stable. In contrast, the most probable flat-rhombic molecular orientations (depicted in the left and right plots of Fig. 2c) correspond to only two hydrogen bonds. This is true even though orientations with three hydrogen bonds (like the one in the middle plot of Fig. 2c) are possible. We begin to observe these orientations with three hydrogen bonds via the entropic exploration of θ and Ï• angles at high temperatures and in the presence of quantum nuclear fluctuations (see Supplementary Note VI). However, these orientations still do not emerge as local probability maxima, suggesting they are not stable or long-lived states. This behavior reinforces the observation that the two-dimensional flat-rhombic ice phase prefers to remain hydrogen bonded only along one dimension, even at finite temperatures.Fig. 3: Molecular orientation and the number of putative hydrogen bonds.a The two-dimensional log probabilities along the Ï• and θ parameters for the flat-rhombic phase at 2 GPa and three different temperatures. b The average number of putative hydrogen bonds associated with each (θ, ϕ) molecular orientation for the flat-rhombic phase under the same conditions. Relevant molecular orientations are depicted in Fig. 2c. We employ the geometric definition of the hydrogen bond from ref. 68. Source data are provided as a Source Data file.The peculiar hydrogen bonding behavior of the flat-rhombic phase hence emphasizes the importance of incorporating dynamical information into the hydrogen bonding definition, as is needed in, e.g., the description of supercritical water38. As shown in Fig. 4a, if we employ a purely geometric definition of a hydrogen bond, counting the number of instantaneous putative hydrogen bonds that exist on average, we observe an apparent increase in the number of hydrogen bonds with temperature. In the same figure, we also observe that nuclear quantum effects apparently exhibit an additional 0.2 putative hydrogen bonds, even at temperatures as low as 20 K. This is the result of the zero-point fluctuations lowering the energy barrier for sampling these additional hydrogen-bonded configurations, as we elaborate upon in Supplementary Note VI. In these configurations, hydrogen bonds form between different chains in the flat-rhombic crystal structure, as they do classically at high temperatures. The increased disorder due to quantum nuclear effects is not surprising, as the effect of quantum nuclear motion on properties of water has often been mapped to a temperature increase39. To characterize the hydrogen bonding motifs in this structure, we introduce the notation NDMA to indicate a water molecule that donates N hydrogen bonds and accepts M hydrogen bonds. Using this notation, we see that this increase in the number of putative hydrogen bonds can be attributed to the increasing prevalence of instantaneous 2D1A motifs, i.e. instantaneous hydrogen bonds between chains.Fig. 4: Temperature dependence of hydrogen bonding.a The average number of putative/geometric hydrogen bonds across the range of temperatures in which the flat-rhombic phase is stable using classical and quantum simulations. The error bars show the standard error of the mean, as computed by block averaging over 10 blocks. b The proportion of water molecules of each hydrogen bonding motif (defined in the text). The 1D2A proportions are perfectly aligned with the 2D1A proportions. c The 2D1A motif’s lifetime remains extremely short across the range of temperatures in which it is observed. The error bars show the standard error of the mean, computed across all instances of hydrogen bonding. Dashed lines serve as a visual guide for the eye. Source data are provided as a Source Data file.This apparent increase emerges from the shortcoming of the standard geometric definition of a hydrogen bond, which cannot distinguish fluctuating molecular orientations from long-lived orientations that correspond to true hydrogen bonds. As shown in Fig. 4c, the lifetime of the 2D1A orientations is extremely short – on the order of femtoseconds. As suggested by Schienbein and Marx38, these fleeting hydrogen bonds are best referred to as either putative or geometric, as they do not last long enough to be considered hydrogen bonds from a traditional perspective. In Supplementary Note II, we compute the number of hydrogen bonds that survive for typical intermolecular oscillations. Within this definition by Schienbein and Marx38, a hydrogen bond is not counted if its lifetime is shorter than the period of typical intermolecular oscillations. The resulting dynamical hydrogen bond count indeed exhibits the expected disorder-induced decrease as we raise temperature.Dielectric nature of the flat-rhombic phaseA consequence of the quasi-one-dimensional hydrogen-bonded structure is that it may possess ferroelectric behavior due to its net dipole moment along the direction of the hydrogen-bonded chains. Ferroelectricity in ice has been conjectured previously in force field simulations40, first-principles calculations of nanoconfined water41, and experiments on confined28 and supported42 films. To investigate the dielectric behavior of flat-rhombic ice, we model the temperature dependence of its dielectric response in Fig. 5. We analyzed the in-plane ε∥ and out-of-plane ε⊥ dielectric constants based on the variance of the system’s polarization. Since our MLP only predicts the potential energy surfaces, a first principles investigation of the dielectric response would be computationally demanding due to numerous single point calculations of the electronic polarization43. Therefore, for a semi-quantitative understanding, we use a simple linear polarization model based on TIP4P charges44.Fig. 5: Temperature dependence of the dielectric response.The classical dielectric constant of the flat-rhombic phase across the range of temperatures considered in this work. The singularity at 380 K corresponds to the transition to the hexatic phase of nanoconfined water53. Source data are provided as a Source Data file.We select the TIP4P water model due to its low computational cost and a semi-quantitative description of the dielectric response of bulk44 and confined45 water. The TIP4P model underpredicts the polarization of aqueous systems as it doesn’t incorporate the electronic polarization of the water molecules. For instance, TIP4P predicts a molecular dipole moment of 2.348 D which is 12.8% lower compared to first principles estimates46, and effectively leads to a  ~ 25% lower dielectric response. To incorporate the lack of the electronic polarization in the TIP4P model, we rescale the calculated dielectric response by a factor of 1.25, as has been done previously45. We report the classical dielectric response, as quantum nuclear effects are expected to only make a small quantitative difference44.We observe low and near-constant in-plane and out-of-plane dielectric constants for the flat-rhombic phase, in agreement with the experiments on nanoconfined water from ref. 47. The low dielectric response of this phase is only an indirect outcome of vdW interactions, as they dictate the thermodynamic stability of the flat rhombic phase. The dielectric nature of flat-rhombic ice remains the same from 0 K up to its phase transition into the paraelectric hexatic (disordered) water phase16. The temperature independence of the flat-rhombic phase’s dielectric behaviors reflects the resilience of the quasi-one-dimensional structure to thermal and quantum nuclear fluctuations. The singularity observed in Fig. 5 indicates the temperature at which the phase transition to the hexatic phase occurs. The hexatic phase exhibits an increased in-plane dielectric constant but a similarly small out-of-plane dielectric constant.Long-ranged spatial ordering and coherent proton dynamicsThus far, the results in this work have considered large simulation boxes with 576 molecules spanning tens of nanometers. However, in simulation boxes containing 144 molecules in a relatively small but experimentally accessible system size on the order of nanometers48,49, we identify unusual concerted dynamics at intermediate temperatures involving the breaking and forming of hydrogen bonds between equivalent structures. In this concerted motion, we observe an exchange between the individual molecular ϕ = ± π/4 states as a collective motion throughout the entire system, in which all the water molecules simultaneously swap the signs of their Ï• angles (see Supplementary Movie 1).To describe this concerted motion, we define a σ order parameter (see Supplementary Note III) to distinguish the two states shown in Fig. 6a. The structural change to the flat-rhombic phase between σ values can be seen in Fig. 6a. When the sign of σ changes, the directions of the dipoles in each row switch in a concerted fashion. During this motion, we also observe the complete rearrangement of the hydrogen bonding network in the lattice. Figure 6a shows labeled water molecules in each σ state, illustrating that the hydrogen bonds in these two states are between different pairs of water molecules. At low temperatures, the system remains frozen in a single σ value. The resulting σ free energy profile contains two minima separated by a large free energy barrier, shown in Fig. 6c. At intermediate temperatures, we begin to see exchanges between σ signs on the order of  ~10 ps, as shown in Fig. 6b. Compared to the low-temperature setting, the intermediate temperature free energy profile exhibits a clear finite free energy barrier between these states. Eventually, at high enough temperatures, the free energy profile along σ exhibits an extremely small free energy barrier, allowing the system to explore molecular configurations freely. While the highest temperatures explored here are greater than the melting temperature of ice, we believe that our trajectories correspond to a metastable state associated with the solid phase. Since the free energy profiles at the intermediate and high-temperature conditions are qualitatively the same, the near-free exploration of molecular configurations in the high-temperature simulations is a result of thermal activation. To confirm this hypothesis, we compute rotational autocorrelation functions at these different temperatures in Supplementary Note IV.Fig. 6: Hydrogen bonded switching behavior.a The two symmetry-related structures of the flat-rhombic phase, with alternating rows of aligned dipoles, discerned by the different signs of σ. A zero value of σ indicates that the molecular dipoles are not aligned row-wise. The side views in each panel show the system as viewed from the right side, where the switch in the row-dependent dipole direction between the two states is clear. b For smaller unit cells with 144 water molecules, we observe coherent switching between these two states at intermediate temperatures, as shown in the σ trajectories for three different temperatures. c The free energy profiles for this σ parameter show that raising the temperature lowers the free energy barrier between these states. The smooth, low-barrier free energy profile at high temperatures indicates that the molecular dipoles lose spatial coherence. Source data are provided as a Source Data file.We also note that quantum nuclear motion makes the system more disordered due to zero-point fluctuations lowering the barrier for the system’s hydrogen bond network to switch between the free energy minima, as shown in Supplementary Fig. 8. Furthermore, since nuclear quantum effects lead to increased proton disorder, the σ parameter gets pushed closer to 0 on average. This causes the minima in the free energy profiles along σ to move closer to 0 in our PIMD simulations as compared to the free energy profiles from our classical simulations.As stated earlier, this concerted motion is not observed in the larger simulation cell size that contains 576 molecules for the timescales we consider: it is only seen in simulations with 144 water molecules. This suggests that the free energy barrier for this concerted motion is size-extensive. As a result, we would expect a larger system size to either exhibit an ordered phase or a phase in which smaller domains that exhibit this behavior coexist. However, crystallites of nanoconfined water molecules containing around 100 water molecules are experimentally accessible, so such concerted motion could be observed in experiments as small ice crystallites forming between graphene sheets48 or high-pressure graphene nanobubbles formed by irradiation49. Previous work has suggested that the simulation setup we employ is still reasonable for studying such encapsulated systems, even without explicit confining atoms48,50. A careful treatment of such encapsulated systems would require a detailed analysis of the water-carbon interactions at the edge of the confined water pockets. In this work, we ignore these edge effects and instead focus on characterizing the behavior of nanoconfined water under idealized atomistically flat nanoconfinement. A careful consideration of the edge effect would be an interesting and relevant topic for future work.
