Controlled Electrochemical Barrier Calculations without Potential Control

The knowledge of electrochemical activation energies under applied potential conditions is a prerequisite for understanding catalytic activity at electrochemical interfaces. Here, we present a new set of methods that can compute electrochemical barriers with accuracy comparable to that of constant potential grand canonical approaches, without the explicit need for a potentiostat. Instead, we Legendre transform a set of constant charge, canonical reaction paths. Additional straightforward approximations offer the possibility to compute electrochemical barriers at a fraction of computational cost and complexity, and the analytical inclusion of geometric response highlights the importance of incorporating electronic as well as the geometric degrees of freedom when evaluating electrochemical barriers.


Proton transfer to Au(111) -atomic structures and potential dependence
To demonstrate our methods we model the proton transfer from a single hydronium ion surrounded by implicit electrolyte to a Au(111) surface.We find that the proton is transferred first to an (energetically unstable) atop position of the neighboring metal atom and subsequently comes to rest in the nearby hollow adsorption site.The structures of initial, transition, and final state (IS, TS, FS) and their dependency on electrode potential U are shown in Fig. S1.While IS and TS respond to changes in the electrode potential, there is vanishing response for the FS.We use the position along the 3N-dimensional path of the reaction from IS to FS as reaction coordinate as implemented in the Atomic Simulation Environment (ASE) 1 Nudged-Elastic-Band (NEB) package.

S-2
All density functional theory (DFT) calculations were performed using the ASE 1 and the DFT code GPAW 2,3 in combination with the Solvated Jellium Model (SJM) 4 for simulating electrochemical interfaces.We employ the BEEF-vdW exchange-correlation functional 5 and a real space grid with a spacing of below 0.17 Å.This results in a 48 × 48 × 192 grid for the studied (hexagonal) 3 × 3 Au(111) slabs with a thickness of 4 atomic layers, the bottom two of which were frozen in the bulk geometry that was optimized resulting in a lattice constant of 4.275 Å.The vacuum layer between the slabs has a thickness of at least 20 Å.The Brillouin zone is sampled using an even-numbered Γ-centered Monkhorst-Pack grid with a maximum distance of the K-points of 0.15 Å −1 in lateral directions and a single K-point in z-direction yielding a 6 × 6 × 1 K-point grid.
For the bulk calculation we use a 18 × 18 × 18 k-point grid.Pulay Mixing is applied (linear mixing coefficient β = 0.05, at maximum with 3 previous densities, and a weight of 50.0 for the metric) and Fermi-Dirac smearing on the occupations with a width of 0.1 eV.
The electrolyte is modelled by an implicit solvent (linear dielectric constant of 78.36, surface tension of 1.1484 • 10 −3 eV/Å 2 ) where the cavity of the quantum-mechanical region is defined by an effective power12potential with the default potential strength (0.18 eV) and a temperature of 298.15 K.We scaled the Van-der-Waals radii of the gold atoms by a factor of 1.08 yielding a capacitance of the clean Au(111) slab of 19 µF/cm 2 and a PZC of 4.7 V on an absolute scale which is in close agreement to the experimental values of ≈ 15 − 30 µF/cm 2 6 and ≈ 4.9 V 7,8 .The jellium counter charge in the electrolyte above the slab used for both potential-controlled grand canonical as well as charged canonical calculations is inserted at a distance of 3 Å to the highest atom and 1 Å below the dipole-correction layer.The deviation of the electrode potential from the target potential in potential-controlled calculations is converged to below 5 mV and to 1 mV for the grand canonical Hessian calculations discussed below.
Structure optimizations are performed with the BFGS algorithm (convergence of forces below 0.03 eV/Å) and transition state searches with the dynamic Nudged-Elastic-Band (dyNEB) method, [9][10][11] both as implemented in ASE.The dyNEBs are first pre-converged then refined using the climbing-image method 12 (forces below 0.05 eV/Å).When used as the central point for Hessian calculations, we converged until the atomic forces at resting states are below 0.01 eV/Å and at transition states below 0.03 eV/Å.Hessian calculations use a displacement of 0.1 Å, and the central difference method between positive and negative displacement directions.We performed the analyses both taking into account all non-metal atoms and all metal atoms of the first layer as well taking only into account the Au atom involved in the reaction and the non-metal atoms yielding essentially identical results.
For all figures showing absolute energies, we choose the grand canonical energy of the initial state at q = 0 as the reference energy.

