Molecular dynamics simulations of metal-electrolyte interfaces under potential control

The interfaces between metal electrodes and liquid electrolytes are prototypical in electrochemistry. That is why it is crucial to have a molecular and dynamical understating of such interfaces for both electrical properties and chemical reactivities under potential control. In this short review, we will categorize different schemes for modelling electrified metal-electrolyte interfaces used in molecular dynamics simulations. Our focus is on the similarities between seemingly different methods and their conceptual connections in terms of relevant electrochemical quantities. Therefore, it can be used as a guideline for developing new methods and building modularized computational protocols for simulating electrified interfaces.


Introduction
The metal-electrolyte interface plays a key role in electrochemical energy storage and conversion.In general, there are two kinds of idealized electrodes: ideal polarizable and ideal non-polarizable.The former as exemplified by the Hg electrode involves only capacitive charging and zero exchange current (density); the latter, taking the example of the Ag/AgCl electrode, allows free passing of faradic current with zero charge transfer resistance [1].Different from the meaning of dipole-moment density in physical chemistry, polarization here refers to the applied potential.In other words, an ideal polarizable electrode can sustain an applied potential and form an electric double layer.Therefore, a realistic computational model of the metal-electrolyte interface under electrochemical conditions should include the effects of an external potential.
Despite that the grand canonical (GC) formulation of density functional theory (DFT) was introduced right after the birth of canonical DFT in the mid 60s [2], its implementation for a slab system which resembles an metalelectrolyte interface in electronic structure calculation came out much later.Notably, Lovozoi et al. [3] discussed strategies of applying GC-DFT to metal slab systems under periodic boundary conditions (PBCs) already in 2001.In particular, they realized that PBCs convoluted the electrostatic potential of systems with a net charge because of the homogeneous compensating background and that an electrostatic correction must be applied to restore the physical one.
Parallel to the development of GC-DFT approaches, attempts were made to use charge transfer from ions to the metal interface to mimic the electrified interfaces with canonical DFT.Early examples were given by Skúlason et al. [4] and Rossmeisl et al. [5] around 2007, where different numbers of hydrogen atoms were added at the the Pt-water interface, to study the system at different electrode potentials.In these setups, the hydrogen atoms turn to solvated protons in the water bilayer and release electrons to the metal surface, which provides a static model of the electric double layer.
Since metal-electrode interfaces involve both electronic and ionic responses (to an external potential), a dynamical picture of such interfaces beyond a static DFT description is clearly needed.The pioneering work in this direction was done by Price and Halley in 1995 [6].They carried out the first Car-Parrinello-type molecular dynamics (MD) simulations of metal-water slab systems with an applied potential difference between the slabs (of 1.36 V), plane-wave basis set and local pseudopotential.
These seminal works that developed independently have influenced later theoretical works that try to describe structural, dynamical and electronic properties of electrified interfaces between metal electrode and electrolyte solution on an equal footing.Therefore, an up-to-date view of this topic is desirable.On this note, we will focus on works based on DFTMD approaches and highlight the recent method developments in this field.Thus, discussions on their applications in electrocatalysis and electrochemical energy storage, can be found in excellent reviews elsewhere [7,8,9,10,11,12,13].Furthermore, we will only include methods based on equilibrium calculations, therefore, non-equilibrium methods such as Green's functions [14] are not discussed here.
In the following, we will first introduce necessary concepts in electrochemistry, which will help us to understand the motivation and the strategies behind different types of methods.Then, we will sort seemingly different methods into three categories and discuss examples in each type of method.Finally, a summary and outlook for future developments is also given.

Fermi level and electrode potential
The work function of a metal surface W M e is defined as the energy required to take an electron from it into vacuum.For bulk metal, it is well-known that the opposite of the metal work function equals to the Fermi level of the metal E M F .The Fermi level is equal to the electrochemical potential μM e , which can be partitioned into the chemical potential µ M e and the contribution from the Galvani potential ϕ M (see Figure 1a for illustration): Then, the corresponding thermodynamic electrode potential will be written as: where the electron charge is defined as −e 0 .
Similarly, the thermodynamic electrode potential of metal in solution equals to the work function of metal in solution (see Figure 1b).
This is a key conclusion from a series of works from Trasatti and elaborated in Section 2.4 in Ref. [15].In the expression above, we emphasize that the work function is well-defined when the surface charge σ is zero as by definition it does not include any contributions from the Volta potential ψ (see Figure 1a and its caption).Subsequently, it means the corresponding Fermi level of the metal electrode in solution can be expressed as: For the electrified metal-electrolyte interface, the double-layer potential ϕ EDL will be built up and this leads to the expression of the Fermi level as: In addition, by definition, the double-layer potential is also equal to the change in the Galvani potential difference S ∆ M ϕ = ϕ M − ϕ S between metal (M) and solution (S), i.e.
As noted by Trasatti, despite "that the experimentally measured potential has nothing to do with the actual electric potential drop across the interface" [16], they change by the same amount when the metal electrode is polarized.Thus, one can adjust the change in the Fermi level of metalelectrolyte interfaces in order to control the double-layer potential.
As shown in Section 3, this idea is exploited in grand canonical DFT or constant Fermi DFTMD, in which the total number of electrons in the system is varied in order to control the change in the Fermi level.

