The ability to simulate electrochemical reactions from first-principles has advanced significantly in recent years. Here, we discuss the atomistic interpretation of electrochemistry at three scales: from the electronic structure to elementary processes to constant-potential reactions. At each scale, we highlight the importance of the grand-canonical nature of the process and show that the grand-canonical energy is the natural thermodynamic state variable, which has the additional benefit of simplifying calculations. We show that atomic forces are the derivative of the grand-potential energy when the potential is fixed. We further examine the meaning of potential at the atomic scale and its link to the chemical potential and discuss the link between charge transfer and potential in several situations.

## INTRODUCTION

Electrochemical transformations are of increasing importance in the transition of society away from fossil fuels, such as for electricity storage or the conversion of renewable electricity into chemicals or fuels. Simultaneously, atomic-scale calculations are increasingly able to describe the electrochemical interface with greater levels of sophistication. Fundamentally, electrochemical transformations occur in systems that are open to the flow of electrons, from an external circuit, which allows the potential to be controlled at the point of reaction. Increasingly, atomistic and electronic-structure methods that simulate the controlled potential are being employed to describe such transformations. In this Perspective, we discuss the thermodynamic implications of switching from standard atomistic simulations, in which the number of electrons is fixed, to grand-canonical simulations, in which the system’s electrical potential is fixed. We will discuss the implications at three different scales: at the electronic structure level, at the atomistic level, and, finally, at the reaction thermodynamics level.

## ELECTRONIC STRUCTURE CALCULATIONS AT CONSTANT POTENTIAL

Electrochemical reactions are chemical transformations in which the electric potential of an electrode is controlled by an external circuit, such as a potentiostat. Here, we examine the concept of potential at the atomistic scale. Introductory textbooks often show the potential as dropping off smoothly from an electrode surface as it is screened by the solvent’s double layer. The electric potential is a quantity that is accessible—and fundamental—in electronic-structure calculations, so we performed some simple simulations to analyze the behavior of the potential near an electrode. In Fig. 1(b), we have plotted some traces, along the *z* axis, of the electrostatic potential throughout a periodic calculation of the type often used to simulate such reactions. (Calculation details are in the section titled Methods.) The first observation is that the potential is not single-valued, but varies dramatically with position in the cell, with strong variations, in particular, near atomic nuclei. This picture is quite different from that often shown in textbooks.

In Fig. 1(c), we instead plot the *xy* average of the potential at each *z* coordinate to smooth out variations in-plane with the electrode surface. Here, we do this at two different charges—that is, with different numbers of electrons in each simulation—and we can see a slight difference in the two potential profiles. At both charges, we still observe large spatial variations in potential. If we examine the difference between these two profiles [Fig. 1(d)], we finally see something that resembles the smooth drop-off near the surface, as is often shown in textbooks. Thus, it is best not to think of the electric potential as a single-valued quantity at an electrode, but rather as a small—but significant—perturbation to the scalar potential.

While we often think of using charge (*q*) to control potential (*ϕ*), this can be difficult to assess in electrochemical simulations since *ϕ* is not single-valued, but rather varies in space. Instead, it is simpler to consider the electron’s *chemical* potential (*μ*_{e}), which is controlled by the number of electrons (*N*_{e}) in the system. The chemical potential—which is defined as the incremental change in the system’s free energy when electrons are added or removed—has a single value for the system and, thus, is unambiguous to define. The potentiostat controls *μ*_{e} by injecting or removing electrons from the system—that is, through current—which in static atomic configurations is manifested as changes in *N*_{e}. Since *μ*_{e} and *N*_{e} are conjugate variables, they cannot be individually varied—that is, if *μ*_{e} is specified, then we have no control over the number of electrons in the system, assuming that the other key variables (such as atomic positions, temperature, or system volume) are fixed.

