Electrochemistry from the atomic scale, in the electronically grand-canonical ensemble

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.

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₂ reduction,^6,16–18 formic acid oxidation,¹⁹ and CO adsorption.²⁰ In particular, these methods allow for the potential-dependent study of electrochemical reaction barriers—something otherwise accessible through computationally expensive extrapolation methods²¹—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)⁸ 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.²⁴ 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 (−μ_eN_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.²⁵ 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).

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¹⁵ sites per cm² 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₂ 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).³² 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 pK_A since the protonation reaction

involves a partial electron whose chemical potential is a function of ϕ, the applied potential. This implies that pK_A should vary with potential as − xe/(k_BT × 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⁻¹, which is very close to what was observed experimentally. The SJ calculations were also in close agreement with the related compound 2MBA,³³ 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,³⁴ 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¹⁵ 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³⁵ 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.³⁹ 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,⁵⁰ 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,¹⁵ 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^ϕ > 4b) and the activation less region for high-magnitude driving forces (ΔE^ϕ < − 4b).

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 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₂ 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μ[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,¹⁵ 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−⟷12H2$⁠, 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.

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,⁸ and implemented in GPAW. Atomistic manipulations were carried out in ASE,⁵⁴ the atomic simulation environment. All calculations used the RPBE⁶⁰ 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⁵⁵ 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,⁵⁶ with parameters of u₀ = 0.180 eV, ɛ_∞ = 78.36, γ = 18.4 × 10⁻³ 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³⁴ 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⁶⁰ exchange–correlation functional. Free energy corrections are included in the harmonic limit at 298 K,⁵² 

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.

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.

The authors have no conflicts to disclose.

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).

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