Grand canonical ensemble approach to electrochemical thermodynamics, kinetics, and model Hamiltonians

The unique feature of electrochemistry is the ability to control reaction thermodynamics and kinetics by the application of electrode potential. Recently, theoretical methods and computational approacheswithin the grand canonicalensemble(GCE) have enabled to explicitly include and control the electrode potential in first principles calculations. In this review, recent advances and future promises of GCE density functional theory and rate theory are discussed. Particular focus is devoted to considering how the GCE methods either by themselves or combined with model Hamiltonians can be used to address intricate phenomena such as solvent/electrolyte effects and nuclear quantum effects to provide a detailed understanding of electrochemical reactions and interfaces.


Introduction
Electrochemical reactions take place at extremely complex electrified solideliquid interfaces.Experimentally, the properties of such electrochemical interfaces and (elementary) reaction thermodynamics and kinetics are controlled by the choice of the electrode material and reaction conditions including the electrode potential (U), temperature (T), and concentrations (c) or equivalently (electro)chemical potentials e m.Computational and theoretical models should explicitly use the same variables but realizing such a level of control continues to pose difficulties [1,2].In recent years, the grand canonical ensemble (GCE) approach [3e6] has been established as a general framework for treating electrochemical thermodynamics and kinetics as an explicit function of U ; T ; e m d here some of these developments and their foundations are reviewed and some outstanding open questions highlighted.
Particular focus is placed on methods to describe electron transfer (ET) and proton-coupled electron transfer (PCET) reactions at electrochemical interfaces, which are crucial in both fundamental and application-oriented electrochemistry [2].It needs to be recognized that ET and PCET by themselves are very complex reactions and exhibit various physical/chemical phenomena such as nuclear tunneling, electronic/vibronic nonadiabaticity, sensitivity to solvent dynamics, and so on.[7,8] Electrochemical ET/PCET is even more complex necessitating the use of advanced theoretical and simulation approaches in collaboration with (spectro) electrochemical experiments to understand how the pH, solvent, potential, and electrode material control electrochemical reactions [2, 9,10].
Currently, most theoretical works in electrochemistry are computational studies using density functional theory (DFT) to study PCET thermodynamics within the computational hydrogen electrode (CHE) method [11] and Sabatier's principle to predict catalyst performance [1,12].Often, the applied computational and theoretical methodologies are highly simplified leaving many fundamental aspects of ET/PCET and electrochemical interfaces unanswerable [2,13].For instance, the importance of tunneling or solvent dynamics has been experimentally realized [10,14], but such effects are often overlooked or cannot even be properly described using basic DFT methods.Such methodological limitations cannot be circumvented by building large materials databases [15,16] only but requires developing more insightful models of electrocatalysis.Although DFT is very useful, I share the concern of Schmickler et al. that "Theoretical treatments only based on Density Functional Theory calculations are not sufficient."[17] Yet, I would change the word "theoretical" with "current computational" d as shown below in Section Thermodynamics from GCEeDFT, GCEe DFT methods provide an exact quantum mechanical treatment of electrochemical thermodynamics.
Besides thermodynamics, advanced theoretical and computational models are needed for treating electrochemical (ET/PCET) reaction rates.Traditionally, such rate theories have been developed and used in the context of model Hamiltonians [18e21] to provide a deeper physical/chemical understanding of electrochemical kinetics.However, as discussed below in Sections GCE rate theory and Combining model Hamiltonians and GCE methods, the recently developed GCE rate theory [5] provides a general approach to complex electrochemical ET/PCET kinetics directly from GCEeDFT.Alternatively, the GCEeDFT and GCE rate theory can be merged with model Hamiltonians and associated rate theories to obtain rich atomiclevel insight into electrochemical systems.By doing so, theoretical and computational efforts will hopefully provide even more chemically and physically insightful descriptions of electrochemistry from first principles.

