MetalWalls: A classical molecular dynamics software dedicated to the simulation of electrochemical systems

1 Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France 2 Sorbonne Université, CNRS, Physico-chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France 3 Department of Mathematical Sciences, Politecnico di Torino, I-10129 Torino, Italy 4 Centre Européen de Calcul Atomique et Moléculaire (CECAM), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland 5 School of Chemistry, University of Leeds, Leeds, UK 6 Department of Materials, University of Oxford, Oxford, UK DOI: 10.21105/joss.02373


Summary
Applied electrochemistry plays a key role in many technologies, such as Li-ion batteries, fuel cells, supercapacitors, solar cells, etc. It is therefore at the core of many research programs all over the world. However, fundamental electrochemical investigations remain scarce. In particular, electrochemistry is among the fields for which the gap between theory and experiment is the largest. From the computational point of view, there is no classical molecular dynamics (MD) software devoted to the simulation of electrochemical systems while other fields such as biochemistry or material science have dedicated tools. MetalWalls, a MD code dedicated to electrochemistry, fills this gap. Its main originality is the inclusion of a series of methods which allow a constant electrical potential to be applied to the electrode materials. It also allows the simulation of bulk liquids or solids using the polarizable ion model and the aspherical ion model. MetalWalls is designed to be used on high-performance computers and it has already been employed in a number of scientific publications. It was for example used to study the charging mechanism of supercapacitors (Merlet et al., 2012), nanoelectrowetting (Choudhuri et al., 2016) and water desalination devices (Simoncelli et al., 2018). A typical snapshot of a simulation cell for the latter is shown on Figure 1. Water molecules are in red and white, Cl − ions in cyan, Na + ions in green and the carbon atoms are colored according to their charge (from light grey to black when the charge increases).

Constant potential electrodes
MetalWalls employs most of the standard methods and algorithms available in typical MD packages. Here we discuss only its main originality, which lies in the possibility of simulating materials at constant electric potential instead of constant charge. This is interesting for simulating conductors where the charges are not fixed on the atoms but can spread in the material and adapt to the external medium. In MetalWalls, this is done by allowing the charge on the material atoms to fluctuate under the condition of constant potential. The method was first introduced by Siepmann and Sprik (Siepmann & Sprik, 1995). It was later adapted to the case of 2D periodic boundary conditions and extended to a Born-Oppenheimer-like approach in the first version of MetalWalls (Reed, Lanning, & Madden, 2007). The method was also discussed in the framework of Mass-Zero (MaZe) constrained molecular dynamics in ). An updated description of the statistical mechanics of the constant potential ensemble can be found in the reference .
Two different models are employed for the charge distributions of the electrode and electrolyte atoms constituting the electrochemical cell. On the one hand, electrode atoms carry a Gaussian charge distribution of amplitude Q i and width η −1 where R i are the electrode atoms positions.
On the other hand, electrolyte atoms have a point charge distribution where r i are the electrolyte atoms positions. The Coulomb potential energy for the full charge distribution ρ = ∑ i ρ i is given (in atomic units) by As it is long-ranged, it cannot be truncated at a cut-off distance. In MetalWalls, it is calculated using the Ewald summation method, whose specific expression depends on the periodicity (2D or 3D) of the system (Reed et al., 2007). Given the above form of the charge density, for fixed atomic positions it can also be expressed as a V c (Q), where Q is a vector containing all the N electrode atom charges (Q 1 , ...Q N ).
The electric potential generated by the charge distribution associated with an electrode atom is ∂Qi . It follows that the constant electric potential condition within the electrode materials is given by where Ψ i is the applied potential difference on the electrode sites. This is equivalent to where the total potential energy V (Q) is defined as and includes the work necessary to charge the electrodes. The charges Q that satisfy the constant potential condition are found by solving the system of N equations given by Eq. (5).
It is also frequent, in constant potential simulations, to enforce the so-called electroneutrality constraint, defined by the equation As the electrolyte ions undergo thermal motion, thereby changing the external potential, these equations must be solved to determine the electrode charges consistent with the constant potential condition. Two methods, described below, have been applied to achieve this.

Born-Oppenheimer approach
A first method to simulate the dynamics of the electrolyte atoms under the constant potential condition is to solve directly the equation for each electrolyte configuration. This method relies on the quadratic form of the Coulomb energy with respect to the electrode charges to solve efficiently the equation for the Q. In particular, V c can be written in the form which yields a set of equations for the charges given by AQ − b = Ψ, where the matrix A depends only on the geometry of the electrodes while the vector b depends also on the positions of the electrolyte atoms (Ψ is a vector containing the N electrode atom potentials (Ψ 1 , …, Ψ N )). The set of equations can then be solved by direct numerical matrix inversion yielding the charges This is very efficient when the electrodes are fixed as the matrix A is constant and the inversion has to be performed only once at the beginning of the simulation. MetalWalls also allows constant pressure P z simulations to be run by treating the electrodes as rigid pistons. In this case, a more efficient way to solve the linear system is by using the conjugate gradient iterative method.
The fulfillment of the electroneutrality constraint in this framework and its repercussions on the statistical mechanics of the constant-potential ensemble are discussed at length in .

Mass-Zero constrained MD (MaZe)
Another approach, based on the method introduced in (Coretti, Bonella, & Ciccotti, 2018), is to consider the motion of the electrolyte particles subject to the set of holonomic constraints given by This approach is based on an extended Lagrangian, where the electrode charges Q are considered additional dynamical variables for the system. Unlike previous extended-Lagrangian approaches à la Car and Parrinello, the mass-zero limit is enforced rigorously in this method and the additional degrees of freedom are properly separated from the physical ones in an adiabatic sense. This yields dynamical equations of motion for the additional variables which can be solved numerically by symplectic and time-reversible algorithms, while keeping the adiabatic separation. In (Bonella, Coretti, Vuilleumier, & Ciccotti, 2020) it has been proven that the dynamics also samples the Born-Oppenheimer distribution exactly. The application of the method to constant potential simulations is described in detail in  together with the possibility of enforcing the electroneutrality constraint.

Finite electric field variant
An extension of the constant potential approach was recently proposed, in which the electrochemical cell is simulated by including one electrode only, using 3D periodic boundary conditions. The electrode is set at constant potential and a finite field is applied to the whole cell. Two electrochemical interfaces form on the opposite sides of the metal, and we have shown that this approach yields similar results to the two electrodes setup for bulk electrode materials (Dufils, Jeanmairet, Rotenberg, Sprik, & Salanne, 2019).