Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems

Electrochemical systems involve many complex interfaces which are difficult to characterize experimentally and to simulate accurately. Different phenomena may occur at the molecular scale, such as electron transfers, ions adsorption, chemical reactions, etc Groß (2018); Fedorov and Kornyshev (2014). In the case of electrochemical energy storage devices, which are often operated at rather large voltages and for the largest number of cycles possible Armand and Tarascon (2008); Simon and Gogotsi (2008), understanding all these processes is compulsory to design new electrode materials and/or electrolytes. From the simulation point of view, several methods can be used depending on the target properties and the scale at which the system needs to be represented. In particular, density functional theory (DFT) calculations allow to determine the thermodynamics of reactions at the close vicinity of the electrode Cheng and Sprik (2012); Cheng et al. (2014); Chan and Norskov (2015). Nevertheless there are many problems for which it is necessary to treat the interface at length scales going beyond the nanometer, for example in order to understand the structure of the adsorbed layer that forms and how it screens the electrode charge Fedorov and Kornyshev (2014); Salanne et al. (2016). In such cases, thermal fluctuations play an important role, but DFT-based molecular dynamics cannot be used because of the large system sizes that are involved (typically thousands of atoms). It is thus necessary to use a more approximate method such as classical molecular dynamics (CMD).

Indeed, CMD is commonly used to study extended systems over long time-scales. As an example, hundreds of thousands of atoms can be simulated for microseconds to study biochemical mechanisms. Yet, in most classical calculations, the molecules/materials at play are electronically insulating, so that they are conveniently represented as collections of fixed charges over the whole simulation. This is of course not adequate to represent electrodes. More complex electrostatic schemes have been derived, among which a constant potential method was proposed by Siepmann and Sprik to account for the delocalized nature of the electronic cloud inside a metallic material and its perturbation by an external charge arising from the electrolyte Siepmann and Sprik (1995). Their method allowed to fix the potential of each electrode, therefore providing the correct framework for studying electrochemical energy storage devices. The method was later extended by Reed et al. to the case of two-dimensional periodic systems, Reed, Lanning, and Madden (2007) opening the way towards systematic studies of complete electrochemical cells Merlet et al. (2012).

The main ingredient of the Siepmann-Sprik method is the use of fluctuating electrode charges, in the same spirit as the charge equilibration Mortier, Ghosh, and Shankar (1986); Rappe and Goddard III (1991) scheme for molecular systems. These charges are thus treated as additional degrees of freedom in a similar way as the induced dipoles involved in polarizable force fields Madden and Wilson (1996) or as the wavefunction coefficients in DFT-based molecular dynamics Vuilleumier (2006). The same algorithms were therefore used as for these much more widely used cases. In their seminal work, Siepmann and Sprik have used the Car-Parrinello algorithm Car and Parrinello (1985), which consists in introducing a fictitious mass to each additional variable and treating it as a full dynamical degree of freedom. Along the years this algorithm was progressively replaced by Born-Oppenheimer-based approaches, in which the degrees of freedom other than the atomic positions are determined either in a self-consistent way or through energy minimization techniques. Recently, a new algorithm was proposed to solve numerically the temporal evolution of polarization degrees of freedom, in the general framework of constrained molecular dynamics. Coretti, Bonella, and Ciccotti (2018) The algorithm was shown to sample correctly the Born-Oppenheimer probability density and to be applicable to systems where the additional degrees of freedom must also satisfy conditions such as specific conservation or orthonormality properties, providing a viable alternative to current algorithms also for evolution based on orbital-free DFT. Bonella et al. (2020) Here we show that this approach, named mass-zero constrained dynamics, is also very efficient for computing the electrode fluctuating charges under constant applied potential conditions, allowing further to enforce the overall electroneutrality of the system. We test the algorithm on a simple capacitor made of graphite electrodes and a pure water electrolyte. It yields instantaneous charges in perfect agreement with other techniques, and accurate properties of the system when long simulations are performed.

In this section, we show how the mass-zero constrained dynamics can be used to simulate the fluctuating charge model. We begin by recalling the main features of the model, highlighting only the points more directly related to the derivation of the new approach. The extended Lagrangian associated with the mass-zero constrained dynamics and the corresponding evolution equations are then introduced. Finally, the most effective algorithm to integrate these equations for fluctuating charges, which differs from the one adopted in previous applications due to the characteristics of the model, is discussed.

The system of interest is an electrochemical cell composed of two metallic electrodes kept at constant potential and separated by an electrolyte. To set the stage, let us recall that the electric potential $V(\bm{r})$ generated by a spatial charge distribution $\varrho(\bm{r})$ is given by the Poisson equation