Thermodynamics from GCE-DFT
DFT is a thermodynamic theory where the (free) energy is determined by equilibrium densities as was shown by Mermin for electronic systems [22].Subsequently, DFT was shown to provide a statistical mechanical description and free energies of classical nuclei as well [23].Later, it was established that nuclear quantum effects are captured with multicomponent DFT, which introduces nonadiabaticity and explicit coupling between the nuclear and electronic densities in DFT [24].These theoretical advances showed that nuclear and electronic densities provide an exact thermodynamic statistical mechanical framework and that the (free) energy is obtained by minimizing the DFT (free) energy functional [25].

Standard DFT uses the Hamiltonian b
H for describing the internal energy E(S, V, N)[r] at fixed entropy S, volume V, and the number of particles N as a functional of the density r.In thermal DFT, a Legendre transformation is performed to exchange S to temperature T leading to the canonical ensemble and Helmholtz free energy AðT ; V ; N Þ½b r as a functional of the thermal density operator b r.A further Legendre transformation on the canonical ensemble exchanges the constant particle number with a fixed chemical potential m d that is the GCE where the grand free energy UðT ; V ; mÞ½b r gce is a functional of the grand canonical density operator b r gce .Crucially, the electron electrochemical potential is directly linked with the absolute electrode potential while the electrolyte/solvent electrochemical potentials account for the concentration and interfacial electrostatic potentials; the GCE naturally provides a fixed electrode potential/concentration description of electrochemical interfaces [4].The (free) energy and the electronic/nuclear densities can be obtained from the EulereLagrange equation: where r(r) is the joint nuclear/electronic density in a given ensemble, F is the free energy functional, and d denotes a functional derivative.For small systems typically treated within DFT, the choice of ensemble is absolutely critical as thermodynamic parameters define the state of the system.For large systems, correct expectation values may be obtained in all ensembles but the fluctuations, variations, and other thermodynamic correlations will depend on the ensemble [26].To correctly mimic electrochemical experiments and to capture correct thermodynamics at a fixed temperature, electrode potential, and electrolyte concentrations, GCE provides the theoretically bestfounded and natural framework as depicted in Figure 1.By self-consistently carrying out the minimization in Eq.
(1) within GCE, the free energy and the corresponding equilibrium densities of electrochemical interfaces are obtained from a single self-consistent calculation d there is no need for sampling the phase space.
All nuclear and electronic quantum effects are by construction included in GCEeDFT yielding an exact thermodynamic formulation of electrochemical systems [4].Exact equations can, however, rarely be solved and hence several approximate computational techniques have been developed [3].A systematic coarse-graining Schematic description of a GCE-DFT simulation at constant electron (e m electrons ) and electrolyte (e m electrolyte ) electrochemical potentials.
[4] of the quantum GCEeDFT has elucidated the connections between different solvent/electrolyte, shown what is neglected in different computational approaches, and set a hierarchy for constructing welldefined GCEeDFT approximations at varying levels of detail.For example, the reacting protons and all electrons can be treated quantum mechanically and selfconsistently within the nuclear-electronic orbital DFT [27], while the electrode comprises explicit classical nuclei, whereas the electrolyte/solvent is treated with implicit GCE-DFT models.
Although the constant electrode potential or electron grand canonical treatment was an outstanding problem until 5 years ago [1], now several practical methods exist for carrying out such calculations using either a single self-consistent calculation [6] or an iterative approach [28].The iterative approach is particularly easy to implement taking only around 20 lines of code in the popular ASE simulation interface [29] to turn ASEsupported DFT calculators to GCEeDFT calculators as shown in Figure 2. Yet, there is much need for more efficient and accurate GCEeDFTebased methods to treat the solvent and electrolyte at fixed chemical potentials; the ability to describe the equilibrium structure and thermodynamics of the solvent/electrolyte from a single calculation will be worth the effort.Dielectric continuum (PoissoneBoltzmann) models present the least detailed and most widely used electrolyte description [30], but inherent approximations [3,4] limit their accuracy and scope.More spatially resolved continuum models based on statistical liquid theories, such as reference interaction site model (RISM) [31] or molecular DFT [32], have been combined with electronic DFTand present highly promising alternatives to achieving a good balance between Modifications needed to turn an ASE-supported [29] DFT-calculator to a fixed electrode potential GCE-DFT calculator.Note that a calculator supporting charged supercell calculations is needed.Tested with the GPAW code [33] with adopted boundary conditions and Poisson-Boltzmann methods [4].
accuracy and computational cost.Yet, these methods are still rather new and will likely require extensive testing and benchmarking to establish their performance.
Besides providing a general, systematically improvable approach to electrochemical thermodynamics, GCEe DFT facilitates benchmarking other approaches.For instance, the widely used CHE [11] neglects the explicit effect of the electrode potential and is applicable to PCET thermodynamics only.Recently, GCEe DFTwithin a continuum description was applied [34] to show that the accuracy of CHE decreases for large potentials, high capacitance, and adsorbates with significant dipole moments.GCEeDFT also enables computing electrosorption valencies to address pure ET reactions and to generalize the CHE-like treatment to decoupled PCET reactions as well [35,36].

