ELECTRODE: An electrochemistry package for atomistic simulations

Constant potential methods (CPM) enable computationally efficient simulations of the solid-liquid interface at conducting electrodes in molecular dynamics (MD). They have been successfully used, for example, to realistically model the behavior of ionic liquids or water-in-salt electrolytes in supercapacitors and batteries. The CPM models conductive electrodes by updating charges of individual electrode atoms according to the applied electric potential and the (time-dependent) local electrolyte structure. Here we present a feature-rich CPM implementation, called ELECTRODE, for the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), which includes a constrained charge method and a thermo-potentiostat. The ELECTRODE package also contains a finite-field approach, multiple corrections for non-periodic boundary conditions of the particle-particle particle-mesh solver, and a Thomas-Fermi model for using non-ideal metals as electrodes. We demonstrate the capabilities of this implementation for a parallel-plate electrical double-layer capacitor, for which we have investigated the charging times with the different implemented methods and found an interesting relationship between water and ionic dipole relaxations. To prove the validity of the one-dimensional correction for the long-range electrostatics, we estimated the vacuum capacitance of two co-axial carbon nanotubes and compared it to structureless cylinders, for which an analytical expression exists. In summary, the ELECTRODE package enables efficient electrochemical simulations using state-of-the-art methods, allowing one to simulate even heterogeneous electrodes. Moreover, it allows unveiling more rigorously how electrode curvature affects the capacitance with the one-dimensional correction.


I. INTRODUCTION
A common approach to treating electrodes in atomistic simulations is to assume them to be uniformly charged walls, either structureless or atomically resolved. In the case of equilibrium electrolytes near planar electrodes at low charge densities, this approach is known to capture the electrochemical properties well. Several studies, however, have emphasized the importance of polarization of the electrodes by the ions and molecules in their vicinity [1][2][3][4][5] . In more realistic electrochemical scenarios, constant potential method (CPM) molecular dynamics (MD) results are often significantly different from those obtained with uniformly charged electrodes 6 .
While the CPM 7,8 is a popular tool for modeling metal electrodes by dynamically updating individual charges on electrode atoms, alternatives such as image charge methods [9][10][11][12] are commonly used to enforce a constant potential for planar electrodes. While one of these methods can handle nonplanar surface by inducing a charge density on the interface between two media most of them are limited to planar electrodes. 12 Those approaches faithfully reproduce the behavior of electrolytes near electrodes, particularly the correlation between thermal fluctuations in the electrolyte near the electrode and the induced-charge polarization of the electrode, while obtaining a realistic picture of the electrical double-layer. A CPM MD is able to capture the temporal response in the build-up and break-down of electric doublelayer and thereby allows realistic capacitor charging and dis-charging curves to be generated in silico 6,[13][14][15] . Interestingly, near highly charged planar electrodes [16][17][18] and nonplanar electrodes (such as curved substrates or nanoporous carbons) 6,[19][20][21][22][23][24][25] , CPM MD and Monte Carlo simulations 26 yield a spatially-specific charge polarization and a non-trivial electrolyte structure.
Here, we present a package for treating electrodes in MD simulations which interfaces with the Large-scale Atomic/-Molecular Massively Parallel Simulator (LAMMPS) 27 . Our ELECTRODE package uses the highly parallelized and efficient computational infrastructure of LAMMPS and allows interaction with many other packages and features already available in LAMMPS. This work builds in part on an earlier work in which we showed how a particle-particle particlemesh (P 3 M)-based calculation makes the electrostatic calculations of a CPM simulation more efficient 28 . In addition to some new enhancements to the CPM, this implementation provides a constrained charge method (CCM) and a thermopotentiostat (TP) 29 . To capture the electronic response of nonideal metals, a Thomas-Fermi (TF) model 30 is included. Both an Ewald and a P 3 M k-space solver are available for various constraints of the periodicity of the systems, such as infinite slabs, cylinders or fully 3D periodic systems. The ELEC-TRODE package also contains the closely-related finite field (FF) 31 and finite displacement (FD) 32 methods, which extend CPM MD with a slab geometry to fully periodic boundary conditions for increased computational efficiency. 33 A list of new features in the ELECTRODE package is presented in Section II including a brief description of their theo-retical background. In Section III, we summarize the concept of CPM MD using data from various CPM MD runs and rationalize it based on the charging times of an electrical doublelayer capacitor. We discuss briefly for which situation each approach is suitable and give in Section IV an overview of future development directions and possible applications.