Procedures
In the following we give a short schematic description of the various methods discussed in the main text.

Potentiostated constant potential calculations
• Perform all calculations at constant potential conditions using a potentiostat.
• For every potential U i of interest perform the following steps: -Relax the initial and final state (IS and FS) to their resting positions.
-Converge a dyNEB transition state search from IS to FS, first without the climbing image algorithm and subsequently with the climbing image algorithm.
-Obtain the grand canonical energy E * i and the calculated excess charge q * i for IS, TS, and FS.If the grand canonical energy is not directly provided by the DFT code, derive it from the canonical energy E * i , the potential U i and the excess charge q * i for every state using Eq. 2 in the main text.

S-4
-The difference in energies E * i between the respective states yields the forward and reverse activation barriers Interpolation / Extrapolation of grand canonical energies from canonical calculations • Perform all calculations at constant charge conditions.
• For every excess charge q i of interest perform the following steps: -Relax the initial and final state (IS and FS) to their resting positions.
-Converge a dyNEB transition state search from IS to FS, first without the climbing image algorithm and subsequently with the climbing image algorithm.
-Obtain the canonical energy E * i and the calculated electrode potential U * i for IS, TS, and FS.
-Derive the grand canonical energy of each state for every state (Eq. 2 in the main text).
• Perform an inter-or extrapolation of E * i in dependence of U * i for each of the states IS, TS, and FS, yielding a function E * (U).If desired, use the slope dE dU (U * i ) = −q i as additional information for inter-/ extrapolation.
• Derive the forward and reverse activation barriers as a function of U using

