Water's ionization constant stops relying on error cancellation: quantum nuclear corrections computed per dissociation, phosphate's six steps land within 0.45

Acid dissociation constants can now be computed directly from atomistic simulation without relying on experimental data for similar compounds, and the same calculation yields a microscopic picture of the dissociation and solvation mechanism. Such calculations previously either predicted constants from molecular structure with data-driven models, or read the constant off a free energy profile with estimators that differ by several units on the same underlying profile. In a box of 128 water molecules, the reversible-work theorem gives a classical value of 17.86, while taking the difference between the free energy minima of the associated and dissociated states gives 22.06, a gap of 4.2 units.

The dissociation constant comes from integrating the potential of mean force between the products, whose asymptotic normalization is fixed by matching to the Coulomb interaction of point charges in a periodic box of the same size. A simulation delivers the free energy profile only up to an additive constant, and the long-ranged Coulomb attraction between the ions keeps the profile from reaching a plateau at any distance accessible in a periodic box of a few hundred molecules. In a box of 128 water molecules, the thermodynamic-limit and periodic-box reference curves differ visibly over the accessible range of separations, showing that finite-size effects are not negligible at this system size.

The nuclear quantum correction is obtained by thermodynamic integration over nuclear mass, which requires converging only the kinetic energy in the associated and dissociated end states rather than a quantum free energy surface along the whole dissociation path. Running enhanced sampling with path integral molecular dynamics multiplies the cost of every step by the number of ring-polymer replicas, and the cost of a quantum free energy surface has been reported to be more than an order of magnitude above the corresponding classical calculation. For water autoionization the correction is 3.54 units, lowering the classical 17.86 to 14.31 against the experimental 13.93, and the integrand can be integrated with five path integral simulations.

The same workflow moves from water to the entire dissociation ladder of phosphoric acid without changing the order-parameter definitions or the analysis, adapting only the metadynamics parameters to the magnitude of the dissociation free energy. Previous studies extracted the dissociation constant and treated nuclear quantum effects inconsistently, making values hard to compare across studies. The six classical constants of phosphoric acid are systematically too large by 2.19 units on average, and after the quantum correction the mean absolute error against experiment drops to 0.45 units, within what is conventionally taken as chemical accuracy.

The workflow is demonstrated in boxes of 128 water molecules with sodium counterions, covering water autoionization and the six dissociations of phosphoric acid, with sampling performed in the isothermal-isobaric ensemble so that the result is a Gibbs free energy of dissociation comparable with experiment. The authors note that moving from water to the phosphate species required no change in the order-parameter definitions or the analysis, only metadynamics parameters adapted to the magnitude of the dissociation free energy, and they expect the workflow to extend to polyelectrolyte simulation, determination of protein residue protonation states, and the acid-base chemistry of surface sites.

The authors state that the main practical limitation is the choice and optimization of the potential: the level of theory to which the model is fine-tuned determines the final accuracy, a generalized-gradient functional is unlikely to give quantitative results, and even the chosen hybrid functional may itself benefit from a cancellation of errors. The deviation observed for the first dissociation constant of phosphoric acid may come from residual error of the reference method, of the model, or from counterions in a small box, and these cannot be disentangled without further reference calculations. The quantum correction varies with species and with potential and cannot be treated as a constant transferable across systems or levels of theory; agreement with experiment in purely classical calculations more likely reflects a compensation of errors than the sufficiency of a classical treatment.