Can machine-learning potentials reliably compute acid-base dissociation constants in water?

Computing pKa from a single solute in a box

The paper asks whether acid-base dissociation constants can be obtained bottom-up from atomistic simulation rather than from data-driven fits to known compounds. The target quantity is the absolute pKa or pKb of a species in water, defined in the dilute limit through concentrations and the standard concentration c° = 1 mol/L, with pKa = −log10 Ka. For a conjugate acid-base pair in water these constants are tied to the solvent autoionization constant by pKw = pKa + pKb, a relation the authors later use as an experiment-free consistency check.

The difficulty is that the reaction breaks and forms covalent bonds, so the potential driving the dynamics must be reactive. Density functional theory is too expensive over the time scales needed for a free energy that can reach 100 kJ/mol, dissociation is a rare event requiring enhanced sampling, and the proton's quantum behavior cannot be ignored. The authors' answer is to replace the electronic-structure bottleneck with a machine-learning interatomic potential (MLIP), a model trained to reproduce a reference quantum-chemistry method at a small fraction of its cost, and to handle the remaining costs with targeted approximations rather than by cutting corners.

The same free energy profile can give different pK values

The paper's central methodological claim concerns how the equilibrium constant is read off a simulation. The authors use the reversible work theorem, in which the potential of mean force W(r) between the two dissociation products is converted into a dissociation constant by integrating the Boltzmann weight of the associated basin out to a cutoff distance rc. Because the integral is a probability, it counts how likely the ion pair is to sit in the associated basin relative to a non-interacting reference, so any dissociation distance that identifies the same basin yields the same constant.

That property does not hold for two estimators common in the literature. Taking the free energy difference between the minima of the associated and dissociated states depends on the collective variable's Jacobian and does not correspond to a well-defined standard state; for water autoionization the authors report that this estimator differs by 4.2 pK units from the reversible-work value evaluated on the same profile. Integrating the Boltzmann weight separately over associated and dissociated regions is closer in spirit but depends strongly on the interval chosen for the dissociated state, and the resulting box-size dependence is hard to correct. The reversible-work route instead concentrates all arbitrariness into one additive constant C.

One constant carries the finite-size correction

The additive constant C fixes the normalization of the pair distribution function, and it matters: an error of kBT ln10 ≈ 5.7 kJ/mol in C translates into one full pK unit. The authors determine C by matching the simulated potential of mean force, at separations beyond the contact ion pair, to the screened Coulomb interaction of two point charges in a periodic box of the same size and shape as the simulation cell, computed with Ewald-type summation. Because the reference is the periodic-cell interaction rather than the infinite-dilution one, the periodic electrostatics and box geometry are absorbed into C at no extra cost.

Two details make this workable. The accessible solid angle at a separation r is clipped by the box faces beyond L/2, so the reference profile must include a Jacobian-related correction term; the authors note that for their dodecahedral boxes this term is handled numerically by placing point charges on a grid. When the products carry a net charge, a neutralizing background charge shifts the Ewald result, and the authors subtract the corresponding Wigner self-energy term. They also report that the calibration still works when the MLIP does not include explicit long-range electrostatics, provided the simulated profile follows the screened Coulomb shape in the matching window. A sensitivity test varying rc from 4 to 8 Å found no species changing by more than one pK unit.

A signed distance puts acid and base branches on one axis

The products of dissociation are not tied to fixed atoms: the excess proton and the proton hole of hydroxide both migrate through the hydrogen-bond network by a Grotthuss-like mechanism. The authors therefore assign each hydrogen to nearby oxygen sites with smooth weights, producing a formal charge per site that is close to +1 for hydronium, −1 for hydroxide, and near zero elsewhere, with differentiable transitions during proton transfer. This smoothness is required because the collective variable is used for biasing.