complemented with appropriate boundary conditions. (Atomic units are used throughout this section.) In our case, these conditions must reflect the fact that the electrodes are conductors and therefore the potential is constant in the space they occupy. Thus

where $\bm{r}$ is a point in the region, $\Omega_{\pm}$, occupied by electrode $\pm$ and $\Psi_{\Omega_{\pm}}$ is the known value of the potential on each electrode. An equivalent representation of the system, better suited for the introduction of the extended Lagrangian required by mass-zero constrained dynamics, is obtained by recalling that the Coulomb energy associated to the spatial charge density $\varrho(\bm{r})$ can be written as

The above equation implies that the electric potential can be obtained as the functional derivative of the Coulomb energy and that, in particular, the boundary condition in eq.(2) can also be expressed as

for $\bm{r}\in\Omega_{\pm}$. A more explicit expression for the Coulomb energy of the system is obtained by specifying the charge distribution of the electrolyte and of the electrodes. In analogy with standard models of charged solutions, the electrolyte charge density is represented as a set of point charges located at the positions of its $N$ atoms. This part of the charge density will be indicated as $\varrho^{q}_{i}(\bm{r})=q_{i}\delta^{3}(\bm{r}-\bm{r}_{i})$ where $q_{i}$ is the fixed charge of particle $i$, $\delta(x)$ is the Dirac delta function and $\bm{r}_{i}$ are the positions of the particles of the electrolyte. The specific feature of the fluctuating charge model, on the other hand, is to replace the continuous physical charge density of the metallic electrodes with a discrete set of $M$ Gaussian charges centered at fixed locations, $\bm{R}_{\alpha}$, in the $\Omega_{\pm}$ regions. The charge distribution assigned to the site $\bm{R}_{\alpha}$ is given by

where $\eta$ is a system-dependent parameter, fixed a priori in the model. The $Q_{\alpha}$, that represent the integrated charge on each electrode site, have to be determined for each configuration of the system to enforce the conditions in eq. (4). In particular, in a molecular dynamics simulation of the capacitor, these variables must be determined at each time step, as they respond to the evolution of the coordinates of the particles in the electrolyte along the trajectory. Substituting the total charge density $\varrho(\bm{r})=\sum_{\alpha=1}^{M}\varrho_{\alpha}^{Q}(\bm{r})+\sum_{i=1}^{N}\varrho_{i}^{q}(\bm{r})$ in eq. (3) and performing the integrals in the second equality, leads to the following explicit expression for the Coulomb energy

where we have used the notation $\bm{r}\equiv\{\bm{r}_{i}\}$, $Q\equiv\{Q_{\alpha}\}$, and where $R_{\alpha\beta}=|\bm{R}_{\alpha}-\bm{R}_{\beta}|$, $r_{ij}=|\bm{r}_{i}-\bm{r}_{j}|$. Note that in this potential the $\bm{R}_{\alpha}$ are time-independent parameters, specified by the geometry of the model. Given the assumed form of the charge density on the electrodes, the condition in eq. (4) can be reformulated in terms of the discrete “variational” parameters $Q_{\alpha}$ as

($\alpha=1,\dots,M$). In this discrete model for the electrode charge distribution, the condition of constant potential is enforced only at sites $\bm{R}_{\alpha}$, i.e. at the center of the Gaussian charge distributions, in contrast with (4) which holds at all positions in the metal. Previous work (Reed, Lanning, and Madden, 2007) has shown that, in spite of the approximate treatment of the charge distribution, this still enables accurate simulations of electrochemical system.

Recently, it was suggested to modify the original fluctuating charge model by forcing the $Q$ variables to respect also the so-called electroneutrality condition

for all accessed configurations of the electrolyte. The condition above mimics the presence of an ideal generator connected to the electrodes and enables a closed description of the system, i.e. one that does not require to take explicitly into account interactions with a charge reservoir. As discussed in previous workScalfi et al. (2020), eq. (8) also ensures that relevant physical properties such as the capacitance depend only on the potential difference between the electrodes and not on their absolute potential. In simulations, this condition has the added benefit to simplify the Ewald treatment of electrostatic interactions. While electroneutrality is only enforced on average by real generators, the more stringent condition (8) has been used in simulations Limmer et al. (2013a); Scalfi et al. (2020) providing accurate results.

The dynamics of the system within the fluctuating charge model is then defined as follows. The charges $Q_{\alpha}$, mimicking the reorganization of the electronic charge density within the electrodes, adapt instantaneously to the configuration of the particles in the electrolyte. The evolution of the electrolyte’s degrees of freedom is then governed by

