Grand canonical rate theory for electrochemical and electrocatalytic systems I: General formulation and proton-coupled electron transfer reactions

. Electrochemical interfaces present a serious challenge for atomistic modelling. Electrochemical thermodynamics are naturally addressed within the grand canonical ensemble (GCE) but the lack of ﬁxed potential rate theory impedes fundamental understanding and computation of electrochemical rate constants. Herein, a generally valid electrochemical rate theory is developed by extending equilibrium canonical rate theory to the GCE. The extension provides a rigorous framework for addressing classical reactions, nuclear tunneling and other quantum eﬀects, non-adiabaticity etc. from a single uniﬁed theoretical framework. The rate expressions can be parametrized directly with self-consistent GCE-DFT methods. These features enable a well-deﬁned ﬁrst principles route to address reaction barriers and prefactors (proton-coupled) electron transfer reactions at ﬁxed potentials. Speciﬁc rate equations are derived for adiabatic classical transition state theory and adiabatic GCE empirical valence bond (GCE-EVB) theory resulting in a Marcus-like expression within GCE. From GCE-EVB general free energy relations for electrochemical systems are derived. The GCE-EVB theory is demonstrated by predicting the PCET rates and transition state geometries for the adiabatic Au-catalyzed acidic Volmer reaction using (constrained) GCE-DFT. The work herein provides the theoretical basis and practical computational approaches to electrochemical rates with numerous applications in physical and computational electrochemistry.


