AB-stacked square-like bilayer ice in graphene nanocapillaries

Water, when constrained between two graphene sheets and under ultrahigh pressure, can manifest dramatic differences from its bulk counterparts such as the van der Waals pressure induced waterto-ice transformation, known as the metastability limit of twodimensional (2D) liquid. Here, we present result of a new crystalline structure of bilayer ice with the AB-stacking order, observed from molecular dynamics simulations of constrained water. This AB-stacked bilayer ice (BL-ABI) is transformed from the puckered monolayer square-like ice (pMSI) under higher lateral pressure in the graphene nanocapillary at ambient temperature. BL-ABI is a proton-ordered ice with square-like pattern. The transition from pMSI to BL-ABI is through crystal-to-amorphous-to-crystal pathway with notable hysteresisloop in the potential energy during the compression/decompression process, reflecting the compression/tensile limit of the 2D monolayer/ bilayer ice. In a superheating process, the BL-ABI transforms into the AB-stacked bilayer amorphous ice with the square-like pattern.


Introduction
The unique strength and orientation of the hydrogen bond and molecular polarity of water molecules endow a rich but still incompletely understood phase diagram of water/ice. [1][2][3][4][5][6] The interplay between the hydrogen-bonding interaction and nanoscale confinement results in additional complex phase behavior of water/ice. [4][5][6][7][8] Understanding the behavior of water/ice in nanoscale confinement is an active area of research with many technological implications such as nanotribology, nanofluidics, and synthetic ion channels. Intriguing phase behaviors discovered theoretically thus far include 2D crystalline/amorphous ice phases, metastable phases in the compression limit of confined liquid water, as well as water transport and adsorption in a narrow channel. [9][10][11] For theoretical studies of ice formation at high pressure, different computational methods such as molecular dynamics (MD) simulations with force fields and density functional theory (DFT) computations have been utilized. [12][13][14][15][16][17][18] Liquid-to-solid, solid-to-solid, and solid-to-liquid-to-solid phase transitions are found in the simulation of water/ice polymorphs and polyamorphous, while some of these transitions are the first-order and others are continuum. 19,20 The metastability limit extension of liquid water at high compression/tensile pressure or high degree of superheating/supercooling provides additional features that manifest themselves in the equilibrium states of water/ice. 2,3,[21][22][23] The effects of confinement and associated van der Waals pressure on the mechanical stability of water/ice are yet to be comprehensively studied. The main focus of this study is the phase transition between puckered monolayer ice and AB-stacked bilayer ice under ultrahigh lateral pressure in graphene nanocapillaries.
Monolayer and bilayer crystalline/amorphous ice are widely investigated in previous studies. Monolayer ice has been seen in the MD simulations with the five-site TIP5P water model. [24][25][26] Puckered monolayer ice can be formed when the width of the silt pore is slightly larger than that giving the flat monolayer ice. 12,27,28 The degree of puckering of the monolayer ice is closely related to the lateral pressure. Numerous bilayer ice polymorphous and polyamorphous structures have been shown in hydrophobic slit pores from computer simulations. [29][30][31][32][33][34][35][36] Most bilayer ices reported from previous studies are the AA-stacked structures at low pressure, while the AB-stacked bilayer ice is little noticed. This is because most AA-stacked bilayer ices are crystalline phases which satisfy the ice rule, whereas most AB-stacked bilayer ices are amorphous. A recent study indicates that monolayer ice can transform into AB-stacked bilayer amorphous ice at certain nanocapillary width and lateral pressure. 12 Similarly, the phase transformation from triangular AA-stacked bilayer ice to trilayer square ice with the ABA-stacking order occurs in relatively narrow graphene channels as well, under ultrahigh lateral pressure. 18 Moreover, some bulk ice phases can be formed at very high pressure, e.g., on the order of TPa (Mbar) magnitude. [37][38][39][40] Since the stacking order of water molecules in graphene nanocapillaries is determined by the pore width and lateral pressure, we explored the possibility of a new phase transition from monolayer ice to AB-stacked bilayer ice, that may occur under higher lateral pressure. Indeed, a novel AB-stacked bilayer crystalline ice is seen in our MD simulations.