A. Constant potential method
In atomistic electrochemical simulations, the system of interest is often a fluid electrolyte confined between two electrodes (cf. inset of Fig. 2). This could serve as an in silico nanoscale model of a capacitor, to optimize some metric, such as energy or power density, by modifying electrolyte composition or electrode structure.
The distinctive feature of CPM MD is the calculation of electrode charges which keep electrodes at a desired electrostatic potential. To achieve this, we first partition the potential energy, U, of an MD simulation into: Here U non−Coul includes all non-Coulombic interactions, U elyt includes all Coulombic interactions between electrolyte particles, and U elec includes all Coulombic interactions involving electrode particles (both with electrolyte particles, and with other electrode particles). While the former are treated with regular force field approaches, the last term is treated somewhat special. U elec , i.e., without electrolyte-electrolyte interactions, is written in terms of an electrode charge vector q comprising all electrode charges as with a matrix A, and vectors b and v; where b depends on the electrolyte positions {r}. The applied potential v has an entry for every electrode atom. The interactions between electrode atoms are represented by A, called elastance matrix due to the analogy between a vacuum capacitor and a spring. If the electrode atoms do not move, A can be pre-computed, allowing significant computational savings. The electrolyte vector b({r}) represents the electrostatic potential on each electrode atom due to the electrolyte atoms. At each step, q is updated to minimize the Coulombic energy contribution U elec , possibly subject to additional constraints. The desired energy-minimizing charge vector q * is straightforward to calculate 34 : Here, the elastance has been inverted to yield A −1 , which is called the capacitance matrix in light of its role in equation (3): the response of the charge vector q * can calculated as the product of the capacitance matrix with the vector of external potentials, analogously to the well-known scalar equation Q = CV linking the capacitance C to the charge Q.
Provided the electrode atom positions and thus the vacuum capacitance are constant, the main computational burden is the calculation of b({r}) at every time step, which is necessary due to the motion of the electrolyte. The primary purpose of ELECTRODE is to compute the electrode-electrolyte interaction in b efficiently and update the electrode charges accordingly. Alternatively, the electrode charges could be obtained with the conjugate gradient method, which solves the minimization problem without a matrix inversion. 35,36 Yet another approach is to treat the electrode charges as additional coordinates and perform mass-zero constrained dynamics for them 37 .
In MD with periodic boundary conditions, the simulation cell ideally is charge neutral. Scalfi et al. 34 showed that this constraint could be imposed by using the symmetric matrix as capacitance matrix instead of A −1 . Non-ideal metallic electrodes have been recently modeled by Scalfi et al. 30 using a semiclassical TF approach. We have implemented this promising approach in our ELECTRODE package, as its implementation is very similar to the selfinteraction correction of the Ewald summation 38 and contains only a single summation over the electrode atoms. An interesting alternative to effectively model a wide range of materials between insulator and ideal metal was proposed by Schlaich et al. 39 and involved using a virtual TF fluid within the electrodes. However, the virtual TF fluid approach appears computationally more expensive. Both models require free parameters, most crucially the TF length, l TF , in Scalfi et al. 30 and (a rather artificial) parameterization of the virtual TF fluid in the approach of Schlaich et al. 39 . Further, assumptions such as atom-centered densities prohibit effects such as quantum spillover and delocalization of the image plane. It should be noted that ELECTRODE provides a more flexible implementation of the TF model, allowing heterogeneous electrodes with different l TF for different atom types. For more general information on CPM MD approaches, the interested reader is referred to the excellent and thorough review of current electrode-electrolyte simulations by Scalfi et al. 40 , the well-written theory part of MetalWalls 41 or the thesis of Gingrich 42 .

