Evidence for Stable Square Ice from Quantum Monte Carlo

Recent experiments on ice formed by water under nanoconfinement provide evidence for a two-dimensional (2D) `square ice' phase. However, the interpretation of the experiments has been questioned and the stability of square ice has become a matter of debate. Partially this is because the simulation approaches employed so far (force fields and density functional theory) struggle to accurately describe the very small energy differences between the relevant phases. Here we report a study of 2D ice using an accurate wave-function based electronic structure approach, namely Diffusion Monte Carlo (DMC). We find that at relatively high pressure square ice is indeed the lowest enthalpy phase examined, supporting the initial experimental claim. Moreover, at lower pressures a `pentagonal ice' phase (not yet observed experimentally) has the lowest enthalpy, and at ambient pressure the `pentagonal ice' phase is degenerate with a `hexagonal ice' phase. Our DMC results also allow us to evaluate the accuracy of various density functional theory exchange correlation functionals and force field models, and in doing so we extend the understanding of how such methodologies perform to challenging 2D structures presenting dangling hydrogen bonds.

Recent transmission electron microscopy (TEM) measurements and classical molecular dynamics simulations report that a new two dimensional (2D) square phase of ice forms. 1 This phase is not part of the bulk ice phase diagram, but it was suggested that it is stabilized under confinement because of lateral pressure (estimated to be in the gigapascal, GPa, regime), arising from the van der Waals (vdW) attraction between the graphene sheets. However, these experiments have been questioned, 2 with it even being suggested that it is sodium chloride contamination and not ice that is responsible for the square symmetry observed. So far, it is not clear under what conditions (if any) square ice is stable.
Theoretical investigations of the stability of confined 2D ice at high lateral pressures can, in principle, help in disentangling this issue and in complementing experimental findings. From a theoretical perspective, the prediction of square 2D ice can be traced back to Nagle's 1970s 'unit model' of ice. 3 However, later atomistic force field (FF) simulations found that 2D ice prefers a buckled rhombic structure. 4,5 More recently, density functional theory (DFT) based investigations [6][7][8][9] have been performed. However these have produced qualitatively different results depending on the precise details of the calculations and, in particular, on the choice of exchangecorrelation (XC) functional. For instance,Chen et al. 7 found hexagonal and pentagonal structures to be stable phases at low pressures ( Fig. 1). They also found that square ice is only stable in the GPa pressure regime and that even at these pressures it is only favored enthalpicly by < 10 meV/H 2 O (1/4 kcal/mol). Meanwhile, DFT results ofCorsetti et al. 6 showed that the square ice structure is more stable than the hexagonal phase at all pressures. These disparate findings raise serious questions about the reliability of the adopted computational approaches (both DFT and FF) applied so far to 2D ice.
Indeed this is exemplary of a broader longstanding issue: the water and ice phase diagram is extremely challenging for any computational approach, because there can be competing phases within an energy range of only a few tens of meV/H 2 O. 10 Achieving this accuracy is often beyond the capabilities of DFT XC functionals and most FF approaches. 11,12 Therefore, a study of 2D ice with a more accurate theoretical method is needed.

Rhombic (R ) Square (S )
Hexagonal (H ) Pentagonal (P ) The development and application of electronic struc-ture approaches with meV accuracy is now a thriving area of condensed matter research (see e.g. Refs. [13][14][15][16] ). Of the various methods available diffusion Monte Carlo (DMC) is particularly attractive. [17][18][19][20][21][22][23] First, DMC has already been shown to offer the requisite accuracy for bulk ice phases by producing results in excellent agreement with experiment . 10,11,24,25 Second, thanks to recent improvements in computational efficiency, 26 it is now possible to obtain converged results on the large unit cells that must be considered when treating 2D ice with a many-body electronic structure approach. With this in mind, herein we report a DMC study of 2D ice. Our DMC calculations reveal that hexagonal and pentagonal phases are indeed the most stable 2D ice phases identified at low pressures. Perhaps of more interest though, DMC calculations at 2 GPa clearly support the existence of the square ice at high lateral pressure. As a further step we then use our DMC reference data to understand how the much more widely used DFT and FF approaches perform for such systems, so as to provide guidance for future studies on 2D ice and interfacial water. From this, several DFT XC functionals which perform well are identified and we also find that commonly used FFs such as SPC/E 27 and TIP4P 28 tend to overstabilize high density phases; helping to explain why such phases have been widely observed in FF studies.
Previous DFT studies suggest that four structureshexagonal, pentagonal, square, and rhombic (see Fig. 1)are the most stable monolayer ice structures, 7 so these are the focus of the current study. In the rhombic and square structures the water molecules are fourfold coordinated, as in bulk ice, but are arranged in the plane. The pentagonal and hexagonal structures have water molecules arranged similarly to a cut in bulk ice, so that there are some water molecules that are not fourfold coordinated and some with dangling hydrogen bonds. Specifically, in the hexagonal structure all water molecules are threefold coordinated and 1/2 of them have one dangling hydrogen bond; in the pentagonal structure 1/3 of the water molecules have one dangling hydrogen bond and are fourfold coordinated, the remaining molecules are threefold coordinated and have no dangling bonds. The casino code 29 has been used for the DMC calculations using Dirac-Fock pseudopotentials 30,31 with the locality approximation. 32 Slater-Jastrow trial wave functions with single Slater determinants were used and the single particle orbitals obtained from DFT-LDA plane wave calculations reexpanded in terms of B-splines. 33 We used the DMC algorithm with the prescriptions of,Zen et al. 26 which allows a DMC time step as large as 0.02 a.u. to be used with negligible time step errors. DMC calculations have been performed on structures optimized with DFT at a range of lateral pressures and under a uniform 2D confining potential fitted to DMC for values for the water-graphene interaction . 34 See Ref. 7 and the supplemental material therein for details of the confining potential used and discussions on the impact of using explicit graphene as the confining material. Structures used for DMC calculations have been obtained by performing geometry optimizations with the optPBE-vdW functional. 35 Additional calculations have also been performed on selected structures obtained from the revPBE-vdW functional, 35  In order to assess the stability of the different 2D ice structures, it is necessary to consider their enthalpy (H) at 0 K. 36 The enthalpy is the sum of the binding energy, E b 37 , the confinement energy, E conf , the pressure volume work (P V ), and the zero point energy (ZPE). To begin, we first discuss the accurate evaluation of E b using DMC. In DMC, and other many-body methods, finite size (FS) errors can be sizeable unless large supercells and/or correction terms are considered. [38][39][40] Since we are concerned with very small energy differences between the various phases, we have carefully addressed this issue with a series of calculations for increasing supercell size. Supercells with 8 up to as many as 192 water molecules were considered and for each system considered DMC simulations were run until the stochastic error of E b was ≤ 3 meV/H 2 O. The results obtained are plotted in Fig. 2 as a function of the inverse number of water molecules in the simulated supercell, 1/N w . From these calculations we have extrapolated E b for the different structures to infinite system size. The extrapolated binding energies are summarized in Table I from where it can be seen that the hexagonal phase has the largest binding energy and the pentagonal phase is only marginally (4 meV/H 2 O) less stable. These calculations also reveal that the FS error leads to an overestimate of E b by more than 10 meV/H 2 O in the smallest cells (N w < 20), and the relative energies between the different structures are fairly insensetive to the size of cell used. We note that the N w → ∞ extrapolations are obtained by assuming that FS errors are proportional to 1/N w . In the S.M. we show that other choices do not alter the extrapolated values of E b . 36 TABLE I. Properties of the various 2D ice structures at 0 pressure and 2 GPa. The values reported are the binding energy E b (meV/H2O) obtained from DMC, the confinement energy E conf (meV/H2O), the lateral area per water A (Å 2 /H2O), the pressure volume work P V = P × A × w (meV/H2O; w is the width of the confinement), the zero point energy (ZP E) (meV/H2O), and the enthalpy H (inclusive of ZP E). Also reported is the enthalpy difference with respect to the most stable structure: ∆pen = (H) − (H) [ With the binding energies obtained, the enthalpies can be further calculated considering E conf , P V , and ZPE. Each of these additional terms has been obtained with DFT and upon putting everything together we find that at vanishing pressure the pentagonal and hexagonal structures have essentially the same enthalpy. The values obtained come within 1 meV/H 2 O of each other which is within the stochastic errors of our DMC simulations (3 meV/H 2 O). Although it is hard to differentiate the pentagonal and the hexagonal structures at vanishing pressure, based on the enthalpies obtained we would expect the higher density pentagonal phase to become more stable than the hexagonal phase at small applied pressure, as first proposed using DFT. 7 The square and rhombic structures are 13 and 42 meV/H 2 O higher in enthalpy, respectively, than the pentagonal structure at zero pressure. Note that the rhombic structure is buckled, thus also the confining energy penalizes it compared to the other structures.
In the 2D ice phase diagram, the question of the stability of square ice is particularly interesting because it is arguably the only experimentally observed 2D ice and disagreements exist between FFs and DFT, and also within DFT itself. 1,2,6,7 We have already seen that the square phase is less stable than the pentagonal and hexagonal phases at low pressures. However, it has been estimated that the pressure in graphene nanocapillaries can be as large as several GPa due to the van der Waals forces pulling the graphene sheets together. 1 Therefore, we carried out DMC calculations on structures at 2 GPa. In Table I we report results for pentagonal, square and rhombic structures; the hexagonal structure becomes unstable at 2 GPa. We find that the square structure is indeed the most stable. Its enthalpy is lower than the pentagonal structure by 28 meV/H 2 O and lower than the rhombic structure by 17 meV/H 2 O. Although this study has focussed on relative enthalpies at 0 K, we have also estimated the relative free energies of the various phases by taking the vibrational contributions to the free energies into account. As shown in the S.M. (Table S.II) at 300 K the square phase remains the most stable phase at 2 GPa. DMC has helped to clarify the relative stabilities of the various 2D ices at ambient and high pressure. We now use these DMC benchmarks to understand how various DFT XC functionals and FF models perform on such structures. This is important to establish as DFT and FFs are widely used to examine 2D ice, interfacial, and confined water. Here we discuss the main indications that come out of these comparisons. In Fig. 3 (left panel) we plot the average error in the binding energy of 2D ice structures as obtained from a variety of approaches. From this we find that several XC functionals, namely HSE-vdW(TS), rPW86-vdW-DF2, optB88-vdW, optB86b-vdW, optPBE-vdW, BLYP-D3 atm , and PBE perform well, yielding an average error of < 25 meV/H 2 O. Apart from binding energies averaged over the four structures, it is also crucial to predict correct relative energies. As shown by the relative energies of the four 2D ice structures ( Fig. 3 right), we find that only the vdW inclusive approaches mentioned above provide satisfactory predictions. The long range part of the vdW force therefore plays an important role in 2D ice. Interestingly, SCAN, a recently developed meta-GGA functional that partially accounts for the medium range vdW force, 42,43 does not perform particularly well for these systems and overbinds all phases. Extending the analysis, in Fig. 4 we plot the energy error of several XC functionals and FFs for both 2D and 3D ice. Results are shown in order of increasing density for both 2D and 3D ice. On comparing 2D and 3D ice we find that whilst vdW is important in 2D it plays a smaller role than it is known to play in 3D. 10,24 This can be seen for example by comparing the slopes of the PBE curves in 2D and 3D, where it can be seen that the PBE curve in 3D is much steeper. We also show in the S.M. that the reduced significance of vdW is due to the lower coordination of water molecules in 2D ice than in 3D. The difference of the vdW contribution in 2D and 3D leads to contrasting performance of some of the vdW functionals for 2D and 3D ice, as shown in Fig. 4. Overall we find that the rPW86-vdW-DF2 functional yields the smallest errors for both 2D and 3D ice structures, and represents a good choice for future studies of these and related systems.
In Fig. 4 we also present results with several FF models which have been widely used to study confined water and ice. 1,4,5,44,45 Our calculations show that the SPC/E and TIP4P models overestimate the binding energy of the high density 2D ice structures more than the low density ones. The difference between the errors on the hexagonal and rhombic structures is about 40 meV/H 2 O for TIP4P and about 60 meV/H 2 O for SPC/E. This error comes from both the Coulomb and the van der Waals components of the potential (Fig. S6). Results of TIP4P/2005 46 and TIP4P/Ice 47 are also reported in the S.M. and typically they show a consistent shift of energy for all structures compared with TIP4P. Therefore, our results help to explain why rhombic and square structures are often seen in 2D ice simulations using these models. We have also performed calculations using the mW model, 48 a widely used coarse grained model of water. Although it significantly underestimates the binding energies of the 2D ice structures by ca. 50 meV/H 2 O, it does very well in reproducing the relative energies of 2D ice. This good performance on the relative energies is also consistent with the fact that a pentagonal ice structure has been observed in mW simulations of confined water. 45 To conclude, our large scale DMC simulations of 2D ice reveal that at ambient pressure the most stable structures are hexagonal and pentagonal phases. In addition, our calculations show that at high pressure 2D square ice is more stable at 0 K than the other phases considered. The data and insight obtained here is important for current understanding of confined 2D ice in general and is also relevant to the TEM measurements of Algara-Siller et al.. 1 We note that although our calculations find that square ice is stable they do not rule out the possibility, as suggested by Zhou et al., that sodium chloride could have been observed in the initial measurements. 2 The influence of finite temperature, growth kinetics, and the finite size of the particles that form have also yet to be investigated in detail. Clearly more work is needed from experiment and theory in order to substantiate the existence of square ice, and to potentially observe the predicted pentagonal ice phase. Finally, from a theoretical point of view, our study reveals the different performance of many widely used DFT functionals and FF models on 2D ice. We find that the role of vdW forces in 2D and 3D ice is different, thus it is important to consider both 2D and 3D ice in order to reach a consistent picture of vdW inclusive XC functionals. We also rationalize observations in previous FF studies by showing that widely used FF models incorrectly stabilize the high density 2D ice phases. and A.M.'s work is also sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant number FA8655-12-1-2099. A.M. is also supported by the Royal Society through a Royal Society Wolfson Research Merit Award. J.G.B acknowledges support by the Alexander von Humboldt foundation within the Feodor-Lynen program. We are also grateful for computational resources to ARCHER, UKCP consortium (EP/ F036884/1), the London Centre for Nanotechnology, UCL Research Computing, and Oak Ridge Leadership Computing Facility (No. DE-AC05-00OR22725).

SUPPLEMENTAL MATERIAL
The material here reported is the following: computational details in Sec. S1; extended data in Sec. S2; additional notes in Sec. S3; and structure files used for DMC calculations are described in Sec S4.

S1. ADDITIONAL DETAILS OF COMPUTATION
A. Evaluated quantities: binding energy, confinement energy and enthalpy The binding energy E b (sometimes called cohesive energy) of a given ice structure is defined as where E H2O is the energy of a water molecule in vacuum and E ice is the energy per water in the corresponding 2D ice structure. The latter is calculated as where E tot,ice is the total energy and n H2O is the number of water molecules in the simulated supercell of 2D ice.
In all our calculations the graphene layers are not simulated with an ab-initio approach, because that would be extremely difficult as the unit cells of graphene and of the 2D ice structures are incommensurable, so huge supercells should be considered (out of the reach for abinitio approaches). However, in our previous paper we did an extensive set of tests which showed that the uniform confining potential and explicit graphene confinement led to similar structures. 7 Two different types of confinement potentials were used to evaluate the confinement energy and to optimize the 2D ice structures. The first one is a Morse potential fitted to the quantum Monte Carlo (QMC) results for the interaction of a water monomer with graphene, 7,34 and is defined as where z is the distance between the oxygen atom and the wall, second one is a Lennard-Jones 9-3 (LJ93) potential, and is defined as where = 21.7 meV and σ = 3.0Å. Assume that the 2D ice is in the xy-plane and that the two walls are at z = 0 and z = w (w is the width of the confinement), such that the waters are found in the region 0 < z < w. At 0 pressure, in order to benchmark force field models, we report results with the LJ93 confining potential with w = 5Å, which is also available in force field programs and has been widely used. As shown in our previous paper, 7 the 5Å LJ93 confinement leads to very similar results to the Morse potential. At 2 GPa, as it is relevant to the experiments in graphene confinement, we have used the Morse potential as the confinement to achieve a better description. At 2 GPa the width of confinement used is 6Å. The impact of the width on the stability of 2D ice has also been discussed in our previous paper. 7 The corresponding confinement energy E conf is defined as: where the index i runs over all the water molecules in the simulated supercell, z i is the z coordinate of the i-th water, and the potential V (·) can either be V Morse or V LJ93 . Note that in both cases E conf is only a function of the structure, and in particular of the distance between the oxygen atoms and the two walls that realize the confinement.
The enthalpy is defined as, where P is the lateral pressure and A is the lateral area per water. Zero point energy is calculated as, and vibrational free energy at finite temperature is calculated as, where ω Γ,i is the ith gamma point vibrational frequency calculated using finite displacement method, and is the reduced Planck constant.

B. DMC
For all the DMC calculations, we used trial wavefunctions of Slater-Jastrow type with single Slater determinants and a Jastrow factor that includes electronnucleus, electron-electron and electron-electron-nucleus interaction terms. The single particle orbitals are obtained from DFT-LDA plane-wave calculations, performed with pwscf, 49 using a plane-wave cutoff of 600 Ry and re-expanded in terms of B-splines 33 with the natural grid spacing a = π/Gmax, where Gmax is the modulus of the longest G vector in the PW expansion. Dirac-Fock pseudopotentials 30,31 and the locality approximation 32 are used. With this wave-function ansatz, the variance of the VMC local energies per molecule is of ∼ 0.3 Ha 2 in all the considered systems, included the isolated molecule. 2D periodicity is used in all calculations of 2D ice structures, whereas for the single molecule open boundary conditions are used.
Based on previous work, we know that one of the most challenging parameters to choose is the DMC time-step τ ; exact results are obtained for τ → 0, but the computational cost is ∝ 1/τ . In previous DMC investigations on water systems E b was evaluated with a τ as small as 0.002 a.u. in order to get rid of undesired time-step errors (see S.I. of Ref. 10). Recently,Zen et al. 26 have introduced some improvements to the DMC algorithm that drastically reduces finite time-step errors. Thus, in this work we used the prescriptions of Ref. 26 and a τ = 0.02 a.u., which guaranties a negligible time-step error. In order to show this, we considered the square 2D ice structure (see Fig. 1 in main paper) to evaluate the extent of time-step errors. The unit cell includes four water molecules, and here we consider a 4 × 4 supercell, for a total of 64 waters in the system. At the VMC level of theory the evaluated binding energy is -0.108(4) eV, that is severely underestimated (by a factor 4) with respect to the DMC evaluations. In Fig. S1 we display the energy of the isolated water molecule, as well as the energy per water in the square structure, as a function of the DMC time-step τ . No sizeable time-step errors are observed in E b for τ ≤ 0.05, and the inset shows that the absolute energies are almost converged for τ ≤ 0.03 a.u. A value of τ = 0.05 would be enough for having no sizeable finite-time-step bias on the E b evaluations on the selected system, however in this work we are going to consider also larger 2D ice structures (up to 192 waters for supercell). Thus, we took for the production calculations a safe value of τ = 0.02 a.u., that is 10 times more efficient than the values that was considered for instance inSantra et al. 10 with the old DMC algorithm.
The second aspect that needs to be carefully checked in order to provide correct benchmark values are finitesize (FS) effects. In contrast to DFT, which is an effective one-particle method that can exploit Bloch's theorem in calculations for extended periodic systems, DMC is a many-body method and the finite size errors can be sizeable unless large supercells and correction terms are considered. [38][39][40] To eliminate the finite size errors, we performed several calculations for increasing supercell sizes up to ca. 200 water molecules and the binding energy is extrapolated to the infinite size, as shown in Fig. 2 of the main manuscript. Larger systems yield smaller FS error, but the computational effort increases with N w roughly linearly. The reason for having E b scaling as N w α , with α ∼ 1, is because in DMC the stochastic error σ is proportional to the inverse of the square root of the computational cost, and, given a target stochastic error σ on the total energy, an evaluation for a system with N electrons has a computational cost ∝ N γ , with γ ∼ 3. Each water has 8 valence electrons (we use pseudo potentials), so N = 8N w . Here, our target stochastic error is on the energy per water, thus we need an error ∝ N w σ, yielding α = γ − 2. However, the actual scaling is slightly worse than linear, because in DMC there is an  equilibration time that gets removed; in small systems this equilibration time is a negligible fraction of the total computation, but it scales as N γ so it becomes relevant as the system size increases. Thanks to the algorithmic improvements of Ref. 26, we can use large time-steps such that also in the huge systems here considered the equilibration time is not a bottleneck.
Here we further discuss the FS effects and corrections using the model periodic Coulomb interaction [50][51][52] or the Kwee, Zhang and Krakauer method. 40 Again, we performed our tests on the square 2D ice structure, this time keeping τ = 0.02 a.u. and considering different sizes of the supercell, namely by considering an increasing number N w of water molecules in the simulated supercell in order to observe the convergence of E b . Fig. S2 shows the results of these tests. Panel a leads to several interesting observations. First, we observe that the black asterisk in the plot, representing a twist average DMC evaluation, is superimposed with the purple square, that is the evaluation at Γ. Thus, the twist average technique is not necessary here, as expected since the system is an insulator and the smallest considered supercell is already pretty large (it comprises 16 water molecules); we do not use twist average for larger supercells. If the Kwee, Zhang and Krakauer 40 (KZK) corrections are added on top of the Ewald evaluations, we observe that the FS errors are overcorrected, probably because of the difference between these systems and the homogeneous electron gas where KZK corrections were tuned. Thus, bare Ewald evaluations are to be preferred here. We also tested the use of the model periodic Coulomb (MPC) interaction [50][51][52] in place of the Ewald one (actually, in the reported MPC calculation the DMC propagation is performed using Ewald, as recommended in CASINO's manual, because there is some evidence that the MPC interaction can distort the XC hole). MPC interaction mitigates the FS effects, although they are not completely removed and some sort of extrapolation to the infinite system is still necessary. We have tested different possible extrapolations for the bare DMC evaluations with Ewald interaction, and they are reported in Fig. S2b. It is not know a priori what the functional form of the FS bias is for a quasi-2D system as the one we are dealing now; for a bulk system the FS error is often expected to be ∝ N w −1 , and for strict 2D systemsCeperley 53 proposed that they are ∝ N w −3/2 , whereasDrummond et al. 54 proposed ∝ N w −5/4 . We are not able to determine which is the best fitting function because the error bars of the DMC evaluation are much larger than the residue of the fitting. This is also because our smallest system is already very large: one or more orders of magnitude larger than what has been considered in other studies where infinite size extrapolations have been performed. 54 Thus, we are much less affected here by the specific form chosen for the fitting function. Indeed, we obtain no significant different in the extrapolation of E b in the three methods; the maximum difference in 2.1 meV, that is small enough for the accuracy that we need for our evaluations. Thus, in the letter we took the easiest choice for the fitting, that is a FS error ∝ N w −1 , upon the bare DMC evaluations with Ewald 2D periodicity.
The structures used for DMC calculations are obtained from relaxation using the optPBE-vdW-DF functional; The hexagonal structure is not considered in the case of 2 GPa pressure because its area per water A is to much larger that for the other structures, thus the P V = P × A × w term in the entropy makes this structure unstable already at much smaller pressures.

C. DFT
DFT calculations were performed using the Vienna Ab-initio Simulation Package (VASP), using hard projector augmented wave (PAW) potentials and a 1000 eV plane wave cutoff. A periodic slab model was used containing a vacum region and the cell dimension perpendicular to the ice layer is 15Å. The cell dimension perpendicular to the slab was fixed during structure optimization, and the lateral cell dimensions were relaxed until the lateral diagonal elements of the stress tensor were converged to the target lateral pressures. K-points were sampled using the Monkhorst-Pack grids, with 4 × 4 × 1 for the hexagonal and pentagonal structures and 6×6×1 for the square and rhombic structures. Various DFT functionals of LDA, 55 GGA, meta-GGA and hybrid type were used. GGA functionals used include the Perdew-Burke-Ernzerhof (PBE), 56 revPBE 57 and Becke-Lee-Yang-Parr (BLYP). 58,59 SCAN is a meta-GGA functional. 42 Hybrid functionals considered are the PBE0 60 and HSE 61,62 functionals. Van der Waals corrections were accounted in several different schemes. The first is the D2 and D3 corrections of Grimme, 63,64 and superscript atm indicates a three-body dispersion term of Axilrod-Teller-Muto type. The second is the corrections of Tkatchenko and Scheffler (+vdW(TS)). 65 The self consitent screening (scs) and many body dispersion (MBD) 66 were also considered with the TS scheme. The third contains the van der Waals functionals of Langreth and Lundqvist, including PBE-vdW-DF, optPBE-vdW-DF, revPBE-vdW-DF, optB88-vdW-DF, optB86b-vdW-DF 35,67-69 and the dW-DF2 functional (rPW86-vdW-DF2). 70

D. Low-cost electronic structure methods
Three low-cost electronic structure methods has been considered. A hybrid functional evaluated in a double-ζ type basis set HSE-3c, 71 a minimal basis set Hartree-Fock method HF-3c, 72 and a density functional tight-binding Hamiltonian with re-parametrized dispersion correction DFTB3-D3. 73 The setups can be found in Ref. 74 48 were carried out using the LAMMPS code. 77 A 10 × 10 × 1 supercell was used and periodic boundary conditions were applied in the two lateral dimensions. The long-range Coulombic interaction is calculated using a particle-particle particle-mesh (PPPM) algorithm with an accuracy of 10 −4 . Energy minimization was performed with fixed cell dimensions and the lateral cell dimensions were tuned manually until the total energy is converged to within 1 meV/H 2 O.

S2. EXTENDED DATA
Benchmark on the high pressure structures is shown in Fig. S3. The effect of geometry relaxation on the enthalpy is shown in Fig. S4. Binding energy values calculated with different methods are reported in Table S1. A comparison between the integrated number of neighbors as a function of the O-O distance for the 2D ice structures and the bulk ice (Ih, II, VIII) is shown in Fig. S5, and in Fig. S6 we show the contribution in the interaction given by waters within a given distance.  is the free energy difference where the vibrational free energy at 300 K is considered.

S3. ADDITIONAL NOTES
Recently, a hierarchy of the low-cost electronic structure methods has been developed in the Grimme group (HSE-3c, 71 HF-3c, 72 and DFTB3-D3 73 ). They allow speed-ups of one, two, and four orders of magnitudes compared to the DFT calculations. We find that 2D ice is still challenging and especially the tight-binding accuracy is not satisfactory. This is probably due to very approximate description of short-range exchange repulsion and induction based on monopole charges. However, the accuracy of these methods is reasonable comparing to many DFT methods and might be useful for screening applications.