Charge transfer and Volta potential
The alternative way to realize the "constant potential" is through charge transfer within the whole system itself (in contrast to the external reservoir) .In the following, we will use two examples to illustrate this point.
Supposing that we have two metal species α and β.When they are separated apart, the electrochemical potential or the Fermi level of each species is: where µ is the chemical potential and χ is the surface potential (see Figure 1a and its caption).Then, we put them in contact, as shown in Fig. 2a.Supposing μβ is more negative as compared to μα , then a charge transfer will take place and the two metal species will become electrified.Now the electrochemical potential of each species is: where the Volta potential ψ is no longer zero.The electrochemical potential of two metal species in contact must be the same.This means the chemical potential, the surface potential and the Volta potential have the following relation: where we have used the relation between the work function, the chemical potential and the surface potential.It is worth noting that in this conceptualization, the change in the surface potential χ upon contact is not considered.
The significance of the above equations is that the system of two metal species in contact is actually under a constant potential, since the work function of each metal species is an intrinsic material's property.This means one could use a counter-electrode or counter-ion which has a different work function to control the Volta potential difference.
Along the same line, the follow-up example is two metals of the same species in contact under an external voltage ∆V (Fig. 2b).Since two metals must be in equilibrium, this implies: In other words, charge transfer happens between two metals of the same species under an external voltage ∆V .In fact, since no field can be sustained within the same species of a metal, this means the two surfaces of the same metal are electrified.
Comparing two scenarios shown in Fig. 2, we can see both systems are under potential control without varying the number of electrons.This means the whole system is canonical while each of two metals in contact is grand canonical.As shown in Section 3, this can be realized through either counterion/pseudo-atom methods or finite-field methods.In grand canonical DFT the number of electrons in the system is allowed to vary over time, which can be used to charge an electrode.The electrode Fermi-level is then forced to become equal to μtarget e .The change in Fermi level corresponds to a build up the double layer potential as shown in Equation 7.During constant Fermi-level DFTMD, the electrochemical potential is set to a target value during SCF in practice.This calculation may be difficult to converge which is why multiple canonical DFT calculations were often used instead to approximate the grand canonical ensemble [17].Another approach by Bonnet et al. to achieve this is to devise a potentiostat scheme where the calculation of the number of electrons is separated from the optimization of the electronic states [18].The potential energy of the system in contact with an electron resevoir at μtarget e is

Potential control with half cell and full cell models
and the derivative with respect to the number of electrons can be interpreted as a fictitious force: Now we can let the electronic charge be a dynamical variable governed by where p ne and m ne are fictitious momentum and mass.The fluctuation of the number of electrons about ⟨n e ⟩ then allows for sampling the grand canonical ensemble to make ⟨μ M e ⟩ = µ target e .This method was used by Bouzid and Pasquarello [19] to study the charged Pt(111)-water interface.With the addition of a hydronium ion to the system, the electrode potential can be aligned to the standard hydrogen electrode (SHE) by considering hydrogen adsorbed to Pt as an intermediate step in the SHE reaction, similar to the proton insertion method that developed for doing the band alignment at the semiconductor-water interfaces [20].This allowed them to relate properties of the interface, such as double layer capacitance, to experiment.
The setup for constant Fermi-level DFT-MD is a half cell.The electrode charge then has to be compensated somehow to make the overall system charge neutral under PBCs.Bonnet et al. used an electronic screening medium (ESM) whereas Bouzid and Pasquarello made use of a homogeneous background charge.The homogenous background charge results in spurious interactions and corrections to the energy and potential have to be included [21].Therefore, the key difference between different methods in the Type I category is about how the electrolyte model (therefore counter-charge) is designed and implemented.Indeed, various implicit solvation models have been develope for this purpose [22,23,24,25].
It is worth noting that simulations of the metal-electrolyte interface with GC-DFT can be done without explicitly referencing the Fermi-level.This point was illustrated by the constant inner potential (CIP) DFT method from Melander et al. [26].If the inner potential in bulk solution is set to zero for both neutral and charged interfaces, ϕ S = 0, then following Equation 6the double layer potential is ϕ M can be controlled by setting a constraint on the average potential in the bulk electrode.

