Converging Divergent Paths: Constant Charge vs. Constant Potential Energetics in Computational Electrochemistry

Using the example of a proton adsorption process, we analyze and compare two prominent modelling approaches in computational electrochemistry at metallic electrodes — electronically canonical, constant-charge and electronically grand-canonical, constant-potential calculations. We first confirm that both methodologies yield consistent results for the differential free energy change in the infinite cell size limit. This validation emphasizes that, fundamentally, both methods are equally valid and precise. In practice, the grand-canonical, constant-potential approach shows superior interpretability and size convergence as it aligns closer to experimental ensembles and exhibits smaller finite-size effects. On the other hand, constant-charge calculations exhibit greater resilience against discrepancies, such as deviations in interfacial capacitance and absolute potential alignment, as their results inherently only depend on the surface charge, and not on the modelled charge vs. potential relation. The present analysis thus offers valuable insights and guidance for selecting the most appropriate ensemble when addressing diverse electrochemical challenges.


Introduction
Electrochemical processes at electrified interfaces play a pivotal role in a plethora of applications, ranging from energy storage systems to catalysis.Assessing the nuances of such processes by means of first-principles-based, atomistic simulations has long been a topic of interest for researchers.][29][30][31] Using the adsorption of a proton at metallic electrodes as a prototypical electrochemical modelling problem, this study assesses the fundamental properties and practical implications of both computational approaches.First we prove that in the thermodynamic limit, both methods offer identical results concerning the differential free energy change.Subsequently, we elaborate in more detail on the advantages of the grand-canonical method in terms of interpretability and size convergence, but also their inferior sensitivity to robustness against certain modelling parameter errors.Consequently, this work offers vital guidelines in terms of validity, accuracy, and applicability of both methods and can aid researchers in selecting the most suitable ensemble for their specific electrochemical modelling requirements.

Canonical, Constant-Charge Description of Proton Adsorption
Within the canonical constant-charge (CC) approach by Nørskov et al. 2,3 a proton adsorption reaction is assessed by evaluating the free energy change of an explicitly modelled electrode surface s while adding a proton-electron pair (effectively a hydrogen atom) consistent with the reaction equation The method thus approximates the derivative dG s dN H and maps the obained results to a free energy change ∆G at electrochemical conditions by including the energy cost to remove the exchanged particles from their thermodynamic baths according to μH + hereby denotes the electrochemical potential of the proton and μe − the electrochemical potential of the electron, which is trivially related to the electrode potential ϕ E via μe − = −eϕ E .The potential dependence ∝ μe − enters thus only in the a posteriori energetic referencing step, eq. ( 3), and not at the electronic structure calculation level, i.e., when evaluating G s .Referencing is typically done against the equilibrium at SHE conditions where .Both methods add essentially dipole-fieldlike energy contributions (see discussion below and section S1 in the SI) and thus yield by and large identical results if N net e − is chosen consistent with the field E. 36 While the so-discussed CHE methods are thus based on canonical energy differences of systems of predefined charge states (e.g.fixed electronic excess charge N net e − ), more recent, constant potential methods approach the proton adsorption process at electrified interfaces differently, by modelling directly the effect of the potential.

Grand Canonical Description of Proton Adsorption
Constant-potential descriptions treat the electronic degree of freedom grand canonically, i.e., they allow the electron number to adapt according to the constant-potential boundary condition.The intellectual difference to the CHE method becomes clearer when we reanalyze the free energy change for the adsorption reaction, eq.(2).Within the constant-potential approach 18,20,21,37,38 energy differences are to be taken for energies G s that were Legendretransformed in the electronic degrees of freedom in order to yield the electronically grandcanonical (GC) energies G s at given electron electrochemical potential μe − , where Here, it is of central importance to recognize that N e − = N e − (N H + , μe − ) is a function of the thermodynamic variables N H + and μe − .The GC energy change of the interfacial system when adding a proton is then obtained as the finite difference which approximates ∂G s ∂N H + μe − at fixed electrode potential and thus allows to evaluate the energy change of the proton adsorption reaction as Note that GC simulations thus simply and exclusively describe a proton adsorption process.
The concomitant number of exchanged electrons is only obtained a posteriori from the simulations at fixed μe − = const, and GC calculations are valid and applicable independent of charge-transfer considerations. 20This is in contrast to the CHE method which assesses the canonical energy difference by a simultaneous exchange of a proton-electron pair H + + e − as expressed by the derivative in the number of hydrogen atoms H as in eq.(3).

