Accurate Thermophysical Properties of Water using Machine-Learned Potentials

Simulating water from first principles remains a significant computational challenge due to the slow dynamics of the underlying system. Although machine-learned interatomic potentials (MLPs) can accelerate these simulations, they often fail to achieve the required level of accuracy for reliable uncertainty quantification. In this study, we use MACE - an equivariant graph neural network architecture that has been trained using an extensive RPBE-D3 database - to predict density isobars, diffusion constants, radial distribution functions, and melting points. Although equivariant MACE models are computationally more expensive than simpler architectures, such as kernel-based potentials (KbPs), their significantly lower total energy errors allow for reliable thermodynamic reweighting with minimal bias. Our results are consistent with those of previous studies using KbPs; however, equivariant models can be validated against the ground-truth density functional theory (DFT) ensemble, providing a critical advantage. These findings establish equivariant MLPs as robust and reliable tools for investigating the thermophysical properties of water with DFT-level accuracy.

💡 Research Summary

This paper addresses the longstanding challenge of obtaining first‑principles accuracy for liquid water simulations while keeping computational costs tractable. The authors employ the equivariant graph neural network architecture MACE (Multi‑Atomic‑Coordinate‑Equivariant) trained on an extensive RPBE‑D3 density‑functional theory (DFT) database. Two variants of MACE are considered: a scalar (L = 0) model and a fully equivariant model with tensor ranks up to L = 2. The latter contains roughly 880 k trainable parameters, about twice as many as the scalar version, and is designed to respect rotational and permutation symmetries, thereby capturing long‑range electrostatic and dipolar interactions that kernel‑based potentials (KbPs) neglect.

The training set comprises 2 370 configurations of 63–64 water molecules, generated via active learning with a previous KbP and recomputed at the RPBE‑D3 level. During training, the authors discovered systematic force‑error growth when the MACE models sampled configurations not represented in the original KbP database. To remedy this, additional structures were generated from preliminary MACE ensembles and incorporated until the error plateaued. A residual mean stress shift (≈50 bar) was also identified and corrected post‑hoc, highlighting the sensitivity of equivariant models to subtle DFT artifacts.

Model evaluation on a hold‑out test set of 128‑molecule configurations shows that equivariant MACE (L = 2) achieves an energy root‑mean‑square error (RMSE) of 0.09 meV per atom and a force RMSE of 8.9 meV Å⁻¹, outperforming the KbP (0.2 meV/atom, 27 meV Å⁻¹) by a factor of five. The scalar MACE (L = 0) performs worse than KbP for total energies on larger cells, indicating that inclusion of higher‑order tensor features is essential for accurate long‑range physics.

Thermodynamic reweighting (Thermodynamic Perturbation Theory) is a central theme. Because the MACE energy errors are sufficiently small, the Boltzmann distribution generated by the MACE‑driven molecular dynamics (MD) overlaps strongly with the true DFT distribution. Consequently, observables measured in the MACE ensemble can be reweighted to the DFT ensemble with negligible bias. In contrast, the larger errors of KbP lead to poor overlap and unreliable reweighting, limiting their utility for rigorous uncertainty quantification.

Using the MACE potentials, the authors compute several key water properties:

1. Density isobars: Parallel‑tempering MD with 128 molecules (≥10 ns per replica) yields density–temperature curves that agree closely with KbP results. The MACE‑predicted density maximum occurs at 288.8 ± 2.4 K, slightly higher than the KbP value (~280 K) but within statistical uncertainty. Finite‑size effects were examined by additional simulations of 1024 molecules, confirming that the 128‑molecule results are robust.

2. Diffusion coefficients and radial distribution functions (RDFs): Both MACE and KbP reproduce experimental diffusion trends, with MACE providing marginally better agreement in the first RDF peak position and height, reflecting its superior description of short‑range hydrogen‑bond networks.

3. Melting point: A solid–liquid interface simulation containing 800 molecules was performed across a range of temperatures. The internal energy discontinuity indicates melting at 283.75 ± 1.0 K for MACE, about 4 K lower than the KbP estimate but consistent with the higher temperature of the density maximum.

Computational performance is discussed in detail. MACE simulations are executed on NVIDIA RTX 6000 GPUs with CUDA equivariance acceleration, achieving roughly a four‑fold slowdown relative to CPU‑only KbP runs. Memory consumption is substantially higher, and the current implementation requires one GPU per replica, limiting the number of parallel tempering images that can be run simultaneously. Nevertheless, the authors argue that the gain in accuracy and the ability to perform reliable reweighting outweigh the increased resource demand, especially on modern GPU clusters.

The paper concludes that equivariant MACE potentials represent a significant step forward for water modeling. Their low energy and force errors enable DFT‑level thermodynamic predictions and rigorous reweighting, while still offering speedups of several orders of magnitude over direct DFT MD. The authors suggest that the approach is readily extensible to other systems where long‑range electrostatics and many‑body polarization are critical, such as electrolytes, ionic liquids, and complex interfaces. Future work will focus on improving GPU scalability, exploring larger cutoff radii, and integrating more diverse DFT functionals to broaden the applicability of equivariant machine‑learned potentials.