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.

Figure 1

The four most relevant 2D ice structures considered as part of this study: hexagonal (H), pentagonal (P), square (S), and rhombic (R). The gray boxes represent the primitive unit cells.

Figure 2

DMC binding energies $E_{\text{b}}$ for the zero pressure structures of free-standing 2D ice in the hexagonal (H), pentagonal (P), square (S), and rhombic (R) structures (see Fig. 1) as a function of the inverse number of water molecules in the simulation cell $1/N_w$. Open symbols correspond to DMC evaluations of $E_{\text{b}}$, and error bars are one standard deviation. For each point we indicate the supercell size in terms of the primitive unit cell. The solid symbols on the yellow region are the extrapolated values $E_{\infty }$ from a fit with the function $E\left(N_w\right)=E_{\infty }-c/N_w$; values and errors of the fit are reported on the top right corner.

Figure 3

The left panel reports the average error (blue bars) for the $E_b$ evaluations, taking DMC as the reference, using a selection of XC functionals. See the SM for information on each method. The right panel reports the energy difference $E_{\text{b}}-E_{\text{b}}^{\text{hex}}$ at zero pressure for each of the considered methods. The benchmark DMC evaluations are shown as color-coded bands, the widths of which represent the stochastic error (i.e., $\pm \sigma$). Color and symbol conventions are the same as in the previous figures. Results are on optPBE-vdW optimized structures at zero pressure. Similar plots for relaxed structures and for 2 GPa are reported in the SM.

Figure 4

Binding energy difference $E-E^{\mathrm{DMC}}$ with respect to DMC for various DFT XC functionals and FF models for the 2D and some bulk (3D) ice structures at zero pressure. For clarity, only a subset of the FF potentials and DFT XC functionals are plotted here; data with other approaches are included in the SM. Bulk ice results are taken from Refs. [11, 25] (DMC and DFT) and Ref. [42] (FF). In both panels the densities of the phases increase from left to right.