B. Simulating an arbitrary number of electrodes
A CPM MD is typically performed with two electrodes, which means there are only two possible values for each of the n components of the potential in Eq. 2. In the ELECTRODE package, an arbitrary number N of electrodes is allowed with every electrode atom belonging to exactly one electrode. We define an electrode-wise indicator vector g α for every electrode α with n entries, which are equal to 1 if the respective electrode particle belongs to that electrode and 0 otherwise. The indicator matrix comprising the indicator vectors of all N electrodes allows us to connect electrode-wise quantities to particle-wise quantities. From hereon, we use tildes for electrode-wise quantities. For instance, we defineṽ as the electrode-wise potential and use it to write the potential v = Gṽ. Likewise, energyminimizing charges q * for a given set of electrode-potentials are C. Simulating electrodes at specified total charge In the CCM, the user sets the electrode-wise total chargẽ q * for each electrode. Such a fixed-charge setup corresponds to an open-circuit configuration 43 . This type of simulation has recently been attempted as a variation of the finite-field method 32 . Therein, it was found that ramping the total charge up or down over time could be considered computational amperometry, and a faster non-equilibrium response was observed.
Working with the capacitance matrix A −1 rather than the symmetrized matrix S (since charge neutrality is explicitly enforced by the appropriate choice ofq * ), we have: q * b defines the total charge each electrode would carry at zero potential, andC is the electrode-wise capacitance matrix. To subsequently estimateq * b , Eq. 7 is solved forṽ which is then applied using the CPM. This results in an energy minimization w.r.t. the charge distribution with a constraint on the total electrode charges. Analogous to how constant volume and constant pressure simulations can be thermodynamically equivalent, CCM and CPM simulations will give the same capacitances under suitable conditions. However, a thorough proof of that assertion is out of scope of this work and will be discussed in an upcoming work.

D. Simulating electrodes with a thermo-potentiostat
Deißenbeck et al. 29 recently presented a thermopotentiostat (TP) that takes into account the fluctuationdissipation relation of electrode charges at a given voltage and temperature in an electronic circuit. They have also provided a TP implementation based on a uniform charge distribution using the scripting capability of LAMMPS 27 . The ELECTRODE package provides an implementation that minimizes the energy with respect to the charge distribution and conforms to the formalism described by Deißenbeck et al. 29 . Our TP approach is currently limited to only two electrodes and instead of a vector of applied potentialsṽ, a potential difference between two electrodes is used. At every time step, the potential difference ∆v(t) between the two electrodes is evaluated to find the new capacitor charge according to Here, k B is the Boltzmann constant, τ v and T v are parameters of the TP and X is a normally distributed random number with a mean of 0 and a standard deviation of 1. The vacuum capacitance C 0 is obtained from the capacitance matrix 34 and the effective potential ∆v(t) is computed from the electrode charges and the electrolyte configuration (cf. Eq. 7). Hence, all quantities required to evaluate Eq. 9 are readily available in the CPM. The obtained capacitor charge ±q(t + ∆t) is applied using the CCM on both electrodes, respectively.

E. Simulations with different periodicity
The Ewald summation commonly assumes periodic boundary conditions in all three directions and has to be modified for systems with slab and one-dimensional periodic geometries. As shown by Smith 44 , a regular 3D Ewald summation for slab-like systems, which are periodic in the xy-plane but confined in z-direction, results in a dipole term M z is the z-component of the dipole of the simulation cell. This dipole term was subsequently used for correcting the infinite boundary artifact of slab-like systems 45 . This is known as the EW3Dc method, which is implemented in many MD codes including LAMMPS 27 . Several authors 38,46-48 have shown that the infinite boundary contribution in slab-like geometries can also be solved in an exact form. This rarely implemented EW2D solver is another cornerstone of the ELEC-TRODE package. Just like slab-like geometries, systems with only one periodic dimension require an appropriate treatment of the longrange electrostatic interactions. As shown by Bródka andŚliwiński 49 , the approach of Smith 44 can be extended for an infinitely extended one-dimensional summation: Here, z is the periodic dimension and M x and M y are the respective components of the total dipole of the unit cell. Contrasting established codes, the ELECTRODE package contains these corrections for one-dimensional periodic systems. And even more crucially, ELECTRODE is the first package (as far as we know) to implement these corrections in combination with a CPM. As an outlook, the electrostatic layer correction 50 in combination with PPPM 51 is also considered for implementation as an alternative to the EW3Dc approach. As demonstrated exemplarily for slab-like two-dimensional periodic systems, the boundary corrections can be easily incorporated into the CPM formalism by splitting the dipole components into their electrode and electrolyte contributions: This way dipole corrections fit into the linear form of the Coulombic energy in Eq. 2 that is used in the CPM and the computational effort for the electrode-electrolyte interaction scales linearly with the number of particles.

