Structural and configurational properties of nanoconfined monolayer ice from first principles

Understanding the structural tendencies of nanoconfined water is of great interest for nanoscience and biology, where nano/micro-sized objects may be separated by very few layers of water. Here we investigate the properties of ice confined to a quasi-2D monolayer by a featureless, chemically neutral potential, in order to characterize its intrinsic behaviour. We use density-functional theory simulations with a non-local van der Waals density functional. An ab initio random structure search reveals all the energetically competitive monolayer configurations to belong to only two of the previously-identified families, characterized by a square or honeycomb hydrogen-bonding network, respectively. We discuss the modified ice rules needed for each network, and propose a simple point dipole 2D lattice model that successfully explains the energetics of the square configurations. All identified stable phases for both networks are found to be non-polar (but with a topologically non-trivial texture for the square) and, hence, non-ferroelectric, in contrast to previous predictions from a five-site empirical force-field model. Our results are in good agreement with very recently reported experimental observations.

of square due to a distortion of the unit cell. These configurations are higher in density than the Archimedean tiling and hexagonal phases, and are therefore expected to be stabilized at high lateral pressure. Furthermore, there appears to be a transition between an almost completely flat rhombic phase, which competes with the low-density phases for small confinement widths, and a rippled rhombic phase, which is instead stable over almost all of the pressure range for larger confinement widths [46][47][48][49] .
Another interesting question is the ferroelectricity of monolayer ice. Both the hexagonal and the flat/rippled rhombic phases have been predicted to be ferroelectric by Zhao et al. 47 . A key prerequisite for this is that the phases are ordered, and that the most stable configuration in each case is one with a net polarization. Such a conclusion is supported for the rippled rhombic phase by the findings of Zangi and Mark 33 ; however, other studies report disordered quasi-square 33,35,41,43,45,46,48 or hexagonal 44 phases.
Very recently, experimental observations of monolayer ice have been presented 50 , showing a strictly square lattice of O ions. The overall picture from all previous computational studies does not provide clear support for this finding, and a physical explanation of the relative stability of the square and rhombic crystal symmetries is yet lacking.
The simulations undertaken so far have made use of well-established empirical force-field models of water: TIP4P 35,46,48 ,TIP5P 33,34,41,45,47,49 , and SPC/E 43,44,50 . Zhao et al. 47,49 have performed some additional checks with all of these models plus others (SPC, TIP3P, TIP4P/2005), as well as with first principles calculations using density-functional theory (DFT) with a semi-local (PBE-GGA 51 ) exchange and correlation (xc) functional. However, the DFT calculations were only carried out for a few examples in order to confirm the structural stability of the phases found by TIP5P; furthermore, these were not tested with the same confining potential, but either free-standing or confined between graphene sheets.
Most of the development and parameterization of the various empirical models has focussed on reproducing characteristics of the bulk phases (even for the bulk, however, there are known to be important discrepancies between them 52 ); their accuracy in a very different and extreme environment such as that considered here is an open question. Therefore, this is an area in which the predictive power of first principles calculations can provide valuable information. This is especially true thanks to the development of fully non-local xc functionals that account for van der Waals (vdW) interactions from first principles 53 ; these have helped to overcome the severe problems previously encountered by (semi-)local functionals in describing the properties of water, and have given a series of very promising results both for liquid water [54][55][56] and ice 57,58 .
In this paper, we investigate the phase diagram of nanoconfined monolayer ice entirely ab initio, using DFT with a non-local vdW xc functional. We do not base ourselves on previously identified phases, rather employing a preliminary ab initio random structure search 59 (AIRSS) procedure to identify promising low-enthalpy structures. The questions we aim to answer are: What are the intrinsically stable phases of monolayer ice over a range of confinement widths and lateral pressures? Are these phases ordered or disordered? What contribution do configurational entropy and vibrational effects have on phase stability? What factors determine the stability of individual configurations? Finally, is the electronic structure of the water molecule altered by confinement?