Introduction
Electrochemical reactions and especially electrocatalysis are at the forefront of current green technologies to mitigate climate change. To realize and utilize the full potential of electrocatalysis, selective and active catalysts are needed for various applications and reactions including e.g. oxygen and hydrogen reduction/evolution reactions, nitrogen reduction to ammonia and CO 2 reduction. [1] These and other electrocatalytic/electrochemical reactions are based on successive proton-coupled electron transfer (PCET), electron transfer (ET), and proton transfer (PT) reactions; the unique aspect of electrochemistry is the ability to directly control PCET, ET, and PT kinetics and thermodynamics by the electrode potential. [2] Besides the catalyst material, electrocatalytic performance is controlled by the electrolyte composition and electrode potential. To translate these to microscopic, computationally treatable quantities, it is the combination of the electrolyte and electron electrochemical potentials which determine and control the (thermodynamic) state of electrochemical systems. Therefore, an atomic-level computational model needs to provide an explicit control and description of these chemical potentials as depicted in Figure 1. In thermodynamics fixing the chemical potentials is achieved through a Legendre transformation from a canonical ensemble to a grand-canonical ensemble (GCE) for both electrons and nuclei. [3] This calls for theoretical and computational methods to treat systems where the particle numbers are allowed to fluctuate and the chemical potentials are fixed.
The theoretical basis for fixed potential electronic structure calculations was developed by Mermin who formulated electronic density functional theory (DFT) within GCE. [4,5]. Later, GCE-DFT has been generalized for treating nuclear species either classically or quantum mechanically [3,6,7,8,9]. The GCE-DFT provides a fully DFT, atomistic approach for computing free energies of electrochemical and electrocatalytic systems at fixed electrode and ionic/nuclear chemical potentials. [3] Importantly, the free energy from a GCE-DFT calculation is in theory exact and unique to a given external potential. In practice, the (exchange-)correlation effects in both quantum and classical systems need to be approximated. The thermodynamic GCE framework has already been adopted by the electronic structure community to model electrocatalytic thermodynamics at fixed electrode [10,11,12,13,14,15,16,17,18,19,20,3] and ion potentials [3,14,12]. Based on the large number of theoretical and computational works utilizing GCE-DFT, the computational framework for thermodynamics within GCE seems generally accepted. The thermodynamic approach has provided fundamental atomic level insight on reactions at complex electrochemical interfaces and enabled computational catalyst screening using free energy relations, Volcano curves, and scaling relations. [1]  However, it has been shown that a purely thermodynamic perspective on electrocatalysis is not sufficient for understanding and predicting activity, selectivity, or catalytic trends. [21,22,23,24] Besides applications in catalysis and material science, electrochemical kinetics are fundamentally important and provide a way to understand complex solvent effects, electron and nuclear tunneling, and non-adiabatic reactions. Ideally both fundamental and applied kinetic computational/theoretical studies should make use of general and self-consistent first principles Hamiltonians within GCE. This has, unfortunately, remained unattainable due to theoretical and methodological difficulties and omissions. [25] Surprisingly, a general GCE rate theory has not yet been established; mending this deficiency is the central goal of the present work.
Before diving to the development of the GCE rate theory, it is worth considering what new and important information can be obtained from a general electrochemical rate theory. First and foremost, the theory needs to accurately capture the intricacies of ET, PT, and PCET reactions as function of the electrode potential. Therefore, a general treatment of electrochemical reaction rates needs to be applicable to 1) both inner-sphere and adiabatic as well as outer-sphere and non-adiabatic reactions, 2) sequential ET/PT or decoupled reactions as well as simultaneous PCET reactions, 3) tunneling of both electrons and nuclei, and 4) be combined with general first principles GCE Hamiltonians. The motivation for including each of these four requirements is discussed next.
First, adiabatic inner-sphere reactions present a large and important class of electrocatalytic reactions as demonstrated by a large body of computational works aiming to evaluate rate constants for this class of reactions [11,12,20,19,21,26,27,28,29,30,31,32]. However important adiabatic reactions are, all electrochemical reactions are certainly not innersphere nor adiabatic. In particular, both vibronic and electronic non-adiabatic effects are frequently encountered in outer-sphere and long-range ET, PT, and PCET reactions. [25,33] Even for electrocatalytic reactions, non-adiabaticity may be present and the importance and contribution of non-adiabaticity may depend on the electrode potential. [34,35] As a concrete example, it has been shown that only the inclusion of vibronic non-adiabaticity in electrochemical hydrogen evolution reaction can explain experimentally observed Tafel slopes and kinetic isotope effects. [34] Second, there are several reactions where the PT and ET are decoupled for kinetic reasons. For example, in alkaline ORR pure ET has been proposed as the rate determining step [36,37,38,39]. Recent experiments of ORR on carbon-based materials conclusively demonstrate that ET is the rate-and potential-determining step. [40,41]. Also solution pH can alter the reaction mechanism and e.g. CO 2 reduction can proceed through simultaneous PCET in acidic and through decoupled PCET (ET-PT) in alkaline solutions [42,43]. In general, decoupled ET and PT are expected to play an important role on weakly bonding electrode surfaces in oxygen, CO 2 , CO, alcohol etc. reduction reactions. [44] In such reaction-catalyst combinations long-range ET/PT may take place warranting the inclusion non-adiabaticity effects.
From an applied perspective, decoupled steps may enable circumvention of thermodynamic scaling relations and lead to identification of novel electrocatalysts. [45] Third, ET, PT, and PCET include the transfer of very light particles and therefore quantum effects may be very important. Especially nuclear tunneling has a long tradition in electrochemistry [46] and experiments have conclusively demonstrated that room-temperature hydrogen tunneling takes place during ORR on Pt, and at low over-potentials tunneling is the prevalent reaction pathway. [47] Tunneling contributions are rarely considered in the field of computational electrocatalysis which is mainly due to tradition and methodological difficulties; the computational electrochemistry community has adopted tools and classical transition state theory (TST) from computational heterogenous catalysis where reactions take place at high temperatures and quantum effects are considered negligible. On the other hand, the theoretical electrochemistry community has traditionally considered ET, PT, and PCET in the non-adiabatic, tunneling framework [36,48,49,50,51,52,53,34,54,55,56,57,58,59,60,61,62]. The computational community has been slow in adopting the language and approaches developed in the theory community which has resulted scarcity of first principles study of tunneling in electrochemical environments.
Fourth, theoretical electrochemistry has a long tradition of using model Hamiltonian formulations to understand reaction kinetics.
Yet, widely different parametrizations for the same reaction can result in similar rates. For instance, differences as larges as ∼ 3-4 eV in reorganization energy and the coupling matrix elements [66,67,35,68] lead to practically identical reaction rates; it is clear that some unphysical error or parameter cancellation takes place. The difficulty of parameter estimation and error cancellation limits the physical/chemical insight obtained from model Hamiltonians. Furthermore, model Hamiltonians are static and (usually) not selfconsistent. Typically, the electrode potential serves to role of changing the Fermi-level in an otherwise static electronic structure. Even when potentialdependent electrostatic interactions and work terms are included [35,51,67,69], most parameters such as the solvent reorganization energy, chemical bonding characterized by Morse potentials, electrode structure, tunneling matrix elements etc. remain unchanged by the electrode potential. As such, it is unlikely that model Hamiltonians can quantitatively capture the complexity of electrochemical reactions. Besides issues related to self-consistency, model Hamiltonians studies of non-adiabatic reactions implicitly rely on the single orbital picture which is highly problematic for first principles Hamiltonians as discussed in the Supporting Information Section 1. Instead, modern fixed potential first-principles methods explicitly incorporate the effect of electrode potentials on the interfacial properties and bonding. Especially the GCE-DFT has proven to provide a well balanced and rigorous description electrochemical interfaces.
However, using general first principles methods for addressing ET/PCET kinetics in general have remained largely elusive thus far.
The above discussion highlights how different reactions and phenomena have been and can be addressed in the theoretical and computational communities.
Computational works utilize highquality ab initio Hamiltonians but rate constants are based on tools derived from heterogeneous catalysis and electrocatalytic reactions have been studied only using classical adiabatic TST theory. These computational studies describe the electrochemical interfaces in a self-consistent way and there is no need for empirical parametrization of the TST rate equation. Thus far, these methods have only given access to the reaction barrier but not the prefactor beyond the TST approximation.
Estimates on importance of the prefactor has relied on perturbative rate theories with model Hamiltonians at the nonadiabatic limit to describe electron/proton tunneling. Other theoretical works extend the Newns-Anderson-Schmickler model Hamiltonian to study both classical adiabatic TST and non-adiabatic tunneling reactions. While both barriers and prefactor have been computed, the models are evaluated using non-self-consistent parametrization. Therefore researchers have been be faced with a difficult choice: Should the study include all the complexity addressed in a self-consistent manner using an ab initio approach but with the restriction of classical TST approximation without general prefactors? Or should the studies include prefactors to reflect non-adiabaticity or tunneling but with a empirically-parametrized model Hamiltonian?
In this work this difficulty is resolved by developing a generally valid electrochemical rate theory which can be directly combined with fixed-potential ab initio methods. This is achieved by deriving a grand canonical rate theory building on Miller's general equilibrium (micro)canonical rate theory [70,71,72]. As Miller's theoretical framework is equally valid for adiabatic and non-adiabatic as well as quantum, semiclassical, and classical rate expressions [73] and can utilize both model or first principles Hamiltonians [33,57,58,59,60,61,74,75] the presented novel GCE extension provides a generally valid electrochemical rate theory; the developed GCE rate theory enables using all canonical rate theories in constant potential simulations. In particular, the work herein provides a unified rate theory for computing reaction barriers as well as the prefactors making the theory applicable to treat adiabatic and non-adiabatic reactions, classical and tunneling reactions, and PT, ET, and PCET on equal footing using GCE-DFT methods.
Besides developing a general and exact GCE rate theory, approximate techniques for adiabatic reactions are developed; non-adiabatic reactions are treated using the same formalism in a future publications. First, for adiabatic ET, PT and PCET reactions a generalized GCE transition state theory (TST) is derived. Second, adiabatic Marcus-like [62] empirical valence bond theories (GCE-EVB) are developed. These lead to well-defined non-linear free energy relationships ideally suited for materials' screening purposes with kinetic information as demonstrated for the acidic Volmer reaction on Au (111) in Section 4. Crucially, the developed rate theories can be seamlessly combined with modern computational methods based on (GCE-)DFT to facilitate self-consistent evaluation of rate constants without experimental parameters. The fixed potential rate theory will expand the type of systems, conditions, and phenomena in electrocatalysis amenable for first principles modelling.
The paper is organized as follows. In Section 2 a general rate theory and TST within GCE are developed. Rest of the paper focuses on ET and PCET kinetics within GCE. Section 3 shows how the adiabatic barrier and rate of ET and PCET reactions are computed using GCE-EVB and free energy perturbation theory to developed a fixed potential version of Marcus theory. Tafel slopes and other useful quantities as extracted from GCE-EVB are analyzed. A simple computational demonstration of the GCE-EVB for Au-catalyzed Volmer reaction is presented in Section 4. Next, additional computational aspects for evaluating the rate constants are discussed in 5. Finally, the advances and results are summarized.