The separation is then defined as a charge-weighted sum of minimum-image distances from a reference point at the acid or base center. An acidic dissociation gives a positive signed distance and a basic dissociation a negative one, so both reactions of an amphiprotic ion such as H2PO4− appear on opposite sides of a single one-dimensional profile, with the undissociated species near zero. For water autoionization, where both products are fluxional, the same idea is written as a dipole-like expression with a periodic circular-mean reference point. The practical payoff is that only a one-dimensional free energy profile needs to be converged, rather than the two-dimensional surface that earlier collective variables required.

Metadynamics, replica exchange, and a mass-integration correction

Sampling proceeds in two stages: metadynamics builds a bias that drives the system along the signed distance, and replica exchange with a fixed bias then produces converged Boltzmann weights. Ten independent runs per species, started from different configurations, give a direct statistical uncertainty estimate. For water, the spread corresponds to a few kJ/mol in the dissociated region and ±0.37 pK units on the classical pKw. The authors also note that the density increases by more than 1% from the associated to the dissociated state, which is why sampling is done at constant pressure.

Nuclear quantum effects are handled separately rather than by running full path-integral simulations along the whole path. The authors compute the quantum correction to the free energy difference by thermodynamic integration over the nuclear mass, which requires converging the quantum kinetic energy only in the two end states. This is the step that changes the headline numbers most: the correction is −3.54 ± 0.07 pK units for water and ranges from −0.11 to −3.05 pK units across the phosphate dissociations. Because the correction varies that much within one acid's ladder, the authors argue it cannot be replaced by a constant shift, even within the same potential.

Quantum corrections bring the phosphate ladder within chemical accuracy

The classical and quantum-corrected values are reported for water and six phosphate dissociations at 300 K. For water, the classical pKw of 17.86 ± 0.37 becomes 14.31 ± 0.37 after a −3.54 correction, against an experimental 13.93. For phosphoric acid, the first pKa goes from 3.00 ± 0.06 to 2.10 ± 0.09 against 2.16; the second from 9.00 ± 0.22 to 6.87 ± 0.23 against 7.19; the third from 15.06 ± 0.28 to 12.02 ± 0.29 against 12.33. Across the set, the mean absolute error relative to experiment drops from 2.19 to 0.45 pK units, which the authors place within the roughly 0.7 pK units corresponding to the conventional 1 kcal/mol chemical-accuracy threshold.

The clearest outlier is the basic dissociation of H2PO4−: 15.22 ± 0.13 classically and 12.85 ± 0.14 after a −2.38 correction, against an experimental 11.77. The authors flag this reaction as the only one in the set whose products, neutral H3PO4 and OH−, do not interact through a charge-charge term, so the Coulomb calibration of C cannot be applied to it. They state that the deviation may reflect the residual error of the reference method, the model, or counterions in a small box, and that these cannot be disentangled without further reference calculations.

Conjugate-pair sums test the protocol without experiment

Because the six phosphate constants are not independent, summing each acid's pKa with its conjugate base's pKb gives three separate estimates of pKw. The quantum values are 14.95 for the H3PO4/H2PO4− pair, 14.05 for H2PO4−/HPO4²−, and 14.20 for HPO4²−/PO4³−, to be compared with the directly computed 14.31. The latter two agree within combined uncertainties; the first deviates by 0.6 pK units. The classical sums, 18.22, 17.60 and 17.35 against a direct 17.86, show comparable deviations without a clear outlier, so the discrepancy for the first pair comes from about 0.4 units in the classical profiles plus about 0.3 units in the quantum corrections.

The authors recommend this check as a routine diagnostic for any polyprotic acid or conjugate pair, on the grounds that agreement with experiment can always arise from compensating errors, whereas a violation of pKw = pKa + pKb cannot. The paper's practical limitation sits with the potential: foundation models make a reactive MLIP obtainable for a new solute with a few thousand reference calculations, but the level of theory used for fine-tuning sets the final accuracy, and the authors consider it unlikely that a GGA functional would give quantitative results. They used revPBE0-D3, which has a good track record for water, while noting that even this choice may benefit from error cancellation. The transferable part is the protocol: moving from water to the six phosphate dissociations required no change in the order-parameter definitions or the analysis, only adjustment of the metadynamics parameters to the size of the dissociation free energy.

Sources used in this answer