F. Simulating electrodes with the finite field method
In the FF method 31 , the potential difference between two electrodes is not directly specified using the applied potential v. Instead, the simulation cell is periodic in the z direction, i.e., without adding the artificial vacuum between the slabs required otherwise. The FF method allows efficient simulations of infinite electrode slabs, since no additional vacuum is required. However, complexly shaped electrodes or electrodes with a one-dimensional periodicity cannot be simulated with the FF method.
A potential difference ∆v 0 (cf. Eq. 8), is created in the FF method by introducing a z-directed electric (polarization) field of magnitude −∆v 0 /L z , creating a discontinuity of ∆v 0 across the periodic z boundary (and thus between the two electrodes on either side of the slab). In this formulation, the electrode Coulombic energy is Here, ζ ζ ζ is a vector containing the normalized z-positions of each electrode atom, namely {z/L z } with an offset for the bottom electrode to make the system symmetric along the z-direction, replicating the conductor-centered supercell in Ref. 31. The energy-minimizing charge q * in this model is which is equivalent to the standard CPM, replacing v with −∆v 0 ζ ζ ζ . by only 1 %, showing a good agreement between the methods at equilibrium.
When a voltage is initially applied with the CPM, the charge induced on both electrodes is very small since the capacitance of the electrode pair in vacuo is small 28,34 . However, the electrode charges induce the formation of a dipole in the electrolyte, which in turn induces additional charge on the electrode. CPM MD thus models the process of charging an electrical double layer capacitor physically correctly, and the charging-discharging curves obtained from CPM MD can be used to fit parameters for equivalent macroscopic electrical circuits 53 . Two charging times τ 1 and τ 2 are obtained by fitting a bi-exponential charging function to the z-component of the electrolyte dipole 54 . M eq z is the extrapolated equilibrium dipole reached at late times. Comparing the individual contributions to the total electrolytic dipole reveals that τ 1 describes relatively fast water dipole relaxations and τ 2 describes charging times related to the ion diffusion. To validate this statement, we show both individual components in Fig. 1. From the mixing parameters c in the panels, it is clear that water dipole relaxation dominates at the beginning, while ion diffusion prevails at later times. τ 1 of the water in Fig. 1b largely corresponds to that of the total dipole at the beginning of the charging. While the ions in Fig. 1c relax mainly on the slower timescale of τ 2 , the water dipole decreases as the ion dipole is slowly built up.   2 compares results using CPM, CCM and TP, the latter with a time constant τ v = 100 fs, to each other and to their uniformly-charged counterparts. We focus on the charging term with the faster timescale τ 1 because of its large contribution to the total dipole. For the uniform variants, the charges are always evenly distributed across the inner layers of the electrodes. The uniform methods are in general very close to their heterogeneous counterparts, which is consistent with previous studies that found only small differences between a heterogeneous and uniform CPM at low voltages for simple planar electrodes 16 . The charging times obtained with CCM appear to be too fast, since τ 1 is about two orders of magnitude smaller than obtained with the potential-based methods, which is a well-known effect 6,55 that could be related to rather unphysical high temperatures and voltages when applying an instantaneous charge with CCM on the electrodes 6 . CCM also predicts two times smaller τ 2 than other methods, while the contribution of the second exponential decay is almost negligible since c = 0.98 in Eq. 15. However, these values should be interpreted with caution, given that the biexponential curve is a poor fit in the case of a CCM (cf. Fig. S2). Using the TP little change in the slower timescale τ 2 is observed compared to the CPM and it only weakly depends on the time constant τ v (cf. Table S1). In contrast, the applied voltage v 0 in a CPM has an effect on the ionic charging times τ 2 and for a smaller applied voltage of 0.3 V, τ 2 drops to 50 ps (cf. Fig. S4).
Using a TF model for representing real metals impacts both charging times and the total dipole. It is interesting to note that while a decrease for τ 1 with increasing l TF in Fig. S3a is observed, interpretation of τ 2 is more complex when using a TF (cf . Table S1). Especially when comparing this to a regular CPM, i.e. l TF = 0, at the same voltage. The quite significant difference between a regular CPM and the TF for small l TF might be an artifact due to the rather thin metal slab model or is due to a complex interplay between charge screening in the metal and the water/ionic relaxation and/or the smaller total dipole obtained with TF model.
To understand the range of the second relaxation time τ 2 ≈ 50 − 100 ps, it is instructive to consider the product RC of the aforementioned areal capacitance C = 2.94 µF cm −2 and the areal electrolyte resistance R. Continuum models for ion dynamics have shown that the ionic relaxation time decently agrees with RC for applied potentials up to around the thermal voltage e/(k B T ) ≈ 25 mV; 56,57 a recent dynamical density functional theory for a dense electrolyte found that the ions relaxed with RC even around 1 V. 58 In a bulk electrolyte at infinite dilution, the areal resistance between two electrodes spaced L apart is given by R = Lρ, where ρ = k B T /(2e 2 Dc b ) is the ionic resistivity, k B T is the thermal energy, e is the elementary charge, D is the ionic diffusion constant, and c b is the salt number density. 59 In our simulations, the plate separation was L = 5 nm and the salinity in the bulk phase was approximately 0.95 M, corresponding to c b ≈ 0.57 nm −3 ; the ion diffusivity D ≈ 1.5 · 10 −9 m 2 /s was obtained from a separate bulk electrolyte simulation. Using these values, we obtained RC = 14 ps; roughly 7 times smaller than the largest fitted τ 2 . This discrepancy must be due partly to our underestimation of ρ, which, at the salinity of our interest, is larger by a factor of about 1.7; 59 accounting for this effect yields a relaxation time of RC = 23 ps. Another cause of the remaining factor 4 discrepancy between the largest fitted τ 2 and predicted ionic relaxation times is the nanoconfinement, which could affect the diffusivity D and, in turn, the areal resistivity R. Finally, the mentioned increase of τ 2 with the applied potential is in line with the potential dependence of the capacitance of the Gouy-Chapman model, though in disagreement with that of the Kilic-Bazant-Ajdari 60 model. Although the analytical estimates of RC times presented here are interesting and provide starting points for further research on the implications of nanoconfinement and finite salt concentration on charging times, these results should not be overinterpreted as the analytical models contain simplifications that may not apply to such nanoscopic systems.

