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, modeling 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.
We recall that the Rα are fixed parameters in the potential. This prescription can be relaxed to include the location of the Gaussian charges to change but this extension of the model, that can be easily incorporated in the mass-zero constrained dynamics, is not considered here.
Note that a different initialization scheme was proposed in Ref. 17. This scheme is more rigorous, but its implementation in MetalWalls turned out to be impractical. Use of the conjugate gradient minimization, on the other hand, is available in the code. This choice of initialization did not cause visible problems and has negligible additional cost.