Quantum-mechanical exploration of the phase diagram of water

Water exhibits rich polymorphism with at least 17 experimentally confirmed ice phases in addition to liquid and vapor. Despite extensive study, uncertainties remain about whether all stable phases have been identified and about the precise locations of some coexistence lines, especially at high pressures where experiments are difficult. Predicting phase stability requires accurate free energies that include thermal and nuclear quantum fluctuations and account for proton disorder in many ice phases. Empirical water models have provided useful insights but involve strong approximations and often fail to reproduce the full qualitative phase behavior. Ab initio electronic-structure methods promise higher fidelity but are too costly for direct free-energy computations across many phases and thermodynamic conditions. This study aims to compute the phase diagram of water from first principles by combining machine-learning potentials, hybrid DFT, and advanced free-energy methods, explicitly treating nuclear quantum effects and proton disorder, and to assess whether any hypothetical ice phases become thermodynamically stable.

Prior work with empirical models (e.g., TIPnP, SPC/E, TIP4P/2005, TIP3P, TIP5P, iAMOEBA) has explored water’s phase behavior, with only a subset reproducing the qualitative phase diagram; some models even predict ice Ih stability only at negative pressures. The MB-pol force field, while accurate for many properties, has not been parameterized for high-pressure phases. Machine-learning potentials (MLPs) have recently enabled near-ab initio accuracy at reduced cost, elucidating effects of van der Waals interactions and reproducing several thermodynamic properties at ambient pressure. Multithermal–multibaric ML simulations have computed phase diagrams in metals (e.g., gallium) and suggested supercritical behavior in high-pressure hydrogen. Proton disorder and partial disorder in several ice pairs (I/XI, III/IX, V/XIII, VI/XV, VII/VIII, XII/XIV) introduce configurational entropy and sampling challenges. The need to incorporate nuclear quantum effects (NQEs) has been recognized, with path-integral methods revealing isotope effects (H2O vs D2O) on melting and phase stability. The accuracy and limitations of DFT for water, including the role of dispersion corrections and hybrid functionals, are well documented, motivating benchmarking across functionals.

• Candidate structures: From 15,859 hypothetical ice structures, 57 phases were pre-screened using a generalized convex hull. Defected and dynamically unstable phases and ice X were removed; ice IV was added, yielding 50 candidate ice structures (including all experimentally known phases) for simulation.
• DFT setup: Energies were computed with CP2K using three hybrid functionals with dispersion corrections: revPBE0-D3, PBE0-D3, and B3LYP-D3. Input settings followed established protocols with these functionals.
• Free-energy workflow: A three-step thermodynamic integration (TI) and perturbation scheme was employed: (i) reference harmonic crystal → fully anharmonic classical system described by a machine-learning potential (MLP); (ii) free-energy perturbation from MLP to DFT; (iii) inclusion of NQEs via path-integral methods by integrating the centroid virial kinetic energy with respect to a fictitious mass. • Step (i) Classical MLP free energies: For each ice phase and state point (T=25–300 K, P=0–10,000 bar), equilibrated NPT simulations with a Parrinello–Rahman barostat determined lattice parameters. A perfect crystal was prepared; its harmonic Helmholtz free energy was computed from the Hessian eigenvalues (center-of-mass motion corrected). TI yielded the fully anharmonic classical chemical potential at the MLP level. Finite-size effects were minimized by using large cells (hundreds to thousands of molecules). Chemical potentials at other (P, T) were obtained via numerical integration of the Gibbs–Duhem relation along isotherms and the Gibbs–Helmholtz analogue along isobars. • Liquid reference: Direct-coexistence simulations between ice I and liquid at 0 bar determined the MLP melting point (Tm ≈ 279.5 K). Equality of chemical potentials at coexistence provided the liquid reference, which was propagated across (P, T) via TI. • Proton disorder: For fully and partially disordered phases (e.g., I, III, V, VI), 5–8 distinct depolarized proton-disorder realizations were generated (Buch algorithm and GenIce). For partially disordered III and V, configurations matching experimental site occupancies were selected. For each realization, chemical potentials were computed independently and averaged; configurational entropy for (partial) disorder was added. • Step (ii) MLP → DFT: Using configurations sampled with the MLP (56–96 molecules per cell), DFT single-point energies were evaluated, and free-energy perturbation computed: Δμ_MLP-DFT = −kBT ln ⟨exp[−(U_DFT − U_MLP)/kBT]⟩. Averages were taken over disorder realizations. This correction accounts for long-range electrostatics and other residual MLP errors. • Step (iii) Nuclear quantum effects: PIMD with 24 beads evaluated centroid virial kinetic energies at scaled masses y = m̃/m = 1/4, 1/2, √2/2, 1. The NQE correction was obtained by numerically integrating with respect to y. H2O and D2O phase behaviors were computed by integrating to hydrogen or deuterium masses, respectively.
• Uncertainty estimation: Repeated reference calculations (5–10 times) at 125 K, 0 bar provided standard deviations. For disordered phases, uncertainty was decomposed into TI sampling versus disorder realization components. Additional uncertainties came from MLP→DFT perturbation (evaluated at 225 K, 3000 bar) and from NQE integration (propagated from kinetic energy statistical errors via numerical quadrature).

