Ab initio simulation of the first-order proton-ordering transition in water ice

Abstract

Proton ordering in water ice is a paradigmatic order-disorder transition in a locally constrained system. The ice rules require exactly two hydrogens close to each oxygen, restricting the disorder to an exponentially large yet strongly correlated manifold of hydrogen-bond configurations. Within this constrained space, meV-scale energy differences drive the transition from disordered ice Ih to ordered ice XI, while distinct configurations are separated by eV-scale barriers. These barriers hinder equilibration in experiments, and efficient sampling of this space with the required energy accuracy has remained a long-standing challenge in simulation. We address this by combining a machine learning interatomic potential with loop updates that preserve the ice rules and continuous updates of atomic coordinates, enabling equilibrium sampling with ab initio accuracy and capturing configurational entropic effects. In systems of up to 360 water molecules, with over 10^6 samples retained per temperature point, the simulations reveal clear first-order transition signatures at 83 K: a negative Binder cumulant, a bimodal potential energy distribution, and a sharp step in the lattice aspect ratio. Nuclear quantum effects are estimated to lower the transition temperature by approximately 20 K, bringing the prediction closer to the experimental value of 72 K.

Figures (9)

• Figure 1: Relative potential energies of hydrogen-bond configurations.a, Potential energies of all $16$ symmetry-inequivalent hydrogen-bond configurations in a cell containing $N=8$ water molecules, computed with DFT-$\mathrm{r^2SCAN}$ and the MACE potential. Energies are referenced to the $Cmc2_1$ configuration, $\Delta E = E_i - E_1$. b, Potential energy distributions at different temperatures for an $N=128$ cell, comparing the proton-ordered $Cmc2_1$ configuration with a random proton-disordered configuration. Atomic structures were visualized using VESTA software momma2011vesta with oxygen in red and hydrogen in pale pink.
• Figure 2: Thermodynamic and structural results of the proton ordering transition.a, Binder cumulant $B$ of the total dipole moment magnitude $M$ for different system sizes $N$ (number of $\mathrm{H_2O}$ molecules). The inset shows the temperature at which $B$ reaches its minimum for each $N$, which is extrapolated to the thermodynamic limit. b, Polarization density magnitude $\langle P \rangle$. c, Orthorhombic lattice anisotropy ratio $\langle b/a\rangle$ of the unit cell. The dashed line marks the ideal hexagonal close-packed value $\sqrt{3}$. d, Potential energy per molecule $\langle E \rangle$ with the classical harmonic contributions $\frac{9}{2}k_{\mathrm{B}}T$ (dashed line) subtracted, isolating the anharmonic and hydrogen-bond configurational contributions. e, Molar heat capacity $C_{\mathcal{P}}$ estimated from the enthalpy fluctuations excluding kinetic contributions, including vibrational and hydrogen-bond configurational contributions.
• Figure 3: Potential energy and polarization distributions with representative configurations. (a) Two-dimensional probability density of polarization magnitude $P$ versus potential energy per molecule $E$ from $10^6$ Monte Carlo samples at $65$, $85$, and $105$ K for an $N=360$ supercell. Black markers highlight two selected configurations, labeled $C^{\bigtriangleup}$ and $C^{\bigtriangledown}$. (b) Structural snapshots corresponding to the marked points. Oxygen (dark teal) and hydrogen (light teal) atoms in the winding loop are highlighted. A loop update flips the hydrogen atoms along the loop, transforming $C^{\bigtriangleup}$ into $C^{\bigtriangledown}$. (c) Probability density function (PDF) of potential energy $E$ at the three temperatures. (d) Monte Carlo trajectories at $65$ K for polarization magnitude $P$ (left axis) and potential energy $E$ (right axis). Dashed vertical lines indicate configurations $C^{\bigtriangleup}$ and $C^{\bigtriangledown}$.
• Figure S1: Comparison between the MACE model and DFT-r$^{2}$SCAN. Results on the test set for (a) energies per water molecule and (b) forces. Dots denote individual data points in the test set, and the dashed red line marks perfect agreement. Insets show the distributions of prediction errors.
• Figure S2: Sampling statistics of the Monte Carlo simulations.a, Number of recorded post-thermalization samples as a function of temperature. b, Total number of Monte Carlo move proposals. c, Ratio of the number of unique visited hydrogen-bond configurations to the total number of recorded samples. The inset shows the same quantity on a logarithmic scale.