Methods based on s single NEB calculations
All of the single NEB based methods follow the same procedure: • Perform all calculations at constant charge conditions at excess charge q 0 .
• Relax the initial and final state (IS and FS) to their resting positions. S-5 • Converge a dyNEB transition state search from IS to FS, first without the climbing image algorithm and subsequently with the climbing image algorithm.
• Obtain the canonical energy E * 0 and the calculated electrode potential U * 0 for IS, TS, and FS.
• Derive the grand canonical energy E * 0 = E * 0 −U * 0 q 0 for every state (Eq. 2 in the main text).
• Derive the respective approximation for the total capacitance C * total associated with each method as described below.
• If using q 0 = 0, i.e.U * 0 = U * PZC , the grand canonical energy E * (U) of each state in dependence of U is given by (cf.Eqs.6 of the main text): If using q 0 ̸ = 0, either compute the main text) and use aboves formula or directly use that in a parabolic approximation the grand canonical energy E * (U) of each state in dependence of U is given by: 13 The different methods -single capacitance approximation (SC), electronic capacitance approximation (EC), electronic + geometric capacitance approximation (EGC), and the electronic + geometric capacitance approximation along the reaction path (Eξ GC) only differ in how the capacitance is derived: Single capacitance approximation (SC) • Use an appropriate guess for C * total , e.g.compute it for a certain state of the system.In order to do so, perform electronic single point calculations at several charges q i for the chosen state and obtain the electrode potentials U * i .
• Derive the capacitance as the slope of a linear fit of q i with respect to U * i : Electronic capacitance approximation (EC) • For each of the considered states IS, TS, and FS of the system perform electronic single point calculations at several charges q i and obtain the electrode potentials U * i .
• For each state, derive the capacitance as the slope of a linear fit of q i with respect to U * i : Electronic + geometric capacitance approximation (EGC) The total capacitance is given by the sum of electronic and the geometric capacitance contributions C * el and C * geom , respectively (cf. Eq. 8 in the main text or Ref. 13): Derive the electronic contribution C * el as described in the EC approximation above.
For the geometric contribution C * geom (using purely canonical calculations) perform the following steps for each of the considered states IS, TS, and FS: • Compute the canonical Hessian H * and the corresponding gradient in electrode potential − → ∇U * of the relevant part of the system (e.g. the metal atoms in the first layer of the electrode and any non-metal species) at the excess charge q 0 .For this, choose the optimized structure of each state as the center structure (from relaxation or from the transition state search).Then, for all of the relevant coordinates r i , i.e. for each of the coordinates of the relevant atoms: -Displace the system along r i by a certain amount d in both positive and negative directions (+), (−).
-Perform an electronic single point calculation for both displacements.
-Obtain the resulting forces F +,− j along each of the relevant coordinates r j of the system and obtain the resulting electrode potentials U +,− caused by the positive and the negative displacement.
-The value H * i, j of the Hessian and the value of ( − → ∇U) * i are given by: • Hessians, by definition, are symmetric matrices.In order to eliminate numerical noise, symmetrize the Hessian by replacing H * i, j with 1 2 (H * i, j + H * j,i ) from the just derived values.
• Compute the grand canonical Hessian H * via Eq. 10 in the main text: using the electronic capacitance derived above for the considered state.
• Calculate the inverse of the grand canonical Hessian H * −1 by matrix inversion of H * .
• Calculate the geometric capacitance contribution by: Similarly, we could use potentiostated grand canonical calculations, directly calculate H * and − → ∇ q * (analogous to H * and − → ∇U * above), and use Eq. 9 in the main text.Note however, that in general, one has to use a very high potential-accuracy for the potentiostat, if one wants to determine small contributions to H * and − → ∇ q * accurately enough.A detailed analysis of Hessians and geometric capacitances for the considered states and for the present system are reported below in Sections 7 S-8 and 9 of the SI.
One might wonder, why we do not provide a simple closed formula for the geometric capacitance derived from canonical properties only.The reason is that the inverse of the grand canonical Hessian expressed in canonical properties yields a rather complicated expression.We provide a closed formula and an extensive derivation and discussion in Ref. 13.
Electronic + geometric (only reaction path) capacitance approximation (Eξ GC) For deriving the capacitance in the Eξ GC approximation along the reaction path at the TS (using purely canonical calculations), first compute the electronic capacitance as described in the EC approximation above.
Then perform the following steps using the images along the reaction energy path obtained from the dyNEB at charge q 0 : • Interpolate the canonical energy E(ξ) along the reaction path ξ using a cubic Hermite spline taking into account the projected forces (e.g. as implemented in the ASE 1 NEB module).
• Obtain the position of the transition state ξ TS along the reaction path ξ .
• Interpolate the electrode potential U(ξ) along the reaction path ξ in a smooth way (e.g. a cubic spline).
• Fit a parabola (aξ 2 + bξ + c) to E(ξ) in the vicinity of the transition state, e.g. in the vicinity of ξ TS .The obtained curvature of the parabola yields H TS ξ ,ξ = 2a.
• Fit a linear function (mξ + n) to U(ξ) in the vicinity of the transition state, e.g. in the vicinity of ξ TS .The obtained slope of the linear fit yields ∂U ∂ ξ • Calculate the corresponding value of the grand canonical Hessian via . yielding H TS ξ ,ξ and ∂ q ∂ ξ (similar as the canonical NEB calculation yields H TS ξ ,ξ and ∂U ∂ ξ ) that we could insert into Eq.15 of the main text.At least for the studied system, both approaches yield essentially identical results, as shown in Fig. S4.Note, however, that in general the normal modes (and hence the direction of the reaction path at the transition state) can slightly deviate in the canonical and the grand canonical ensemble since the Hessians differ (cf.Ref. 13 and Section 6 below for the math and Section 7, where we analyze in more detail normal modes and Hessians for the TS.).