• Phase stability set: Starting from 50 candidate ices, only experimentally known phases (Ih/I, II, III, V, VI, VII, IX, XI, XII, XIII, XIV, XV) and the liquid were competitive; none of the additional hypothetical phases became thermodynamically stable in the explored range (T=25–300 K, P=0–10,000 bar), suggesting completeness of the experimental phase diagram in this region.
• Agreement with experiment: The ab initio phase diagram (MLP corrected to DFT with NQEs) qualitatively matches experiment, especially for P ≤ 8000 bar. The I–liquid coexistence line slope is captured; higher-pressure regions show small discrepancies (e.g., overly steep V–VI boundary).
• Sensitivity to proton disorder: Many phases are within a few meV in chemical potential. Typical per-phase uncertainty is ≤ ~2.5 meV, which can shift coexistence lines noticeably. Ice III stability is highly sensitive: it lies within ~2 meV of the stable phase near ~250 K and ~3000 bar; small shifts within estimated uncertainties render ice III stable over the experimentally correct region.
• MLP→DFT corrections: Free-energy perturbation corrections to chemical potentials range from ~7.9 meV (e.g., VI at 125 K, 10,000 bar) to ~10.8 meV (e.g., XV at 275 K, 5000 bar). Proton-ordered phases are over-stabilized by the MLP; DFT corrections reduce their stability. Differences in corrections between ordered–disordered pairs are roughly constant across relevant conditions: ΔΔμ(V − XIII) ≈ 5 meV; ΔΔμ(VI − XV) ≈ 11 meV.
• Nuclear quantum effects: NQEs stabilize denser phases (VI > V > II > liquid > I in stabilization order). Relative to ice II, NQE corrections to chemical potentials are on the order of ~5–8 meV in magnitude. Including NQEs significantly shifts high-pressure boundaries and yields distinct H2O vs D2O phase behavior; the H2O melting point of ice I is ~4 K lower than D2O, and D2O’s I–liquid line closely follows the classical (no NQEs) line.
• Melting points at 1 bar: revPBE0-D3 predicts a melting point near 274 K with an uncertainty ~6 K (consistent with a corrected 274(2) K from prior work at the same level). For alternative functionals, B3LYP-D3 gives 274(6) K and PBE0-D3 gives 268(6) K.
• Overall: The combined ML–DFT–PIMD workflow enables first-principles prediction of water’s phase diagram, providing a stringent thermodynamic benchmark for electronic-structure methods.

The study addresses whether ab initio methods, when coupled with machine-learning surrogates and rigorous free-energy techniques, can predict the complex phase behavior of water, including the role of proton disorder and NQEs. The results show strong qualitative agreement with experimental phase boundaries at low to moderate pressures and demonstrate that small differences in chemical potential (few meV) due to disorder or methodological choices can significantly affect phase stability, explaining some discrepancies. The approach highlights both the capability and limits of hybrid DFT for water: while multiple functionals yield similar phase diagrams, residual differences remain and can alter specific coexistence lines (e.g., V–VI). The explicit inclusion of NQEs is essential to capture isotope effects and high-pressure behavior. The findings validate the use of MLPs trained on liquid data for solid phases when corrected to DFT, but caution against relying on uncorrected MLP predictions, especially for proton-ordered phases. Overall, the workflow provides a robust, thermodynamics-based benchmark for assessing and improving electronic-structure methods and interatomic potentials.

This work delivers a first-principles phase diagram of water by integrating a machine-learning potential with hybrid DFT corrections and nuclear quantum effects via PIMD, explicitly treating proton disorder. It finds no new thermodynamically stable ice phases among 50 candidates in the studied range and reproduces key experimental features, particularly below ~8000 bar. The method demonstrates that accurate phase behavior can be predicted for a polymorphic system from quantum mechanics, and that phase diagrams provide a sensitive benchmark for electronic-structure accuracy. Future research directions include: extending the workflow to other electronic-structure methods (e.g., random-phase approximation, double-hybrid DFT) to benchmark functionals; exploring negative-pressure regions and the low-density clathrate regime; investigating the ice VII–liquid transition at higher pressures; applying the framework to nucleation thermodynamics and to assess (meta)stability of novel materials; and leveraging experimental phase diagrams to refine and correct interatomic potentials.

• Proton disorder sampling: Equilibration of proton disorder is infeasible on simulation timescales; results depend on the limited set of disorder realizations and assumed experimental site occupancies, introducing meV-level uncertainties that can shift phase boundaries.
• DFT accuracy: Even hybrid DFT has known limitations for water; functional choice (revPBE0-D3, PBE0-D3, B3LYP-D3) affects boundaries (e.g., V–VI slope). Finite basis sets and energy cutoffs may contribute additional errors.
• MLP biases: The MLP trained on liquid configurations tends to over-stabilize proton-ordered phases without DFT correction, necessitating perturbative corrections for quantitative accuracy.
• Finite-size and statistical errors: Despite large cells and repeated calculations, residual finite-size, sampling, and integration errors persist (typical chemical potential uncertainty up to ~2.5 meV), particularly for disordered phases.
• Scope limits: The study covers 25–300 K and 0–10,000 bar; behavior outside this range (e.g., very high pressures, negative pressures) and kinetics (nucleation, transformations) were not resolved here.