Results and discussion
The simulation system is similar with that used in two recent studies, 7,12 in which two water reservoirs containing 1000 molecules, each connected by a relatively long capillary formed by two parallel graphene sheets. The lateral pressure forces water molecules from reservoirs into the graphene nanocapillary. The distance h between the center planes of each graphene wall is varied from 6.5 to 7.5 Å, allowing one or two monolayers of water accommodated within the graphene channel. Periodic boundary conditions are imposed in all three directions. The MD simulations are performed in the isothermal-isobaric (NP zz T) ensemble, in which the temperature (T) and lateral pressure (P zz ) are controlled by the Nosé-Hoover thermostat and barostat, respectively. The total number of water molecules N in the nanocapillary can change, depending on P zz , T and h. The total potential of interaction is taken as the sum of interactions among TIP4P/2005 water molecules, and the external Lennard-Jones potential of the interaction between water molecules and the graphene wall. 12,41 The potential energy per molecule, the oxygen-oxygen radial distribution function [g O-O (r)], the oxygen-hydrogen radial distribution function [g O-H (r)], and the mean-squared displacement (MSD) are computed in the simulations to characterize the phase transition and phase behavior during the compression and superheating processes.
In the monolayer square-like ice, the flat hydrogen-bonding network becomes puckered when the pore width is slightly enlarged. 27 Under ultrahigh lateral pressure, the puckered monolayer square-like ice can transform into bilayer ice, 12 akin to the phase transition from triangular bilayer-AA stacking ice to square trilayer ABA-stacking ice. 18 In Fig. 1, the new bilayer crystalline ice with the AB-stacking pattern is demonstrated. This new ice is transformed from the puckered monolayer square-like ice (pMSI) under lateral pressure P zz > 4.5 GPa and at ambient temperature in our MD simulations. The red and blue colored water molecules in Fig. 1(a) denote the molecules in two different layers, respectively. From the top view, water molecules in the second layer (red color) are located at the center sites among the lattice sites in the first layer (blue color). Water molecules in BL-ABI form a regular rhombic hydrogen-bonding network [longer green lines in Fig. 1(a) represent hydrogen bonds]. Here, each oxygen atom (larger sphere) is connected with four hydrogen bonds (two shorter and two longer green lines in Fig. 1(a)). Thus, the ice rule is satisfied over the entire network except a few point defects. 28 This rhombic hydrogenbonding network is a bilayer network rather than a puckered monolayer network. The side view of BL-ABI in Fig. 1(b) demonstrates that the hydrogen bonds in BL-ABI can be classified into two types, that is, in-plane and interlayer, respectively. Note that, the interlayer hydrogen bonds display an ordered zigzag pattern from the side view [ Fig. 1(b)], indicating that some hydrogen bonds are nearly parallel to the graphene surface while others are tilted with respect to the graphene surface. Interestingly, the hydrogens [white balls in Fig. 1(b)] in BL-ABI have a regular and inerratic arrangement. In the two layers, every four water molecules form a square-like (slightly rhomboidal) arrangement with in-plane hydrogen bonds. Each layer of the BL-ABI [ Fig. 1(c)] resembles the monolayer ice reported by Zangi and Mark. 24 The formation of the AB-stacked ice is mainly due to the relatively small distance h between the center planes of the graphene, resulting in higher density with respect to the AA-stacked structure. The oxygen-oxygen radial distribution function in Fig. 1(e) illustrates the first sharp peak in g O-O (r), located at 2.79 AE 0.06 Å, which corresponds to the first neighboring distance between oxygen atoms in BL-ABI. The well-separated peaks and valleys indicate that the BL-ABI is a crystalline phase.
The slightly rhomboidal bilayer ice exhibits a square-like structural feature although the water molecules in each layer exhibit a near-rhombic-network arrangement. Different to the rhombic and triangular configurations of BL-AB and BL-AAI structures reported previously, 12   neighboring water molecules (two red and two blue molecules) in the BL-ABI form a square-like pattern. To gain more insight into the structural characteristics of BL-ABI, a randomly selected section is shown in Fig. 2, including 25 water molecules. The center of this section is relocated as the origin of coordinates in Fig. 2(c). Then the positions of the oxygen atoms are plotted on this coordinate plane. The red and blue circles correspond to the oxygen atoms in Fig. 2(a). The red and blue dotted lines indicate that oxygen atoms in every layer form a square-like (slightly rhomboidal) configuration. The green boxes reflect that four neighboring oxygen atoms in the BL-ABI form a square-like pattern. The green dotted lines denote that four neighboring oxygen atoms in two layers (two red and two blue symbols) can also form a rhombic unit, corresponding to the rhombic hydrogen-bonding net unit in Fig. 1(a), and the side view of this rhombic unit displays the in-plane and zigzag interlayer hydrogen-bonds formed in BL-ABI. In Fig. 2(c), the length of the red and blue dotted lines between two oxygen atoms is B2.74 Å, the length of the green solid-line between the red and blue symbols on the coordinate plane is B2.00 AE 0.08 Å (in 2D view plane). In the green-dotted rhombic unit (considering the interlayer spacing), the lengths of two sides are B2.735 Å and B2.84 Å, respectively. The two characteristic lengths (B2.74 Å and 2.84 Å) correspond to the location of the first high peak Fig. 1(e). The interlayer spacing of the two layers is B1.96 Å.
Two independent MD simulations of compression and decompression processes (at h = 7.0 Å, T = 300 K, 1.0 GPa r P zz r 6.0 GPa) are illustrated in Fig. 3. In the compression process (red line), the lateral pressure is increased in a step of 0.1 GPa per nanosecond, and the potential energy first increases gradually and then suddenly jumps by about 1.2 kcal mol À1 at B4.1 GPa. This large energy increase is due to the phase transition from a highly puckered monolayer square-like ice (hpMSI) to the BL-ABI in the graphene nanocapillary. Before the formation of hpMSI, the monolayer liquid water first forms pMSI when the lateral pressure is beyond the compression-limit (P zz > 0.6 GPa), 12 which then forms hpMSI with increasing lateral pressure (P zz o 3.8 GPa). In the reversed decompression process (blue line), the potential energy first decreases gradually and then suddenly drops to a lower value at B2.2 GPa, overlapping with the path of the compression process. Here, the large hysteresis loop manifests that the transformation between pMSI and BL-ABI is a strong first-order phase transition. 12 In particular, in the range of 2.3-4.2 GPa of the lateral pressure, the compression/ decompression alters the ice layer-number, giving the compression/ tension limit of the 2D monolayer/bilayer ice. Moreover, the structural transformation from pMSI to BL-ABI exhibits typical structural evolution of the distribution of hydrogen atoms, that is, in pMSI the distribution of hydrogen atoms is along the rows and columns, while in BL-ABI it shows an ordered arrangement.
To achieve a better understanding of the phase transition from pMSI to BL-ABI, the simulation temperature is controlled by a Berendsen thermostat (T = 300 K) for 25 nanoseconds in the NP zz T ensemble with a given lateral pressure. Six typical snapshots under different lateral pressures are demonstrated in Fig. 4, corresponding to the lateral pressure P zz = 1.0, 2.0, 3.0, 4.0, 5.0, and 6.0 GPa, respectively. Overall, for h = 7.0 Å, three different ice phases can be formed in the graphene confinement. As discussed above, the liquid is transformed into pMSI ( Fig. 4(a)) when lateral pressure is beyond the water/ice compression-limit. With increasing lateral pressure, pMSI turns into hpMSI ( Fig. 4(b and c)) gradually. Finally, hpMSI transforms into BL-ABI (Fig. 4(e and f)) under ultrahigh lateral pressure, during which, an amorphous phase (Fig. 4(d)) fors when the potential energy exhibits a sharp jump. When the lateral pressure is B1.0 GPa, pMSI shows a typical hydrogen-bonding square-network (green lines in Fig. 4(a)). In Fig. 4(b and c), hpMSI displays more clearly rectangular and rhombic characteristics due to the high lateral pressure. Another two ice phases are the AB-stacked amorphous and crystalline ices with the same square pattern. The main difference among the two phases is the proton arrangement. In Fig. 4(d), the AB-stacked amorphous ice is a proton-disordered phase. The AB-stacked crystalline ice in   is similar with that in the bilayer structures predicted by hard or soft sphere models. 47,48 From the structural evolution in Fig. 4, the phase transition between pMSI and BL-ABI goes through the crystal-to-amorphousto-crystal transformation. Note that this monolayer to bilayer phase transition differs from the bilayer to trilayer phase translation 12,18 in that the transformation from bilayer triangular AA-stacking ice to ABA-stacked trilayer square ice is through the solid-to-liquid-to-solid process. 18 To characterize this novel crystal-to-amorphous-to-crystal phase transition, in Fig. 5, the oxygen density profiles, computed MSDs, g O-O (r)s and g O-H (r)s at different lateral pressure are displayed, respectively. The red and blue lines denote hpMSI at P zz = 2.0 and 3.0 GPa, respectively. The magenta lines represent the AB-stacked amorphous ice with square pattern at P zz = 4.0 GPa. The olive and green lines denote the BL-ABI at P zz = 5.0 and 6.0 GPa. The peak value of oxygen density profile in Fig. 5(a) for hpMSI is only half that for BL-ABI, indicating that the AB-stacked bilayer ice is a very-high-density phase. The MSD results in Fig. 5(b) show that the diffusion coefficient of hpMSI is B3.0 Â 10 À9 cm 2 s À1 , while that of BL-ABI is even smaller. The radial distribution function can also demonstrate the structural difference among the three ice phases. The locations of the first sharp-peak in g O-O (r)s and g O-H (r)s for the three ices are at 2.79 AE 0.06 Å and 1.75 AE 0.03 Å, respectively. While the locations of the second sharp-peak are quite different: the second sharp-peak in g O-O (r) for hpMSI is located at 4.32 AE 0.03 Å, while that for BL-ABI is located at 4.88 AE 0.02 Å (Fig. 5(c)). The magenta g O-O (r) indicates that the AB-stacked bilayer amorphous ice with square pattern is an oxygen-ordered phase. In Fig. 5(d), the peak values in g O-H (r)s of hpMSI are greater  than those of BL-ABI. The second sharp-peaks in g O-H (r)s for hpMSI and BL-ABI are located at 3.08 AE 0.02 Å and 2.84 AE 0.02 Å, respectively. The peak values in magenta g O-H (r) are slightly smaller than those in olive and green g O-H (r)s, which manifests that the AB-stacked bilayer amorphous ice with square pattern is a proton-disordered phase. Thus, the phase transition from pMSI to BL-ABI is a crystal-to-amorphous-to-crystal transformation.
Ice formation and ice structures in nanoscale confinement depend strongly on lateral pressure P zz and pore width h. 12 In addition to the pMSI to BL-ABI transition at h = 7.0 Å, we have also performed extensive MD simulations over a range of h (6.5-7.5 Å) and P zz (0.0-6.0 GPa) to explore the compressionlimit (phase) diagram in the h-P zz plane. The marked symbols in Fig. 6 denote the compression limit of the corresponding phases. The compression-limit curves (solid lines) are approximately drawn by connecting the compression-limit points in Fig. 6. The dashed lines are schematic phase boundary lines. As illustrated in Fig. 6, the 2D liquid is transformed to pMSI when the lateral pressure is beyond the compression limit (black line with square symbols). At much higher lateral pressures, the pMSI is transformed to BL-ABI in a relatively narrow channel (h o 7.2 Å), or to BL-AB in a relatively wide channel (h > 7.3 Å). The difference between the formation of BL-ABI and BL-AB is mainly due to the nanocapillary width h. When the channel width is slightly wider (h > 7.3 Å), water molecules are easily transformed into the BL-AB phase. We overlap the compression-limit lines of the amorphous (blue line) and BL-AB (olive line) at the point (7.2 Å, 3.2 GPa). But the blue line is due to a vertical path (hpMSI-toamorphous), and the olive line is due to a horizontal path (increasing the pore width). The two paths are different. The region between the blue and red lines denotes the amorphous phase during the crystal-to-amorphous-to-crystal phase transition. Oxygen atoms in BL-ABI and the amorphous phase exhibit square configurations, while those in BL-AB show rhombic arrangements. 12 Hydrogen atoms in pMSI and hpMSI are along the rows and columns, whereas those in the BL-ABI exhibit regular arrays.
Lastly, the superheating of BL-ABI in a graphene nanocapillary (h = 7.0 Å, P zz = 5.0 GPa) is investigated. The simulation starts with the ice structure at a given initial lateral pressure and temperature (300 K). The system is then heated to the corresponding temperature in NP zz T ensemble for 20 nanoseconds. The potential energy, radial distribution function, and density profiles are computed to reveal the melting properties of BL-ABI. In the superheating process, the oxygen atoms in the channel under high temperature still maintain the AB-stacking order, while the hydrogen atoms are disordered.
The lower-right inset in Fig. 7(a) shows that the superheated structure is similar to the AB-stacked amorphous bilayer ice shown in Fig. 4(d), where the oxygen arrangement is square-like instead of the rhombic configuration in BL-AB. The potential energy per molecule curve in Fig. 7(a) is not fully linear with increasing system temperature. If the temperature goes beyond 500 K, the potential energy increases gradually. The values of the first sharp-peak in g O-H (r)s ( Fig. 7(b)) decrease gradually with increasing temperature, indicating that the proton arrangement is not as regular as that in BL-ABI. When the temperature T Z 350 K, the shape of g O-H (r) is similar to the magenta line in Fig. 5(d). The density profiles of oxygen and hydrogen atoms in the channel are illustrated in Fig. 7(c) and (d), respectively. The peak values of oxygen density profiles decrease gradually with increasing temperature. The hydrogen density profiles in Fig. 7(d) exhibit additional results. The superheated phase is the proton-disordered structure.
In conclusion, a new bilayer crystalline ice is found in MD simulations for water constrained in graphene nanocapillaries under high lateral pressure. A compression limit (phase) diagram for the pMSI, BL-ABI, and BL-AB phases is essentially for the nanocapillary width versus lateral pressure (h-P zz ) plane (6.5 Å r h r 7.5 Å, P zz r 6.0 GPa) at T = 300 K. BL-ABI is an AB-stacked crystalline bilayer ice with regular and inerratic proton arrangement. It is transformed from puckered monolayer square-like ice under higher lateral pressure and ambient temperature. Water molecules in BL-ABI demonstrate obvious square pattern from the top view. The phase transition from pMSI to BL-ABI is through the crystal-to-amorphous-to-crystal transformation, as well as the first-order phase transition with sudden change in potential energy and with hysteresis-loop in compression/decompression process. The sudden jump and drop in potential energy (per molecule) curves during compression/decompression process mark the structural evolution among monolayer/bilayer ice, simultaneously, giving the compression/tension limit of the 2D monolayer/bilayer ice. The intermediate amorphous phase between hpMSI and BL-ABI is an AB-stacked bilayer protondisordered ice while the oxygen configuration is similar to that in BL-ABI. Under superheating, BL-ABI transforms into AB-stacked bilayer amorphous ice with a square-like pattern.