Relation between both Descriptions
In order to understand the relation between ∆G and ∆G one needs to analyze how the partial derivative ∂G s ∂N H + μe − can be reexpressed in terms of derivatives of the canonical Gibbs energy G s (N H + , N e − ).Using the Legendre transform, eq. ( 4), one can write where the total derivative (first term) signifies that one has to consider the total change of G s (N H + , N e − ) at constant μe − which arises on the one hand from the change of N H + but as well from the adaption in N e − = N e − (N H + , μe − ) at constant potential conditions.The second term represents the energy cost to take out the respective total number of electrons from the bath which derives from the second term in the Legendre transform in eq. ( 4).
Without loss of generality, the change in N e − at metallic electrodes can be decomposed into where the 1 arises from the exchange of one electron to compensate for the +1 charge of the proton at the interface, and the second term on the right hand side derives from the consecutive adaption of the double layer charge N net e − to keep the electrode potential fixed.In the same logic, one can decompose the first term in eq. ( 7), dG s dN H + μe − , into a first part at constant double layer charge (N net e − = const.)and a second part that accounts for the change in N net e − .Evidently, as we exchange identically one proton and one electron in the former contribution, it is identical to the considered energy difference in the CHE method, and will thus be written for convenience (and consistent with the considerations above) as the total derivative in the number of hydrogen atoms (see section S2 in the SI).We thus obtain from the identity ∂G s ∂N e − = μe − at equilibrium.Insertion of eqs.( 8) and ( 9) into eq.( 7) then yields By inserting this back into eq.( 6) and comparison with eq. ( 3) we then find ∆G = ∆G when canonical energy differences at constant double-layer charge are used and whenever the relation N net e − (μ e − ) is identical in the canonical and grand-canonical treatment.Hence, for metallic electrodes, GC methods are expected to yield identical results with the (E-field-augmented) CC CHE method in the thermodynamic limit, i.e., whenever energy differences become differentials.Importantly though, we stress that this applies only if the proton state under study is identical and well defined (stable) in both ensembles.As we showed recently, one can in principle stabilize states in the canonical ensemble that are not stable against infinitesimal geometric perturbations and thus not observable in the GC ensemble. 39

Advantages and Disadvantages of Both Methods
While thus ∆G = ∆G numerically for calculations in the infinite cell size limit, both methods have advantages and disadvantages in terms of interpretability, size convergence and sensitivity to computational errors.
Interpretability The meaning of ∆G is arguably more evident within the GC picture: Considering eq. ( 6), ∆G simply measures the difference between ∂G s ∂N H + μe − and the proton electrochemical potential in the bulk solution reservoir μH + .Note that one can rewrite this partial derivative as using eq.( 4) and the identity ∂G s ∂N e − = μe − at equilibrium.The right hand side of eq. ( 11) is the standard expression of the electrochemical potential of the proton as it is the partial derivative of the canonical energy G s at fixed total number of all other components of the system, in particular also of the electrons.Hence, we find that the GC energetics directly assesses the proton electrochemical potential μs H + on the surface at applied potential conditions where Hence eq. ( 6) becomes