B. Co-axial cylindrical capacitor
As a sanity check of our approach for systems which are periodic in just one dimension, we study the capacitance of two co-axial carbon nanotubes of radii R 1 and R 2 , with R 1 < R 2 (cf. inset of Fig. 3). The vacuum capacitance C 0 of the coaxial carbon nanotubes can be calculated from the electrodewise capacitanceC 34 . At large radii, the atomic structure of the tubes should have a negligible effect, and thus the capacitance should approach that of structureless cylinders. The analytical line capacitance for a given ratio of the radii is C 0 /L = 2πε 0 / ln(R 2 /R 1 ); in which ε 0 is the vacuum permittivity and L is the length of the simulation box in the periodic dimension. As shown in Fig. 3  large compared to the bond length d CC between carbon atoms. In the CPM electrode atoms are assigned a Gaussian charge distribution ρ i (r) = q i η 2 /π 3/2 exp −η 2 (r − R i ) 2 at their position R i with the reciprocal charge width η. In agreement with Serva et al. 61 increased capacitances are observed for larger Gaussian width (i.e. smaller η) in Fig. 3. However, the impact is almost negligible.

IV. CONCLUSIONS
We presented the ELECTRODE package as an efficient implementation of the constant potential method (CPM) and closely related methods for the popular LAMMPS simulation environment. Initially, the main goal was to bundle many different approaches to electrochemical simulations into one package and ensure that they are handled in the most computationally efficient way. However, we also found interesting relationships between the two relevant charge-time contributions, i.e., water dipole relaxation and ion diffusion. We also implemented several new features, such as the EW2D summation and a correction for systems periodic in just one dimension, whose capabilities and full potential have been scarcely explored and which also work independently of the CPM. Recent improvements to the CPM such as the finite field (FF) method and a Thomas-Fermi (TF) model were included and compared for consistency to results found in the literature. Remarkably, using the TF model with varying TF lengths has a suprising and complex impact on the water and ionic relaxation times.
These results demonstrate that the ELECTRODE package can efficiently simulate electrified interfaces, including unusual systems such as infinitely long charged nanotubes. For a capacitor composed of co-axial carbon nanotubes, the vacuum capacitance agrees well within the limit of the analytical result of a structureless cylindrical capacitor and enables the investigation of curvature-dependent effects 24,57 more rigorously in the future by avoiding interactions between the nanotubes through the periodic images. Moreover, the charging process of a plate capacitor with an aqueous NaCl electrolyte in between illustrates vividly the differences between the range of methods introduced here and which are used to estimate the electrode charges. Interestingly, in these simulations, it was observed that the water dipole initially responds very quickly to the applied potential but then slowly drops off as the ionic dipole slowly builds up, as if the water dipoles were shielded from the ions.
While the package is in a stable state, the development is ongoing and will include in the future features like a conjugate gradient solver or compatibility to TIP4P water models.