Absolute one-dimensional grand canonical PES along the reaction path
In the following we shortly want to outline how we created Fig. 1 of the main text which illustrates the one-dimensional grand canonical PES along the reaction path.
• Perform a constant charge NEB calculation at q = 0.
• Freeze the obtained geometries, the NEB images, along the reaction path and perform several charged single point calculations for every image.Essentially evaluate the NEB without further geometric optimization at various charges q i = −1.111,−0.999, ..., 0.999, 1.111.For all images determine their positions ξ j along the reaction path ξ .Obtain the canonical energies E i , the forces F ∥ i projected onto the reaction path, and the electrode potentials U i for the various q i .
• In order to take into account the information encoded in the projected forces, for every q i fit the canonical energy along the reaction path using a cubic Hermite spline taking into account the projected forces as implemented in the ASE 1 NEB module.
• For every q i evaluate the just derived canonical energy fit on a fine grid along the reaction path.Note that incorporating the projected forces at the various excess charges q i leads to a S-10 slight shift of the transition state.Essentially it is the canonical shift that we derive in Ref.
13 (considering only one spatial dimension).Note, that, hence, we take into account only the geometric shift in one dimension, neglecting a possible shift perpendicular to the q = 0 reaction path that might results from other q i ̸ = 0.The limitation to the purely one-dimensional case is why we consider this more a qualitative illustration.
• Perform a two-dimensional spline interpolation of E in dependence of q and ξ .This essentially yields the canonical PES E(q, ξ ).
• Perform a two-dimensional spline interpolation of U in dependence of q and ξ yielding U(q, ξ ).
• Define a function returning the grand canonical energy as E (U, ξ ) = E(q(U, ξ ), ξ )−q(U, ξ )U.This is the grand canonical PES as plotted in Fig. 1 of the main text.
• For the constant potential lines, evaluate E (U, ξ ) at constant U.The transition state position is where E reaches is maximum.
• For the constant charge lines evaluate E (U(q, ξ ), ξ )) at constant q.Since we are at constant charge conditions along this path, the transition state position is where where E(q, ξ ) reaches its maximum (i.e.not where E (U(q, ξ ), ξ )) reaches its maximum).While this might seem counter-intuitive at the first glance, it illustrates that a stationary point of the grand canonical PES at constant potential conditions is a stationary point in the canonical PES at constant charge conditions.

