Water Phase Diagram from a General-Purpose Atomic Cluster Expansion Potential

TL;DR

The paper addresses predicting water's phase diagram using a data-driven interatomic potential. It trains an Atomic Cluster Expansion (ACE) potential on revPBE-D3 DFT data and uses OPES-enabled biased coexistence to obtain melting points, followed by Gibbs–Duhem integration to extend boundaries across pressure and temperature. The results reproduce the main ice polymorphs Ih, II, V, VI, VII and liquid water, with ice III metastable and systematic shifts relative to experiment, and are benchmarked against MB-pol, DeepMD-SCAN, and NNP-revPBE0-D3. This work demonstrates that ACE can accurately and efficiently capture complex phase behavior from first-principles data and provides a platform for future extensions to electrolytes and interfaces, employing the Clapeyron relation $dP/dT = Δh/(T Δv)$ to propagate coexistence lines.

Abstract

Water's phase diagram remains one of the most intricate and challenging benchmarks in molecular modeling. In this study, we compute the phase diagram of water using an Atomic Cluster Expansion (ACE) potential trained on density-functional theory (DFT) calculations based on the revPBE-D3 exchange and correlation functional. We compute solid-liquid chemical potential differences and melting points using biased coexistence simulations with the On-the-Fly Probability Enhanced Sampling (OPES) method. Starting from these points, we trace coexistence lines using Gibbs-Duhem integration. This combination of methods allows us to consistently map pressure-temperature phase boundaries and reconstruct the full phase diagram between approximately 100-500 K and 0-4 GPa. The stability regions of the main ice polymorphs (Ih, II, V, VI, and VII) are reproduced in close agreement with experiments. As in earlier studies based on DFT, ice III is metastable and there are systematic shifts of coexistence lines with respect to experimental results. Our results demonstrate the capability of our general-purpose ACE potential to capture the complex phase behavior of water across wide thermodynamic conditions.

Figures (6)

• Figure 1: Schematic overview of the computational workflow employed in this study, combining ab initio data, Atomic Cluster Expansion (ACE) training, biased coexistence simulations with OPES, and Gibbs–Duhem integration to construct the phase diagram of water.
• Figure 2: Representative OPES biased-coexistence simulations used to determine solid-liquid free-energy differences for ice Ih polymorph at 0.1 GPa. See text for details. (Top row) Time series of the collective variable $N_{\mathrm{ice}}$ at selected temperatures, showing repeated and reversible growth and melting events of an ice layer. (Middle row) Biased probability distributions of $N_{\mathrm{ice}}$. (Bottom row) Corresponding reweighted free-energy profiles $\Delta G(N_{\mathrm{ice}})$ together with linear fits used to extract the chemical-potential difference $\Delta\mu = \mu_{\mathrm{ice}} - \mu_{\mathrm{liq}}$.
• Figure 3: ACE potential training. Panels a and b show the convergence of the root-mean-square error (RMSE) of energies and forces, respectively, as a function of the number of basis functions for different number of optimization steps. Panels c and d show parity plots between the DFT reference data and the ACE predictions for energies and forces, respectively. Insets show the distribution of errors.
• Figure 4: Chemical-potential differences for ice polymorphs relative to the liquid as a function of temperature at fixed pressure. The horizontal dashed line marks $\Delta\mu = 0$; vertical dotted lines indicate the melting point.
• Figure 5: Gibbs--Duhem integration of phase coexistence lines. Each curve represents a solid--liquid or solid--solid coexistence boundary obtained by integrating from coexistence points obtained using OPES simulations. Colored lines correspond to coexistence lines between the phases indicated in the legend (L for liquid, and ice polymorphs Ih, II, III, V, VI, and VII), while stars denote triple points (TPs) identified from the intersections of coexistence lines. The inset highlights the region where we tested the robustness of our calculation by computing the Ih-II-V triple points using different routes.