  The use of accelerated sampling methods such as metadynamics has shown
significant advantage in calculations that involve infrequent events, which
would otherwise require sampling a prohibitive number of configurations to
determine, e.g., 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 $K_a$ 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$_3$O$^+$ cation and the acetate anion results in a
substantially more efficient sampling both accelerating the dynamics and
virtually eliminating the computational box size dependence.