Including the reference state
We want to introduce one additional state, that we denote as the 'reference state' RS, which essentially is identical to the initial state, however with the H + that takes part in the reaction S-11 considered in the bulk electrolyte and the corresponding electron e -in the external electron bath.a For simplicity, in the main text of this work we focused on states of same composition (IS, TS, and FS), i.e. we do not consider a change in the number of atoms.As shown in previous work [16][17][18][19] , we can incorporate the effect of compositional changes by accounting for changes in the number of (possibly charged) adsorbates N a (i.e.non-substrate atoms) and similarly for changes in the number of (neutral) substrate atoms N s via the corresponding (electro-)chemical potentials μa and µ s .We will in the following, however, assume a constant number of substrate atoms and hence ignore their contribution.For a derivation including the energetics of substrate atoms see e.g.Refs.16 and 17.
Before we continue, note that the absolute excess charge q is not uniquely defined, i.e. in principle we could assign q to the total number of electrons of the simulated system.However, it is more sensible in an electrochemical context to consider q as the net excess charge of the system.
Essentially this means, we assign q = 0 to the charge-neutral state of the surface, where there is no counter charge present in the (continuum) electrolyte.With this in mind, if we add a charged species, e.g. an isolated ion with charge q a , to the system, we similarly have to add a compensating electronic counter charge −q a in order to obtain an overall charge neutral interface again, i.e. the absolute electronic excess charge q abs = q − N a q a has to compensate any charges of explicit ions and the charge of the electrolyte.The energy of N a species of type a together with their respective compensating electrons in the respective bulk reservoirs is hence given by: Note that this way of counting electrons does not couple ionic and electronic degrees of freedom since the net excess charge q is still an independent variable.Now, we want to compare the energy of the RS state to the energies of IS, TS, and FS.To do so, we can add the contribution E a (U) from Eq. 2 to the grand canonical energy of the reference state a The RS is sometimes also called the 'true initial state' in contrast to the 'pseudo initial state' that we just denote as initial state. 14,15nstead of computing the electrochemical potential of a proton H + in the bulk electrolyte, we make use of the Computational Hydrogen Electrode approach 20 where the energy of (H + + e -) in their respective reservoirs at an arbitrary absolute potential U (Eq. 2) and at pH = 0 is given as: 16 E H + + e -(U ) = 1 2 µ(H 2 (g)) + 4.44eV − eU , leveraging that the Standard Hydrogen Electrode (SHE) scale is related to the absolute electrode potential via: 21 U SHE = 4.44V .
We derive µ(H 2 (g)) by performing a vacuum calculation and computing the free energy in gas phase at a temperature of 300 K and a pressure of 1 bar via the ideal gas approximation as implemented in the thermochemistry module of ASE 1 .Note, that in this work we only discuss the results at pH = 0, however, we could also include pH effects by considering the Reversible Hydrogen Electrode potential scale (RHE) similar to the procedure in Ref. 16. Note, that instead of adding E H + + e -(U ) to the grand canonical energy of the reference state, we could similarly subtract it from IS, TS, and FS.Again, the specific choice is arbitrary and -as long as only considering relative energies -irrelevant.For the absolute energetics we show the additive referencing in Fig. S2 and the subtractive referencing in Fig. S3, respectively.Note, that we could equally consider the 'excess energy' of every state, i.e. the energy cost of creating the entire system from suitable reservoirs (see Ref. 17 for more details).S-15 6 Derivation of the geometric capacitance and the grand canonical Hessian In this section we summarize the derivation of the central equations Eqs.9-11 of the main text, which are extensively explained, derived, and discussed in Ref. 13.
We consider the grand canonical energy E a function of electrode potential U and system geometry − → r of the system: E (U, − → r ).We denote the potential dependent 3N-dimensional atomic coordinates of a specific geometrically stationary point, e.g.IS, TS, or FS, as − → r * (U), which can be interpreted as a smooth path in the potential window where the stationary point is energetically stable.The path is defined by the stationarity criterion − → 13 To increase clarity, for partial derivatives we explicitly indicate quantities held constant, e.g. by writing for the spatial derivative of a quantity X in direction of coordinate i at constant potential U. Note, that − → r * (U) is only a function of one variable, the electrode potential U, so partial and total derivative are identical.
The total derivative of E with respect to electrode potential U for a stationary point is given by the excess charge q: using the chain rule, the stationarity condition yielding the potential dependent change of the path − → r * (U).
By differentiating Eq. 3 a second time while inserting Eq. 4) we can derive the total capacitance of the system as the negative total curvature of the grand canonical energy of a stationary point with respect to electrode potential: 13 where we used our definition of the purely electronic capacitance ) and in the last step Schwarz's theorem: By defining the geometric capacitance as C * geom = C * total −C * el and by incorporating the dependence on the potential-dependent positions of the stationary state − → r * (U), i.e. by rewriting a quantity X that depends on (U, − → r * (U)) as a new quantity "X * (U)", Eq. 5 reads in matrix notation: 13 which is essentially Eq. 9 in the main text.
As written in the main text and discussed in Ref. 13, it is possible to express the grand canonical Hessian H in terms of the canonical Hessian H by starting with: Using the equivalence of grand canonical and canonical forces ), which directly follows from the relation between grand canonical energy E and canonical energy E given by the Legendre transform 13,15 this yields: 13