Type II: counter-ion/pseudo-atom methods
In a full cell setup with two electrodes of the same metal, constant Fermilevel DFT cannot be used to induce a potential bias since the system only has one Fermi level.From the example in Figure 3a) however, a system with two different metals in contact will have a potential bias across the cell stemming from the charge transfer induced by the equalization of their Fermi levels.A similar situation can also happen when a metal is put into contact with a charged semiconductor/insulator.If the Fermi level of the metal is in the band gap of the charged semiconductor/insulator, then any excess charge is going to transfer to the metal.
Surendralal et al. [27] exploited this idea to control the potential of the anodic Mg/water interface using doped (pseudoatom) Ne as a counterelectrode.Ne was chosen for its large band gap and its valence band maximum/conduction band minimum relative to water.To generate a surplus or shortage of electrons in the system, the number of valence electrons and proton charge of each of the n Ne Ne atoms is changed by a fraction q/n Ne .This way the total system stays charge neutral while each electrode obtains a charge q.To maintain a constant potential in the cell the counter-electrode charge can be altered during the dynamics.Similarly, in the work by Le et al. [28] the Pt electrode is charged by adding Na/F atoms close to its surface.Khatib et al. [29] used the same strategy in their ion imbalance method but with Na/Cl atoms instead.In these setups, the electrode charge is constant (corresponding to the ion charge) while the double layer potential ϕ EDL fluctuates.Although the double-layer can be modelled effectively in this manner, it is worth noting that the ion distribution is not fully in equilibrium as there is no real driving force to keep these ions next to the surfaces.

Type III: finite-field methods
Recently, Dufils et al. [30] showed that constant potential simulation with the Siepmann-Sprik (SS) model which is commonly used in together with 2D PBCs can be also realized with 3D PBCs and the finite-field methods [31].
The basis of the SS model [32,33] is to allow the electrode charges to fluctuate in response to the external potential, which mimics the physical process of an ideal polarizable electrode.Although the original model was built for simulating metal-water interfaces under potential bias, its applications to electrochemical interfaces with electrolyte solution are straightforward.
In the SS model, each response charge of the electrode atoms follows a Gaussian distribution of magnitude c i centered on the position of the electrode atom R i where ζ i is an adjustable parameter related to the Gaussian width.
The original model at zero applied potential can be rewritten as follows using the chemical potential equalization ansatz [34] where U 0 corresponds to the energy of electrode atoms in absence of an external potential (field).The term U q 0 −∆ν corresponds to the electrodeelectrolyte interaction (so electrostatic interactions between the atomic charges of electrolyte atoms and the base charges q 0 of electrode atoms plus their van der Waals interactions).η is the hardness kernel, describing the interaction between response charges and ∆ν is the potential generated by the electrolyte at the electrode atom sites.This energy is minimized with respect to the response charge c at each MD time-step under the constraint of charge neutrality.
With this model, a full-cell setup can be realized with two electrodes kept at constant potentials.If one electrode is set at zero potential we can enforce a constant potential difference ∆ψ between the electrodes by altering the model Hamiltonian as where the sum is over all charges belonging to one electrode.This requires an electrolyte centered cell and the use of 2D PBCs, which makes the potential calculation computationally more expensive than the Ewald summation with 3D PBCs.
Instead, with finite-field methods and 3D PBCs, the electrode charges are now coupled to an applied electric field E through where is the polarization in the direction of the electric field.In this scenario, one electrode with a single inner potential was employed instead and we can use an electrode centered cell.The cell potential bias now comes from the different charges at the two electrode surfaces, in analogy with the example in Figure 2b.Knowing the constant E field directly leads to constant potential simulation under 3D PBCs, applying a potentiostat on top of constant D simulations [35,36] may seem to be a detour.Nevertheless, interesting attempts have been made to control the potential difference using the constant electric displacement D Hamiltonian, which generates a constant ⟨E⟩ ensemble instead [37].