Ensemble considerations
The GCE is open and the system exchanges matter with its surroundings.
The thermodynamics of GCE are well understood [76] and we have recently shown that both electrochemical thermodynamic quantities for both classical and quantum particles can obtained rigorously from GCE multi-component DFT [3]. GCE provides a rigorous and natural way to compute all thermodynamic expectation values at fixed electrode potentials by including the electrode potential explicitly in the ab initio Hamiltonian. This is also the case for rate constants and fixed potential rate constants are GCE expectation values of canonical rate constants as shown below.
To address the GCE rate constants, one needs to consider the dynamics of open (quantum) systems which is still an active area of research. [77,78]. The treatment of open system dynamics directly affects the GCE rate theory. First, GCE phase space volume is not globally conserved and the Liouville theorem does not hold in general and computed ensemble properties will depend on time if the system is not in equilibrium or is non-stationary. [79,78,80]. For the present work it is important that equilibrium and short-time properties are unique and time-independent in the GCE [79,81]. At other times the expectation values depend sensitively on the coupling between the system with the particle reservoir and introduces the reservoir time scales. [80] As a result, time-dependent quantities such as particle fluxes and correlation functions entering the general flux formulation of rate theory (see below) would require extensive sampling and careful computation. [78,80] To avoid the treatment of explicitly time-dependent quantities, the GCE theory developed herein only utilizes equilibrium and instantaneous quantities. Therefore, non-equilibrium processes cannot be treated using the approach taken here. Neglecting the bath time-scale and coupling also means that electron transfer kinetics from the electron "bath" to system (see Figure 1) are assumed fast, a condition satisfied by well-conducting electrodes. Neither of the the above restrictions on treating the bath coupling and time scale are expected to greatly affect the use or validity of the developed GCE rate theory in electrochemical and electrocatalytic systems.
A related consideration from on the treatment of open systems is particle conservation. If a quantum system is characterized by particle conserving operators (Ĥ Hamiltonian,Ŝ entropy, andN particle number), even time-dependent observables are obtained as ensemble weighted (p n ) expecta- Note, that changes between states with different number of particles are not included in the propagator when both the propa-gatorÛ and the operatorÔ are particle conserving. [82] Hence, even explicit propagation of the wave function does not allow sudden jumps in particle numbers. Therefore, in the extension of (micro)canonical rate theory to the GCE, only particle conserving reactions are considered. Then, all equilibrium quantities are always well-defined but jumps between states with unequal number of particles are suppressed. While this is not an issue for adiabatic reactions with smooth changes in the number of particles, the prefactors entering e.g. non-adiabatic rate constants need to be formulated so that particle conservation is respected. Therefore, all rate expressions derived herein will only utilize particle conserving operators.

General grand canonical rate theory
After establishing the particle conserving and equilibrium nature of the rate constants, the GCE rate constants can be formulated. To allow various types of reactions to be described, the exact equilibrium canonical rate expression due to Miller [70,71,72,83] is adopted: where Q I is the canonical partition function of the initial state, and β = (k B T ) −1 .
The first expression is written in terms of transition probability at a given energy P (E). The second expression utilizes a canonical flux-side correlation function constrains the trajectories to start from the dividing surface,q is the initial flux along the reaction coordinate, and h[f (q t )] is the side function which includes the dynamic information whether a trajectory is reactive or not.
Based on the discussion in Section 2.1 on the dynamics of open systems, only the t → 0 + and t → ∞ should be considered for the fluxside correlation function in the equilibrium rate expressions. Depending on the choice of P (E) or H and h[f ] non-adiabatic and adiabatic (nuclear) quantum effects are included in the rate. [84,85,86,87]. It is noteworthy that P (E) and C f s are computed using only particle conserving operators [71] and the conditions discussed above are satisfied when (1) is used as the starting point for formulating GCE rate constants.
To compute reaction rates at fixed potentials a straight-forward, yet novel, extension of the canonical rate theory to the GCE is made: where Ξ I = exp[βµN ]Q I is the initial state grand partition function and k(T, V, N ) was introduced in (1). Above, N is the number of species (nuclear or electronic) in the system and C µ f s is the GCE flux-side function.
The previous equation shows that all canonical rate equations can be applied to electrochemistry within GCE approach and that fixed potential electrochemical rate constants are GCE averaged canonical rates constants.
The above equations are completely general and various flavors of rate theories can be extracted by invoking different Hamiltonians and transition probabilities, but they are somewhat cumbersome for computational purposes. Indeed, it would be convenient if the GCE rates could be directly evaluated without explicitly summing over different particle numbers.
One way to achieve this is to make the transition state theory (TST) assumption [72,71,70] but generalized to GCE herein. In canonical TST, the instantaneous lim t→0+ C f s (t) is considered corresponding to the assumption that there are no recrossings of the dividing surface.
Both quantum/classical and adiabatic/non-adiabatic TSTs are written as [88,89,90,91] and the exact rate is recovered by introducing a correction where κ(t) is the time-dependent transmission coefficient which at long times is κ = k(T, V, N )/k T ST (T, V, N ). [92] Inserting this equation in (2) results in the most general grand canonical rate constant. Significantly simplified rate constants are obtained when focusing on classical nuclei and using TST. As derived in the SI Section 2, for classical nuclei the TST result is [71,72]: where P cl (E) denotes transition probability for classical nuclei but the electrons are of course quantum mechanical [75,93] with details given in [72] and the SI Section 2. The previous equation shows that the structure of GCE-TST and canonical TST are similar which is true for open system in general if memory effects are neglected [94]. To obtain the GCE rate constant without invoking the TST approximation, one can use the transmission coefficient κ to write where it is assumed that an effective transition probability κ µ can be used.
To complete the derivation for the classical GCE rate constant, the rate is expressed in terms of grand energies with the definition Ω i = − ln(Ξ i )/β and ∆Ω ‡ = Ω ‡ − Ω I for the GCE barrier. Above the only new assumption besides grand canonical equilibrium distribution and TST, is that the flux out of the transition state κ µ can be treated as an expectation value and separated from the barrier. For large enough systems and small variations in the particle number this is a justified assumption.
The above development establishes the general fixed chemical potential rate theory. For classical, adiabatic reactions the rate constants in GCE are essentially the same as in the canonical ensemble. Within TST approximation the rate constant is determined by the grand free energy barrier and effective prefactor. The transmission coefficient needs to be approximated but this depends on the case at hand; examples for the adiabatic and non-adiabatic harmonic GCE-TSTs expression valid for fully open system are derived in Supporting Information section 3.
A more thorough treatment on the theory and computation of non-adiabatic and tunneling corrections will be presented in forthcoming work.