In practice, to control the potential in atomistic simulations—whether we think in terms of *μ*_{e} or *ϕ*—the most straightforward approach is to use grand-canonical simulations in which *N*_{e} is a free variable used to set the potential to a target value. Grand-canonical density functional theory (GC-DFT)^{1–12} methods have emerged as tools to study potential-dependent reactions and processes at user-specified applied potentials. These methods add or subtract electrons to the simulation cell in order to achieve a desired potential and can therefore control the potential to any arbitrary precision. Since the electron's chemical potential is equivalent to the Fermi level against a reasonable reference state, this quantity can be accessed through the work function. Grand-canonical methods have been used to study the potential-dependence of a multitude of interfacial reactions, including hydrogen evolution,^{13–15} CO_{2} reduction,^{6,16–18} formic acid oxidation,^{19} and CO adsorption.^{20} In particular, these methods allow for the potential-dependent study of electrochemical reaction barriers—something otherwise accessible through computationally expensive extrapolation methods^{21}—but at relatively low computational expense and with precise control over the potential at which the barrier is evaluated.

Here, we use the solvated jellium (SJ)^{8} method to calculate reaction energies and barriers at the constant applied potential, as illustrated in Fig. 2. This electronically grand-canonical method is implemented in the GPAW^{22,23} electronic structure code and allows for computationally efficient calculations at the constant applied potential. The countercharge (“jellium”) is spatially constrained in the implicit solvent, and a solvent exclusion region is included in the innermost region of the interface to avoid spurious interactions between the implicit and explicit solvents. In our view, the primary purpose of the implicit solvent is to screen the large fields that would otherwise be created by the charge imbalance; what little energetic effects of solvation that occur on the top side of the water tend to cancel out between stages of a reaction. Importantly, we find an approximately linear relationship between *N*_{e} and *μ*_{e} (shown in Fig. 2) at a fixed geometry, which leads to rapid convergence to a specified *μ*_{e}, often with a computational penalty of only ∼30% over the course of a trajectory. In this work, we control the potential to within a threshold of 0.01 V for structural optimizations and 0.001 V for normal mode analyses. Full computational details are included in the section titled Methods.

A traditional electronic structure calculation takes place at fixed nuclear coordinates {*x*_{i}} and with a fixed number *N*_{e} of electrons—typically, *N*_{e} matches the positive charge of the nuclei (or pseudopotentials) in the system. Other information input to the calculation includes the electronic temperature *T*_{e}. These fixed quantities are the natural variables of the Helmholtz energy, *A*_{pot}, and thus, the potential energy is calculated in this thermodynamic potential, whose fundamental equation is given as

where the electronic entropy *S*_{e} is the conjugate variables of *T*_{e} and *F*_{i} is the component of the force vector associated with nuclear coordinate *x*_{i}. (Additional terms may be present if the system volume changes.) Here, we are using the subscript “pot” on *A*_{pot} to clearly distinguish that this energy only includes the system at a static nuclear configuration (that is, the potential energy surface), as opposed to the system free energy including nuclear motion, which we will examine later in this work.

In typical electronic structure methods, such as DFT, the force can be acquired at minimal additional cost by employing the Hellmann–Feynman theorem.^{24} By differentiating Eq. (1), we find the relationship between the force and potential energy,

In atomistic calculations, *A*_{pot} as a function of {*x*_{i}} is typically referred to as the “potential energy surface” (PES), and we get the expected relation that the force on any atom is the negative gradient of the PES at a constant number of electrons, electronic temperature, and position of all other nuclei. This relationship is the key to many atomistic methods. In stationary-point searches, the gradient of the potential energy surface is used to quickly converge to a local minimum or saddle point in *A*_{pot}. The consistency of the force with the potential energy surface is also required to conserve energy in molecular dynamics simulations.

In electrochemical systems, the number of electrons, *N*_{e}, is a free variable that is used to control the electronic chemical potential, *μ*_{e}. Thermodynamics allows us to transform our thermodynamic potential by invoking a Legendre transform,

