Insights about accelrated dynamics calculations of the acid <i>pK_a</i> beyond biasing the coordination number collective variable

Knowing the dissociation constant K_a of a weak acid is fundamental to understand many chemical and physical processes. Theoretical calculation of the K_a remains a challenge. The use of accelerated methods such as metadynamics has shown promise for numerical computation of K_a with explict solvent molecules. In the most common approach, a single collective variable (CV) representing the coordination number of the proton donor group yields results that are in reasonable agreement with experiments, however there remains questions as to whether the configuration space of the deprotonated state was sufficiently sampled given the small typical simulation box. We studied the deprotonation of acetic acid using the ReaxFF in simulation boxes of varying sizes and observed significant size dependence of ΔG when biasing with a single CV representing the coordination number of the acetate oxygen atoms. However, biasing with an additional CV describing the H₃O⁺ and acetate anion distance results in a virtual elimination of the size dependence. The improvement is due to accelerated sampling of the deprotonated configuration state.