Results
Ab initio random structure search. In order to identify configurations of interest in an unbiased way, we begin our investigation by performing an AIRSS on the confined geometry. For a given value of the confinement width d and the 2D lateral pressure P l (see definitions in the Methods), we generate many trial unit cells, which are then relaxed. We use 100 unit cells of four molecules and 200 unit cells of eight molecules; the length of the two lattice vectors in the x-y plane and the angle between them are set randomly (within a sensible range), and then the water molecules are placed inside the confinement region with random positions and orientations. No symmetry is imposed on the system. We do this at each point on a regular grid on the P l -d phase diagram, for d from 5 to 10 Å in steps of 1 Å, and P l from 0.01 to 1000 GPa⋅ Å logarithmically in order of magnitude steps. For each point, we rank the 300 relaxed configurations by enthalpy/molecule (H = U + P l A).
The first important result is that stable monolayer configurations are only recovered for d ≤ 7 Å and P l ≤ 10 GPa⋅ Å. Therefore, we focus on this reduced part of the phase diagram, and leave the analysis of multilayer ices for the future 75 .
It is noticeable that, for all grid points in the monolayer region, the lowest-enthalpy structure found by our search is a square configuration (as explained previously, this refers to the network topology based on hydrogen bonding, and not necessarily the symmetry classification of the ionic structure, discussed below). Furthermore, for the points at d = 6 Å and d = 7 Å, the majority (75%) of structures within 50 meV of the lowest one belong to this square network (SN) family; the remaining configurations are SNs with bonding defects (20%), and a small number of bilayer configurations (5%) in the 40-50 meV range.
Instead, for the points at d = 5 Å, more variety is found within the 50 meV range: the number of non-SN configurations ranges from ~20% of the total at P l = 10 GPa⋅ Å up to ~90% at P l = 0.01 GPa⋅ Å. These can be classified into four types: honeycomb network (HN) configurations, amorphous configurations (made up of irregular combinations of squares, pentagons, hexagons, and some larger voids), SNs with bonding defects, and sub-2D networks (1D chains and ribbons, 0D square islands). The latter type can also be understood as SNs with systematic bonding defects, as the O ions are found to form quasi-square lattices. The HN and amorphous configurations are found in the higher energy range (20-50 meV above the SN ground state).
Configurational energetics and entropy. The AIRSS seems to indicate quite definitely the square as the most stable bonding network over the entire monolayer range. However, as with the bulk phases of ice, this does not correspond to a single ionic arrangement, even neglecting the possibility of bonding defects. The same is true of the HN. More importantly, such an analysis cannot tell us about the crucial entropic contribution. We therefore turn to a systematic investigation of the possible configurations for the two networks.
The square network. Figure 1(a) shows the six possible configurations for a four-molecule unit cell that are non-equivalent by symmetry. Many examples for all of these six are found by the AIRSS. As mentioned previously, each configuration exists in distinct flat and rippled arrangements, with a first-order phase transition in d between the two 47 . This is recovered naturally from the AIRSS, in agreement with previous studies [46][47][48][49] .
It is important to note that each molecule must donate one proton along a row of the grid and another along a column. Linear configurations (i.e., both protons along a single row or column) carry a significant energy penalty (they are found in some of the higher-energy defective cells). There is therefore both a decoupling of the ordering along rows and columns, and a long-range order imposed along a single row or column. As noted by Koga and Tanaka 35 , the number of allowed configurations for N molecules is 2 N 2 , which is not sufficient for a non-zero macroscopic entropy.
An important consequence of this is that, unless they are all energetically degenerate, all SN configurations correspond to separate phases, only one of them thermodynamically stable for any given set of external conditions. No transition from an ordered to a disordered phase on the same lattice is possible, unlike what happens in bulk between ice XI and ice Ih. Figure 2 shows the energetics of the different SN configurations. Surprisingly, the most stable is Ab/Cd, which is the only one with no net polarization (see Fig. 1(a)). This is confirmed for both the flat and rippled arrangements by the symmetry analysis of the relaxed crystal structures, discussed below. Conversely, the configuration with the largest polarization (Aa/Aa) is the highest in energy. The energetic ordering of configurations is constant for all d and P l . Therefore, our simulations predict the stable SN phases to be ordered and non-polar, and, hence, non-ferroelectric.
Despite going against previous predictions, the energetic ordering of configurations is readily explained by a simple electrostatic model of identical point dipoles on an ideal 2D periodic square grid (Fig. 3). We have calculated the energy of the six configurations using this model, perfectly reproducing the DFT ordering and giving excellent estimates of relative energy differences, also shown in the figure.
In this model it is straightforward to see how every water molecule is represented as a dipole. However, it is also possible to approach the problem from a different limit, that of an idealized arrangement of symmetrized hydrogen bonds. In this case, dipoles can be obtained by considering the shift of each proton away from the bond centre, resulting in two smaller dipoles per water molecule. It can be shown that the two limits (protons close to O ions/protons close to bond centres) are perfectly equivalent; furthermore, a point charge model in which the protons are free to move between the two limits shows that the ratios in energy differences between configurations are constant at all points to a very good approximation (Supplementary Note 1). This goes some way towards explaining the robustness of such a simple model.
A close analysis of the dipole model (Supplementary Note 1) reveals that there is only a single feature accounting for practically all of the energy difference between configurations: the orientation of protons in neighbouring rows or columns of the SN grid (it is trivial to show by symmetry that energy differences are decoupled between Scientific RepoRts | 6:18651 | DOI: 10.1038/srep18651 rows and columns). The energetically favourable arrangement is for the two rows/columns to have opposite proton orientation; if instead it is the same, there is a relative energy penalty of: where n is the number of lattice spacings separating the two rows/columns on the lattice, and the energy scale is p 2 /(4πε 0 s 3 ), determined by the magnitude of the molecular dipole p, and the lattice spacing s. For nearest  neighbours (n = 1), this gives a value of 7.28 × 10 −2 p 2 /(4πε 0 s 3 ). For second-nearest neighbours (n = 2), the value is already three orders of magnitude smaller; therefore, contributions from n ≥ 2 can be safely neglected, and only nearest-neighbour interactions need to be considered. The schematic configuration diagrams from Fig. 1(a) can now be used to accurately estimate energy differences given by the dipole model (and, hence, by DFT), simply by counting the number of neighbouring rows and columns in the unit cell with the same orientation: four for Aa/Aa, three for AaAb and AaBb, two for Ab/Ab and Aa/Cc, and none for Ab/Cd. More generally, the average number of parallel rows/columns per molecule (P number) is a useful measure for classifying all allowed SN configurations independently of unit cell size. Its value can range between zero and four (i.e., the number of neighbours/molecule). Therefore, Ab/Cd must be the lowest-energy configuration for any size unit cell, and Aa/Aa the highest; it would be appropriate to rename them according to their P number. A pair of parallel rows/columns can be seen as a defect of the P = 0 ground state, with a formation energy of E(1). Such a defect behaves similarly to a domain wall, which can propagate through the lattice but only be annihilated by another equivalent defect.
It is interesting to note that the P = 0 ground-state configuration has two unique features: (a) it provides the best possible screening of electrostatic interactions, leading to the fastest convergence of the dipole model energy with system size; and (b) it is the only one out of the six four-molecule configurations to exhibit topological charges, with the dipoles arranged into two square sublattices of vortices and antivortices (equal and opposite charge). All other configurations, instead, have infinite field lines for the coarse-grained vector field of dipoles. In general, sources and sinks are prohibited, as well as net electric charges.
We can therefore summarize the expanded ice rules governing the arrangement of water molecules in a SN monolayer as follows: • The four-fold coordination of each molecule is still valid.
• The presence of two short and two long OH distances for each O ion is also still valid, but there is an additional restriction forbidding linear molecular configurations (variants of this rule have already been noted by several authors 33,35 . This restriction prevents proton disorder. • The energetic ordering of the allowed configurations depends linearly on the number of neighbouring rows/ columns with the protons along both being oriented in the same direction. The honeycomb network. Graphene-like honeycomb monolayer ice cannot obey the bulk ice rules; the number of bonds/molecule ensures that half the molecules must donate two protons but receive only one, while the other half must donate one and receive two. This leads to N/2 dangling protons, which can be either above or below the monolayer plane (as shown in Fig. 4), and N/2 dangling lone pairs. Figure 1(b) shows the ten non-equivalent configurations for a four-molecule unit cell (note that each of the five configurations shown has two possible placements of the dangling protons, both in the same direction or in opposite directions). It is therefore reasonable to expect an energy penalty, but also a non-zero entropy, and, hence, a single disordered phase which is stabilized with temperature respect to the ordered SN phases. A naive estimation of the configurational entropy using Pauling's approach 15 gives W = 2 N/2 × (9/2) N/2 = 3 N , where the first term accounts for the up/down position of the dangling protons, and the second for satisfying the ice rules for the in-plane bonds, and so S = k B ln3. This is almost three times that of ice Ih. Given the enthalpy difference between the SN and HN configurations shown in Fig. 2, this would be sufficient to obtain a transition between the two at room temperature.
However, this is not the case. Indeed, contrary to the SN, the AIRSS finds only a subset of the allowed HN configurations. We test this systematically by building all ten unit cells by hand and relaxing them. The results confirm the fact that not all configurations are stable (details are given in Fig. 1(b)). The unstable ones either spontaneously relax to one of the stable ones, or remain trapped in a high-energy metastable state. Furthermore, the available data suggests that the instability arises from placing two in-plane or two out-of-plane molecules next to each other on the lattice; the stable configurations can then only be those for which the two types of molecules are segregated into the two triangular sublattices (i.e., the a and b sites). Interestingly, this results in neighbour pairs similar to the favourable configuration reported for polyhedral water clusters 60 . A further confirmation is obtained by testing two eight-molecule unit cells, one which obeys this restriction (aAaA/aAaA) and one which does not (abBA/BAab). The latter relaxes to the former, which is stable regardless of the placement of the dangling protons. We note that the energies of all the stable four-and eight-molecule configurations are almost degenerate, within 5 meV (Fig. 2).
Interestingly, the stable configurations are not the ones expected by considering the possible arrangements of layers in bulk hexagonal ice. This is because in the bulk the dangling protons are locked into hydrogen bonds, while in the monolayer they are free to rotate. This gives rise to a mechanism by which two neighbouring molecules, one in and one out of plane, can both rotate in concert to swap the proton in the hydrogen bond connecting them. This is therefore a localized relaxation pathway for moving between configurations, which can be expected to have a small or even no energy barrier.
The reduction in the number of allowed configurations has a considerable impact on our previous estimate of the entropic contribution, which we can now revise to W = 2 N/2 × (9/8) N/2 = (3/2) N , and so S = k B ln 3/2, equivalent to ice Ih. Consequently, the phase transition is shifted upwards to significantly higher temperatures (≃ 860 K), clearly no longer relevant for the solid phase.
Our conclusion, therefore, is that the HN phase is only ever metastable (we note that this analysis has been carried out at the most favourable point on the phase diagram for HN configurations). It is predicted to be disordered, and, hence, non-polar and non-ferroelectric. The lowest-enthalpy configuration, which should constitute an ordered phase at very low temperatures, cannot be identified with certainty; from Fig. 2, it seems that relaxation effects in larger cell sizes might be crucial in stabilizing some configurations over others.
The ice rules governing the HN monolayer are starkly different both to the bulk and the SN: Scientific RepoRts | 6:18651 | DOI: 10.1038/srep18651 • Each molecule must be three-fold coordinated, with half the molecules featuring out-of-plane dangling protons.
• Two in-plane or two out-of-plane molecules cannot be neighbours in the network.
• In-plane proton disorder is allowed, following the 'two short/two long' rule for OH distances.
• Proton disorder for the dangling protons (above/below the plane) is allowed with no restrictions, and is decoupled from the in-plane disorder.
Vibrational effects. We have checked for any vibrational effects by calculating the phonon frequencies of the lowest-enthalpy four-molecule configurations for the two networks at d = 5 Å, P l = 0.01 GPa⋅ Å (flat Ab/Cd for the SN, and aAaA for the HN). We use the finite displacement method with a (5 × 5) supercell and an ionic displacement of 0.02 Å. The phonon DOS is shown in Supplementary Fig. 3. The zero-point energy of the two monolayer configurations is almost identical, 111 ± 1 meV/molecule in both cases. This is ~15% smaller than the value reported for bulk ice at a similar level of theory 58 . However, the vibrational free energy calculated from the phonon DOS decreases slightly faster for the SN as the temperature increases (e.g., a relative gain of 10 meV at T = 430 K). This further increases the stability of the SN phases at finite temperatures, thus reinforcing our previous conclusion on the metastability of the HN phase.
We have also checked that vibrational effects are very similar between SN configurations, and that switching from light to heavy water, despite lowering the zero-point energy by 26%, is an entirely negligible effect when considering the relative stability of different configurations and phases.
Crystal structure. We now examine in detail the ionic positioning within the unit cell, and the shape of the cell itself. We do so for a fine grid of points in a range of confinement widths from 5 to 6 Å and lateral pressures from 0.01 to 10 GPa⋅ Å. Figure 5 shows how the enthalpy of the SN and HN phases depends on the confinement width for near-zero lateral pressure. The phase transition from the flat SN phase for small d to the rippled SN phase for large d can be observed at ~5.5 Å. The discontinuity is shown clearly by the contribution to the total enthalpy from the external wall potential; when viewed with respect to the well depth at the centre of the confinement region (also plotted), this gives a measure of the amplitude of the monolayer in z. While the HN follows this baseline value quite closely, thereby remaining fairly flat for all d, the SN switches discontinuously between an almost perfectly flat arrangement and a rippled one. The ripples are seen in the z-position of the O ions for alternate rows of molecules, as shown in Figs 1(a) and 4. There is a clear gain in internal energy from buckling (related to the improved alignment of the hydrogen bonds), which competes with the penalty from the confining potential; hence, the transition when increasing d.
In terms of Bravais lattices, the flat SN phase is indeed perfectly square within the relaxation tolerance, as should be expected from the symmetry of the bonding arrangement; the layer group is p4b2 (p4/mbm for the O sublattice). Once ripples are introduced, the symmetry is reduced; the rippled SN phase is therefore slightly rectangular in practice, with an a/b ratio of 0.9 at d = 5.5 Å, P l = 0.01 GPa⋅ Å, and decreasing with both increased d and P l ; the layer group is now p2 1 /b11 (equivalent for the O sublattice). Nearly all configurations from Fig. 1(a) follow a similar pattern; the only two exceptions, Aa/Aa and AaAb, are the ones in which there is a non-zero net polarization which is not parallel to one of the bonding directions. In these cases, the unit cell switches from rhombic (flat) to oblique (rippled). The distortion in the unit cell angle with respect to the square/rectangular configurations is on the order of 1-10%. For Aa/Aa, the layer group of the flat arrangement is p211 (cmmm for the O sublattice), and that of the rippled arrangement is p11a (p112/a for the O sublattice).
The HN phase is close to hexagonal in symmetry for all d and P l , and for all stable configurations from Fig. 1(b). Distortions are small, on the order of 1% in the unit cell angle, and 1-10% in the a/b ratio. Since the phase is disordered, only the symmetry of the O sublattice averaged over all configurations is relevant; the layer group is therefore p6/mmm. Figure 6 shows how the enthalpy and area/molecule vary with the confinement width and lateral pressure for the SN and HN phases (note that the area of the flat SN phase is almost completely insensitive to d, and that there is a discontinuous transition to the rippled SN phase, which can be seen in the large gap between curves for the SN in Fig. 6(b)). An interesting phenomenon is observed: the metastable HN phase becomes unstable above a critical pressure, and spontaneously relaxes to a SN configuration. The critical pressure decreases with d. This explains why the AIRSS only finds HN configurations for small values of d and P l .
There are three possible strains which transform a hexagonal cell to a square one (Fig. 7). Using the cell angle as the order parameter, we follow this transformation for a representative confinement width, and show the effect on the enthalpy caused by varying the lateral pressure. The behaviour is that of a first-order phase transition; therefore, what is observed is the crossing of the phase stability limit. This also suggests the existence of a phase transition from the SN to the HN, although the critical pressure is in the negative regime.
It is also interesting to note that the distortion can be seen to have an intermediate stage, that of molecules arranged on a square or quasi-square lattice but maintaining the honeycomb bonding network, as illustrated in Fig. 7. This is because half of the molecules in the HN configuration need to rotate to bring their dangling proton in plane and complete the SN bonding. This intermediate configuration can itself sometimes be metastable (especially for small values of d), and, indeed, accounts for some of the previously mentioned low-lying defective SN configurations found by the AIRSS. It is therefore reasonable to assume that this mechanism will be important for introducing a degree of disorder in the SN phases at moderately high temperatures, and as a way of propagating changes between different configurations through the lattice.  Electronic structure. In order to examine how the electronic charge is distributed around the molecules, and how this differs with respect to bulk water, we turn to a maximally-localized Wannier function 61 (MLWF) analysis of the system. The MLWFs are obtained from our DFT calculations by interfacing 62 with the Wannier90 63 post-processor code. Figure 8 shows the typical MLWFs for the water molecule. The molecular bonding is very strong, so we cannot expect a qualitative change in this picture; however, subtle changes in the electronic structure will be reflected in the position of the MLWF centres and their spreads. We analyze the monolayer SN and HN phases, as well as bulk ice XI and an isolated water molecule. We note that the monolayer results are very similar for the different configurations of Fig. 1, and also for the flat and rippled SN arrangements.
In several respects the monolayer MLWFs are found not to differ significantly from the bulk crystal. This is the case for the spreads: 0.49 ± 0.01 Å 2 for the OH bonding orbitals, and 0.58 ± 0.01 Å 2 for the lone pairs. (The values are somewhat lower for the isolated molecule: 0.47 and 0.52 Å 2 , respectively.) The distances between the WF centres and the O ion are also very similar, including for the isolated molecule: 0.50 ± 0.02 Å for the OH bonding orbitals, and 0.32 ± 0.02 Å for the lone pairs.
The main difference is found for the HN monolayer when considering the angle between pairs of MLWFs and the O ion. This is shown in Fig. 8. Interestingly, the values obtained for the SN phases are very similar to those of bulk ice. Within the level of detail offered by the MLWFs, therefore, the electronic structure of the water molecules can be considered not to be altered significantly between these phases. The HN monolayer, however, is more complicated. First of all, there is a noticeable difference between the out-of-plane and in-plane molecules. The MLWFs of the former are positioned similarly to those of an isolated molecule, as can be seen in the figure (this is true for all of the molecule's MLWFs, not only the one associated with the dangling proton). The MLWFs of the latter are positioned differently from all other phases. In particular, the angle between the two lone-pair orbitals is reduced; this can be understood to favour the formation of a hydrogen bond pointing between the two, as is required for the quasi-planar sp 2 -like bonding of the honeycomb lattice.
The magnitude of the molecular dipole moment is found to be quite insensitive to this angle: it is 3.0 ± 0.04 D for all molecules in the SN and HN phases, compared to 3.4 ± 0.06 D for ice XI and 2.0 ± 0.06 D for the isolated molecule (additional calculations show that there is a consistent overestimation of 5-10% due to the limitations of the basis set). It is also important to note that the changes in the angle between lone-pair orbitals are not correlated with the HOH angle. Therefore, we expect the HN monolayer to pose a particularly interesting challenge for force-field modelling.

Discussion
The picture offered by first principles simulation differs significantly from the previous review of empirical force-field results 49 , mainly based on data obtained with the TIP5P model. The overall phase diagram is greatly simplified by ruling out the existence of any stable phase other than the SN ones. We note that the previously reported 41,47 Archimedean tiling phase could not be recovered from our AIRSS due to the limitation in unit cell size; therefore, we have checked this separately by relaxing the ionic configuration reported by Zhao et al. 47 , finding it to be high in enthalpy (44 meV/molecule above the ground state in Fig. 2, at the same level of the metastable HN configurations), thus not affecting our conclusions. As with the HN configurations, the double hydrogen bond arrangement between pairs of atoms predicted by TIP5P is found to be unstable, with our relaxation resulting in one dangling proton per two molecules.
The SN and HN phases are also found to differ from these previous results in terms of their proton ordering, crystal symmetry and net polarization. It is particularly noticeable that our prediction for the stable flat SN phase is in perfect agreement with recent experimental observations 50 , both in terms of the square symmetry of the O sublattice and the value of the lattice parameter (2.82 Å, compared with 2.83 ± 0.03 Å from experiment). The first principles approach, therefore, successfully untangles the contradictory predictions from different force-fields models; furthermore, the point dipole 2D lattice model provides a simple physical explanation as to why the Ab/Cd configuration is the most stable.
It is surprising that previous force-field studies do not all recover the correct ordering of SN configurations (with some models placing the most polar and rhombic one, Aa/Aa, as the most stable instead of the least stable). This is unexpected because the striking agreement between our simple dipole model and DFT suggests that the ordering is entirely dominated by electrostatics, which all force-field models should have no problem in reproducing. Indeed, we have confirmed this to be the case by comparing the electrostatic interaction of different configurations of water dimers placed on the SN for DFT and a variety of force-field models (TIP3P, TIP4P, TIP5P). The value for DFT is approximated by a seven-site point charge model, using the three ionic positions and the four MLWF centres. Both the DFT and force-field models are in excellent agreement with the dipole model.
However, previous studies using TIP5P show that the relaxation of different configurations is not consistent with DFT. In particular, the polar Aa/Aa configuration is severely distorted from the idealized square lattice, with a unit cell angle of ~70%. Such a distortion significantly lowers the electrostatic energy of Aa/Aa with respect to Ab/Cd, reversing the ordering of configurations. For DFT, instead, the rhombic distortion is much smaller, and so the dipole model remains accurate.
In general, our results suggest that TIP5P suffers from an over-emphasis on the tetrahedral bonding, resulting in overly distorted SN configurations, as well as the incorrect bonding for HN configurations. Results for other force-field models are much more limited, although there is some evidence that three-and four-site models correctly order the SN configurations 47,49,50 , which is consistent with their reduced emphasis on the tetrahedral bonding geometry.
Another important point to make is that the results presented here, in particular the relative stability of the SN and HN phases, will depend on the detailed form of the confining potential. We have purposefully chosen one of the simplest possible confinement schemes in order to investigate the intrinsic tendencies of monolayer ice, without bias from a particular chemical environment. Therefore, our results offer a comprehensive characterization of the most relevant possibilities for water to crystallize in a monolayer, which will then be subjected to different such biases in different realistic scenarios.
In particular, several reasonable modifications to the confining potential might favour the HN over the SN, e.g., adding an attractive term to the H ions, and/or introducing a periodic modulation to the walls with a quasi-commensurate hexagonal lattice, typical of many surfaces. Such modifications can be expected to decrease the energy penalty of the HN, and shift the first-order phase transition with the SN towards positive pressures. Indeed, various specific cases have been shown to favour a HN configuration [18][19][20]43,44 . It is however interesting to note that the only available experimental results show the intrinsic tendency towards the SN phase to overcome the hexagonal confining potential provided by graphene sheets 50 . The formation of monolayer ice requires very small, sub-nm confinement regions; it is therefore reasonable to question the possibility for experimental synthesis, or for natural formation in biological or geological structures. A simple thermodynamic argument, combined with predictions from DFT, suggests that only moderately high pressures are necessary to cause the insertion of ice at equilibrium between sheets with such narrow spacing. We reach this conclusion by comparing the chemical potential of water molecules in bulk ice XI under isotropic pressure and our most stable nanoconfined phases; the two are equilibrated at a value of the external pressure of ~1 GPa for d = 6 Å (increasing to ~2 GPa for d = 5 Å). The value will depend both on the confinement width and the form of the confining potential, which defines the well depth (as shown in Fig. 5) with respect to the zero of potential outside the confinement region. A strongly hydrophilic or hydrophobic surface will shift the well depth, and, hence, the external pressure required for insertion.

Methods
Simulation setup. All the simulations are performed with the SIESTA 64 DFT code, with norm-conserving pseudopotentials in Troullier-Martins form 65 and a basis set of finite-range numerical atomic orbitals 66 (NAOs). The pseudopotentials are the same as those described in a previous study of liquid water and ice 67 . We employ a variationally-obtained 68 double-ζ polarized NAO basis, also used in previous studies 56,67,69 . This basis has been extensively tested against both larger NAO bases and plane-wave calculations 56,67 (in these studies, it is referred to as (P) dζ + p); it has shown a high degree of transferability, and accuracies comparable to plane-wave cutoffs of ~850 eV for energy differences and ionic forces, and ~1500 eV for absolute pressures. SIESTA represents the electronic density on a real-space grid, for which we use a cutoff of 2040 eV (150 Ry).
For xc, we use the vdW-DF functional by Dion et al. 53 , but substituting the revPBE exchange energy with PBE, as previous studies have shown this to be among the best possible descriptions currently available for water from first principles 55,56,69 .
The confining walls are described by a classical Lennard-Jones (L-J) 9-3 potential with parameters σ = 2.5 Å and ε = 13.0 meV (1.25 kJ/mol), applied only to the O ions; the separation between the origin of the two walls along z is referred to as d. This potential is structureless in the x-y direction, and was originally developed to approximate a hydrophobic paraffin-like surface 28 . We chose it due to its simplicity, and because it is the same one used in many previous studies, thus allowing for a direct comparison of results 34,35,41,[46][47][48][49] . Tests on our system show that the results are not significantly affected by the details of the confining potential, in agreement with the study of Kumar et al. 37 .
The simulation cells are periodic in all directions; we add a buffer region of 11 Å in z, and use a dipole correction for slab geometries to prevent artifacts. We use a Monkhorst-Pack (MP) k-point sampling grid 70 in the x-y plane only, with a length cutoff of 10.6 Å (20 a 0 ) as defined by the scheme of Moreno and Soler 71 . Convergence tests were performed for the k-point sampling, size of the buffer region, and real-space cutoff, to ensure an accuracy within 1 meV/water molecule.

Structural relaxations.
Relaxations are performed at particular values of the 2D lateral pressure (defined below), with no restrictions on either the ionic positions or the two vectors defining the 2D unit cell (the third vector in z is held constant and perpendicular to these). The convergence tolerances used are 1 meV/Å for the maximum ionic force, and 0.1 MPa for the maximum error on a component of the stress tensor. A robust and reliable relaxation procedure is achieved by cycling repeatedly through 20 steps of conjugate-gradient minimization, 20 steps of a modified Broyden optimization 72 , and 20 steps of the FIRE algorithm 73 , until reaching convergence.
Lateral pressure definition. The 2D lateral pressure P l which we refer to throughout the paper is given in dimensions of [pressure]⋅ [length] (GPa⋅ Å). This is because, for the ground-state calculations performed here, the conventional components of the stress tensor are ill-defined due to the arbitrariness in setting a length for the z dimension of the unit cell. An approach based on the virial 34 is also not feasible, as we do not perform dynamic simulations. Therefore, we use a (2 × 2) 2D stress tensor with components σ σ = l ij z ij 2D 3D , where {i, j} = {x, y} and l z is the length of the unit cell vector in z used for the simulation. The 2D lateral pressure is then P l = − ∂E/∂A, where E is the ground-state energy of the system, and A is the area of the unit cell in x-y. A simple estimate of the 3D lateral pressure can be obtained by using the confinement width as a reasonable approximate z-width for the system: .