The protonation state of molecules and surfaces is pivotal in various disciplines,
including (electro‐)catalysis, geochemistry, biochemistry, and pharmaceutics. Accurately and
efficiently determining acidity constants is critical yet challenging, particularly when
explicitly considering the electronic structure, thermal fluctuations, anharmonic vibrations, and
solvation effects. In this research, we employ thermodynamic integration accelerated by com‐
mittee Neural Network potentials, training a single machine learning
model that accurately describes the relevant protonated, deprotonated, and intermediate states.
We investigate two deprotonation reactions at the BiVO4 (010)‐water
interface, a promising candidate for efficient photocatalytic water splitting. Our results
illustrate the convergence of the required ensemble averages over simulation time and
of the final acidity constant as a function of the Kirkwood coupling parameter.
We demonstrate that simulation times on the order of nanoseconds are required for statistical
convergence. This time scale is currently unachievable with explicit ab‐initio molecular dynamics
simulations at the hybrid DFT level of theory. In contrast, our machine learning workflow only
requires a few hundred DFT single point calculations for training and testing.
Exploiting the extended time scales accessible, we furthermore asses
the effect of commonly applied bias potentials.
Thus, our study significantly advances calculating free energy differences with
ab‐initio accuracy.