Semi-grand canonical ensemble
The above development is valid when both nuclear and electronic subsystems are open. A significant simplification results if one assumes that the reaction rate does not explicitly depend on the number of some nuclei in the system. In a typical first principles calculation this simplification is often exploited when the system can be divided to two subsystem: 1) classical electrolyte species consisting of nuclei and electrons and 2) electrode + reactants treated either classically or quantum mechanically. Typically the number of nuclei constituting the electrode and reactant are fixed while the electrolyte and electron chemical potential are fixed.
Fixing only the electron and electrolyte chemical potentials defines a semi-grand canonical ensemble used for deriving the thermodynamics of electrocatalytic systems within GCE-DFT [3]. In this treatment is often utilized in e.g. Poisson-Boltzmann type models where the electrolyte is at a fixed chemical potential but the energetics do not explicitly depend on the number of electrolyte species. Then, summation over the number of electrode/reactant nuclei or the electrolyte species is not needed.
Herein the semi-GCE is applied to derive rate constants as a function of the electrode potential. From now on, I assume that the reaction rates depend explicitly only on the number and/or chemical potential of electrons in the system. Then, the state of the system is determined by T , V , number of nuclei of the electrode+reactant N N , chemical potential of the electrolyte, chemical potential of the electrons µ n , and number of electrons in the system N unless explicitly specified otherwise. Electroneutrality is maintained by the electrolyte. A widely utilized harmonic TST rate for constant number of nuclei and constant electrochemical potentials are derived in section 3 of the Supporting Information.

Adiabatic barriers and rates from GCE-EVB
To compute the GCE-TST rate at a given electrode potential, the grand energy barrier in (6) needs to be obtained. For electronically adiabatic reactions methods like the constant-potential [20] nudged elastic band [95] can be used. An alternative method for computing the grand energy barrier is to formulate a Marcus-like [62] approach or empirical valence bond (EVB) theory [96] within GCE. Such models are commonly utilized in electron [62] and proton transfer theories. [65,96,97,98,99] Here, the treatment is based novel development of GCE diabatic states and the extension of the canonical thermodynamic perturbation theory to the GCE to facilitate derivation of a GCE-EVB rate theory (see SI sections 3 and 4). The GCE-EVB theory provides a theoretically welljustified and computationally affordable way to obtain fixed potential barriers at various electrode potentials; the adiabatic barrier needs to be explicitly computed only at a single electrode potential while barriers at other potentials can be obtained using well-defined extrapolation of (17). The utility of the GCE-EVB theory is demonstrated in Section 4.
In canonical EVB and Marcus theories use diabatic states, effective wave functions and free energies [62].
This can be extended to GCE by using two fixed potential, diabatic ground state surfaces which represent a GCE-statistical mixture of states with probabilities given by the density operator in GCE [3]. Importantly, the diabatic states obtained using the GCE density operator naturally include many-body effects of the coupled electrodereactant-solvent system and the complexity of the electrochemical interface is explicitly included in the model. Also, there is no need to decompose the rate constants to orbital dependent quantities(see Section 1 in the Supporting Information for additional discussion). Then, two grand canonical diabatic allelectron wave functions are used to form an effective diabatic GCE Hamiltonian. This is analogous to molecular Marcus theory utilizing a canonical diabatic Hamiltonian containing an initial (oxidized) I and final(reduced) molecule F .
Following the treatment in the Supporting Information Section 3, a diabatic 2 × 2 grand canonical Hamiltonian in (7) can be formed from two diabatic GCE states. The resulting form is analogous to the canonical EVB methods [96], electron [62], proton [98,99] and proton-coupled electron [65] theories. The present form is, however, crucially different from its predecessors; based on the approach developed in this work, all quantities are defined and computed at fixed electrode potentials. In the basis two GCE diabatic states the GCE Hamiltonian is as derived in Supporting Information Section 3. Here the diagonal elements are the grand energies of the initial (II) and final (FF) systems. The off-diagonal elements account for the interaction and mixing between the initial and final states. They can be computed as GCE expectation values of contributions from different N-electron states ψ N i as Ω I . This is rather straightforward for ET reactions using e.g. constrained DFT discussed in Sections 4 and 5. For PT and PCET reactions computing these matrix elements would require computing the vibronic matrix elements using e.g. the semiclassical approach of Georgievskii and Stuchebrukhov [100] This direct computation is particularly useful for non-adiabatic rate constants which are investigated in future work.
For adiabatic reactions, the direct calculation can be replaced by the diagonalization of the 2 × 2 diabatic Hamiltonian in Eq. (7). This diagonalizating produces the adiabatic ground and first excited states as As shown below, the diabatic states cross (Ω II = Ω F F ) at the transition state. This makes it possible to compute the coupling matrix element as the difference between the diabatic states and the adiabatic states. For the ground state one has Ω IF = Ω II − Ω − ad . The adiabatic ground transition state grand free energy can be computed using e.g. NEB calculations and the coupling matrix element is simply the difference between the diabatic transition state grand energy and the adiabatic one as shown below in Eq. (18).
Finally, the (diabatic) grand canonical states correspond to a single electron density which is guaranteed by the Hohenberg-Kohn-Mermin [4,3] to be unique for a given electrode potential.
By definition, GCE diabatic states are unique ground states. Such diabatic statas also include the interaction and exchange between all the electrons in the system and for adiabatic, ground state reactions there is no need to include addition excited states despite the continuum of (single-electron) states of the electrode. In principle it is possible to add other, possibly excited states as basis states but here the focus in on treating adiabatic reaction and excited states beyond the first excited state are neglected. If a general quantum mechanical Hamiltonian is used, bond breaking is naturally included in the GCE-EVB model. The only disambiguity is the definition of diabatic states. In practice the GCE diabatic energies, (Ω II and Ω F F ), can be computed directly by applying using e.g. cDFT [101,102,103] with fixed potential DFT as discussed in Section 5 and shown in Section 4.