where $m_{i}$ is the mass of the electrolyte particle $i$ and $\hat{Q}=\{\hat{Q}_{1},\dots,\hat{Q}_{M}\}$ indicates the values of the electrode charges that satisfy eq. (7) and (8). As mentioned in the Introduction, this dynamical model is formally identical to that of other adiabatically separated systems such as classical polarizable models or first principle molecular dynamics. Consequently, algorithms such as Car-Parrinello dynamics or minimization schemes — such as conjugate gradient — typically adopted in Born-Oppenheimer propagation have been employed also to solve the evolution of the fluctuating charge model. In the following, we show how to adapt the recently proposed mass-zero constrained dynamics to this problem. In this approach an extended dynamical system is defined that leads to the same evolution equations as above for the electrolyte degrees of freedom. On the other hand, the constrained evolution of a set of auxiliary dynamical variables is used to satisfy exactly the conditions in the model. Lagrangian mechanics offers the most suitable framework to discuss the method in detail. We begin by observing that the constraints (7) can be trivially interpreted as the requirement that the function $\mathcal{W}(\bm{r},Q)\equiv U_{c}(\bm{r},Q)-\sum_{\alpha=1}^{M}\Psi_{\alpha}Q_{\alpha}$ (where we indicated with $\Psi_{\alpha}$ the potential at site $\bm{R}_{\alpha}$) is at a minimum with respect to $Q$. Eq. (8), however, implies that not all variations of the $Q_{\alpha}$ are independent, and this must be accounted for when stating the minimum condition. This is conveniently done employing the formalism of Lagrange multipliers. To that end, we introduce the auxiliary function

where $\nu$ is the Lagrange multiplier associated to the electroneutrality condition. It is easy to see, looking at the definition of $\mathcal{W}(\bm{r},Q)$ and of $f(Q)$, that $\nu$ actually corresponds to a potential shift as already noted in previous work Scalfi et al. (2020). The minimum solution $\hat{Q}$ with $Q$ satisfying also eq. (8) is then given by the stationary point $\{\hat{Q},\hat{\nu}\}$ of $\mathcal{U}(\bm{r},Q,\nu)$. This leads to the $M+1$ conditions

As usual in the Lagrange multiplier scheme, the last equation above is in fact the electroneutrality condition but now obtained as a result of an optimization problem in the space that includes the Lagrange multiplier $\nu$ as a variable. Proceeding in analogy with the derivations of the mass-zero constrained dynamics presented in ref. 17; 18, we then extend the space of the dynamical degrees of freedom to include the variables $Q$ and $\nu$ and define the Lagrangian of the extended system as

In the equations above, we have introduced two (arbitrary) masses, $\mu_{Q}$ and $\mu_{\nu}$, associated with the variables $Q$ and $\nu$, respectively. To set the stage for the homogeneous limit that we will take on these quantities (see discussion below eq. (13)), we set $\frac{\mu_{Q}}{\mu_{v}}\equiv\kappa$, with $\kappa$ a constant of appropriate physical dimensions. Conditions (11) are then interpreted as a set of $M+1$ holonomic constraints to obtain the evolution equations

where the $\lambda_{\beta}$ are the (undetermined) Lagrange multipliers that enforce the constraints. The first derivatives of $\mathcal{U}(\bm{r},Q,\nu)$ in the second and third line of the equations above are null due to the conditions imposed, eq. (11). The equations for the $Q$ and $\nu$ variables can then be reorganized by dividing the non-vanishing terms (i.e. the constraint forces) by the masses of these degrees of freedom. In the third equation above, this means dividing by $\frac{\mu_{Q}}{\kappa}$. At this stage, the key step of the method is performed taking the limit $\mu_{Q}\rightarrow 0$ (see also ref. 21). In order for the auxiliary variables to have a finite acceleration in this limit, the ratio $\lambda_{\beta}/\mu_{Q}$ must remain finite implying that the $\lambda_{\beta}$ are proportional to the masses. In the limit, the constraint forces acting on the physical degrees of freedom vanish and the dynamical system takes the form

with $\gamma_{\beta}=\lim_{\mu_{Q}\rightarrow 0}\frac{\lambda_{\beta}}{\mu_{Q}}$ finite. In the first equation above, we used the fact that, given eq. (10) and the definition $\mathcal{W}(\bm{r},Q)\equiv U_{c}(\bm{r},Q)-\sum_{\alpha=1}^{M}\Psi_{\alpha}Q_{\alpha}$, we have $\bm{\nabla}_{\bm{r}_{i}}\mathcal{U}(\bm{r},Q,\nu)=\bm{\nabla}_{\bm{r}_{i}}U_{c}(\bm{r},Q)$.

