Advantages of Multidimensional Biasing in Accelerated Dynamics: Application to the Calculation of the Acid pK_a for Acetic Acid

The use of accelerated sampling methods such as metadynamics has shown a significant advantage in calculations that involve infrequent events, which would otherwise require sampling a prohibitive number of configurations to determine the difference in free energies between two or more chemically distinct states such as in the calculation of acid dissociation constants K_a. In this case, the most common method is to bias the system via a single collective variable (CV) representing the coordination number of the proton donor group, which yields results in reasonable agreement with experiments. Here we study the deprotonation of acetic acid using the reactive force field ReaxFF and observe a significant dependence of Ka on the simulation box size when biasing only the coordination number CV, which is due to incomplete sampling of the deprotonated state for small simulation systems and inefficient sampling for larger ones. Incorporating a second CV representing the distance between the H₃O⁺ cation and the acetate anion results in substantially more efficient sampling, both accelerating the dynamics and virtually eliminating the computational box size dependence.