2D ice from first principles: structures and phase transitions

Despite relevance to disparate areas such as cloud microphysics and tribology, major gaps in the understanding of the structures and phase transitions of low-dimensional water ice remain. Here we report a first principles study of confined 2D ice as a function of pressure. We find that at ambient pressure hexagonal and pentagonal monolayer structures are the two lowest enthalpy phases identified. Upon mild compression the pentagonal structure becomes the most stable and persists up to \textit{ca.} 2 GPa at which point square and rhombic phases are stable. The square phase agrees with recent experimental observations of square ice confined within graphene sheets. We also find a double layer AA stacked square ice phase, which clarifies the difference between experimental observations and earlier force field simulations. This work provides a fresh perspective on 2D confined ice, highlighting the sensitivity of the structures observed to both the confining pressure and width.

1 arXiv:1508.03743v1 [cond-mat.mtrl-sci] 15 Aug 2015 Confined and interfacial water-ice is ubiquitous in nature, playing an important role in a wide range of areas such as rock fracture, friction, and nanofluidics [1][2][3][4]. As a result of a delicate balance of forces (hydrogen bonding, van der Waals, and interaction with the confining material or substrate) confined and interfacial water forms a rich variety of structures [5][6][7][8][9][10]. Almost every specific system examined has revealed a different structure such as a 2D overlayer built from heptagons and pentagons on a platinum surface or the square ice observed within layers of graphene [8,10]. This shows that in contrast to bulk ice the phase behavior of 2D ice is much less well understood.
From an experimental perspective a full exploration of the phase diagram of 2D ice has not been achieved yet. However recent experiments revealed the exciting possibility of exploring 2D ice structures at specific conditions by trapping water within layered materials [2,10,11].
For example by confining water between layers of graphene it is possible to create so called nanocapillaries in which water experiences a pressure in the GPa regime due to the van der Waals forces pulling the sheets together [10]. Using transmission electron microscopy (TEM) square ice structures from a single up to a few layers were observed in such graphene nanocapillaries [10], and in the multilayer regime square ice layers located directly on top of one another in an AA stacking arrangement. Force field simulations performed as part of the same study of confined water in graphene layers reproduced the square monolayer ice but failed to explain the AA stacking [10]. Indeed prior to this recent study there was already a long tradition of computer simulation studies of nanoconfined water, mostly involving classical force field approaches [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. Such work has been incredibly valuable and provided considerable insight into the structures and phase transitions of monolayer and multilayer ice.
The 2D ice structures predicted include hexagonal, pentagonal, quasicrystalline, hexatic, and orthogonal phases as well as various amorphous structures. Although the relative stability of the structures depends on the particular force field used, for monolayer ice an orthogonal phase has been widely predicted to be stable across a broad pressure regime [13,14,19,23].
Given the sensitivity of force field studies of confined ice to the potential used, the aspiration to understand the observation of monolayer square ice, and the unexplained AA stacking order of square ice layers, a systematic study with an electronic structure method such as density functional theory (DFT) is highly desirable. DFT studies of ice do not come without their own sensitivity to the exchange-correlation functional used, however, they have proved to be very useful in predicting and understanding structures of adsorbed water and bulk ice.
Here we report a systematic study of 2D phases of water-ice from first principles based on an unbiased exploration of monolayer ice structures. The much greater computational cost of DFT compared to force fields means that we cannot map the entire phase diagram of monolayer confined water. Instead we focus on the phase transitions as a function of lateral pressure and confinement width at zero Kelvin. The stable structures identified at different pressures include an hexagonal structure, a pentagonal Cairo tiling (CT) structure, a flat square structure, and a buckled rhombic structure. The observation of a flat square structure is consistent with the recent experimental observation. However, the sequence of low energy phases identified differs significantly from that predicted in force field based studies [13,14,19,21,23] and a recent DFT report [26]. Interestingly the sequence of structures observed depends sensitively on the confinement width used, suggesting that it should be possible to tune the monolayer ice structures produced in experiments by e.g. varying the confining material. We also propose a double layer square ice phase with interlayer hydrogen bonds, which further explains the AA stacking order of multi-layer ice recently observed [10].
In order to explore ice structures in an unbiased manner, we used the ab initio random structure search (AIRSS) technique [29,30], an approach which has previously predicted new ice, 2D, and interfacial structures [31,32]. Structures from ambient up to a lateral pressure of 10 GPa were considered. The 2D confinement was introduced via a Morse potential fit to quantum Monte Carlo results for the binding of a water monomer to graphene [33].
By tuning the confinement width, we are not only able to study the general properties of water under flat and smooth confinement but also to compare with the recent experiments of ice in graphene confinement [10]. Our electronic structure calculations were carried out using the Vienna Ab-initio Simulation Package (VASP) [34] with the DFT+vdW approach [35] in conjunction with the Perdew-Burke-Ernzerhof exchange-correlation functional [36].
Tests with other exchange-correlation functionals show that whilst the transition pressures between the various phases depend to some extent on the choice of exchange-correlation functional, the overall conclusions do not change. See the supporting information (SI) for these results as well as further computational details [37].
From a preliminary set of calculations we established that the optimal separation between graphene sheets in which a monolayer of water is sandwiched is somewhere between 6.0 and 6.5Å (Fig. S2) [37]. With this in mind we first report results for water within smooth confining potentials that are either 6.0 or 6.5Å wide. With such confinement two phases have been identified at ambient pressure, which have exceedingly similar enthalpies. These are an hexagonal monolayer structure, resembling an hexagonal bilayer, and a Cairo tiled structure built exclusively from water pentagons. With the particular exchange-correlation functional used, the less dense hexagonal structure is marginally more stable than the pentagonal structure by 5 meV/H 2 O. This difference drops to just 2 meV/H 2 O when harmonic zero point energy effects are taken into account (Table S2, Fig. S9). Tests with other exchangecorrelation functionals generally concur that the energy difference between the two phases is tiny, with a vanishingly small preference for the hexagonal structure (Table S2) [37].
Anharmonic zero point energy differences between the two phases and finite temperature effects could also easily exceed the energy difference [38], suggesting that both the hexagonal and the CT phases could be observed at ambient pressures. and 6.5Å (b) confinement. ∆H is the relative enthalpy with respect to the CT phase. Enthalpy H = E water +E confinement +P×A×w, where E water is the total energy per water molecule, E confinement is the energy (per water molecule) in the confinement potential, P is the lateral pressure, A is the lateral area per water molecule, w is the width of confinement.
The hexagonal structure has p6mm wall-paper group symmetry if only oxygen atoms are considered and is built exclusively of six-membered rings (Fig. 1a). Water molecules are three-fold coordinated, with half of them having one OH bond directed out of the monolayer (a so-called "dangling OH"). The average O-O separation in this monolayer hexagonal phase is about 2.72Å which is similar to the bulk O-O separation in ice I [39]. However the confined hexagonal structure identified here is quite flat with the vertical separation between oxygen atoms < 0.3Å, much smaller than the 0.9Å buckling within an hexagonal layer in bulk ice I [40]. Given that bulk ice I is built from hexagonal layers and double layer hexagonal structures have been observed frequently in force field simulations [16][17][18], it is not surprising that a low enthalpy hexagonal structure should be identified. However monolayer hexagonal ice has yet to be observed experimentally and in force field simulations it has only been found when an hexagonally patterned substrate has been used as a template [21].
The pentagonal structure identified has a wall-paper group symmetry of p4gm (Fig. 1b).

5
The unit cell has twelve water molecules, four with dangling OH bonds. One third of the water molecules are four-fold coordinated and the rest are three fold coordinated. The higher average coordination and smaller ring size of the CT phase renders the density of this phase higher than that of the hexagonal phase (Fig. S4, S5) [37]. Therefore, upon increasing the pressure the CT phase clearly becomes more stable than the hexagonal phase. After searching for the lowest enthalpy structures at finite pressure we find that the CT structure is clearly the most stable in a broad range of pressures all the way up to ca. 2 GPa (Fig.   2). The stability regime of the monolayer CT phase and the small energy difference at ambient pressure suggest that five-membered rings are more important in 2D ice than in 3D ice where they only appear in the 0.2-0.7 GPa range (in ices III, V, IX, XIII) [41]. There is some precedent for the pentagonal based structure being proposed here. First and foremost, it is a monolayer version of the double layer confined ice structure identified in Ref. 17 on the basis of simulations with a coarse grained model of water. In addition a 1D pentagonal ice structure has been observed with scanning tunnelling microscopy on a metal surface, although the structure of the 1D pentagonal chain is very different from 2D pentagonal ice [6]. The only purely 2D structure similar to our prediction, that we are aware of, is a recently proposed allotrope of carbon [42].
Previous force field studies have suggested that a rhombic phase is the most stable at ambient pressures [13,14]. Our own force field studies with either TIP4P/2005 [43] or SPC/E [44] indeed find that the rhombic structure has the lowest enthalpy (Fig. S8) [37].
However with DFT it is considerably less stable than the hexagonal and CT structures.
Similarly a recent DFT study concluded that a square structure is more stable than any hexagonal structure [26]. Here we find that our most stable square structure has a higher enthalpy by 43 meV/H 2 O than the hexagonal phase at ambient pressure. In the SI we trace this difference to the different computational setups [37]. As shown in the SI we are confident that the hexagonal and CT monolayer ice structures are indeed more stable than any square ice structure at the low pressure limit. are distributed on two orthogonal antiferroelectric sublattices (Fig. 1c). The HB network of the b-RH phase is similar to f-SQ but it is buckled and has a higher lateral density (Fig. 1d). The relative stability of the f-SQ and b-RH structures in the 2-4 GPa regime depends sensitively on the confinement width ( Fig. 2). At 6.0Å f-SQ is more stable while at 6.5Å b-RH has the lower enthalpy. Beyond 4 GPa the b-RH phase is more stable than any other structure identified. Several other metastable structures belonging to the b-RH family with different hydrogen bond ordering have also been observed. However, since a more delicate discussion of hydrogen ordering is beyond the scope of this study, we only show one of the most stable members of the b-RH family at the pressure and confinement conditions considered (Fig. 1d).
We also identified a second metastable square phase which we dub "b-SQ" because of its buckled "basketweave"-like pattern of HBs (Fig. 1e). Its lattice structure resembles a 2D projection of bulk ice VIII but it is unique in that the two sublayers are hydrogen bonding with each other. The energy of the b-SQ phase is higher than the most stable phases identified, however because of its unique hydrogen bonding arrangement and as it might be possible to observe it in systems where the substrate has a square lattice we feel it is worth reporting.
Beyond the phase behavior of monolayer ice confined within the 6.0Å to 6.5Å regime, we also explored a broader range of confinement widths at different pressures. These additional calculations support the validity of the conclusions reached but also show that there is scope for altering the relative enthalpies of the various phases by tuning the confinement.
The phase diagram for monolayer water with respect to lateral pressure (0-10 GPa) and confinement width (5-8Å) that emerges from these calculations is shown schematically in For larger confinement widths the b-RH phase is generally favored.
The f-SQ phase is found to be stable in the 2-4 GPa regime at 6.0Å confinement. The structure of the f-SQ phase and the approximate pressure at which it appears are consistent with the recent experimental observation of a square ice phase in graphene nanocapillaries [10]. This lends some support to the predictions made here and suggests that it might also be possible to observe the other structures predicted, by for example controlling the density of water inside the nanocapillaries or with another 2D material with lower interlayer adhesion. In addition to TEM which has already been used, techniques such as scanning tunneling microscopy and atomic force microscopy might be able to further substantiate 2D ice lattices and hydrogen ordering in the future [45,46]. In Fig. S6 we also show that different 2D ice monolayer structures have quite different vibrational properties [37]. For example the hexagonal and the CT phase have vibrational modes around 3700 cm −1 due to dangling OH groups and lower frequency stretching modes < 3000 cm −1 arising from strong HBs. The stretching regions also have quite different total widths for the different phases and the bending mode of the b-SQ phase is softer than the other phases by ca. 120 cm −1 .
Therefore it should also be possible to discriminate one phase from another with vibrational spectroscopy.
Multi-layer square ice was also observed experimentally with a structure in which oxygen atoms are located directly on top of each other, suggesting an AA stacking square phase [10].
However, the force field simulations performed as part of the experimental study obtained an AB stacking arrangement. The difference was attributed to a possible inaccuracy of the force field models used. In order to clarify this difference between experiment and theory, we calculated the enthalpy for different confinement widths at 2 GPa for a number of AA and AB stacking structures (Fig. 4a). In agreement with the force field simulations we find AB stacking is preferred for the f-SQ phase. However the other square structure identified in this study, b-SQ, is more stable with AA stacking, which agrees with experiment. The main interaction between the layers in these structures is van der Waals, with the HB networks within the layers remaining unperturbed. We find, however, that these stacked van der Waals bonded structures are considerably less stable than double layer structures with interlayer hydrogen bonds, the most stable of which we dub "dl-SQ". The structure of dl-SQ is shown in Fig. 4b,c. It is AA stacked and as can be seen from Fig. 4a it is ca. 80 meV/H 2 O more stable than any of the stacked van der Waals bonded double layer structures. Thus we conclude that the AA stacked double layer square structure observed in experiments is unlikely to be a van der Waals bonded layered structure but rather one stabilised by interlayer hydrogen bonds.
In summary, monolayer ice phases and their phase transitions in confinement have been studied with DFT and a random structure search approach. At ambient pressure we have predicted hexagonal and pentagonal Cairo tiled structure, which are similar in enthalpy 9 and more stable than other structures. The CT structure becomes more stable than the hexagonal structure when lateral pressure is applied. Upon increasing the pressure to above about 2 GPa high density square and rhombic phases are observed. Looking forward a complete description of the temperature dependent phase diagram of 2D water is desirable.
Experimentally it would be interesting to explore a broader range of temperatures, water densities, and confining materials. From the computational perspective, it would be desirable to explore 2D ice at higher temperature with electronic structure methods, such as ab initio molecular dynamics. In addition, in light of the small ring sizes and relatively short intermolecular separations, 2D ice also provides further opportunities for investigating collective proton quantum dynamics [47][48][49].