where Φ_{pot} may be referred to as the grand-potential energy. This transformation allows us to re-express the fundamental equation with natural variables *μ*_{e}, *T*_{e}, and {*x*_{i}},

Therefore, each force can equivalently be expressed as

That is, the forces are simultaneously the negative gradient of the potential energy (*A*_{pot}) at constant charge and the grand-potential energy (Φ_{pot}) at constant potential.

The implications of this are significant for atomistic simulations. That is, all of our standard algorithms—such as for structure optimization, saddle point searches, or molecular dynamics—will work unmodified as long as the transformed energy Φ_{pot} is used when constant-potential simulations are employed. As we will see later, the same transformation (−*μ*_{e}*N*_{e}) is required for thermodynamically consistent reaction energies. Thus, Φ_{pot} can be considered the natural energy associated with the grand-canonical ensemble—and it is the energy returned to the user by default in the SJ method.

In Fig. 3, we illustrate the force-equivalence using the results of constant-charge and constant-potential DFT calculations. In these calculations, we start from a given charge and potential condition (the black dot) and then calculate *A*_{pot} along a constant-charge path and Φ_{pot} along a constant-potential path. Here, we clearly see that the two curves share a slope at the black dot—where the charge–potential conditions coincide—and we further see that this slope matches that of the force calculated analytically by the Hellmann–Feynman theorem. Thus, although there are two different potential energy surfaces (*A*_{pot} and Φ_{pot}), the force at any point is valid for either path. This also implies that stationary points on the potential energy surface—local minima and saddle points—are stationary points in either ensemble.^{25} That is, the atomic coordinates that describe a local minimum or saddle point found in the canonical ensemble (*A*_{pot}) are equivalently a local minimum or saddle point in the grand-canonical ensemble (Φ_{pot}).

## ELEMENTARY REACTIONS AT CONSTANT POTENTIAL

We will next turn to the topic of an elementary reaction, which takes place at an electrode surface. As the simplest example, consider a proton reacting to form an adsorbed hydrogen atom in the first step of the hydrogen evolution reaction. Conventionally, reactions such as this are written formally with an integer electron transfer, i.e., H^{+} + e^{−} + $*$ → $H*$, where $*$ represents an adsorption site on the surface of the electrode. However, if the reaction takes place at a fixed potential (expressed in terms of *μ*_{e}), then we have no control over the number of electrons (*N*_{e}) in either the initial or final state: these are the free variables that set *μ*_{e}. That is, we should properly write

where *x* is not specified *a priori*, but is instead inferred *a posteriori* based on the value of *N*_{e} in potential-controlled calculations of the bare electrode $*$ and the hydrogen-adsorbed electrode ($H*$). The ^{(1−x)+} superscript refers not to the charge localized on the final state, but to the net charge difference between the two states. This non-integer electron transfer has gone by various names in the previous literature, including “electrosorption valency” and “partial charge transfer”.^{11,26–31}

To be clear, this does not imply that electrons are not quantized—they certainly come in integer quantities. If *x* is found to be 0.9 electrons in the above reaction, we would physically expect the potentiostat to hold the potential (or work function) constant by injecting, on average, nine electrons every ten times this reaction occurred to an electrode that might have 10^{15} sites per cm^{2} of surface. However, if we examine a complete reaction—i.e., a full catalytic cycle—then we certainly should expect the total number of electrons transferred to be an integer. For example, if the above step is followed by a second protonation to liberate H_{2} gas, leaving the surface in its initial state, then exactly two electrons must be involved. This is also what we infer from grand-canonical electronic structure calculations since the surface will be in an identical state before and after the reaction.

Can this non-integer electron transfer be measured experimentally? Recently, a team (including authors of this work) showed strong agreement between grand-canonical electronic structure calculations and spectroscopic experiments that involved the protonation of a surface-tethered 4-mercaptobenzoic acid (4MBA).^{32} SJ calculations suggested that about −0.1 electrons should accompany this protonation at the constant potential. Experimentally, this is manifested by a potential-dependence of p*K*_{A} since the protonation reaction