MD simulation
We perform classical MD simulations of few-layer water/ice constrained in graphene nanocapillaries, using the LAMMPS program. 42 Fig. 6 Compression limit (phase) diagram of water confined in graphene nanocapillaries at ranges of 6.5 Å r h r 7.5 Å and 0.0 GPa r P zz r 6.0 GPa. The marked symbols denote the compression-limit of corresponding ice. The black and olive square symbols are the same as those in ref. 12. The simulation model in this current study is the same as that in ref. 12 and 18 (an illustration of the model system is displayed in Fig. S1, ESI †). The length and width of the graphene nanocapillary are fixed at 42.60 Å and 36.89 Å, respectively, in all simulations. Periodic boundary conditions are imposed in all three directions. A time step of 1.0 fs is used for the velocity-Verlet integrator. The total potential of interaction is taken as the sum of interactions P f water (r i , r j ) among water molecules, and the external potential P f wall (r i , r j ) of the interaction between water molecules and the graphene wall, where r i stands for the coordinate of water molecule i, f water is the extended four point charge (TIP4P/ 2005) potential of water, f wall is the Lennard-Jones potential for water-graphene interaction. 12,41 The pairwise interactions between any two water molecules are described by the TIP4P/2005 model, 43,44 including the long-ranged Coulomb potential and the short-ranged Lennard-Jones 12-6 potential between the interaction sites. The long-range Coulombic interaction is calculated using a particle-particle particle-mesh (PPPM) algorithm with an accuracy of 10 À4 . TIP4P/2005 water model can be considered as a reliable model of water, offering good performance over a wide range of conditions. 45 The carbon-water potentials are based on parameters for carbon, oxygen, and hydrogen obtained through the Lorentz-Berthelot mixing rule. 46 The force-field parameters used in the simulations are the same as those in ref. 12.