∆G = μs
where ∆G simply and exclusively describes the difference of the proton electrochemical potential on the surface and in the solution (cf.Ref. 12).In contrast to the canonical case, cf.eq. ( 3), the electronic degrees of freedom only play a bystander role as they are fully considered and integrated into the electrode potential dependence of μs H + by means of the Legendre transform.These facts imply on the one hand that ∆G remains meaningful independent on the charge state of the surface and simply describes the proton energetics equally in a chemisorbed state, in interfacial water or along an electrochemical reaction path.In contrast, the imposed charge balance in reaction equation ( 2) that underlies the simulated, simultaneous exchange of protons and electrons within the CC CHE method, makes the association of ∆G with the proton electrochemical potential non-obvious.Nonetheless, our rigorous finding that ∆G = ∆G proves the correctness of CC CHE calculations that assess such states (with partial electron transfer 40 ) by infinite cell size extrapolations of canonical CHE calculations at constant double layer charge N net e − . 9 Size Convergence Such cell size extrapolations of canonical calculations are always necessary, when the considered adsorbate or adsorbate configuration induces a significant dipole in the system, as the native, canonical DFT energy obtained in finite-size periodic boundary condition (PBC) supercells then exhibits significant, unphysical interactions with the periodic images.The significant magnitude of these finite-size errors in practically employed supercells was first recognized by analyzing the energy difference E DFT (H + = ads) between a desorbed and an adsorbed proton, as relevant in the study of proton adsorption in canonical, CC studies in a finite simulation cell size of area A. [7][8][9]41,42 Along such a reaction path a dramatic change in the system's work function is observed that is associated with a change ∆ A ϕ 0 of the potential drop across the inner double layer.
It was found 7,8,41 that ∆ A E DFT exhibits a significant cell size dependence with ∆ A E DFT ∼ ∆ ∞ E DFT − 0.5e∆ A ϕ 0 where ∆ ∞ E DFT describes the energy difference between the desorbed and adsorbed state in the infinite cell size limit and −0.5e∆A ϕ 0 the approximate cell size dependence of the energy difference.The latter term is explained by the energetics of (dis)charging a plate capacitor, created by the proton's charge when in the solvent and the according counter electron on the surface which can be considered as a (partial) charging of the double layer by ∆N net e − ∼ 1.For a metallic electrode of area A and an inner layer capacitance C this partial charging of the capacitor originates in the potential change ∆ A ϕ 0 = −e/(CA) when the proton is in the desorbed state and the size-dependent scaling of ∆ A E DFT is thus explained by the change in the capacitively stored energy which amounts to While these insights and results were obtained based on an original internal energy referencing scheme, 7,8,41 the size convergence of standard, canonical CC energy differences replicate the (approximately) linear dependence of ∆ A E DFT on ∆ A ϕ 0 though with the magnitude of the prefactor being slightly smaller (∼ 0.3e − 0.4e, 9 see Fig. 1).
The capacitor expression, eq. ( 15), not only provides a basis for extrapolating ∆ A E DFT to the infinite cell size limit ∆ ∞ E DFT where ∆ A ϕ 0 A→∞ = 0, it also allows to estimate the considerable magnitude of finite size effects in CC calculations: For a prototypical, e.g.It is important to note here, that while ∆ A E DFT in the respective studies refers to the energy difference between a desorbed and an adsorbed proton, the finite-size scaling is largely related only to the artificial PBC interactions of the desorbed proton state with its strong dipole. 9,28The adsorbed proton state typically only induces a minor work function change 21,43 and thus exhibits intrinsically only small PBC interactions.
In constant-potential simulations, any such finite-size errors are absent, as here any adsorbate dipole is by construction identically canceled by the counter-dipole induced by the adaption of electronic surface and electrolyte double-layer charges (N net e − ) in order to keep the electrode potential constant.As a result, the constant-potential boundary conditions naturally compute the energetics at ∆ A ϕ 0 = 0, and hence the PBC error ∆ A E cap = 0 at any Figure 1: Observed cell size dependence for the energy difference ∆ A E DFT between a desorbed and an adsorbed proton.Datapoints taken from the literature (dots, Chan2015 9 ) describe the interface solvated by a static, full water layer.Results from own calculations in mixed implicit-explicit solvent environments (crosses and stars) employ only a single, static water molecule (crosses) or a static water cluster (stars), respectively.The lower panel reports the averaged size dependence as derived from linear fits through the data in the upper panel.See section S3 in the SI for more information.cell size.Any remaining size effects are only due to higher-order multipole terms.This fact can be seen most clearly by considering that the GC energy difference ∆ A E DFT between a desorbed and an adsorbed proton in cell size A in a constant-potential setup is given in lowest orders in the potential by the canonical energy difference ∆ A E DFT PZC at the PZC plus the energy change ∆ A G s exc,DL (U ) due to adaption of the double-layer charge.
As derived in Ref. 21 1 with the applied potential relative to the PZC of the clean slab.While the first term in eq. ( 16) represents a dipole-field-like interaction term 21,29,30 for the desorbed proton which is largely size invariant (A∆ A ϕ 0 ∼ const., see also section S1 in the SI and the discussion below), the second term in eq. ( 16) identically cancels the PBC error term in eq. ( 14) that is naturally present in the canonical PZC energetics ∆ A E DFT PZC .These results underline that constant-potential calculations are intrinsically free of the significant, capacitor-like PBC interaction errors that are present in canonical CC calculations.
As a result, while the partial differential in eq. ( 6) can be directly evaluated by finite differences of GC energies in moderate cell sizes (L ∼ 10 Å), more care has to be taken when approximating the derivative in eq. ( 3) from finite differences of CC energies at such cell sizes.
Sensitivity to Modelling Errors While the discussion up to now only revealed disadvantages of canonical CC setups for studying the energetics at electrified metal-water interfaces, they do as well exhibit a variety of non-obvious advantages.There is first, of course, their more straightforward implementation in actual simulations, e.g.molecular dynamics setups.
Instead of setting the potential μe − explicitly, which necessitates (the non-trivial) implementation of a potentiostat and concomitant mitigation of existing convergence issues, 23,25,26 canonical CC calculations only necessitate a method to vary the surface charges N net e − .This can be achieved either using a continuum electrolyte model to screen the electronic excess 1 Cf.eq. ( 14) in Ref.21, where states α and β there are associated with the desorbed and the adsorbed proton state, here.In principle, U in eq.16 is the applied potential relative to the PZC of the state β, i.e., the slab with adsorbed proton.However, U corresponds approximately to the applied potential relative to the PZC of the clean slab as the adsorbed proton has no (significant) dipole. 21As a result, we directly associate U with the applied potential relative to the PZC of the clean slab for reasons of consistency with the consecutive discussions.
charges on the surface 35,38,44,45 or by explicitly introducing real or artificial ions. 6,7,36The induced energy change in such CC calculations can be described in general by a dipole-field-like energy contribution 3,21,27,29,30,33,35,36,44,46,47 according to where ϵ 0 A is the average interfacial field and D = −ϵ 0 A∆ A ϕ 0 the adsorbate-induced surface dipole as deduced from the observed potential change ∆ A ϕ 0 in a cell of area A (see also section S1 in the SI).In case of an identical charge vs. potential relation as in the GC calculation with −eN net e − = ACU , one obtains ∆ A G s DF = ACU ∆ A ϕ 0 which is cellsize-independent and equivalent to the first term in eq. ( 16).These equivalent first order potential dependencies reflect nothing but the above prove of the equivalence of CC CHE and GC energetics in case that the charge vs. potential relations are identical and in case that PBC errors are negligible -aka in the infinite cell size limit.
Note, however, that the CC energetics in eq. ( 17) depends only on N net e − in agreement with a dipole-field energy that only depends on the interfacial field (∝ N net e − /A) and the adsorbate dipole (∝ A∆ A ϕ 0 ).In other words, ∆ A G s DF does not explicitly depend on U or C -in contrast to eq. ( 16) -and is thus rather insensitive to the counter charging method employed and the simulation-intrinsic capacitance C. It is thus largely irrelevant how the variation in N net e − (or E) is achieved in practice, which makes CC energy differences robust against variations in the counter-charge model and errors in the simulation-intrinsic interfacial capacitance, as demonstrated by Gauthier et al. in their assessment of the energetics of the proton adsorption problem. 27Here, consistent energetics was observed when evaluating CC energy differences as a function of N net e − but not when evaluating the energetics as a function of the system-intrinsic electrode potential U .Evidently, this robustness of CC energy differences to approximate in eq. ( 3) derives from the significant energy (and thus error) cancellation when subtracting canonical energies G s for N net e − = const.and thus in an environment of identical electrolyte composition.This fact has important implications, in particular, as the simulation of realistic capacitance curves C(U ) and hence N net e − (U ) is difficult in contemporary atomistic simulations: Typical implicit-model-derived capacitances are only approximately correct, 18,48,49  In contrast, such errors can be (largely) mitigated by evaluating CC calculations not for the simulation-intrinsic N net e − (U ) relationship, but for an alternative, more realistic N net e − (U ) relation.This is possible as the electrode potential enters the CC energetics not at the electronic structure calculation level but only as a post-processing step when transforming the dipole-field energy, eq. ( 17), into a function of U , i.e., via −eN net e − (U ) = ACU .Similarly, one can derive from one calculation at given N net e − , the energetics for a variety of alternative interfacial conditions -e.g. for different ion concentrations or ion types -by assessing their impact on C and using an appropriately altered −eN net e − (U ) = ACU relationship as e.g.done in Refs.35,50,51.It is important to understand that using N net e − (U ) relationships that are different from the simulation-intrinsic ones is not merely an uncontrolled approximation but essentially encodes nothing but the principle of locality: While the potential U and the capacitance C are global properties that certainly affect the energetics as a function of U , the local energetics of adsorbates are affected (largely) by the local interactions.These are in good approximation given by the adsorbate-induced surface dipole D and the effective electrostatic field E in the vicinity of the surface 29 which is essentially given by the electronic excess charge N net e − on the surface.While the structure of the diffuse double layer certainly affects C and hence the constant-potential energetics, the principle of locality suggests that it should not affect the local energetics.
In summary, canonical CC calculations are natively less sensitive to the specifics of the double layer structure than GC calculations and are hence less prone to the errors introduced by their approximate description in actual computations.As a result, the CC approach allows to explicitly compute the stability of adsorbates under practically convenient conditions (e.g.high electrolyte concentration) and then to infer their potential dependence for other conditions in a post-processing step, by simply assuming an altered N net e − (U ) relationship. 35onetheless, it is important to keep in mind, that changes in the electrolyte might not only alter the overall charge vs. potential relation.They could for instance also lead to more specific interactions of adsorbates with (partially solvated) electrolyte ions, 52-54 which would then not be captured by a simple rescaling of N net e − (U ) in a post-processing analysis.As a final remark, while yielding correct answers for the differential energetics, any energetic change from the self-consistent adaption of interfacial fields due to adsorbate-induced work function changes 55 is evidently also not included in the CC energetics whenever generic charging relations N net e − (U ) relative to the PZC of the clean surface are used.In contrast, these effects are included naturally in GC evaluations of the interfacial energetics at finite coverage and potential 18,50 which is in particular relevant for understanding the (non-differential) interface energetics at high coverages with adsorbates that strongly affect the PZC.