Computation of diabatic GCE energy surfaces and barriers
An approach often used in molecular simulations for constructing the diabatic free energy curves is to sample the diabatic potentials along a suitable reaction coordinate. For canonical ET, PT, and PCET reactions the reaction coordinate is the energy gap between the two diabatic states as shown by Zusman [104] and Warshel [105]: ∆E gap (R) = E F (R) − E I (R). [76,106] From the sampled energy gap, free energy curves are obtained as A(R) = −k B T ln(p(E gap (R))) + c. If the distribution is Gaussian p(E gap (R)) = c exp −(∆E gap − ∆E gap ) 2 /2σ 2 , the resulting free energy curves a parabolic. The diabatic barrier in EVB or Marcus theory is then obtained from the intersection of the initial and final diabatic curves [106,107,108,109].
Within GCE, the energy gap is simply N ). As shown in the SI section 4, the gap distributions can be formulated and computed by generalizing Zwanzig's [110] canonical free energy perturbation theory to the GCE. This provides a rigorous way to derive the reaction barrier in terms of diabatic states and energies as presented in the Supporting Information Section 4. The reaction energy rate can be computed from the initial-final state energy gap distribution functions using a well-known formula [105,111,112,113,114,115,116] where g i (∆E) is the free energy curve in state i as a function of the energy gap, p I (∆E ‡ ) is the gap distribution at the transitions state, and κ denotes an effective prefactor. The reaction rate is determined by the energy gap distribution function p I (∆E) = δ(∆E(R) − ∆E) I from equation (S30) of the Supporting information.
While the approach is general and valid for complex reactions, assuming that E gap (R; µ) is Gaussian leads to a closed form equation. In this case the GCE-diabatic states are parabolic and the Marcus barrier in GCE is given by (13). As shown in the Sections 4 of the SI, the (Gaussian) gap distribution may be derived using a second order cumulant expansion resulting in where ∆E I is the energy gap expectation value in the initial state obtained from equation (S27) in the Supporting Information and σ I = (∆E I − ∆E I ) 2 I is the gap variance. The Marcus relation is then obtain after standard manipulations [106,112] yielding where σ 2 I = σ 2 F = 2k B T Λ = k B T ( ∆E I − ∆E F ), Λ is the fixed potential reorganization energy and ∆Ω = ( ∆E I + ∆E F )/2 is the reaction grand energy as depicted in Figure 2. These gap identities are valid for symmetric reactions and have been previously established well for the canonical ensemble [112] and generalized here to the GCE. In practice, the reorganization energy is computed as an average of the reorganization energies which are differences of the diabatic free energy at the final geometries shown in Figure 2. Finally, the GCE-EVB rate equation using the above assumptions results in an expression analogous to Marcus equation The energy barrier of (13) is the diabatic energy barrier. The adiabatic barrier is estimated from (7) using the methods discussed in Section 3.2 below. One caveat to keep in mind is the more involved computation of κ within the GCE as discussed in Section 2 and forthcoming work for non-adiabatic reactions. The above result may safely be used when κ ≈ 1 for all particle numbers meaning that the reaction is always fully adiabatic and classical.

Implications of the canonical GCE-EVB rate theory
If the diabatic grand energy surfaces are symmetric and quadratic they have the same curvature and reorganization energy. In this case, the diabatic grand energy barrier is estimated from (13). The assumption on equal curvature can be relaxed [117] (see also SI section 5). One easy approach to realize this is to utilize an asymmetry parameter α as as [118] in terms of the reorganization energies for both the initial and final states Λ I and Λ F , respectively. The transition state is located at the crossing point With these definitions the asymmetric diabatic Marcus barrier and rate become When α as → 0, the regular Marcus barrier and crossing point are obtained. In Figure 3 the effect of asymmetry and reaction energy to the reaction barrier and location of the transition state are compared. It can be seen that both the barrier heights and its location are affected by the asymmetry and reaction energy. While the Marcus-like equation results in a diabatic barrier, the adiabatic reaction barrier can be extracted from the diabatic barrier by diagonalizing (7). The adiabatic barrier can also be obtained from (13) using the Hwang-Åqvist-Warshel adiabaticity correction [119,120] ∆Ω ‡ ad,EV B = where Ω IF is the off-diagonal matrix of the GCE-EVB Hamiltonian in (7). If the Condon approximation is used, the above equation is greatly simplified as Ω IF ≈ Ω IF (x ‡ ) ≈ Ω IF (x I ) becomes a geometryindependent constant.
Next changes in the adiabatic GCE-EVB barrier as function of the parameters is analyzed. From the schematics shown in Figures 2 and 3, one can observe that changes of the minima along the reaction coordinate correspond to horizontal displacements of the diabatic states and changes in Λ. Vertical changes correspond to changes in the reaction grand energy ∆Ω. In general, the reorganization energy of innersphere reactions taking place on or near the electrode surface may depend on the electrode potential and investigations along this direction are on their way.
Under equilibrium conditions, ∆Ω = 0 and the corresponding reorganization energy Λ 0 , the adiabatic barrier is (18) which leads to Λ 0 = 4(∆Ω 0, ‡ ad,EV B + Ω IF ) ≈ 4∆Ω 0, ‡ dia assuming that Ω IF << Λ 0 (this is the case for e.g. the Au-catalyzed Volmer reaction in Section 4). At the equilibrium point, the overpotential is, η = ∆Ω = 0. Assuming for a moment that Λ ≈ Λ 0 and replacing the solution for Λ 0 in (17) gives the diabatic barrier as Inserting (19) in (17) This result is well-known in the canonical EVB and Marcus theories. However, in GCE this relation is valid only when constant reorganization energy is assumed. In general, the driving force can be manipulated easily by changing the electrode potential which is in turn directly related to the absolute electron electrochemical potential as E M (abs) ∼ −μ n . [3,121,122] An experimentally meaningful approach is to study −∂r(T, V,μ n )/∂μ n as done in a Tafel analysis, for example. Tafel analysis can also be understood in a more general context of Brønsted-Evans-Polanyi (BEP) and other free energy relations measuring the change of reaction rate when the reaction energy is changed [123,124,125], as both Tafel and BEP analyses measure the reaction rate as a function the reaction driving force -Tafel analysis focuses on the overpotential and BEP on the free energy. These two quantities are linked by |∆η| = |∆μ n | = |∂∆Ω/∂n|. Defining the rate constant as a function of the electrode potential E as k(E) = k(E = 0)A(E) exp(−βαE) in terms of the prefactor A and the Tafel-BEP coefficient α and the Tafel-BEP coefficient is [2,123,124] where constant α and prefactor are assumed. Within GCE-EVB α is obtained in terms of the reorganization and reaction energies where the first term measures how the rate changes as a function of the reaction energy, γ denotes a BEP coefficient and ∆Ω denotes the grand energy change as a function of the over-potential. The second contribution is novel and unique to the GCE formulation. It measures the sensitivity of the rate to changes in the reorganization energies as a functional of the potential. This unconventional contribution can be observed in e.g. the Volmer reaction treated in [67] with a model Hamiltonian and will be discussed in more detail in a future publication.
To facilitate understand the BEP term, one For macroscopic systems, chemical reactions have N F = N I while simple electrochemical steps have N F = N I ± 1. For chemical reactions ∆Ω = ∆A and the variation ∆Ω is small. Within the computational hydrogen electrode (CHE) concept [126] the reaction energy ∆Ω ≈ ∆A 0 ∓ η for PCET steps with ∆A 0 computed without any bias potential. Then, α CHE = γ for PCET steps and zero otherwise. In general such a simple relationship does not hold in general and models such as GCE-DFT can be used for computing ∆Ω explicitly. Thus far, ∆Ω has been reported in only few studies [20,21,127]. In these works and in Section 4, ∆Ω is found to exhibit a roughly linear dependence on the applied potential.
Next the BEP γ of Eq (22) is analyzed. Using the diabatic barriers in (19) (obtained assuming constant reorganization energy and constant prefactor), one obtains which results in α = −∆Ω (1/2 + ∆Ω/2Λ 0 ). It is seen that γ is not a simple constant but depends linearly on the reaction driving force. If the reorganization energy is small the dependence on the reaction grand energy becomes more pronounced as demonstrated for the Au-Volmer reaction in Section 4. In general, non-linearity of the grand energy barrier has two contributions: non-linearity of the diabatic barrier and the potential-dependent reorganization energy. For macroscopic systems non-linearity is established by including the quadratic part of the diabatic barrier in the model. Lately [29,31,20,21] this has been observed computationally and it is pleasing that the GCE-EVB picture seems qualitatively correct.
To summarize, the generalized BEP-Tafel relationships have been derived from a microscopic perspective using grand canonical rate theory. Both variation in the reaction energy barrier and the transition state location as a function of the potential can be predicted using just a few parameters. This is demonstrated for the acidic Volmer reaction in Section 4. The general form of the BEP-Tafel relation is given in (22). For small over-potentials, the rate is expected to depend linearly on the applied potential. For larger overpotentials non-linear dependence is predicted when the reorganization energy is small. The GCE perspective also predicts a novel potential-dependent reorganization energy which is supported by model Hamiltonian calculations [67] and will be addressed carefully in future publications.