SUPPLEMENTARY MATERIAL
The supplemental material provides more background on the thermo-potentiostat (TP), boundary corrections and the FF method. Further, a description of the package interface to LAMMPS and an overview of the implemented classes are given. More details on the simulations are provided, including plots of trajectories and of charging times as a function of the voltage, the TF length and the time constant of the TP.

DATA AVAILABILITY STATEMENT
The ELECTRODE package has been merged into the official release of LAMMPS which is available under github. com/lammps/lammps/tree/release.

A. Interface
The ELECTRODE package is designed as an add-on to the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) and integrates seamlessly with most molecular dynamics (MD) simulations carried out with that package. The user must first select one of two long-range Coulombic solvers ("kspace styles" in LAMMPS terminology) included in ELEC-TRODE, which implement either an Ewald summation 1 (ewald/electrode) or a mesh-based summation 2,3 (pppm/electrode) providing additional methods to compute long-range contributions for the matrix A and vector b required in the constant potential method (CPM). The user can then include the following command in the simulation script: fix <ID> <group> electrode/<mode> <value> <eta> followed by optional keywords. Mandatory inputs are as follows: • ID is used for identification of the fix. The instantaneous applied potential for each electrode is available via f_<ID>[<x>].
• group specifies the LAMMPS group of particles to which the potential or charge value will be applied via CPM or constrained charge method (CCM), respectively; additional groups could be specified with the couple keyword.
• mode chooses which of the following variants is performed: conp: constant potential conq: constrained charge thermo: thermo-potentiostat (TP) 4 • eta specifies the reciprocal width of Gaussian smearing applied to the electrode charges.
The following optional keywords are available: • couple <group-x> <value> allows additional groups to be specified as an electrode. This keyword can be used multiple times to specify more than two electrodes.
• etypes <types> allows users to specify atom types exclusive to the electrode, so that LAMMPS can provide optimized neighborlists for faster calculations.
• symm [on/off] allows users to enforce charge neutrality for the electrodes via matrix preprojection (cf. Eq. 4).
• ffield [on/off] allows users to turn on or off the finite-field implementation 5,6 . ELEC-TRODE automatically creates the electric field required to polarize the simulation box.
• [write/read]_[mat/inv] <path> specifies a file path for input or output of the elastance or capacitance matrix.
• temp <T-v> <tau-v> <seed> specifies the parameters T v , τ v and the seed for the random number X in Eq. 9. This keyword is only available with mode thermo, i.e. the TP.
In addition, parameters for the Thomas-Fermi (TF) model can be set via: fix_modify <ID> tf <type> <length> <volume> with the parameters length and volume corresponding to the TF length, l TF , and the reciprocal particle density, d −1 , as defined in Ref. 7, respectively. Parameters need to be specified for each atom type in the electrode. Further, the pppm/electrode long-range solver can be configured to calculate its part of the elastance matrix in two steps via: The twostep option can be faster for large electrodes and moderate mesh sizes but requires more memory 8 .
The atomic charges of the relevant electrode groups will then be dynamically updated over time, resulting in a MD simulation with CPM, CCM or TP dynamics as selected. The computed electrode charges are accessible using standard LAMMPS commands.

B. Implementation
In contrast to previous versions, our implementation separates the class for charge updates from the computation of the elastance matrix and the electrolyte vector b. The FixElectrode class is derived from the LAMMPS Fix class and is used to determine and set electrode charges. As shown in Fig. S1, the computation of the short-range, long-range and self interaction as well as the boundary correction are distributed across four additional classes.
The ElectrodeMatrix and ElectrodeVector classes handle the computation of the elastance matrix and electrolyte vector, respectively, with each class containing its own code for the short-range Coulombic interactions. For long-range interactions and boundary corrections, the two classes in turn delegate computations to implementations of the abstract ElectrodeKSpace class. The two implementations inherit from the LAMMPS KSpace classes for the Ewald summation and particle-particle particle-mesh (P 3 M) method, respectively. For the functions of ElectrodeKSpace, additional Fourier transforms for the electrolyte charges and the Green's function at the position of electrode particles have been implemented 8 . A BoundaryCorrection class was also added, which encapsulates the existing LAMMPS implementation of the EW3Dc method 9,10 and adds new implementations for a 2D Ewald summation 11,12 and a correction intended for 1D-periodic systems 9,13 .
As developers, we have found that an object-oriented programming approach has offered significant advantages during our work in writing the ELECTRODE package. Splitting the code across distinct classes has significantly clarified information flows between the different parts of the ELECTRODE package, leading to easier maintenance and optimization.

