  Classical molecular dynamics simulations have recently become a standard tool
for the study of electrochemical systems. State-of-the-art approaches represent
the electrodes as perfect conductors, modelling their responses to the charge
distribution of electrolytes via the so-called fluctuating charge model. These
fluctuating charges are additional degrees of freedom that, in a
Born-Oppenheimer spirit, adapt instantaneously to changes in the environment to
keep each electrode at a constant potential. Here we show that this model can
be treated in the framework of constrained molecular dynamics, leading to a
symplectic and time-reversible algorithm for the evolution of all the degrees
of freedom of the system. The computational cost and the accuracy of the new
method are similar to current alternative implementations of the model. The
advantage lies in the accuracy and long term stability guaranteed by the formal
properties of the algorithm and in the possibility to systematically introduce
additional kinematic conditions of arbitrary number and form. We illustrate the
performance of the constrained dynamics approach by enforcing the
electroneutrality of the electrodes in a simple capacitor consisting of two
graphite electrodes separated by a slab of liquid water.