Application of the GCE-EVB theory to the Au-catalyzed Volmer reaction
Here the first demonstration of the GCE rate theory is provided. I consider the acidic Volmer reaction i.e. proton discharge which is arguably the simplest and yet relevant electrocatalytic reaction for hydrogen production and other electrocatalytic reactions. Similarly, gold can be considered as the simplest electrode material. Yet, the Volmer reaction, even on gold, is not fully understood [128] and the reaction is considered to exhibit nuclear quantum effects and even vibronic non-adiabaticity. [34,35] As a first application of the theory and methodology derived and developed in this work, I consider mostly an adiabatic and classical model for the acidic Volmer reaction -quantum effects and a non-adiabaticity are studied separately in forthcoming publications. The results are discussed in the GCE-EVB framework of Section 3.
The Volmer reaction is modelled as a single hydronium ion on a 3x3x5 Au(111) surface as shown in Figure 5. The needed free energies were computed using GCE-DFT as implemented in GPAW [129] within the surface-jellium approach [20] with a continuum solvent model for water [130]. This approach gives all the thermodynamic quantities at a fixed electrode potential. The potential-dependent minimum energy pathways are computed using a nine image nudged elastic band [95] (NEB) discretization. Geometries and NEB pathways were considered converged when the maximum force was below 0.05 eV/Å. Constant-potential diabatic states and reorganization energies were computed using constrained DFT (cDFT) as implemented in GPAW [131]. As in the canonical case, the constraining potential in GCE-cDFT is introduced as an external potential to the GCE-DFT giving where Ω[n(r); T, V, µ n ] DF T is the GCE-DFT energy functional [3] and n(r) is the electron density. w σ i (r) is the weight function which defines how the charge is to be partitioned, i.e. the regions where charge is to be localized, N c is the desired number of excess electrons within the constrained region, and V c is the Lagrange multiplier enforcing the charge/spin localization. The introduction of constraining terms in Eq. (24) leads to a new effective potential defined as Thus, the cDFT potential is just the sum of the usual KS potential and the constraining potential which is also used in the self-consistent calculation. The constraint is further enforced by introducing the convergence criteria The optimize grand canonical free energy under the specified constraint is [131] To fulfill the cDFT constraints and to perform fix potential calculations, the approach in Fig. 4 has been utilized. An example of GPAW scripts used for performing GCE-cDFT calculations is given in the SI section 6.
The charge states to define diabatic states are chosen as N c = +1 state for the hydronium (H 3 O + ) and neutral N c = 0 for the final water and adsorbed hydrogen (H 2 O +H * ). The reorganization energies are computed from Eq. (12) at the equilibrium potential which resulted in Λ F = 2.1 eV for the initial state geometry, and Λ I = 3.2 eV for the final state geometry. The average Λ = 2.65 eV is used in calculating the Marcus barrier of Eq. (13). These cDFT computed reorganization values are in good agreement with the . values used for model Hamiltonian parametrizations of Huang [68] and Santos [66] but much larger than the one used by Lam [67]. The reaction asymmetry from (14) is 0.2 meaning that the transition state geometry along the reorganization coordinate is closer to the initial state. Analysis of the TST geometries at different potentials shows that that reorganization energy is best presented either by the Au-O distance or the dihedral angle between the H 2 O and surface.
In Figure 5 the adiabatic reaction barrier as obtained from NEB calculation is plotted as a function of the reaction energy corresponding to different electrode potentials. First, one observes that the barriers are very small for all considered electrode potentials. This is in line with explicit water DFT results [20] at all electrode potentials and reaction energies. The figure also shows the adiabatic TST location as a function of the reaction energy from both NEB-DFT and extrapolation using (16a). The extrapolation reproduces the TST geometries very well, and captures the trends in the TST location. For comparison, explicit solvent calculations exhibit a similar trend in the TST position as a function of the potential [20] as the one found here using an implicit solvent. This example demonstrates the Au-O distance is good reorganization coordinate and that the TST location is effectively captured by Eq. (16a). Unlike DFT-NEB calculations, the GCE-EVB requires just one NEB and two reorganization energy calculations to capture the TST geometry for a range of reaction energies and electrode potentials. The stars correspond to TST geometries predicted using equation (14a). Below:The fixed potential Au(111) Volmer barrier as a function of the reaction energy. NEB[*] refers to NEB calculations of the present paper with an implicit solvent while NEB [20] are from [20] with explicit, ice-like solvent. Also the barrier and reorganization energy used in the model Hamiltonian work of [67] are shown and extrapolated using (20).
The results in Figure 5 also show the GCE-EVB barrier as obtained from Eq. (20) which is evaluated using the cDFT computed average reaorganization energy. The adiabatic equilibrium barrier (∆Ω 0, ‡ ad in Eq. (20)) is computed from a NEB calculation close to the equilibrium potential. As the results in Figure 5 show, the extrapolation with Eq. (20) provides a very accurate way for computing the adiabatic potentialdependent energy barrier. It is also observed that the estimate for the reorganization energy used in the model Hamiltonian work of [67] is very small (∼ 0.3 eV) and cannot be used for predicting barriers using GCE-EVB.
Next, the validity of the adiabatic assumption is tested using the coupling matrix elements of Eq. (7). The coupling matrix element is obtained using Eq. (8) at the equilibrium potential transition state and results in Ω IF = 0.37eV . The coupling constant maybe used for estimating transmission coefficient κ µ of Eq. (6) the adiabaticity from the the Landau-Zener factor (P LZ ) for PCET reaction [132,35] adopted to the GCE.
and κ µ = 1 is a signature of an electronically adiabatic reaction whereas non-adiabatic reactions have κ µ << 1. Eq. (28a) is valid for reactions in the normal region. [132,35] ν n is the vibrational frequency along the reaction coordinate. [133] As discussed above, the reaction coordinate n the Au-catalyzed Volmer is the water reorientation and ν n = 1/τ L ≈ 0.5 × 10 9 s −1 computed from the water reorientation time τ L ≈ 2 ps [134]. The expression for γ in Eq. (28c) is valid for quadratic free energy surfaces as derived in Eq. S.15 of the Supporting Information. Evaluating the transmission coefficient at room temperature using the computed average reorganization energy, coupling matrix element, and water reorganization time gives κ µ = 1.0 -this shows that the reaction is adiabatic and justifies the treatment in this section.
Besides enabling a reliable prediction of reaction barriers and TST geometries, the results demonstrate the first GCE-cDFT calculations. Besides showing that GCE-cDFT is technically possible, the results show that ab initio computed diabatic states offer new insight to electrocatalytic reactions. In particular, the results provide a proof-of-principle that GCE-EVB can be used to accurately estimate barriers using just a single NEB calculations and a few cDFT calculations with an expense similar to a standard DFT calculation.