Using the definition of the canonical Hessian
(q, − → r ), the relation of canonical energy, excess charge and electrode potential S-18 this reads: 13 where we used matrix notation and the outer product notation − → ∇U − → ∇ q T and for readability skipped the dependencies in the last line.
Finally, using that we can express the gradient of the electrode potential − → ∇U in terms of using the triple product: 13 which is Eq. 10 of the main text.

Potential dependence of the geometric capacitance
The geometric contribution to the total capacitance principle is a potential dependent quantity (see Ref. 13 and derivation in the previous section).However, we find S-20 it to very good approximation to be constant in the considered potential range as shown in Fig. S5, where we compare the values computed using: • the procedure explained in the main text, based on canonical quantities only (gc from canonical, dark blue circles in Fig. S5), • the direct computation using a grand-canonical, potentiostated Hessian calculation with a potential tolerance of 1 mV and highly accurate converged electronic SCFs (work function convergence to below 0.1 mV) (gc high accuracy, orange diamonds), • as well as a direct computation using grand-canonical, potentiostated Hessian calculation with default electronic SCF convergence criteria (gc, light green diamonds).
The method based on canonical calculations can reproduce the values from the potentiostated grand canonical calculations at comparable accuracy.
As an additional note, it is most important to include the non-metal atoms and only the metal atoms that take part in the reaction in the Hessian calculations (top panel vs bottom panel in Fig. S5).
Note, that the potential dependence might be more significant for more complex systems.that we can also rewrite using Schwarz's Theorem Eq. 6 in terms of − → ∇ q: where we used the matrix notation and dropped the dependence of H and − → ∇ q on (U 0 , − → r * (U 0 )).
Note, that the relevant energetic contributions, i.e. the geometric capacitance contributions C * geom (U) (Eq.7), result solely from the shifts ∆r * ∥∇q in direction of To validate our results, and to illustrate the geometric shift, we hence compute the shift ∆r * ∥∇U for the transition state structure and compare it to the observed shift obtained from actual NEB calculations.Fig. S6 shows the results obtained using the Hessian from the high and normal accuracy potentiostatted calculations (orange and green lines) as well as the ones obtained from S-22 deriving the relevant quantities from purely canonical calculations (blue lines) in comparison to the geometric shifts along − → ∇U obtained from canonical NEB calculations (black dots) using the geometries of the calculation at U 0 = 4.06 V as reference point.
The agreement is remarkable -independent of the underlying Hessian calculations.Note, that the noticeable scatter of the NEB data points is caused by numerical noise originating from the limited accuracy of structure optimizations / NEB calculations.We can roughly estimate the noise δ r d of a structure optimization with a certain force threshold δ F (for the TS structure here: For this we compute the curvature H d of the PES by evaluating the grand canonical Hessian H along this direction: We indicate this estimated accuracy by black error bars around the NEB data points in Fig. S6 confirming the very good agreement between analytically derived and computed geometric shifts.
9 Potential dependence of normal modes, Hessian eigenvalues, and vibrational energy contributions In the main text, we focus on the grand canonical potential energy surface E and discuss its relation to the canonical energy E. Considering that we just showed in the previous sections that the grand canonical Hessian H in principle is a potential dependent property and that it might deviate from the canonical Hessian H, the interested reader might ask, how large the potential dependence is, how large the deviations between grand canonical and canonical Hessian are, and whether they affect vibrational contributions to the energy, i.e. zero point energy ZPE and entropic contributions.
While a more extensive analysis and discussion will certainly be part of future work, we touch on this subject briefly in the following.Grand canonical Hessians H are computed from canonical quantities via Eq. 10 and Eq.11 of the main text (gc from canonical, light green squares), from high accuracy potentiostated grand canonical Hessian calculations as described above (gc high accuracy, orange diamonds), and from normal accuracy potentiostated grand canonical Hessian calculations (gc, blue diamonds).While the difference between canonical and grand canonical eigenvalues is significant, the three different approaches to compute the grand canonical Hessians yield the same result.
All quantities show a similar potential dependence.
but only the vibrational modes orthogonal to it.
While we observe hardly any difference between the canonical and the grand-canonical results (blue vs. green/orange), we can observe a potential-dependence for the IS of around 60 meV/V.
However, for the kinetically more relevant states -the FS, the TS and the RS -the potential dependence is almost vanishing.As a result we find that the vibrational free energy contributions at FS, TS and RS can be considered additive constants as mostly done in electrochemical simulations.
Note that the vibrational free energy contributions at TS and FS are nearly identical thus suggesting the their inclusion would not alter their relative energetic difference aka the relevant electrochemical barrier.Note also that the difference in ∆F between the FS (water molecule & proton at the interface) and the RS (water molecule without proton at the interface) is approximately 0.12 eV, which is very close to the typical estimates of the vibrational free energy correction of (only) an adsorbed H, e.g.
as estimated by considering a zero point energy ∆E ZPE H = 0.17 eV in Ref. 22.However, for more The IS exhibits a slight potential dependence, for the other states the contributions are to good approximation constant.The drop of RS and FS at the highest potentials (close to U PZC , q = 0 e) originates from a weak rotation-like normal mode of the remaining water molecule, that at other potentials is suppressed by interaction with the charged surface.For the other states this mode is as well suppressed by interaction with the substrate.
detailed analyses and possibly more complex systems one might require a detailed analysis of the potential dependence.

Figure S1 :
Figure S1: Atomic structures of IS, TS, and FS (left, middle, and right panels) and their dependency on electrode potential U (inserts) indicated by gray overlays from two different perspectives (top and bottom panels).Also given as inserts are the respective excess charges q.

Figure S2 :
FigureS2: Grand canonical energetics of RS, IS, TS, and FS as a function of absolute electrode potential derived using the various methods discussed in the main text.For the reference state, we add the energy E H + + e -(U ) of H + + e -in their respective bulk reservoirs.Left: Inter-/Extrapolated grand canonical energies derived from eight constant charge climbing image NEBs for q = −0.555,−0, 444, ..., 0.222e (from left to right) in comparison with potentiostat results (diamonds).Cubic Hermite spline interpolation as solid line, quadratic extrapolation as dashed lines.Right: Single NEB calculation based approaches -SC, EC, and EGC approximation (dotted, dash-dotted, and dashed lines) -in comparison to the extrapolated (light solid lines) and reference potentiostat results (diamonds).The electronic SCF calculations used for determining the electronic capacitance of the EC approach are shown as crosses.Note, that the reference state RS without the additional adsorbate related term E H + + e -(U ) had its PZC roughly around the PZC of the FS.

Figure S3 :
Figure S3: Same as Fig. S2, but subtracting the energy E H + + e -(U ) of H + + e -in their respective bulk reservoirs from IS, TS, and FS instead of adding it to the RS.

Figure S4 :
Figure S4: Comparison of kinetic activation barriers as a function of the absolute electrode potential.Forward barrier with respect to the RS (top panel) and with respect to IS (middle panel), and reverse barrier with respect to FS (bottom panel).Multiple-NEB calculation based methods (blue solid lines): interpolation (dark), extrapolation (light), both based on parabolic fits.Single-NEB methods (orange lines): EGC (electronic + geometric, dashed), EC (electronic only, dash-dotted), and SC (single capacitance, dotted).Reference potentiostat data is shown as dark blue diamonds.Note the strong influence of geometric effects, i.e. the difference between EGC (dashed) and EC approximation (dash-dotted) which are to large extent originating from contributions along the reaction path ξ (Eξ GC, bright green lines).Also shown here is the Eξ GC approximation derived from a potentiostated grand canonical NEB calculation (gc NEB) as dark green lines which yields the same result.The PZCs of the relevant states U * PZC are drawn as vertical lines.
− → r ) = −q(U, − → r ) defines the excess charge of the system.Furthermore, since the stationarity is valid along the entire path − → r * , i.e. − → F * = 0 = const.and hence d dU − → F (U, − → r * (U)) = 0, we can write using the grand canonical Hessian H with H