involves a partial electron whose chemical potential is a function of *ϕ*, the applied potential. This implies that p*K*_{A} should vary with potential as − *xe*/(*k*_{B}*T* × ln 10), where *e* is the (positive) charge of an electron and *k*_{B} is the Boltzmann constant. This gives a slope of about −1.7 V^{−1}, which is very close to what was observed experimentally. The SJ calculations were also in close agreement with the related compound 2MBA,^{33} in which the charged group on the MBA is closer to the surface, resulting in a larger charge transfer of about −0.6 electrons upon protonation.

We next turn our attention to reaction barriers, calculated at a specified potential. Consider the initial state of the proton-discharge reaction [Eq. (3)]: the unit cell contains two subsystems, both of which are charged; the acidic water layer is electron-deficient and the electrode surface may have a net negative or positive charge depending on the potential. In the course of this elementary reaction, the proton discharges from the acidic water layer, adsorbing onto the surface, forming an adsorbed hydrogen atom. We show the reaction pathway for this reaction in Fig. 4. As we noted elsewhere,^{34} these calculations are adiabatic and do not introduce discontinuities in the charge transfer profiles; the electron does not “jump” in integer increments. Instead, the electron transfer occurs smoothly across the interior images in the NEB calculation. We recently reported^{15} that while the electron transfer process has a finite width along the reaction coordinate, the inflection point of the electron transfer profile tends to coincide with the saddle point over a wide range of potentials at least for this simple reaction. Here, 0.63 electrons are injected to keep the applied potential constant.

Intuitively, one might expect the magnitude of charge transfer to be close to, but not precisely one, keeping in mind the discussion in the previous paragraphs. However, the charge transfer magnitude is significantly, and consistently, less than unity; grand-canonical methods and Bader charge analyses^{35} typically reveal charge separations on the order of 0.5–0.7 electrons.^{8,15,36–39} This fractional charge transfer was long thought to be an artifact of self-interaction^{40–47} in Kohn–Sham DFT^{48,49}—the interaction of an electron (or basis function thereof) with itself via the electron density in the Hartree potential—but Chen *et al.*^{39} have recently suggested that this is a consequence of hybridization between solvated protons and the electron density of the electrode. Chen *et al.* reported that charge transfer is fractional at different levels of theory—including hybrid methods with fractional and full exact exchange—as long as the solvated proton is in close vicinity to the electrode. However, they reported that the charge transfer approaches unity when the proton is moved farther away from the electrode. The extent of charge transfer is thus dependent on the spatial separation of the two subsystems, and the energetics of hybridized states need not be identical to those of the bulk.

Here, we investigate the effect of the applied potential on the magnitude of charge transfer. Figure 5 shows the reaction pathway and electron-transfer profile for the Volmer reaction on Pt(111)/1ML over a potential window of 1 V. The magnitude of charge transfer decreases by 0.12 electrons when the potential is swept from *ϕ* = 0.2 V_{SHE} to *ϕ* = −0.8 V_{SHE}; this corresponds to an 18% reduction in charge transfer. Hence, the change in reaction energy of interfacial processes is not strictly linear with respect to the applied potential, but exhibits some curvature over wide potential ranges.

The origin of this phenomenon is shown in Fig. 6, where the extent of electron transfer is plotted relative to the Pt–H^{+} distance. The data in Fig. 6 cover the range −0.8 ≤ *ϕ* ≤ 0.2 V_{SHE}. First, we note that these potentials are more negative than the experimentally determined potential of zero charge of 0.3 V_{NHE},^{50} and a hexagonal H-down water geometry is more favorable than an H-up geometry. At *ϕ* = 0.3 V_{SHE}, we find that the water reorients into a mixed H-down/H-up geometry. As the electrode becomes more negatively charged, electrostatic forces cause the water bilayer to approach the surface. When the distance between the two subsystems decreases, the degree of overlap increases, and the increased overlap between the electron density of the electrode and the solvated proton of the initial state decreases the magnitude of charge transfer over the course of the reaction. We should carefully note that the magnitude of such an effect can be expected to have been amplified by the periodic nature of the simulation cells employed and the lack of buffering ions; larger-size simulations will be necessary to explore this effect in more detail.