C. Comparison with other implementations
Another open source implementation of the CPM is found in MetalWalls. 14 In contrast to ELECTRODE it is a stand-alone MD code and is dedicated to electrochemical systems. Met-alWalls offers features such as the conjugate gradient and the mass-zero 15 method for the energy minimization, dipole interactions and the TIP4P water model, which are not yet available in ELECTRODE. Further, MetalWalls offers tools for post processing that are particularly useful for electrochemical simulations. Fortunately, many of these scripts can also be used for LAMMPS trajectories.
On the other hand, ELECTRODE is integrated into LAMMPS which supports a wide range of force fields and MD features which may not be available in MetalWalls, including manybody and machine-learning interatomic potentials, interfaces with quantum calculation software, multi-replica simulations and non-equilibrium MD. ELECTRODE uses self-contained Coulombic solvers to determine electric potentials, including requesting an appropriate neighbor list from LAMMPS. Thus, ELECTRODE will work with any force field that includes long-range Coulombic interactions, including force field combinations enabled under the LAMMPS keyword 'hybrid'. (The usual caveat applies: just because an MD simulation can run does not ensure it is physically valid.) ELECTRODE reuses many data structures native to LAMMPS (such as neighbor lists for nearby interactions and P 3 M implementations), thus benefiting from the high scaling efficiency that has been documented in LAMMPS publications. Additionally, ELECTRODE can use the efficient P 3 M method for long-range interactions required for a CPM. 8 A notable omission from ELECTRODE as compared to MetalWalls is that the charge update of ELECTRODE is not synchronized with other features which attempt to minimize energy midstep, such as Drude model polarizable force-fields and charge equilibration of mobile particles in ReaxFF and other reactive force-fields. MetalWalls includes an implementation of polarizable ion models represented by point dipoles, where electrode charges and point dipoles are simultaneously updated in a synchronized energy minimization (using conjugate gradients). By contrast, in an ELECTRODE simulation with (for example) Drude dipoles, electrode charges will not be reoptimized after Drude dipoles have been optimized (or vice versa), so in theory electrode charges may be out of sync with Drude dipole configurations, and simulation results will change based on the order in which commands are entered in the script. In practice, stable integration requires that both electrode charges and Drude dipoles do not change abruptly over the integration timestep.
We therefore think that the absence of synchronized minimization is not a major obstacle to using ELECTRODE with polarizable force fields, but look forward to verifying this in future with detailed simulations.
Another available implementation in LAMMPS is the 'USER-CONP' code by Wang et al. 16 This code appears to be no longer maintained and is poorly optimized compared to other implementations. A follow-on version entitled 'USER-CONP2' 17 was written and maintained by one author of this paper, and most of its features have now been incorporated into the ELECTRODE package. An implementation of CPM MD for the OpenMM software 18 can also be found at https://github.com/jmcdaniel43/Fixed_Voltage_OpenMM, and may be of interest to users with access to graphical processing units (GPUs) for computational acceleration.

II. SIMULATION DETAILS
A capacitor with an aqueous NaCl electrolyte between two gold electrodes has been adapted from an example provided in the MetalWalls repository 14  Examples of the temporal evolution and the fitting are given in Fig. S2 for the CCM, CPM and the TP. A summary of the fit parameters for the various charging runs with different methods is given in Table S1. Fig. S3 features a CPM setup with the TF approach and various l TF . While not affecting τ 2 , increasing l TF leads to a faster water dipole reorientation (cf. Fig. S3a) and an up to 35% smaller equilibrium dipole compared to a pure CPM with l TF = 0 Å due to the screening of the charge in the metal electrode. In Fig. S3b different potentiostat couplings for the TP are investigated. It is expected that the TP alters charging times because it includes a dissipation term which leads to a delayed response in terms of a weaker coupling to the external (electron) reservoir via τ v = R 0 C 0 . 4 Increasing τ v leads thus to larger charging constants for the water dipole reorientation τ 1 , but does not affect the ionic diffusion relaxation time constant τ 2 (cf. Table S1).  Interestingly, the voltage here only slightly affects τ 1 , but strongly changes the ionic relaxation time τ 2 at voltages below 1 V.