Eq. (14) defines the mass-zero constrained dynamics for the fluctuating charge model with electroneutrality constraint. The evolution equations for the physical degrees of freedom, $\bm{r}$, do not depend directly on the constraints — consistent with eq. (9). The constrained evolution of the auxiliary variables, $\{Q,\nu\}$, however, guarantees that the stationary conditions for $\mathcal{W}(\bm{r},Q)$ are satisfied together with the electroneutrality constraint. The zero mass limit for the auxiliary variables is rigorously enforced by this dynamical system, leading to full adiabatic separation of the physical and auxiliary motions and exact sampling of the Born-Oppenheimer probability density of the system Bonella et al. (2020). The homogeneity parameter $\kappa$ introduced to take the limit of zero masses is a free parameter in the derivation of these equations of motions. To fix this parameter, one can note that it has the dimension $[\kappa]=\frac{[\nu]^{2}}{[Q]^{2}}$, that is the inverse of the square of a capacitance. The numerical value of the parameter can then be estimated from the typical capacitance for the system, e.g. the capacitance of the empty cell, $C_{\text{empty}}$. Exploratory calculations, however, showed that the numerical stability of the algorithm is insensitive to $\kappa$ in a very wide range around this value. An alternative point of view on $\kappa$ is to consider it as a scaling parameter for the additional Lagrange multiplier $\nu$ by introducing a scaled variable $\tilde{\nu}=\frac{\nu}{\sqrt{\kappa}}$. This new variable then has the dimension of a charge and it is natural to set $\mu_{\tilde{\nu}}=\mu_{Q}$. The dynamical equations resulting from this alternative perspective are strictly equivalent to Eq. (14).

We conclude this section noting that extending the set of auxiliary variables to include $\nu$ is not the only way to enforce the electroneutrality condition. For this particular case, it is in fact trivial to identify a set of generalized coordinates constraining the system to the correct hypersurface by reducing the dimensionality of the $Q$ to $M-1$ variables and setting, for example, $Q_{1}=-\sum^{M}_{\alpha=2}Q_{\alpha}$. The method proposed in this section, however, is completely general and can be applied also in cases where the definition of appropriate generalised coordinates is problematic. The only potential limitation of the method is numerical, stemming from the computational cost of enlarging the space of auxiliary variables.

We simulate a system consisting of a slab of liquid water between two graphite electrodes, with the two interfaces perpendicular to the $z$ direction of the box, as shown in Figure 1. Despite its apparent simplicity, this system is very representative of current simulation works, which aim at characterizing the hydration of metal surfaces, both from the structural Willard et al. (2009); Limmer et al. (2013b) and the dynamical Willard et al. (2013) points of view. The extension to more complex systems, such as dilute aqueous solutions Limmer and Willard (2015); Simoncelli et al. (2018), ionic liquids Merlet et al. (2012); Haskins and Lawson (2016); Park and McDaniel (2020) or superconcentrated electrolytes Li et al. (2018), is trivial and does not require further technical development.

Constant applied potential molecular dynamics simulations, which are increasingly used to study the interface between metallic electrodes and many different types of electrolytes, generally rely on a “Born-Oppenheimer” approximation. The models allow the electrode charges to fluctuate by enforcing the local potential to be equal to a value which is fixed for the whole electrode. In this work we have showed that it is possible to reformulate the problem, by treating these charges as dynamical variables with zero mass. The constant potential condition is recovered by using the constrained molecular dynamics formalism. This method is as efficient as the existing alternatives, but it has the advantage to allow additional conditions to be introduced in a natural way. The method also produces an algorithm which is time-reversible and symplectic for all the degrees of freedom. The new algorithm was validated by simulating a capacitor made of two graphite electrodes and a pure water simulation, with an applied potential of 1 V between the electrodes and an overall electroneutrality constraint. We have shown that the computed charges are essentially indistinguishable from those obtained via the Born-Oppenheimer approach for a given atomic configuration. We have then shown that the simulations yield structural properties of the fluid which are in perfect agreement with previous works. The results presented in this work pave the way to more general applications. In future work, in fact, the mass-zero constrained framework will be used to study systems consisting of constant potential electrodes combined with a polarizable fluid, where all the additional dynamic variables would be propagated through similar equations of motions. We also expect this method to allow fixing the charges on the two electrodes independently, thus opening the way towards comparison with more complex electrochemical experiments.