In a recent study,^{15} we demonstrated that the potential-dependence of electrochemical reaction barriers can accurately be described with a curved scaling model developed from parabola intersections, as shown in Fig. 7 (top). In this model, the sole free parameter *b* corresponds to the reaction barrier at thermoneutral conditions, that is, the potential that makes Δ*E* = 0. The functional form of the scaling is then just the intersection of these parabolas, given by

This model not only captures the reaction barriers at intermediate driving forces but also predicts the asymptotic approach to activationless and barrierless discharge for high and low magnitude driving forces, respectively. Importantly, it also captures the curved nature between these limits and fits all these criteria with only a single adjustable parameter. The scaling is shown in Fig. 7 for the Volmer reaction on Pt(111)/1ML. The reaction barrier approaches the barrierless region for low-magnitude applied potentials (Δ*E*^{ϕ} > 4*b*) and the activation less region for high-magnitude driving forces (Δ*E*^{ϕ} < − 4*b*).

## REACTION THERMODYNAMICS AT CONSTANT POTENTIAL

Next, we turn our attention to reaction thermodynamics in which we focus on the calculation of Δ*G* for overall reactions, elementary steps, and barriers. We emphasize that from the standpoint of reactions, electrons are no different from other reactants or products—if they are reactants, they need to be supplied to the system, and if they are products, they must be removed, both at a specified free energy. However, unlike other reactants and products, electrons need not be supplied in integer quantities, as discussed above. See Fig. 8, showing an example of the reduction of gaseous CO to adsorbed CHO; in this case, CO, a proton, and *x* electrons must be supplied to the system. The value of *x* can be inferred from constant-potential electronic structure calculations.

### The decoupled computational electrode and the thermodynamics of barriers

The computational hydrogen electrode (CHE)^{51,52} is often used in elementary reactions in order to access the free energy of a proton/electron pair from that of H_{2} gas. In the CHE scheme, when the chemical potential of a proton/electron pair is encountered (as *μ*[H^{+}] + *μ*[e^{+}]), it is substituted with the thermodynamically equivalent chemical potential of hydrogen gas (as $12\mu [H2]$), which is true at the reference potential and pH and can be thermodynamically adjusted to other potentials and pH’s. The CHE scheme is thus designed for the net energy change across elementary steps involving coupled integer proton/electron transfers and is silent on the treatment of reaction barriers.

When the electron is non-integer or not directly coupled with a proton, such a relation can still be employed in which the proton and electron are decoupled.^{15,53} We refer to this as the “decoupled computational electrode” (DCE) model,^{15} which decouples the proton and electron chemical potential while enforcing thermodynamic consistency. Here, the energy of *any* state along a reaction trajectory, be it a thermodynamically stable state or a saddle point, is referenced to the correct reference state—a proton solvated in the bulk of the electrolyte—and enforces thermodynamic consistency. For example, consider the reaction of a proton with an electrocatalytic surface in order to create the barrier of the Volmer reaction,

(where the final state bears a net charge of (1 − *x*) + relative to the initial state; we omit this for brevity of nomenclature). The change in Gibbs free energy *G* associated with this partial reaction is given as

In the DCE scheme (like in the CHE scheme), the chemical potential of the proton is replaced by the equilibrium relation of the desired reference electrode. If we employ the reversible hydrogen electrode (RHE) scale ($H++e\u2212\u27f712H2$, which is in equilibrium at 0 V_{RHE} at any pH and temperature), then

A simple substitution of this into the previous equation leads to

