  Multistate empirical valence bond (EVB) models provide an accurate
description of the energetics of proton transfer and solvation in complex
molecular systems and can be efficiently used in molecular dynamics computer
simulations. Within such models, the location of the moving protonic charge can
be specified by the so called center of charge, defined as a weighted average
over the diabatic states of the EVB model. In this paper, we use first order
perturbation theory to calculate the molecular forces that arise if a bias
potential is applied to the center of charge. Such bias potentials are often
necessary when molecular dynamics simulations are used to determine free
energies related to proton transfer and not all relevant proton positions are
sampled with sufficient frequency during the available computing time. The
force expressions we derive are easy to evaluate and do not create any
significant computational cost compared with unbiased EVB-simulations. As an
illustration of the method, we study proton transfer in a small liquid water
droplet consisting of 128 water molecules plus an excess proton. Contrary to
predictions of continuum electrostatics but in agreement with previous computer
simulations of similar systems, we observe that the excess proton is
predominantly located at the surface of the droplet. Using the formalism
developed in this paper, we calculate the reversible work required to carry the
protonic charge from the droplet surface to its core finding a value of roughly
4 $k_{\rm B}T$.