Conclusions
Canonical CC and GC constant-potential calculations have emerged as the two prevalent approaches for first-principles atomic-scale simulations of metal electrodes in electrochemistry and -catalysis.Primarily due to practitioner-level approximations diverging results led initially to a separation into two schools, each with considerable scepticism about the other approach.With time, this situation improved, with general considerations on the convergence of different ensembles in the thermodynamic limit, empirical calculations 27 and low-order Taylor-expansions of canonical and GC energetics around the PZC 21,[29][30][31] providing increasing evidence that consistent results can be expected for the energetics of single adsorbates within CC and constant-potential calculations.The rigorous proof provided here finally provides a solid basis for a differentiated assessment of the equivalence and equal validity of both methods.
While fundamentally providing equivalent results in general, a remaining notable difference is the possibility of observing states in CC calculations that are unstable in the GC ensemble. 395][26] Nonetheless, they are more prone to finite-size errors than GC simulations in contemporary supercell setups.The latter point needs to be carefully addressed to obtain consistent differential energetics where ∆G = ∆G holds.
static interfacial water models miss the appropriate, dynamic solvent sampling and all-explicit, ab initio simulations are only possible for very high electrolyte concentrations (typically of the order of 1 M ∼ 1 ion pair per 55 water molecules).With ∆ A G s exc,DL (U ) ∝ ACU • ∆ A ϕ 0 depending explicitly on U and C in GC setups, any imperfect replication of the true charge vs. potential relation leads to errors in the GC energetics.

1 2
H 2 (g) ⇋ H + (aq)+e − (m) and thus μSHE H + + μSHE e − = 1 2 µ H 2 (g) .It allows thereby straightforward referencing of G s against µ H 2 (g) within the so-termed CHE method.Conventionally, one approximates dG s dN H ∼ dG s PZC e − ), i.e., in supercell calculations without electronic excess charges.If ∆G shows a more complex behavior, e.g. a dependence on pH, potential, or electrolyte composition one is bound to make dG s dN H a function of the respective variables.The CHE framework then posits that such effects can be computed by setting up the appropriate interface conditions, e.g. with included interfacial electrostatic fields E or electronic surface excess charges −eN net e − .Such E-field-augmented 32 CHE calculations then evaluate canonical energy differences at constant, finite fields3,33,34 dG s dN H PZC from the energy difference of DFT total energy calculations whose composition differs by exactly one H atom (= H + + e − ) under potential of zero charge (PZC) conditions (μ