where we have used *μ*^{0V}[e^{−}] − *μ*^{ϕ}[e^{−}] = *eϕ*. The free energy of the initial and final states, at potential *ϕ*, is constructed from the potential energy calculations (*A*_{pot} or Φ_{pot} to create *G*^{ϕ} of Ω^{ϕ} ≡ *G*^{ϕ} − *N*_{e}*μ*_{e}). If the grand free energy Ω^{ϕ} is used, then the equation simplifies to

Thus, if grand free energies Ω^{ϕ} are employed, which include the electronic correction, the DCE formulation is nearly identical to the traditional CHE approach, and becomes identical in the limit of integer electron transfer. This free energy difference is apparent in Fig. 4 as the difference between the state marked *G*^{‡} and that marked *G*_{IS}.

### Methods

The electronic structure calculations were undertaken in the open-source code GPAW,^{22,23} employing the solvated jellium method as described by us in an earlier publication,^{8} and implemented in GPAW. Atomistic manipulations were carried out in ASE,^{54} the atomic simulation environment. All calculations used the RPBE^{60} exchange--correlation functional of Hammer, Hansen, and Nørskov, with the exception of those to create the potential traces of Fig. 1.

In simulations to create Fig. 1, a 3 × 3 × 3-atom periodic slab of Pt was used with a hexagonal water overlayer, as shown. The PBE^{55} exchange–correlation functional was used with Fermi–Dirac smearing of the electronic occupations at 0.1 eV and 4 × 4 × 1 *k*-point sampling. The solvent exclusion zone used van der Waals radii for all atoms, and the implicit solvent was only present on the explicit-solvent side of the slab, with a dipole correction present in the direction normal to the slab, to isolate the electronic perturbations to the solvated side of the slab. The implicit solvent was the simplified continuum solvent model of Held and Walter,^{56} with parameters of *u*_{0} = 0.180 eV, *ɛ*_{∞} = 78.36, *γ* = 18.4 × 10^{−3} Pa ⋅ m, and *T* = 298.15 K being used, as described in their works and the GPAW documentation.

Saddle points were calculated with the scaled and dynamic^{34} implementation of the nudged elastic band (NEB)^{57–59} method. Thermodynamically stable intermediates and saddle points were optimized to within 0.03 eV/Å, and the convergence scaling of chains-of-states is based on a displacement metric with a scaling factor, *α*_{s} = 1; full details are provided in Ref. 34. Structures used to deduce subtle trends, e.g., in Fig. 6, were optimized until the forces acting on all unconstrained atoms are below 0.01 eV/Å. Platinum electrodes were modeled as 3 × 3 × 3 Pt(111) surfaces with full or fractional hydrogen coverage in the fcc sites. The grid spacing is 0.15–0.18 Å, and the systems were sampled on a 4 × 4 × 1 *k*-point mesh. The bottom layer was fixed at the bulk-optimized lattice constant of 3.99 Å with the RPBE^{60} exchange–correlation functional. Free energy corrections are included in the harmonic limit at 298 K,^{52}

where Δ*E*^{ϕ} is the potential energy at the constant potential extrapolated to 0 K and subsequent terms come from partition functions obtained in normal-mode analyses.

## ACKNOWLEDGMENTS

The authors thank Alexandar Zaslavsky, Xin Lian, and Matthew Shinkar, all of Brown University, for valuable technical discussions and acknowledge the support from the U.S. National Science Foundation under Award No. 1553365.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Per Lindgren**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Georg Kastlunger**: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Andrew A. Peterson**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

_{2}: New theoretical insights from an improved electrochemical model

_{2}using grand canonical potential kinetics

_{2}evolution on Pt

_{2}over Ag using density functional theory and transport models

_{2}reduction and hydrogen evolution

_{2}electrocatalytic reaction: When details are key

_{2}

^{+}, He

_{2}

^{+}, LiH

^{+}, and Ne

_{2}

^{+}