GCE rate theory
Most (GCEe)DFT studies in electrocatalysis focus on electrochemical PCET thermodynamics while only a handful of works have used GCEeDFT to address electrochemical reaction rates as an explicit function of the electrode potential [2].The potential-dependent rates have been computed using transition state theory (TST) within GCE but the formal proof of GCEeTST and its limitations were only recently established [5] to justify the use of GCEeTST but also to identify ways to go beyond TST.The GCE rate theory [5] is an extension to Miller's general canonical rate theory [37] and has the following form: (3) where the adiabatic potential-dependent reaction barrier DU y ad ;EVB is obtained from the equilibrium adiabatic barrier DU 0;y ad ;EVB , potential-dependent reaction energy DU, and the diabatic barrier, that is, the reorganization energy (L 0 ) at the equilibrium potential.Here all the parameters were obtained directly from GCEeDFT without any fitting and are explicitly potential-dependent.
[5] With this, the barrier height and position as a function of the electrode potential are predicted with a few GCEeDFT calculations.

Combining model Hamiltonians and GCE methods
Various theoretical and computational schemes have been put forward to understand the PCET/ET kinetics and thermodynamics as a function of the electrode potential.Traditional theories use model Hamiltonians to describe the system properties using various effective interactions while more recently DFT-based computational methods have been applied to electrochemical reactions without invoking effective interactions besides those contained in the approximated exchange-correlation functional and the computational setup.Model Hamiltonians are transparent in their physical/chemical interactions and phenomena but the parameters defining these interactions can be ambiguous and difficult to obtain, whereas current (GCEe) DFT methods offer a self-consistent treatment of the reaction environment but are often phenomenologically simple [1].Merging model Hamiltonians with the accurate atomistic (GCEe)DFT description presents a promising way to characterizing and understanding complex phenomena at electrochemical interfaces from the atomic level.Past success for such approaches in heterogeneous catalysis is provided by the powerful dband model [38], which coupled the NewnseAnderson model Hamiltonian [39] and DFT to understand and predict reactivity from electronic structure properties.Electrocatalytic reactions are, however, more complex and will require effects beyond the electronic properties and consideration of finer interactions.
A case in point is given by the alkaline hydrogen evolution reaction (HER) and the Volmer step.McCrum and Koper [40] used a very elegant combination of experiments and DFT calculations to demonstrate that OH adsorption energy is a good descriptor for the alkaline HER.However, the success of the OH binding descriptor cannot conclusively rule out the importance of potential zero charge, water reorganization, that is, interfacial rigidity, or nuclear quantum effects, which have also been used to explain slow alkaline HER kinetics [10,14] beyond adsorption or electronic interactions [41].To quantify the role of potential zero charge and solvent reorganization in alkaline HER, Huang et al. [42] have extended the Schmicklere NewnseAnderson model Hamiltonian [43] to include bond breaking, electrostatic interactions, and solide liquid interface 'rigidity' as shown in Figure 4.A competing approach by Hammes-Schiffer et al.
[44] uses a similar Hamiltonian with a reduced expression for the double-layer effects but treats the transferring proton quantum mechanically to address electronic/ vibronic nonadiabaticity, nuclear tunneling, and solvent dynamics.In both models, the electronic Hamiltonian is used to account for electronic interactions with 3 a as the energy reactive EVB/diabatic state a on H 2 O, k as an index for a quasi-free single electron EVB/diabatic state on the electrode, n i as orbital i occupation, and the interaction term V.Both models yield free energy surfaces (FES) in terms of the electronic interactions, solvent reorganization coordinate q and reorganization energy l, the electrostatic interactions f, and bond forming/breaking (U) Both works treat electronic interactions similarly but they majorly differ in the treatment of the transferring proton and the double layer.As a result, the conclusions on the origin of the sluggish alkaline Volmer reaction are also different: Huang et al. used a classical proton and attributed the slow kinetics to the rigidity of the reaction environment and OeH bond strength, whereas Hammes-Schiffer et al. demonstrated that nuclear tunneling, vibronic nonadiabaticity, and solvent dynamics significantly affect the rate.Other differences can be found in the used parameters: Huang's classical model used a high solvent reorganization energy (w3.5 eV) and high electronic interaction parameter (2 eV), whereas a much smaller reorganization energy (0.35 eV) and electronic interaction (0.9 eV) were used in the quantum picture by Hammes-Schiffer.For the same reaction on metallic electrodes in water, such large differences seem unlikely and some parameter cancellation might be taking place.To resolve these issues, self-consistent (GCEe)DFT and rate theory could be used to resolve the contributions from electronic interactions, solvent, and electrostatic potentials on the HER kinetics.
Although all the complications from quantum effects, interfacial rigidity, OeH bond strength, and so on, can be addressed from GCEeDFT and GCE rate theory directly, physically insightful and hopefully more predictive theoretical models could be obtained by merging GCEeDFT and model Hamiltonians.This would combine the accuracy and atomic level details from DFT with the rich phenomena and interpretability of model Hamiltonians; the goal is theory and computation could contribute both insight and accurate numbers to increase chemical understanding [45].The current state of the art for this approach has been recently demonstrated for Zn deposition [46].
In this scheme, the interaction parameter V can and has obtained by fitting the so-called chemisorption functions to reproduce the DFT-computed density of states [47,42].However, accurate electronic and, for example, electrostatic interactions are already directly available from DFT without any fitting.Addressing solvent effects with DFT is significantly more demanding and here the combination of model Hamiltonians and (GCEe)DFT seems particularly promising.Even though both fast and slow solvent interactions are available directly from implicit solvent GCEeDFT, model Hamiltonians further use the solvent reorganization reaction coordinate to compute the FES to provide deeper insight on solvation effects.The solvent interactions are coined in the reorganization energy l, which has thus far been treated as a free fitting parameter [18], computed using the Born solvation model [42], or using a combination of bulk thermodynamics and molecular simulations [46].Albeit l is an excited state quantity and cannot be obtained from standard DFT, it can be evaluated using constrained (GCEe) DFT [48,5] as has been demonstrated for PCET [5] and ET reactions [49].The solvent description can be systematically improved using a classical DFT description [50] and accounting for nonequilibrium solvent effects using molecular DFT [50] or continuum [51] solvent models.Combining these implicit solvent models with constrained GCEeDFT will hopefully reduce the ambiguity or empiricism of the all-important reorganization energy parameter.
Even more subtle but important interactions, where merging DFT and model Hamiltonians will likely be important, are pH and electrolyte effects [41,10].Huang's model Hamiltonian [42] moves in this direction by including pH and cation effects as effective parameters affecting bond-breaking energies.Such interactions could be evaluated quantitatively with GCEe DFT and translated to effective parameters used in model Hamiltonians or used directly as variables in GCEeDFT computed barriers.This can be realized using the recent GCEeEVB theory [5], which facilitates GCEeDFT evaluation of effective ion interactions on diabatic/EVB states used in Eq. (5).Again, this would make the model Hamiltonian parametrization more transparent and accurate.
Although knowledge of the effective interactions and FES is enough for adiabatic TSTreaction rates, this is not enough for ET/PCET reactions, which often require accounting for various complications discussed earlier [52,2].These are not small corrections to TST rates but can influence the rate by orders of magnitude, which certainly warrants closer theoretical and computational attention [8].For instance, signatures of nuclear quantum tunneling have been recently observed on a metal surface at room temperature for alkaline HER [14] where the proton tunneling is expected to be vibronically nonadiabatic on gold, whereas other electrodes could also have contributions from solvent dynamics as shown in Figure 4 [44].While electronic nonadiabaticity has not been demonstrated for inner-sphere electrocatalytic reactions, outer-sphere ET has been traditionally discussed in the framework of nonadiabatic ET [17].
In all these cases, the computation of ET/PCET rates requires approaches beyond TST to evaluate the preexponential factor k. Thus far practically all studies addressing k in electrochemical systems have used a model Hamiltonian description.However, as discussed in Section GCE rate theory, various computational and theoretical methods to evaluate k for electronically adiabatic reactions already exist as the GCE rate theory allows using previously developed approaches at electrochemical conditions.For electronically nonadiabatic reactions, the situation is conceptually and computationally much more difficult.Model Hamiltonians typically describe electronic nonadiabaticity using Fermi Golden rule transition probability between quasi-free single electron orbitals of Eq. ( 5) where f ð 3 k Þ is the FermieDirac distribution and V ak is the single-orbital coupling constant.Breaking the transition probability to single-orbital matrix elements is very problematic for first principles methods: single-electron orbitals are not unique as they can be (de)localized by a unitary transform without affecting (thermodynamic) observables [5,53].Therefore, prefactors resembling Eq. ( 7) cannot be directly evaluated from DFT-computed orbitals as the nonunitary character of V ak causes arbitrariness when addressing electronically nonadiabatic reactions; translating the quasi-free electron orbital picture from model Hamiltonians to DFT-orbitals is difficult and in principle, there is no need to even include orbitals in DFT.Whilst practically all DFT implementations use orbitals, these are not quasi-free and the index k is not the k-vector used in solid-state chemistry.Furthermore, GCEeDFT already uses the FermieDirac distribution to construct the electronic density (operator) [4,5] and including it again in the rate calculation in Eq. ( 7) is not theoretically well justified.Using unitary invariant many-body wave functions as done for molecular systems [54] leads to different problems as a continuum of electronic wave functions needs to be considered.

Conclusions and outlook
In this work, I have given an overview of GCE approaches for addressing electrochemical thermodynamics and kinetics from first principles.Besides showing the formal background for these methods, I have discussed some applications and future directions where these methods could advance.I have also discussed how GCEeDFT and GCE rate theory could be merged with model Hamiltonian methods to combine the strengths of both approaches to understand and control electrocatalytic reactions at complex electrochemical interfaces.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Figure 4 Left:
Figure 4 These wave functions are not expected to obey FermieDirac statistics, which is valid for single-electron orbitals only.A straightforward solution is using the LandaueZener transition probability in Eq. (2) as proposed in Ref.[5], but this has not been tested in practice and is likely computationally too demanding.Therefore, there is ample room for clever and computationally manageable approaches combining model Hamiltonian and GCEeDFT methods to treat electronic nonadiabaticity.