Potential of metal versus potential of electrolyte
So far, the focus of our discussions is on the potential control of metal.However, in an electrochemical cell with an active redox couple electrons can transfer between electrodes and solution.If we attempt to control the potential by just fixing the electrochemical potential of the metal electrode, μM e , there is no guarantee that the electrode potential is in equilibrium with the redox couple wihtin the time scale of simulation.For the chemical reaction otherwise the system is under an overpotential.This means the electrified interface requires the counter-ions to form a concentration gradient in equilibrium.Therefore, to simulate an electrochemical cell under constant potential conditions, it is crucial to also take into account the potential of the electrolyte.This is indeed very challenging as it requires a grand canonical ensemble for the ions and the corresponding equilibration times goes far beyond the reach of standard AIMD simulations.As a matter of fact, a number of selected works listed in the Table 1 only used water instead of electrolyte solution when building a solid-liquid interfacial system.Therefore, strictly speaking, these setups simulated a (nano)capacitor rather than an electrochemical interface.
Despite being challenging, potential control of the electrolyte can be achieved in various ways.With the help of liquid state theories, a molecular description of the electrolyte can be introduced, e.g. in joint DFT (JDFT) [38,39] where the solution is described by classical DFT [40] or in conjunction with the reference interaction site method (RISM) [41,42].With these descriptions of the electrolyte, the effect of an ion reservoir can be achieved, comparable to the electron reservoir in Type I methods.In Type II & III methods, explicit ions and electrolyte solution can be introduced and one can explore the classical MD simulation to speed-up the equilibration of the electrolyte solution and bring counter-ions to their equilibrium positions next to the electrified interfaces.Similar to the case of metal in Type II & III methods where two sides of electrode are under grand canonical condition but the whole piece of metal is still canonical, the same principle also applies to the electrolyte solution for its potential control.

Summary and outlook
In this short review, we have sorted computational methods for potential control in MD simulations of the metal-electrolyte interface into three categories.In general, the potential control can be achieved either through a grand canonical treatment of a half-cell model or a canonical treatment of the full-cell model.In both cases, charge transfer is induced by equalization of the Fermi level of a system by the exchange of electrons with an external reservoir or via redistribution of electrons within the system.
Work is still in progress on applying the finite-field methods to DFTMD simulations of the metal-electrolyte interface.In particular, it is interesting to see how the constant electric displacement method, that has been shown to be useful to compute Helmholtz capacitance of protonic double layers at metal oxide-electrolyte interface [43,44], can be also applied to the metalelectrolyte interfaces.
Finally, it is impossible not to mention machine-learning potentials (MLP) when discussing MD simulations nowadays.Indeed, a number of implementations toward the goal of modelling electrochemical systems using machinelearning accelerated atomistic simulations have emerged [45,46].Their further developments and the full harnessing of the scalability of MLPs are expected to play a crucial role in the future of DFTMD simulations of electrified interfaces beyond slab geometry and short dynamics.

Declarations of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.programme (grant agreement No. 949012).We thank M. Sprik for critically reading the manuscript and helpful suggestions.
The name of ion imbalance method was coined in this work.A proof-of-concept study using the finite-field molecular dynamics simulations for potential control with the Siepmann-Sprik-type model.

Figure 1 :
Figure 1: Metal interfaces with vacuum and solution.a) Metal-vacuum interface.The electrochemical potential μe consists of the chemical potential μ and a contribution from the Galvani potential ϕ.The Galvani potential is composed of the Volta potential ψ due to net surface charge and the surface (dipole) potential χ.With zero surface charge, ψ = 0. b) Metal-solution-vacuum interfaces.The work function M|S e is defined as the energy required to take an electron from the metal through the solution and into vacuum.The potential of zero charge U M|S pzc can be calculated from this work function if the reference is changed from vacuum to standard hydrogen electrode (SHE).

Figure 2 :
Figure 2: Generating the Volta potential difference due to charge transfer.a) Two metals of different species in contact in which the metal α has a smaller work function.The Volta potential difference equals to the opposite of the work function difference.b) Two metals of the same species in contact under an applied voltage ∆V .The Volta potential difference equals to the applied voltage.The definition of symbols can be found in the Text.

Figure 3 :
Figure 3: Three types of methods to impose the potential control for molecular dynamics simulations of metal-electrolyte interfaces.Type I: grand canonical methods; Type II: Pseuduo-atom methods; Type III: Finite-field methods.The dashed region is not part of the simulation box.

Table 1 :
Examples of molecular dynamics simulations of metal-electrolyte systems under potential control.