Discussion
The distinct advantage of the formalism and theory developed in this paper is that all rate equations can be readily evaluated with GCE-DFT or other first principles approaches. The presented formalism enables the treatment electrochemical and electrocatalytic thermodynamics and kinetics in terms of the prefactor and barriers in the same self-consistent framework -the GCE-DFT. Therefore, the same DFT-based tools can be used to address inner-sphere and outer-sphere kinetics and thermodynamics instead of modifying or changing the theoretical and computational framework for different reaction steps [36].
By construction the rate constants include the interplay between the electronic structure, solvent, electrode potential etc.
All quantities can be computed using self-consistent DFT energies and "wave functions" to include exchange and correlation effects between all the electrons in the system. The Fermi-Dirac distribution is fulfilled at the DFT level and, therefore, there is no need to integrate over the filled/empty orbitals weighted by the Fermi-Dirac distribution in the rate expression as done is traditional single-orbital descriptions (see SI Section 1 and below). Also, the Kohn-Sham-Mermin theorem [3] guarantees that GCE-DFT and GCE-EVB states are unique to a given electrode potential and that the GCE-EVB diabatic inlcude that interactions between all electrons in the metal and the reactants.
The electrode potential is self-consistently treated and all free energies and prefactors depend explicitly on the potential. This is in contrast with traditional model Hamiltonian treatments where the electrode potential rigidly shifts the Fermi-level without modifying any interactions or prefactors [64,69] or modifies only the electrostatic interactions [34,67]. Also, separate computation of work terms [67,69,68,135] is not needed because all relevant interaction can be directly included in the general Hamiltonian. Evaluation of chemisorption functions entering adiabatic Newns-Anderson-based models [63,66,67,64,68,135] is also avoided. Therefore, the current models are free of approximate treatment of semi-elliptic DOSs [67,66,135] or fitting the chemisorption functions to a computed DOS [66,135].
As the developed rate theory utilizes general ab initio Hamiltonians, bond formation/breaking are naturally included. This is again in contrast with model Hamiltonians which require approximate potential-independent terms to describe changes in atomic bonding [136,137,135].
Instead, as demonstrated herein, ET, PT, or PCET and bond rupture/formation are naturally captured with GCE-DFT. Bond formation s is also captured by diabatic models using cDFT as demonstrated herein for the Volmer reaction and previously for ET [138], PCET [139] and general chemical reactions [140,141].
As all necessary terms can be computed from GCE-(c)DFT, adoption and evaluation of the rate expressions is straight-forward (but potentially laborious). While applicability and usefulness of combined DFT and GCE-EVB was demonstrated for the Volmer reaction, it is worth discussing the additional computational requirements in some detail. First, the simulation of charged systems is needed to sample the electrode potential. Electroneutrality can be enforced using some variant of the Poisson-Boltzmann model, for details see [3]. Fixed potential calculations can be ac-complished within a single SCF cycle [10], or iteratively [142,20]. Second, the solvent effects should be included in the model. In traditional TST-based models for adiabatic reactions the main solvent contribution is thermodynamic and stems from (de)stabilization of different structures. GCE-EVB models need to involve the solvent as the reaction barrier is directly related to the solvent/environment reorganization energy and neglecting the solvent contributions will most likely lead to incorrect results.
Given a software capable of handling charged systems and performing constant potential calculations, adiabatic TST rate constants can be readily evaluated. As shown for the Volmer reaction, reaction barriers and adiabatic prefactors at constant potential are obtainable using e.g.
the NEB [95] method. Evaluation of GCE-EVB rate constants requires additional software capabilities for constructing charge/spin localized diabatic states and to evaluate the electronic coupling between these states either through direct calculation or diagonalization of the diabatic Hamiltonian of Eq. (7). Also the reorganization energy, which is an excited state quantity, needs to be computed. One widely implemented and available tool for evaluating the additional parameters is the cDFT methodology [101,102,103] which is implemented in several DFT codes [131,143,144,145,146,147,148,149,150,151,152]. The extension of cDFT to GCE was, for the first time, demonstrated and applied in this work. GCE-EVB should be accompanied with a constant potential simulations to compute fixed potential reaction and reorganization energies.
The general framework can also be extended to treating reactions beyond adiabatic reactions to include e.g. non-adiabatic effects, nuclear tunneling, and solvent-dynamics controlled reactions -treatment of these effects is under current study and will be published separately. Non-adiabatic effects are expected for e.g. outer-sphere ET reactions and several PCET reactions [33]. Several PT and PCET reactions are also likely include adiabatic or non-adiabatic nuclear tunnelling effects. Also, solvent dynamics should be included as these are likely to become increasingly important or even dominant when the reaction is adiabatic and the reaction barrier becomes very small or vanishes. [153] Under such conditions the reorganization will be the slowest process and the reaction prefactor should reflect this. Last, welldefined interpolation [154,155,156] between adiabatic -solvent dynamic -non-adiabatic should be developed or adapted to the fixed-potential rate theory.
Finally, a spectacular feature of canonical Marcus and EVB theory is the observation of an inverted region i.e. the rate constant starts to decline as the reaction becomes more exothermic. However, the inverted region has not been observed for electrochemical reactions even at large over-potentials. The grand canonical Marcus rate of (13) seems to predict an inverted region for highly exothermic conditions and warrants additional discussion. First, the Marcus-like expression is based on linear response or second order cumulant treatment [157] which leads to quadratic free energies along the reaction coordinate -the prediction of the inverted region is a direct consequence of the linear response assumption. To go beyond the quadratic free energy surfaces higher order cumulants can be added (see SI section 5) to modify the existence of the inverted region. [158] For instance, the inverted region is not predicted for Morse potentials. [159] Secondly, the inverted region is very sensitive to tunneling effects, excited states, and solvent effects. [160] For example, nuclear quantum effects are needed to achieve accurate ET rates in the inverted region ET [161] while excited proton vibronic states dominate PCET reaction rate [159]. Even if the inherent approximations on quadratic, classical, and adiabatic nature are acceptable, in GCE-EVB the reorgnization energy is potential-dependent and the inverted region is suppressed if ∆Ω ≈ 0. Also the reorganization energy should saturate at large overpotentials.
The prediction of non-existing inverted region is inherent to several non-adiabatic single orbital approaches. In these treatments the inverted region is avoided by taking the manifold of single-electron states into account and integrating over the Fermi-Dirac weighted transition rates, as discussed in Supporting Information Section 1. This has resulted in e.g. Marcus-Hush-Chidsey [162], Dogonadze-Levich-Kuztnetsov [49,48], Soudackov-Hammes-Schiffer [53] models of ET and PCET. On the other hand, adiabatic rate computed using the Newns-Anderson-Schmickler Hamiltonian [63,64] does not the include the orbitalto-orbital contributions separately as the redox orbital is coupled with all levels on the metal by definition. In the adiabatic Newns-Anderson-Schmickler model, the free energy surface becomes as single well and the reactant and product become indistinguishable. As the vanishing barrier cannot be treated [63,154], the barrier is simply extrapolated to zero in the inverted region and the rate is controlled by the prefactor.
Based on the above discussion, the inverted region is more complicated than the normal region and requires careful consideration of excited states and tunneling effect, for example. The vanishing inverted region is for non-adiabatic single-orbital models and not for adiabatic models where the exchange between all metallic states and the redox molecule are treated. In the GCE-EVB picture the inverted region is inherent to the quadratic grand energy surfaces and simply extrapolation of the barrier to zero removes the inverted region. However, this is not fully satisfactory and more studies are needed to study non-quadratic free energy surfaces, tunneling, and excited states in the inverted region.

