A 2026 benchmark study (O'Neill et al., arXiv:2605.28798) showed that published first-principles simulations of ambient liquid water, all using the same density functional revPBE-D3, disagree by more than 20% in the self-diffusion coefficient and about 10% in density. By training machine-learning interatomic potentials on tightly converged DFT data from six independent community codes, the study traced the scatter to basis-set incompleteness, pseudopotential inconsistencies, and insufficient statistical sampling, and established consensus reference values (density 0.930±0.017 g cm−30.930\pm 0.017\ \mathrm{g\,cm^{-3}}, self-diffusion (2.75±0.15)×10−5 cm2 s−1(2.75\pm 0.15)\times 10^{-5}\ \mathrm{cm^{2}\,s^{-1}} with classical nuclei). It was restricted to a single generalized-gradient functional, and its authors expect the identified numerical-error mechanisms to be even more severe for hybrids, double hybrids, and wavefunction-based theories. The open problem is to extend this reproducible-benchmark protocol to exchange-correlation functionals on higher rungs of Jacob's ladder for liquid water: determining the converged numerical settings each code implementation needs, code-independent reference values for the density and self-diffusion coefficient, and whether the identified error mechanisms account for each functional's literature scatter. Adjacent 2026 studies (a single-pipeline convergence reassessment by Beck and Marsalek; an 11-functional accuracy benchmark by Dào et al.; a single-pipeline dynamics comparison across functionals by Montero de Hijes et al.) address neighboring questions but leave this cross-code extension open.

Background

Density functional theory (DFT) is the standard first-principles method for simulating condensed phases: the interacting-electron problem is recast in terms of the electron density, and all many-body effects are collected into an exchange-correlation (XC) functional that must be approximated. Approximate XC functionals are organized on "Jacob's ladder": rung 1 uses only the local density (LDA); rung 2 adds the density gradient (generalized-gradient approximations, GGAs, such as PBE and revPBE, a revised PBE); rung 3 adds the kinetic-energy density (meta-GGAs, such as SCAN); rung 4 mixes a fraction of exact Hartree-Fock exchange into a lower-rung functional (hybrids, such as PBE0, B3LYP, and the range-separated ωB97X-rV); rung 5 adds wavefunction-based correlation (double hybrids, or the random-phase approximation, RPA). GGA-level water simulations are usually augmented with the D3 dispersion correction, an empirical damped atom-pair attraction available in subtly different variants (zero-damping vs. Becke-Johnson damping, with or without the three-body Axilrod-Teller-Muto term). In ab initio molecular dynamics (AIMD), interatomic forces come from DFT at every time step, so liquid water can be simulated without empirical force fields. But the DFT equations are solved numerically, and each major community code represents the orbitals differently — plane waves with pseudopotentials or the projector-augmented-wave (PAW) method in VASP, Quantum ESPRESSO, and CASTEP; Gaussian-type orbitals with Goedecker-Teter-Hutter pseudopotentials in CP2K; all-electron Gaussian orbitals in PySCF; numeric atom-centered orbitals in FHI-aims. A pseudopotential replaces the inert core electrons with an effective potential; its choice, together with the incompleteness of the finite orbital basis, produces code-specific numerical errors on top of the XC approximation itself. For solid-state properties at 0 K, community verification campaigns — the Δ-code project (Lejaeghere et al.), the precision-verification workflows of Bosoni et al., and beyond-DFT efforts such as the fixed-node diffusion Monte Carlo code comparison of Della Pia et al. — established that independent codes can agree to high accuracy. Liquid properties had no comparable treatment: AIMD needs 10510^5-10610^6 force evaluations, which discourages tight numerical settings, and thermodynamic observables additionally demand statistical convergence of the sampling. O'Neill et al. closed this gap for one case. Surveying the literature on ambient liquid water with revPBE-D3, they found a spread of 0.860–0.991 g cm−30.860\text{--}0.991\ \mathrm{g\,cm^{-3}} in density and (1.82–2.31)×10−5 cm2 s−1(1.82\text{--}2.31)\times 10^{-5}\ \mathrm{cm^{2}\,s^{-1}} in self-diffusion — variations comparable to the differences between functionals, so numerical error can masquerade as functional error. They then trained long-range machine-learning interatomic potentials (MLIPs; here MACE message-passing graph neural networks augmented with learned local split charges for electrostatics) on tightly converged DFT data from six codes (CASTEP, CP2K, FHI-aims, PySCF, Quantum ESPRESSO, VASP), enabling exhaustive sampling, and achieved cross-code consensus: density 0.930±0.017 g cm−30.930\pm 0.017\ \mathrm{g\,cm^{-3}} and diffusion (2.75±0.15)×10−5 cm2 s−1(2.75\pm 0.15)\times 10^{-5}\ \mathrm{cm^{2}\,s^{-1}} with classical nuclei in a 64-molecule cell, the diffusion corrected for finite-size effects by the Yeh-Hummer scheme (D(L)=D∞−ζ kBT/(6πηL)D(L)=D_\infty-\zeta\,k_BT/(6\pi\eta L) with ζ=2.837297\zeta=2.837297 for a cubic cell and the experimental shear viscosity η\eta). Path-integral molecular dynamics (PIMD, a technique that makes nuclei quantum by representing each as a ring polymer; here with 32 beads) gives density 0.927 g cm−30.927\ \mathrm{g\,cm^{-3}} and diffusion (2.549±0.07)×10−5 cm2 s−1(2.549\pm 0.07)\times 10^{-5}\ \mathrm{cm^{2}\,s^{-1}}, closer to the experimental 2.41×10−5 cm2 s−12.41\times 10^{-5}\ \mathrm{cm^{2}\,s^{-1}}. They attributed the literature scatter to basis-set incompleteness, pseudopotential inconsistencies, and insufficient sampling. Three adjacent 2026 studies frame what remains. Beck and Marsalek used MLIPs to compare a popular, pragmatic revPBE0-D3 setup against a highly converged one, finding its widely reported experimental agreement degrades under convergence, while converged ωB97X-rV improves agreement and a common MP2 (second-order Møller-Plesset) triple-ζ setup is under-converged — but within a single computational pipeline. Dào et al. benchmarked 11 functionals spanning the top four rungs of Jacob's ladder for ambient water with MLIPs, mapping accuracy trends rather than numerical reproducibility. Montero de Hijes et al. trained MACE models for SCAN, RPBE-D3, revPBE-D3, revPBE0-D3, and two PBE0-D3 variants on previously reported ab initio datasets and compared dynamical properties across functionals — again within one pipeline. Hilpert and Kresse showed that MLIP architecture itself matters, with equivariant MACE models validating reliably against the DFT ensemble for RPBE-D3 water where simpler kernel-based potentials may not. What no study has done is establish whether the cross-code consensus protocol extends to functionals beyond revPBE-D3, where the source paper itself predicts the numerical-error mechanisms to be "even more severe".

Problem Statement

Extend the tightly converged, MLIP-accelerated, multi-code reproducibility protocol established for revPBE-D3 liquid water to exchange-correlation functionals on higher rungs of Jacob's ladder, for liquid water at ambient conditions with classical nuclei. For a representative set of functionals beyond GGA — including at least one meta-GGA (e.g. SCAN), at least one global or range-separated hybrid (e.g. PBE0, revPBE0, B3LYP, or ωB97X-rV), and at least one beyond-hybrid method (a double hybrid or the random-phase approximation), each with its dispersion treatment, where applicable, fixed across codes and stated explicitly for that functional (its standard choice where one exists): (i) determine, for each DFT code implementation used, the numerical settings — basis-set or plane-wave-cutoff convergence and pseudopotential/PAW consistency — required to reach quantitatively converged energies and forces; (ii) establish reference values of the equilibrium density and the self-diffusion coefficient, with quantified statistical and cross-code uncertainties, demonstrating agreement across independent code implementations (the revPBE-D3 study achieved consensus across six); and (iii) determine, for each functional studied, whether the mechanisms identified for revPBE-D3 — basis-set incompleteness, pseudopotential inconsistency, and insufficient statistical sampling — account for the scatter among that functional's published literature values, or whether additional error mechanisms appear at higher rungs. The single scientific objective uniting these deliverables is a numerically trustworthy, code-independent reference layer for the functionals used to judge the accuracy of first-principles simulation of liquid water, so that functional rankings rest on converged comparisons rather than numerical artifacts.

References

  1. Niamh ONeill, Benjamin X. Shi, William J. Baldwin, Albert P. Bartok, Chris J. Pickard, Angelos Michaelides, Gabor Csanyi, Timothy C. Berkelbach, "How reproducible are first-principles simulations of liquid water?", arXiv:2605.28798 (2026), DOI: 10.48550/arxiv.2605.28798, https://arxiv.org/abs/2605.28798
  2. Hubert Beck, Ondrej Marsalek, "More Converged, Less Accurate? Reassessing Standard Choices for Ab Initio Water Using Machine Learning Potentials", The Journal of Physical Chemistry B 130, 7215-7226 (2026), DOI: 10.1021/acs.jpcb.6c02917 (preprint arXiv:2603.20451, https://arxiv.org/abs/2603.20451), https://doi.org/10.1021/acs.jpcb.6c02917
  3. Nam H. Dào, Daniel S. King, Frederick Stein, Xiaoyu Wang, Dongjin Kim, Gerit Brandenburg, Jürg Hutter, Bingqing Cheng, "Systematic Trends in Water Properties Across Jacob's Ladder Density Functionals", ChemRxiv preprint (2026), DOI: 10.26434/chemrxiv.15000644/v1, https://doi.org/10.26434/chemrxiv.15000644/v1
  4. Tobias Hilpert, Georg Kresse, "Accurate thermophysical properties of water using machine-learned potentials", arXiv:2601.21103 (2026), DOI: 10.48550/arXiv.2601.21103, https://arxiv.org/abs/2601.21103
  5. P. Montero de Hijes, L. Neubeck, G. Kresse, C. Dellago, "Dynamical properties of ab initio water from machine-learning potentials", arXiv:2606.18163 (2026), DOI: 10.48550/arXiv.2606.18163, https://arxiv.org/abs/2606.18163
  6. K. Lejaeghere et al. (69 authors), "Reproducibility in density functional theory calculations of solids", Science 351, aad3000 (2016), DOI: 10.1126/science.aad3000, https://doi.org/10.1126/science.aad3000
  7. E. Bosoni et al. (45 authors), "How to verify the precision of density-functional-theory implementations via reproducible and universal workflows", Nature Reviews Physics 6, 45-58 (2024), DOI: 10.1038/s42254-023-00655-3, https://doi.org/10.1038/s42254-023-00655-3
  8. I. Yeh, G. Hummer, "System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions", The Journal of Physical Chemistry B 108, 15873-15879 (2004), DOI: 10.1021/jp0477147, https://doi.org/10.1021/jp0477147
  9. F. Della Pia et al. (33 authors), "Reproducibility of fixed-node diffusion Monte Carlo across diverse community codes: The case of water-methane dimer", The Journal of Chemical Physics 163, 104110 (2025), DOI: 10.1063/5.0272974, https://doi.org/10.1063/5.0272974