Water's unique properties and complex polymorphism, exhibiting seventeen experimentally confirmed ice polymorphs and several theoretically predicted ones, make its phase diagram a challenging yet crucial area of research. While experimental and theoretical investigations have advanced our understanding, uncertainties remain regarding the completeness of the known stable phases and the precise characteristics of their coexistence curves. High-pressure experiments pose significant difficulties, highlighting the increasing importance of computational simulations. However, accurately computing the thermodynamic stabilities of various water phases is computationally demanding, requiring the consideration of quantum thermal fluctuations and, for proton-disordered phases, configurational entropy in free-energy calculations. Empirical potentials, while offering some insights, often entail approximations, such as neglecting nuclear quantum effects (NQEs) or failing to accurately describe high-pressure behavior. Ab initio methods offer a promising alternative but face significant computational cost. This study bridges this gap by combining machine learning potentials (MLPs) with ab initio calculations and advanced free-energy techniques to construct a more accurate phase diagram.

Previous studies have utilized empirical potentials such as TIPnP, SPC/E, and MB-pol to model water's phase behavior. These models, however, often suffer from limitations, particularly in representing NQEs and high-pressure conditions. While some models like TIP4P-type and iAMOEBA qualitatively reproduce the phase diagram, others fail to accurately capture the stability of ice Ih. The use of MLPs offers a way to reduce the computational cost of ab initio calculations. Earlier work has demonstrated the application of MLPs to study the influence of van der Waals corrections on liquid water properties and to reproduce thermodynamic properties of solid and liquid water at ambient pressure. Recently, similar frameworks have been applied to study the phase diagrams of other materials, such as gallium and hydrogen. The challenge of accurately representing proton disorder in ice phases has also been highlighted, as the choice of proton configuration significantly impacts calculated enthalpies and free energies.

This research employs a multi-step approach to compute the water phase diagram using three hybrid DFT functionals (revPBE0-D3, PBE0-D3, B3LYP-D3). The study considers 50 putative ice crystal structures, including all experimentally known ice phases. To overcome the computational cost of ab initio MD simulations, a previously developed MLP is used as a surrogate model during free-energy calculations. This MLP, trained on liquid water data using the revPBE0 functional, enables efficient exploration of the phase diagram for a wide range of ice structures. The workflow involves three thermodynamic integration steps (TI). The first step uses the MLP to compute the chemical potentials of each phase through classical free-energy calculations. This employs the TI method and a fluctuating box with a Parrinello-Rahman-like barostat to determine equilibrium lattice parameters and perform thermodynamic integration to the MLP potential. Large system sizes (hundreds to thousands of water molecules) minimize finite-size effects. The chemical potentials are then obtained at other pressures and temperatures by numerically integrating the Gibbs-Duhem and Gibbs-Helmholtz relations. To determine the chemical potential of the liquid phase, direct-coexistence simulations are conducted. The second step involves a free-energy perturbation to promote the MLP results to the DFT level, correcting for residual errors, particularly concerning long-range electrostatics not fully captured by the MLP. The equation used is Δμ_MLP-DFT(P, T) = −k_BT ln <exp(−(U_DFT − U_MLP)/k_BT)>_P,T,HML. The third step incorporates NQEs through path-integral molecular dynamics (PIMD) simulations, integrating the quantum centroid virial kinetic energy with respect to the fictitious ‘atomic’ mass to account for quantum effects. This utilizes PIMD simulations for various scaled masses and numerical integration. Error estimations are performed for each step, considering contributions from TI, MLP-DFT corrections, and NQEs, separately analyzing uncertainties for proton-ordered and disordered phases. The phase diagram is then constructed by analyzing the computed chemical potentials.

The chemical potentials calculated using the MLP initially identified 12 ice phases (Ih, Ic, II, III, V, VI, VII, IX, XI, XII, XIII, XV) and liquid water as potentially stable. All are experimentally known phases, suggesting the experimental phase diagram may be complete within the studied region. The treatment of proton disorder is addressed by averaging chemical potentials from multiple configurations for each disordered phase. The MLP results are refined by promoting them to the DFT level and incorporating NQEs. This correction significantly impacts the relative stability of proton-ordered and disordered phases. Analysis of the chemical potentials yields phase diagrams, demonstrating improved agreement with experimental data when corrections for DFT and NQEs are included. The diagrams show the effects of different DFT functionals, indicating that even small changes in relative chemical potential can drastically affect the phase boundaries. Ice III's stability is particularly sensitive to the specific proton disorder configurations and calculation uncertainties. The study also compares phase diagrams with and without accounting for NQEs, showing that NQEs stabilize higher-density phases. The effect of NQEs is quantified by comparing H2O, D2O, and a classical model (no NQEs), revealing significant differences at higher pressures. The melting point of ice I is about 4K lower for H2O than for D2O, a finding consistent with previous simulations.

The computed phase diagrams show good, although not perfect, agreement with experimental data. Discrepancies may arise from the sensitivity to small changes in chemical potentials, particularly for proton-disordered phases, and limitations of hybrid DFT functionals in accurately modeling water. The very close free energies of many phases highlight their metastability and potential for experimental observation under specific conditions. The study demonstrates the significance of considering NQEs, particularly in differentiating the phase behavior of H2O and D2O at higher pressures. The MLP, even without DFT corrections, provides a remarkably good description of the phase diagram, but careful interpretation is needed due to potential overestimation of proton-ordered phase stability. The work establishes a rigorous benchmark for quantum-mechanical methods and offers a pathway for future functional optimization based on accurate phase diagram predictions.

This study presents a comprehensive computational exploration of the water phase diagram, achieving good agreement with experimental findings. The methodology, combining MLPs with ab initio calculations and accounting for NQEs and proton disorder, demonstrates the feasibility of predicting phase diagrams from first principles for complex systems. The results underscore the sensitivity of the phase diagram to even subtle changes in free energy, particularly for partially ordered phases. Future research should focus on exploring less-understood regions of the phase diagram, such as negative pressures and high-pressure transitions, and applying the methodology to other materials to improve interatomic potentials and assess the performance of diverse electronic structure methods.

The study's accuracy is inherently linked to the limitations of the hybrid DFT functionals used and the approximations inherent in the MLP approach. The uncertainty in the chemical potential calculations, particularly for proton-disordered phases, might affect the precise location of phase boundaries. The computational cost of the method still limits its applicability to very large systems or extremely long simulation times needed for studying certain processes (e.g., ice nucleation). Finally, the assumption of the experimental proton disorder as correct might introduce some inaccuracies in the free energy calculation for the partially ordered ice phases.