Figure S5 :
Figure S5: Geometric capacitance contribution C * geom for IS, TS, FS, and RS in dependence of electrode potential U. Shown are values based on canonical Hessian calculations (dark blue circles) as discussed in the main text and based on direct grand canonical, potentiostated Hessian calculations with normal and very accurate electronic SCF accuracy (light green and orange diamonds) as explained in the text.In the top panel we show the results when including the entire first Au layer and the non-metal atoms in the Hessian calculations and in the bottom panel we show the results when only including the Au atom to which the proton is transferred and the non-metal atoms in the Hessian calculations.

8
Potential dependence of stationary point geometries The essence of above derivation of the geometric capacitance contribution originates from slight changes in the stationary state geometry as a response to a change in electrode potential ∆U = U −U 0 .Certainly the energetic effects are most interesting, however, we can also evaluate the geometric changes themselves.By integrating Eq. 4 from an electrode potential U 0 to U = U 0 + ∆U, assuming a negligible change of H and ∂ − → F ∂U r , we obtain the absolute shift ∆ − → r * (∆U) caused by the potential change ∆U accurate to linear order:

Figure S6 :
Figure S6: Computed and analytically derived geometric shift ∆r * ∥∇U of the transition state along − → ∇U for the TS with respect to the geometry at U 0 = 4.06 V. Shown are the NEB reference points (black) from canonical NEB calculations and the analytical results using Eq. 10 based on high and normal accuracy potentiostatted Hessian calculations (orange and green lines) and derived from purely canonical Hessian calculation data (blue line).The estimated accuracy of the geometries, as discussed in the text, is indicated by black error bars.
which in our case yields H d = −1.3eV/Å 2 and use common error propagation considering only the dimension along − → d : δ F d = ∂ F d ∂ r d δ r d = |H d | δ r d .The spatial accuracy of the TS geometry δ r d along − → ∇U, hence is given by δ r d = δ F d |H d | = 0.03 eV/Å 1.3 eV/Å 2 = 0.02Å.