Conclusions
In this work a new theoretical formulation for computing electrochemical and electrocatalytic rate constants at a fixed potential is developed. The rate expressions are obtained by extending the universally valid and exact canonical rate theory [70,71,72] to the grand canonical, fixed potential ensemble. General conditions and limitations for the fixed potential rate theory are developed. It is shown that all rate theories developed within the canonical ensemble can be transferred to the GCE; electrochemical rate constants are "just" GCE-weighted canonical rate constants. This is conceptually important because the fixed-potential rate theory enables treating all potential-driven reactions within a single formalism instead of relying on separate theories for different cases.
This provides a unified framework for computing and understanding both the barrier and prefactor from a single formalism towards modelling beyond adiabatic inner-sphere reactions to include e.g. non-adiabaticity and tunneling. Herein, specific rate expressions are derived for typical electrocatalytic adiabatic ET, PT, and PCET reactions. Fixed potential rate expression have been derived for i) general electrocatalytic reactions with ((5)) and without ((2)) the TST approximation, ii) electronically adiabatic ET, PT and PCET reactions using a grand canonical Marcus-like GCE-EVB theory in (13).
The GCE-EVB formalism enables computing the grand energy barrier in terms of fixed potential reorganization energy and the reaction grand energy in analogy with the canonical EVB or Marcus theory. GCE-EVB can explain and predict the electrocatalytic "Marcus-like" behavior in energy barriers and TST geometries as a function of the thermodynamic driving force. GCE-EVB enables also the computation of non-linear energy relationships and Tafel slopes and general BEP-Tafel relations in goo accuracy using just few DFT parameters; this has shown for the Volmer reaction using fixed potential (constrained) DFT developed here. Also, quantitatively accurate barriers and TST geometries can be predicted using a few selfconsistent GCE-DFT calculations. These features are expected to make the GCE-EVB approach particularly suitable for electrocatalyst screening studies.
The developed theory can be directly combined with modern, solid-state ab initio methods to capture the complexity of the electrochemical interface at constant potential in a self-consistent manner. In this sense, the model is fully ab initio and all parameters can be directly computed. A set of widely implemented DFT-based tools suffices to compute all the needed parameters in a self-consistent manner. This should enable the computational community adopt the theoretical framework and to progress from a thermodynamics-based description of electrocatalysis to addressing also electrocatalytic kinetics under experimentally realistic conditions.
The advances herein enable further development of theory and computational methods to address e.g. tunnelling pathways and non-adiabaticity in electrochemical systems from first principles. Understanding and controlling (non-adiabatic) tunneling can open up new reaction pathways to avoid constraining scaling relations [163,164,165] encountered for adiabatic PCET reactions. Besides applications, the advanced rate theories will improve fundamental understanding of electrochemical kinetics in e.g. de-coupled and nonadiabatic ET, PT, and PCET. This contributions are especially important for weakly-binding catalysts, but neglected in computational studies thus far.

Acknowledgements
I acknowledge support by the Alfred Kordelin Foundation and the Academy of Finland (Project No. 307853). Computational resources were provided by CSC -IT CENTER FOR SCIENCE LTD. I also wish to thank the anonymous reviewer for constructive criticism and for providing comments to improve the manuscript.

Declaration of interest
Declarations of interest: none 9. References