Figure S8 :
FigureS8: Eigenvalues of the pure, i.e. not mass weighted, canonical and grand canonical Hessians H and H for their normal modes along the reaction path at the transition state for various electrode potentials U. Canonical values are shown as blue circles.Grand canonical Hessians H are computed from canonical quantities via Eq. 10 and Eq.11 of the main text (gc from canonical, light green squares), from high accuracy potentiostated grand canonical Hessian calculations as described above (gc high accuracy, orange diamonds), and from normal accuracy potentiostated grand canonical Hessian calculations (gc, blue diamonds).While the difference between canonical and grand canonical eigenvalues is significant, the three different approaches to compute the grand canonical Hessians yield the same result.All quantities show a similar potential dependence.

Figure S9 :
FigureS9: Helmholtz free energy contribution ∆F in dependence of electrode potential U for RS, IS, TS, and FS at T = 300 K considering the vibrations of the oxygen and the proton involved in the reaction.The difference between canonical (blue) and grand canonical results (high accuracy orange and normal accuracy green) is negligible.The IS exhibits a slight potential dependence, for the other states the contributions are to good approximation constant.The drop of RS and FS at the highest potentials (close to U PZC , q = 0 e) originates from a weak rotation-like normal mode of the remaining water molecule, that at other potentials is suppressed by interaction with the charged surface.For the other states this mode is as well suppressed by interaction with the substrate.