Grand-canonical approach to density functional theory of electrocatalytic systems: Thermodynamics of solid-liquid interfaces at constant ion and electrode potentials

Properties of solid-liquid interfaces are of immense importance for electrocatalytic and electrochemical systems but modelling such interfaces at the atomic level presents a serious challenge and approaches beyond standard methodologies are needed. An atomistic computational scheme needs treat at least part of the system quantum mechanically to include adsorption and reactions while the entire system is in thermal equilibrium. The experimentally relevant macroscopic control variables are temperature, electrode potential, choice of the solvent and ions and these need to be explicitly included in the computational model as well; this calls for an thermodynamic ensemble with ﬁxed ion and electrode potentials. In this work a general framework within density functional theory with ﬁxed electron and ion chemical potentials in the grand canonical ensemble is established for modelling electrocatalytic and electrochemical interfaces. Starting from a fully quantum mechanical description of nuclei and electrons, a systematic coarse-graining is employed to establish various computational schemes including i) the combination of classical and electronic density functional theories within the grand canonical ensemble and ii) on the simplest level a chemically and physically sound way to obtain the (modiﬁed) Poisson-Boltzmann (mPB) implicit solvent model. The detailed and rigorous derivation clearly establishes which approximations are needed for coarse-graining as well as highlights which details and interactions are omitted in vein of computational feasibility. The transparent approximations also allow removing some the constraints and coarse-graining if needed. We implement various mPB models in the GPAW code and test their capabilities to model capacitance of electrochemical interfaces as well as study di ↵ erent approaches for modelling partly periodic charged systems. Our rigorous and well-deﬁned DFT coarse-graining scheme to continuum electrolytes highlights the inadequacy of current linear dielectric models for treating properties of the electrochemical interface.


I. INTRODUCTION
Electrochemical reactions take place at the interface between electronic and ionic conductors (electrolyte). Together these two macroscopic conductors form an electrode and electrochemical experiments probe the properties of this interface under the influence of voltage and current between electrodes. A typical current-voltage response obtained from an electrochemical experiment is often di cult to interpret from an atomistic perspective and modeling is at the core of electrochemical analysis. Besides the reaction under study, the outcome of electrochemical experiments is controlled by the temperature, choice electrolyte and the electrode potential which determine the thermodynamic state of the system. 1 The experimental setup sets the stage for computational approach where macroscopic variables like the temperature and chemical potentials of the electrolyte and the electrode need to be controlled similarly as in experiments. For macroscopic systems, the most natural a) Electronic mail: karoliina.honkala@jyu.fi b) Electronic mail: marko.m.melander@jyu.fi thermodynamic ensembles at thermodynamic limit are either the Gibbs G(N, T, P ) or the Helmholtz A(N, T, V ). The large number of species ensures that the, e.g., overall chemical potential of the electrolyte or electrode does not change during a reaction. For a small number of atoms as typically treated in a microscopic, atomistic simulation the grand-canonical ensemble ⌦(µ, T, V ) is a natural choice to allow treatment of a system at fixed chemical potentials with fluctuating particle numbers. Combining the experimental variables with the most widely electronic structure method, density functional theory (DFT), motivates the development of grand canonical DFT (GC-DFT) to enable calculations at fixed chemical potentials and temperatures as done in experiments.
GC-DFT o↵ers a way to study electrochemical microscopic systems in thermodynamic equilibrium characterized by long time and length scales. If some parts of a system are prone to charge transfer events and chemical reactions, short time and small length-scales need to be included in the model. Thus, a general electrochemical modelling approach requires a quantum mechanical treatment of the electrode and reactants combined with a statistical treatment to capture changes in the interfacial charge distribution from the formation of the double-layer and liquid structure, for example. As the potential distribution controls reaction kinetics and thermodynamics, the electrolyte and its e↵ect on the reacting system cannot be ignored. The mixture of several time-and length scales and the need for a grand canonical description present serious challenges for the atomistic model.
Need for constant ion and electrode potentials (grand canonical) rather than constant particle number treatment (canonical) are problematic for atomic-scale modeling and require advanced simulation schemes. 2 In electronic GC-DFT, the electrode potential is controlled by the Fermi-level of the electrode which equals chemical potential of electrons while the electrolyte chemical potential depends on the electrolyte solution and its concentration. Thus, a well-defined DFT model of an electrochemical interface calls for a GC ensemble with a fluctuating number of electrons and ions at a given temperature rather than the commonly applied canonical DFT where the particle number is fixed but the chemical potentials are allowed to fluctuate.
To obtain ensemble averages from atomistic simulations, the most straight-forward solution is to adopt a statistical treatment of the system by using (ab initio or QM/MM) molecular dynamics (MD) where the reactants, electrode and the electrolyte are treated explicitly. To reach an equilibrium state, both the ionic and electronic degrees of freedom need to be sampled extensively. While this approach is well established [3][4][5][6] and justified, reaching an equilibrium requires long simulations on large systems making MD applicable to small systems and limited amount of reactions.
Another formally exact approach to model the ensemble e↵ects, is to use density functional theory where free energy of the entire system is defined in terms of the nuclear and electronic densities. After coarse-graining, classical parts can be fully described by atomic density distributions of di↵erent species. [7][8][9] Classical DFT is naturally formulated in the GC ensemble where the chemical potentials are fixed by the (pure) bulk liquid. By combining GC electronic and classical DFT, a unified DFT theory can be devised; previously this approach has been used to connect classical DFT with canonical DFT 5,10-12 but a more generalized grand canonical DFT framework for both electrons and ions simultaneously has not yet been established; such a general framework is formulated in this work.
Like the classical DFT, also electronic GC-DFT can be realized. 13,14 The most elegant practical approach to electronic GC-DFT is to fix the Fermi-levels and allow the fluctuation of particles during an SCF calculation. While stable algorithms for a constant electron potential (or Fermi-level) calculations have been developed recently 15,16 , usually the fluctuating number of electrons usually causes significant di culties in practical calculations due to poor convergence. Another, and equivalent way, to obtain grand canonical energies is to perform calculations at several points using a fixed number of electrons and interpolate to the desired Fermi-level. 2,[17][18][19][20][21][22][23] .
Besides fixing the chemical potential of the electrons, also the chemical potential due to the solvent and ions needs to remain constant 18,20 which can be included in either the quantum mechanically or classically.
Changing the Fermi-level corresponds to changes in the electron chemical potential which is obtained by altering the charge state of the electrode. This presents a challenge for actual computations because the electrocatalytic system are usually partially periodic and periodic systems need to be charge neutral to avoid divergence of the electrostatic interactions. Typically a homogeneous back-ground charge is introduced and various correction schemes have been devised 17,24,25 to relate charged and charge-neutral systems to a common reference. An elegant solution is provided by joint DFT 10,15,22 and (modified) Poisson-Boltzmann (mPB) implicit solvation models [26][27][28] where the charge neutrality can be maintained by the ionic distribution in the double-layer region. However, in non-linear mPB models the charge neutrality constraint is not automatically fulfilled. 29 Therefore, in practice even the mPB models as applied to charged periodic systems need to either enforce charge neutrality by using Lagrange multipliers 29 , use simple background charges 30 , or to allow charged systems in partially periodic system by using modified mPB boundary conditions 19 mimicking neutralization by image charges. These di↵erent ways of simulation charged periodic system are investigated and compared in this work.
In this paper, we present a formal way to obtain free energies of solid-liquid interfaces at fixed electrode and electrolyte chemical potentials by combining nuclear and electronic DFT within the grand canonical ensemble. The theory is turned into practical computational schemes by gradually coarse-graining the initially fully quantum mechanical description of both nuclei and electrons to i) classical nuclei and ii) eventually a continuum description of the solvent and ions. Besides a general framework, the detailed derivation of the electrochemical GC-DFT makes transparent approximations for particular computational models. Within mPB models we consider a general reference electrode to allow modelling of electrochemistry at fixed absolute electrode potentials. We also review and develop practical computational schemes for modelling charged periodic systems. Finally, several standard and advanced mPB models are implemented in the GPAW software and tested towards double layer capacitances and interfacial potential distributions. In the future, the developed general general framework will be extended to study electrochemical reactions and kinetics, work in this direction is already being made. Modelling of electrochemical systems with quantum chemical methods is commonly performed in the electronic canonical ensemble where the number of electrons and nuclei is fixed and the system is held at a constant temperature and volume. Thermodynamic quantities in the canonical or Helmholtz ensemble are given by 3 where the density operator iŝ and p i is the probability of an electronic state with energy E i corresponding to state | i i with a constant number of particles N i HereĤ el = 1 2 n v(r n ) =T n +V nn +v is the electronic Hamiltonian and i are the electronic wave functions. Using these definitions the electronic Helmholtz energy is written using the density operator 13,14 Minimizing the electronic free energy leads to 14 where is the universal Hohenberg-Kohn functional and v(r) is the external potential. T is the kinetic energy from Tr h⇢T i . Similarly, V nn is the total interaction energy and S is the entropy .

B. Joint DFT
An extension of the canonical electronic DFT is the joint DFT (JDFT) 10,22 which is a combination of nuclear density functional theory 7,31,32 and canonical electronic DFT 33,34 . The combination of these two is based on Mermin's extension of electronic DFT to finite temperature 13 and multicomponent non-Born-Oppenheimer approach 31,32 to incorporate nuclear density in DFT. The application of JDFT is to establish an approach in which a fixed number of electrons interacting with an environment formed by a fixed number of nuclei is described. JDFT 10,22 is coined in a doubleminimization of the free energy with respect to the total electronic density n(r) and the nuclear density N (R), to yield the Helmholtz free energy for the total system as where V (R) is an auxiliary external potential conjugate to N (R) (see Refs. 31 and 32 for details). n is the usual electron density and v is the external potential.
Here we use to the definition of nuclear densities as defined by Capitani 31 where the nuclear density and the conjugate potential are one body objects. A more complete definition 32 would be to use N-body nuclear objects in the body-fixed frame to remove the ambiguities related to translationally invariant Hamiltonians. For our purpose, both should give identical results and we have chosen to use the simpler one body formulation of Capitani. The functional becomes being the density operator for a combined nuclear and electronic system with fixed number of nuclei and electrons. v(r) is a true external potential acting on the electrons while V (R) acts on the nuclei. Note that here the total Hamiltonian is used and H tot =T n +V nn +T N +V Nn +V NN + P NV N + P nv n with N for nuclei and n for electrons. Furthermore, the functional F [n(r), N(R)] is not universal because it depends on the electron and nuclear densities 32 . Also the external potentials are not generated by the ions as in usual DFT but rather by a genuine external field. When Kohn-Shamification is performed for molecular systems 32 ,V nn =V Nn = 0, for the nuclei the auxiliary potential V (R) =V NN , the e↵ective potential is

C. Grand canonical JDFT for electrochemical systems
While the JDFT approach o↵ers an exact picture of the Helmholtz canonical free energy, typical systems treated within (electronic) DFT are very small due to the computational burden of such calculations. Even if the environment formed by the nuclei is treated classically, the systems are far from the thermodynamic limit. Therefore, in a microscopic simulation, the thermodynamic state characterized by macroscopic variables (such the chemical potentials and temperature) of the system can change due to chemical reactions or fluctuations of the environment 35 . To describe equilibrium quantities of a small system in contact with an environment it is more natural to fix the chemical potential of the environment rather than the particle number of the system. For a large enough systems, both choices naturally equivalent.
Fixing the chemical potential corresponds to a grand canonical ensemble. Performing a Legendre transformation in terms of the chemical potential on the canonical free energy gives the grand canonical Hamiltonian 3 ⌦ =Ĥ tot TŜ Nn X iμ iNi (9) so that the grand canonical free energy is ⌦ = ⌦[T, V, {μ i }] specified by the (electro)chemical potentialsμ i and particle number operatorN i . Now the particle number can fluctuate which leads to the following density operator 3⇢ We also define the partition function ⌅ = Tr In the above equation, | i i is the total wave function of both the electrons and nuclei so that the particle number operatorsN i corresponds to electrons or the nuclear identities as specified below.
In the considered electrocatalytic framework as applied to a typical DFT calculation, the number of atoms at the electrode, number of solvent molecules as well as the number of reactants are fixed and we do not consider the DFT system to be open to these species. We therefore require, that the only the chemical potentials of ions are fixed by the bulk solution and that the chemical potential of electrons is fixed by the electrode. These constraints lead to a partial grand canonical DFT in terms of the ions and electrons with fixed chemical potentials.
Before advancing further it needs to be noted that the particle number operators for the ionic nuclei correspond to only nuclei, not electrons attached to the ion. Still, one can simplify the treatment of the nuclear chemical potentials, as only a mean ion chemical potential is needed. This is because experimentally the concentration of ions in the bulk solution needs to respect electroneutrality and, thus, the concentration of a nucleus corresponding to an ion cannot be changed independently. 36 Due to this restriction, chemical potentials in two phases (here the bulk b and interface ⌘) are bound by the constraint for cationic nuclei (+) and anionic nuclei ( ). For a binary electrolyte it holds that where v i is the stoichiometric coe cient,μ i = µ 0 i + q i (r) with q i as the charge and (r) as the Galvani (inner) potential. of the two phase are di↵erent leading to imbalance in concentrations. In the liquid phase we can take b = 0 as will be discussed in the Section II G. Following Guggenheim 36 , the mean activity coe cient for ⌘ phase is The mean absolute activity is defined in terms of the concentration c and standard chemical potential µ 0 = ln 0 where 0 depends on the solvent, temperature, ions etc. (later it will realized that 0 is related to the thermal de Broglie wave length). The chemical potential of nuclei for the ions forming the electrolyte is For this reason, the relevant chemical potential of the ions in the electrolyte is the bulk chemical chemical potential µ b ± = µ N ± where N reminds that at the point the ionic chemical potential is only for the nuclei of the ions, not ions with a nucleus and electrons related to it.
The grand canonical free energy for an electrochemical systems with fixed ionic nuclear ( ± ) and electron ( n )chemical potentials is where T k b p i ln p i = T S as p i ln p i = S is the von Neumann entropy 37 . Now we turn to the minimization of the partial grand canonical energy. This is done following Levy's constrained search 14,38 with the grand canonical Hamiltonian. Recalling that our aim is to describe an electrode treated at the quantum level in contact with a classical electrolyte, we will split the total nuclear density from all the nuclei N in the system to well defined contributions of nuclei comprising the solvent S, nuclei of ions the solution I, nuclei in the electrode M , nuclei of possible reactants R. Also the total electron density n is split to contributions from the electrolyte ions n ± (r), electrode n e (r), and the solvent molecules n s (r). Furthermore, the positions of the atoms constituting the electrode and the reactants are fixed and considered explicitly. Thus, we choose that the external potential v(r) is defined by the electrode and reactant nuclei and given by The first definition corresponds to an e↵ective nuclear potential due to a quantum nuclei a which is approximated by the classical point like nuclei having the charge Z a . By specifying the external potential this way, the corresponding quantity needs to be removed fromĤ tot . By this choice part of the system is to be considered fixed providing a static external potential to the electrons and rest of the nuclei. This is reminiscent of a Born-Huang (BH) approximation for fixed nuclei with positionsR M,R and momentaP M,R (in the absence of non-adiabatic effects the Born-Oppenheimer, BO, approximation is sufficient). Thus, the total wave function is approximated as a product of the fixed (not yet classical) nuclei of the reactants/electrode and the mixed wavefunction of ionic nuclei, solvent nuclei and the electrons of the entire system: where the nuclear-electronic wave functions are written using the BH expansion to factor the wave function with as the nuclear wave function for the reactants/electrode and for the combined electronic and solvent/ion nuclear wave function. As a result, the electrons and solvent/ion nuclei move on the potential energy surface generated by the nuclei of the electrode and reactants. While the BH separation has been assumed, the adiabatic theorem (see for example section 5.1 in Ref. 39) has not been invoked and therefore nonadiabatic effects are still present. Last, either diabatic or adiabatic ion/electron wave function can be utilized in the above expansion; the latter is more suitable for general chemical reactions while the former may prove to be convenient for the study of electron and proton transfer reactions, for example.
To further simplify the situation, the nuclei of the electrode and reactant are assumed classical. Such a division is obtained by splitting the total partition function as a restricted partition function⌅ 40 with fixed classical nuclei with momentumP and subsequently sampling the phase space of the classical nuclei. Thus, the total quan-tum grand canonical ensemble is approximated as a system where the open quantum system interacts with the fixed, classical nuclei. For this system the grand partition function is where the quantum mechanical nuclear kinetic energy has been replaced by its classical equivalent and is the the thermal de Broglie wave length unique for nuclear type. The grand free energy from above is for classical nuclei of the electrode and reactants interacting with the quantum system of electrons, solvent and ionic nuclei. The classical mechanics for reactant and electrode nuclei is implied by the use classical partition function for these degrees of freedom. The free energy contributions from the classical part can be obtained by thermal sampling of the R, M degrees of freedom while free energy of the quantum part is obtained directly from density functional theory. It's also worth noting that electronic nonadiabatic e↵ects are still included and can be computed schemes from mixed quantum-classical dynamics 41 but then one does not integrate out the momenta.
The formal quantum/classical separation leads to an external potential which is defined by the last quantity in Eq. (15). If some nuclei of the reactants or electrode need to treated quantum mechanically (e.g. in tunnelling calculations), these selected nuclei can be modelled using the nuclear-electron orbital DFT (NEO-DFT) 42 . Also the zero point energies for the reactant and electrode nuclei can approximated by e.g. solving the nuclear wave functions on the GC-BO potential energy surface.
After the BH separation with fixed classical nuclei, the minimum grand canonical energy for the quantum system of electrons, ions and the solvent reads where The above yields the exact grand canonical energy for a grand canonically treated quantum subsystem interacting with classical nuclei setting the external potential. Together equations (17) and (18) form general GC-DFT expression for electrochemical systems at fixed electron and ion chemical potentials.

D. Classical/Electronic GC-DFT for electrocatalytic systems
Of course the exact functional in Eq. (19) is highly complex and needs further simplification. A conceptual inconvenience of the above formulation is that the electrons are completely smeared throughout the entire system and the electron chemical potential acts on all electrons. Chemically and physically it is reasonable that deleting/inserting an electron from/to di↵erent parts of the system should a↵ect the energy di↵erently. For example, removing a electron from the electrode or from an ionic nuclei are expected to impact the energy in different way. Similarly inserting ions at the surface, in the double layer, or in the bulk have a di↵erent impact on the energy. On the contrary, a globally constrained electron chemical potential prevents controlling the electron potential independently for the ions, solvent and the electrode which would be convenient for a practical computational scheme; without this partition the charge of the ions cannot be accurately specific according to chemical intuition nor can the surface charge of an electrode be specified, for example. This is especially problematic if the system is further coarse grained to treat ions and solvent molecules as well-defined entities.
The above considerations call for a way to identify how changing concentration or particle number at a specific position changes the energy. Such a spatially resolved energy dependence can be achieved by introducing position dependent chemical potential which measures how a change in the local concentration c i (r) changes the grand energy. A position dependent chemical potential also allows defining regions with di↵erent electronphilicility or ionophilicility. 43 Accordingly, the position dependent chemical potential is with Here the integration volume V is a point in solvent volume for continuum ion models such the Poisson-Boltzmann model, the e↵ective ion volume in explicit particle simulations, or the electrons at the electrode surface in contact the bulk, for example. Such a formal definition allows local partition of the nuclear and electronic densities to chemically meaningful quantities such as molecules as done in molecular grand-canonical ensemble theory 44 . We utilize the spatial dependence of the chemical potential to partition the total electron density to to such as electrons at the electrode and ions with a specified charge confined in the solution.
To allow the electrochemical potential of electrons independently of ion chemical potential, the electron density is partitioned to well defined contributions from the electrode and reactants as well as from the ions and solvent. This separation is achieved by splitting the total nuclear and electronic densities to the electrode/reactant and ion/solvent contributions. We also demand that the partitioned densities sum up to the the densities as done in subsystem DFT 45 : n(r) = n ± (r) + n e (r) + n s (r) and Because the chemical potentials are spatially defined, we can confine them to specific regions. We use this freedom to partition the total volume to quantum V q (comprising of the electrode and reactants) and classical V c volumes (comprising of the ions and solvent). While this not necessary, the partition to "quantum" and "classical" regions reflects the choice of treating the electrons from the electrode and reactants as quantum mechanically and ions/solvent classically with confined electron densities and helps if the two subsystems are to be treated separately. Also with this choice, one can define an ion chemical potential µ ± which corresponds to electrons and an ionic nuclei bound together. Ionic nuclei and electrons in equilibrium form an ion with a specific charge Z + or Z for cations and anions, respectively. As an example of the ion chemical potential µ, for cations in equilibrium we have: such that R Vion dr(µ n (r) + µ N + (r)) = R Vion drµ + (r) = µ + and R Vion dr(n(r) + N + (r)) = R Vion dr⇢ + (r) = Z + . This allows one to specify ions with a well-defined charge. A similar splitting can be performed for the solvent molecules as well. Finally, we introduce these split densities and volumes with the spatial chemical Eq. (18) to obtain with the ionic charge density ⇢ ± (r). Next, taking inspiration from the JDFT approach 10 the electron density due to the solvent and ions is integrated out. This leads to where the solvent charge density ⇢ S has been introduced and The last two equations yield the grand free energy of a microscopic electrode-electrolyte interface where the ion and electron chemical potentials can be explicitly controlled independently and where the solvent and ions have well-defined charges needed for classical simulations of the electrolyte.
The goal of the minimization was to obtain separation of electrons belonging to the electrode and the electrolyte in order to manipulate the electrode potential and ion chemical potentials separately. Furthermore, the solvent and ions were constructed so that the ion and solvent molecules have well-defined charges. This allows for example the treatment of ionic liquids where the solvent molecules can carry net charges. Care was also taken to partition the constraints due to chemical potentials in chemically meaningful entities.
It is important to stress that integrating out the solvent and ion electrons prevents exchange of electrons between the solvent/ions with the electrode/reactants and prevents the formation of a covalent bond between the ions/solvent with the quantum region; only electrostatic, polarization and dispersion-like interactions can be captured with the presented partitioning. Note that a similar situation is encountered whenever e↵ective (atomic) charges are predefined for ions or solvent molecules. Prominent examples of such are classical DFT within the reference interactions site method (RISM) approach 5 discussed below and detailed in the supporting information of Ref. 5 or in the force fields used in classical molecular dynamics simulations. This has direct implications on the applicability mixed classical/quantum systems. In the current setting, the most important restriction is that specific adsorption of the solvent/ions cannot be captured with confidence.

E. Classical Continuum/Electronic GC-DFT for electrocatalytic systems
The next (still formidable) challenge is to find the free energy functionals for the quantum subsystem, classical electrolyte, and their interactions. While the quantum part is available from electronic DFT, the electrolyte density can be obtained from the classical theories such as classical DFT 7 . Currently, the most sophisticated approaches combine electronic DFT with RISM 5,11,12 . The combination of DFT and RISM yield a highly complicated set of coupled equations to be solved selfconsistently. Yet, even in the DFT-RISM approach the solvent and ions are treated using "equilibrium classical statistical mechanics applied to thermodynamics of surface phenomena. Quantum mechanical approaches to surface phenomena are excluded ... Therefore, all aspects related to chemisorption and catalytic processes at surfaces are not considered" as discussed on page 883 of Ref. 46.
To obtain a computationally feasible method, we simplify the situation and for a moment consider the electrode and electrolyte interacting only via a static field created by the other component. The setup is similar (but simpler) to frozen-density embedding where one system is kept frozen while the other one is relaxed in the field of the frozen subsystem. In the simplest case, the electrolyte will feel electrostatic potential set by the electrons and nuclei in the electrode. Within the classical GC-DFT framework the free energy of the electrolyte is written as 47 where V ext is the external potential due to ions,solvent as well as a genuine external potential (in our case from the quantum subsystem) and is the thermal de Broglie wave length. The first term accounts for the Coulomb interactions of the pure component, the second is the ideal gas entropy while F ex accounts for excess free energy due to additional interparticle interactions. Also note that F ex is universal functional of the classical densities.
We adopt the classical DFT and further assume that the solvent can be treated as a (non-homogeneous) dielectric ✏(r). Then following Sluckin 9 , the free energy of the electrolyte can be written as where = V ext (r) + V ions is the e↵ective electrostatic due to the external potential V ext (r) induced by the charge density of the electrode, and charge density of the ionic distribution V ions . Presented in this form, it would be possible to sample the ionic distribution in a dielectric solvent in the presence of a external potential formed by the electrode. We will, however, look for the minimum of the free energy functional. Then, without specifying the model for ion distribution, di↵erentiating the grand energy against the ionic charge density results from which the ion distribution is obtained as a function of the ion chemical potential µ i . As the free energy, also the chemical potential potential can be split to ideal and excess contributions such that Varying the grand potential against the electrostatic potential gives the Poisson equation for the ions The standard model for the ion-density is the Poisson-Boltzmann model to be discussed in Section III A. To go beyond the mean-field Poisson-Boltzmann treatment of the electrolyte, one must include ion specific interaction terms. For example, a relatively simple way to incorporate finite ion-size e↵ects is the Bikerman-Freise equation 48 . A computationally feasible way to to include nonelectrostatic interactions between ions raising from hydration interactions, is to use an (asymmetric) Yukawa potential leading coupled Poisson-Helmholtz-Boltzmann equations. 49 . Even simpler way to account for the finitesize ions and change in the dielectric medium due to hydration is incorporate a concentration dependence in dielectric constant yielding ✏[⇢ ± (r)] in the Bikerman model. 50 Now consider the electrode side where the chemical potential of the electrons is to be fixed. This conforms to the traditional grand canonical electronic density functional of Mermin modified by the presence of ions and the solvent. The electronic system will feel the electrostatic potential from the ion density and from the nuclei. There will also be a universal chemical interaction G e [n] functional akin to the one described above for the ions. For the electrons the grand free energy is 14 where ⌧ is used to define temperature to avoid confusion with the kinetic energy T . Note that we have invoked the BO and adiabatic approximations and terms coupling the nuclear motion to the electronic density have been removed. Minimizing the electronic grand energy against the wave function yields the modified Kohn- and the resulting electron density is given by Also for the electron system we consider the electrostatic potential e = V H + V ext + which can be solved from the Poisson equation When the two components are optimized simultaneously as required by the general grand canonical free energy functional, it becomes clear that the electrostatic potential needs to be solved self-consistently for both subsystems simultaneously and therefore is solved for the entire system from the Poisson equation Since we are not allowing the number of nuclei or electrons from the solvent, nuclei from the electrode or reactant to fluctuate, electroneutrality sets a constraint that the electron density from the quantum system must always be compensated by the ionic charge density: . This aspect is extensively discussed in Section II H below. Now we can write to total grand canonical free energy functional for an electronic system embedded in a continuum electrolyte: The GC-DFT approach introduced in the previous section puts implicit solvent models on a firm theoretical basis. The generalized Poisson equation and the derived potential are the working equations in the common implicit solvation models, such as the one implemented in GPAW 51 . However, electrochemical systems add yet another layer of complexity since the solvent contains mobile ions which from a double-layer at the electrochemical interface.
To include the electrolyte and the double-layer, an ionic density is added to the generalized Poisson equation 22,27 . To obtain a self-consistent scheme for the inclusion of the solvation e↵ects, Eq. 35 needs to be minimized with respect to and then with respect to the solute electron density n(r) as well as the ionic density. Before carrying out the minimization, we need to decide which interactions to include in the model. As we wish to test the performance of linear dielectric models (no explicit dependence on e.g. the electrostatic potential), in this work we make the following choices: 1. It has been shown in Refs. 50 and 52 that dielectric decrement due to ion solvation leads to saltspecific interactions and improves the description of the double-layer in classical simulations. The simplest way to include the include solvated ion effects is a linear dependence between the dielectric constant and the ion density: where ✏ 0 is the dielectric constant of the pure solvent and ↵ ± is an ion specific parameter which are available for a range of ions. The model also changes the dielectric constant of the bulk solution as observed for several simple electrolytes (see Ref. 53 and references therein). The ionic concentrations c i and charges q are related to ionic charge densities via 2. Mean field chemical interactions between the electrolyte and quantum system need to be included. These interactions between the solvent and the quantum subsystem include cavitation, dispersion and (Pauli) repulsion in the following form 54 where rS is the solvent accessible area of the quantum system and is an e↵ective surface tension term to be fitted to experiments.
3. Finite size of the ions is included to account for steric e↵ects and avoid unphysically large ion concentrations near the electrode surface which can also be seen as formation of the Stern layer due to specific adsorption. These are entropic e↵ects due to charge organization and can be accounted for using the Bikerman-Freise approach 48 (also often cited as Borukhov's method 55 ): where it's assumed that all ions have the same effective size a 3 .

4.
In recent studies 26,56 electronic DFT has been combined with mPB which includes a separate Stern exclusion layer. This simple approach has been shown to lead to improved description di↵erential capacitance 56 , radial distribution functions, and free energy of solvation in electrolytes 26 . The exclusion width is linked to the ions hydration number 26 and therefore its size a. Within SCMVD, the exclusion layer is simply: c ex i = (g(r) + f (a))c i and has also been implemented within GPAW in this work using experimentally defined ion radii a.
These e↵ects are included either in the dielectric constant or the free energy in Eq. (35). The additional free energy is written as G[n(r), ⇢ ± ] = G chem T S ions . Under these assumptions the minimization of Eq. (35) with respect to the ionic concentration under the boundary conditions c ± (bulk) = c b and (bulk) = 0 (see also for the discussion in Sec. II G related to the boundary conditions) gives the ionic distributions as 50,52 c ± (r) = C ± (r) 1 + 2c b a 3 + a 3 ± (C + (r) + C (r)) (40a) and c b = exp ⇥ 1 µ ± ⇤ /a 3 .Note that the last quantity of equation (40b) is the dielectric decrement which is related to the water dipole rearrangement in the Kralj-Iglic model 57 , in a microfield model Garish and Promislow 53 , and in Booth's 58 dielectric model where the water dipoles orientation is proportional to ⇠ exp[r ]. The ion decrement a↵ects both the bulk dielectric constant as well as the double layer where the (counter) ion concentration is high.
Because additional "chemical terms" are independent of electron density, minimization of Eq. (35) against the electron density leads the KSM equation presented in Eq. (31). The only modification to this equation is that in practice no separate Hartree, external or electrostatic potential terms are needed as the total electrostatic potential is obtained from the Poisson-Boltzmann equation resulting from the minimization of (35) against the total charge density where ⇢ ± [ ](r) are the spatial density of concentration of ions in the dielectric given by Eq. (37). Finally, the free energy and e↵ective potential for the coupled KSMelectrolyte system are obtained after a partial integration of the electrostatic terms |r (r)| 2 μ n (42b) and the electron and ion densities are given by Eq. (32) and (37), respectively. The above set of equations o↵ers a feasible approach to include the electrolyte and the presence of a double-layer in the simulations. From a computational point of view, inclusion of the ionic density can also simplify the treatment of charged periodic systems; in the best case scenario the electrolyte neutralizes the simulation cell and the use of artificial back-ground charge can be avoided. However, the absolute neutralization is achieved automatically for linear PB models. 29 For other mPB models, the charge neutrality needs to be enforced as discussed in detail in Section.II H. Coupled with the possibility of controlling the Fermi-level via the grand canonical electronic density functional mounts to a powerful machinery for modeling electrocatalytic reactions at the solid-electrolyte interface under a constant electrode potential and the ionic double-layer.

G. Electrode potential
Upon deriving the Poisson-Boltzmann distributed ion concentration, a boundary condition for electrostatic potential was chosen as to approach zero in the bulk liquid phase. This provides a very convenient reference electrode corresponding to electrons solvated in the electrolyte. To appreciate this, one can consider the definition of the absolute electrode potential by Trasatti 59 where K is constant depending on the absolute reference choice and E M (red) is the reduced absolute potential where M S is the Galvani i.e. electric potential difference between the electrode and bulk liquid, and µ M n is the chemical potential of electrons i.e. the Fermi-level. One possible choice for the reference K is an electron interacting electrically but not chemically with the environment corresponding to K = µ S n . 60 . While this choice for the reference cannot be realized experimentally, this is exactly produced by the PB-model for a given model electrolyte and choice of PB boundary conditions.
Making the approach even more transparent, we consider the electrode potential under equilibrium conditions using the solvated electron reference which yields whereμ i n is the electrochemical potential of electrons in phase i. From above choice of PB-boundary conditions we know that S ! 0 in the fluid, and that PB-model accounts only for the electrostatic interactions from which µ S n = 0 follows. Therefore the absolute electrode potential within the PB-solvation model is given the two right-most form of Eq. 45 i.e. the Fermi-energy E F .
In practice, the absolute electrode potential needs to be presented on an experimentally relevant reference scale, e.g. against the standard hydrogen electrode (SHE). Whileμ S n = 0 in the PB-model, this choice depends on the actual implicit solvation model used for the electrolyte and, therefore, a connection between the model fluid and an experimental reference electrode needs to established. While several alternatives exist for the absolute potential reference 60 , the most viable reference is an electron at rest in vacuum close to the surface of the solution -this absolute potential of SHE varies between 4.2-4.8 V. One therefore needs to convert E(abs) M from the PB-solvent reference to this vacuum reference. For this we refer again to the absolute potentials by Trasatti 59 and all relevant quantities are summarized in Fig. 2.
The first approach is to compute a correction due to going from the PB fluid to the vacuum without an internal reference. A practical way to connect the PB fluid with a real solvent compute the potential of zero charge (PZC) using PB and compare it to experimental results for the same electrode and the same electrolyte. Equating the absolute PZCs of the electrodes in PB-solvent with the solvent reference and real electrolyte with a vacuum reference, gives the constant o↵-set: P ZC (abs) E P ZC (abs) (46) This o↵-set can be used to convert the PB-scale to SHEscale by .
Another plausible approach which directly gives the E(abs) P B referenced directly against the vacuum reference 59,60 was introduced and applied by Otani 61 .
Here an asymmetric surface with vacuum on the other and solvent on the other side can be used to compute contact (Volta) potentials of both the electrode and solvent as well as the surface potentials. These can be used to convert the PB-solvent reference to a vacuum reference without any experimental input on PZC.

H. Charged electrodes within continuum models
When modelling periodic systems, the simulation cell needs to be charge neutral to avoid the divergence of electrostatic contributions. In partially periodic systems it is also possible to mimic image charges by changing the boundary conditions of the Poisson equation. The ionic distribution obtained by solving the mPB equation provides a physically and chemically sound method to neutralize the cell making mPB models highly attractive. However, in practice, the neutralization from mPBs is incomplete and significant charges can remain in the simulation cell even after the mPB distributed ions are added. 29 The traditional way of neutralizing the cell is to add a homogeneous background charge, which will create spurious interactions and incorrect energetics unless corrections are applied. 17 To avoid the use of homogeneous back-ground charges, two alternatives are presented.

The Lagrange multiplier approach
Modifying the ion distribution provides a chemically justified approach to the overall charge neutrality. The overall neutrality condition is R dr⇢(r) = R n(r) + ⇢ ± N (r)dr = R n df t (r) + ⇢ ± dr = 0. For the models considered in this work, the charge neutrality condition has the form where g(r) is an exclusion (cavity) function determined in Section III.
In general, the ion concentration can also be written in terms of the ion chemical potential which is obtained by minimizing the grand potential with respect to the ion density resulting Eqs. (28). For example, in the GC- The cavity factor g(r) can be interpreted as the e↵ective pair-density functional which depends on the total interactions between the DFT subsystem and the liquid (see again Section III). Keeping this analogy, we can write the total number of ions to fulfill the neutrality constraint as The µ P B ± term accounts for assumptions of the interactions used for deriving the ion distributions in the (modified) PB models. It can also be seen as a Lagrange multiplier to enforce the bulk concentrations as ! 0. The additional term µ 0 ± (r) can be considered as a Lagrange multiplier to enforce charge neutrality 29 or interpreted as excess chemical potential corresponding to the interaction between the ions and the DFT subsystem. It is stressed that this additional chemical potential can be an ion-dependent local function a↵ecting either the e↵ective potential, cavity function or the chemical potential of the PB ions.
The simplest choice for the neutralizing Lagrange multiploer would be µ 0 ± (r) ⇡ µ 0 and the multiplier would be scaling factor µ 0 = ln[ Q df t /Q ± ]. A second option is µ 0 ± (r) ⇡ µ 0 ± = µ 0 ± , the simplest choice being µ 0 + = µ 0 . Then, the neutrality constraint results in which was also derived and applied in Ref. 29 and subsequent work.
The physicochemical implications of the neutralization conditions above can be understood by studying the behaviour of charge neutralized PB-equation deep in the liquid and to meet the conditions set by µ P B ± . The ionic charge density should be zero in order to respect charge neutrality in bulk electrolytes. Also, the electrostatic potential approaches zero by choice of boundary conditions and by choice g equals unity far from the DFT region. Then ion concentrations for position independent excess chemical potential follow If µ 0 ± = µ 0 , the ion density from above approaches zero correctly. However, the ion concentrations do not yield the correct bulk values and instead c i = exp µ 0 c b is obtained. For the asymmetric constraint µ 0 + = µ 0 , the above results in equation will give ⇢ ± = 2 sinh[µ 0 ] and c ± = exp |µ 0 ⌥ | . Thus the bulk electrolyte is not reached even far away from the electrode and situations is not ideal for the electrochemical setup.
Instead, such a situation is encountered when the electrode and nearby electrolyte solution would be isolated from the bulk liquid by a semi-permeable membrane which allows exchange of ions but not the electrode or reactant nuclei. 62,63 This corresponds to an electrochemical equilibrium of ions across the membrane such that the ion and potential gradient is formed; this is the Donnan equilibrium where a potential di↵erence Donnan = (RT /zF ) ln[c in /c out ] builds over the membrane 64  and (x) is arbitrary. As an example, if one chooses µ 0 = 0 and imposes the charge neutrality condition, S = Donnan and the reference electrode discussed in Sec. II G should be S + Donnan = (1) = 0 and Donnan would the boundary condition for the PB equation. If S = 0 is taken as the reference, the bulk ion densities are not obtained and the free energies would be a↵ected. Thus, enforcing position independent Lagrange multiplier for the neutrality condition destroys either the reference electrode in the liquid, the ion density constraint of bulk charge neutrality, or both.
In Ref. 65 it was concluded that the charge-neutrality condition leads to breaking the spatial independence charge distribution. In other words, only a positiondependent multiplier can satisfy the correct ion density and concentration as well as electrostatic potential limits far from the DFT region. Choosing a local, ion independent constraint, µ 0 ± (r) ⇡ µ 0 (r), will e↵ectively result in an ionic cavity function which depends on the charge state of the DFT subsystem. This would likely allow the ions to penetrate deeper in the DFT region. Thus, ions would reside within the solvent cavity resembling specific adsorption which is outside the applicability of PB-like approaches as highlighted while deriving the continuum GC-DFT approach. Also large ion concentrations would be observed close to the nuclei. These considerations apply also to ionic exclusion models presented in Refs. 56 and 26. Still, such ion-independent spatial correction approaches were suggested in Ref. 66 to account for specific adsorption and in Ref. 65 to derive a charge-conserving Poisson-Boltzmann equation by including a "charge-trap region" near the electrode surface. Bound ions can also be considered as an e↵ective way to renormalize the surface charge by ion adsorption to enforce global charge neutrality. 67 Phenomenologically, the above methods where the neutralizing charge is confined near the electrode can be seen as a Stern layer -type correction with tightly bound ions. However, the underlying mPB models are derived for classical nuclei as discussed in Sec.II E, and no real bonds can form between classical and quantum parts of the system. Thus, the strong chemisorption is not within the assumptions included in mPB models and the applicability of this neutralization approach is questionable. Also for a heterogeneous system, the neutralizing charge should not be homogeneous or lie in the low dielectric region but rather reside in the regions of high dielectric to correctly describe the gradient of the electrostatic potential correctly. 68 In the above methods the ion density mimicking strongly bound ions in the Stern layer would result in ions inside the cavity (g(r) < 1, ✏(r) < ✏(bulk)) where dielectric constant is low. Thus, the artificial Stern layer might have an adverse e↵ect on the description of the interfacial potential distribution.
Capturing the interfacial potential drop is important for an electrochemical experiment and therefore ion distribution close the interface should remain una↵ected by enforcing the neutrality. Far away from the electrode, the ion concentration is small and the potential profile is rather flat making it safer to place the neutralizing ions in high dielectric region. A simple form for an augmenting potential fulfilling these properties is where is the activity coe cients accounting for finitesize e↵ects, interactions beyond electrostatics etc. and which can be used to arrive at di↵erent modified PB models. Setting = 1 corresponds to an ideal ionic liquid which can be treated with the Gouy-Chapmann model. Using these definitions, we obtain To arrive at the final expression, we write ion distribution as When the volume of the ion density is large, the correction term becomes small. The function can be interpreted as either locally increasing the concentration or charge of the ions. Also, the function defining the excess chemical potential leading to neutralizing ion density ✓(r) needs to approach zero far away from the DFT subsystem in order to yield bulk ion density and concentrations. Probably, the simplest choice is to make ✓(r) a smooth step function to equal zero at the boundaries but this choice is rather arbitrary.
The neutralization condition can also be extended to the solvent mimicking solvated ions leading to a uniform neutralizing charge outside the solute cavity residing only in the high dielectric region 68 . Such an assumption seems nonphysical as the neutralizing "ions" would not feel the electrostatic potential of the solute. Surprisingly, such an approach has been successfully applied to the computation of potential of mean force profiles which are directly linked to the interaction of the solvent and ions with the surface. Comparison of classical molecular dynamics with homogeneously distributed charge in the high dielectric regions mimicking charging of water molecules, has been shown to yield potential of mean force profiles comparable to explicit ion simulations 68 .
Given that such a simple homogeneous neutralization successfully mimics explicit ions, a simple way of simulating charged cell is obtained by homogeneous neutralization of the solvent region. Indeed, in very recent independent work by Pettersson 30 such an approach, called solvated-jellium model (SJM), was developed, utilized, and implemented in the GPAW code used also in this work. SJM provides neutralization and also establishes a convenient reference electrode as detailed in the present work but the ad hoc nature of the jellium counter charge requires testing and comparison to other methods.

Modified boundary condition approach
An alternative way to study a charged partially periodic cells is to modify the the underlying electrostatic problem by selecting convenient boundary conditions for the PB equation. A common computational setup in model systems used for studying the double-layer behaviour is a parallel-plate capacitor where both plates are placed at z = ±L carrying opposite but otherwise equal charge densities. Also, in the absence of ion-specific effects, the double-layer densities are identical but carry opposite charges. Thus, the overall system is naturally charge neutral.
In a practical DFT calculation the presence of two electrodes is not feasible due to the increased computational cost nor desirable as one is usually only interested in halfcell reactions. A convenient work-around is to treat the total system using the method of images and to use the image as a neutralizing charge. By imposing a condition that S (z = 0) = 0 between the plates, a reference electrode can be restored. For simple model systems where analytic solutions are available the electrostatic potential from the image charge method (in the relevant regions) is equal to placing a perfect conductor at z = 0 as an imaginary electrode. 69 The system polarizes the imaginary electrode and the charge polarization restores charge neutrality of the simulation cell. Also the electrostatic potential at the imaginary electrode is S = 0.
The imaginary electrode method is equal to the e↵ective screening method (ESM) pioneered by Otani and Sugino 19 to remove periodic boundary conditions in one direction in otherwise periodic plane-wave calculations. The essential outcome from the ESM is a derivation of analytic boundary conditions for the Poisson equation in charged systems. If an asymmetric simulation cell is used, as usually done in the presence of adsorbates or reactants, the imaginary electrode needs to placed on the "active side" of the real electrode to obtain the correct zero electrostatic potential reference on the correct side. To internally reference against the vacuum potential, the gradient of the electrostatic potential is taken to approach zero on the inactive side (see discussion at the end of Sec. II G). Also care needs to be taken that the electron/ion density is vanishingly small at the ESM boundary and that enough vacuum is inserted at the cell boundary.
In real-space approaches, such as the grid-base in GPAW 70,71 utilized in the this work, the unit cell can be truly non-periodic and zero electrostatic boundary conditions in non-periodic directions are commonly employed when solving the Poisson equation. In asymmetric systems the dipole correction 72 is applied to give asymmetric electrostatic potentials leading the di↵erent work functions on di↵erent sides of the slab. In the original dipole correction scheme, a correcting potential is introduced where m = R 1 1 dz 0 ⇢ av (z 0 )(z 0 z 0 ) is the surface dipole density referenced to z 0 which in the present case is chosen to be on the active side of the system to satisfy either (z 0 = 0) = 0 or (z 0 = z m ) = 0 , respectively. Though undocumented, GPAW uses a modification of the original dipole method and a potential V dip (r) ⇠ m erf(r) is used instead of a saw-tooth potential of the original dipole correction. In the presence of mPB ions, the only modification needed is the inclusion of ion densities in the total charge density with ion density restricted to exist only on the active side.
Ideally the both the electrostatic potential and its gradient should be zero on the active side of the electrode to obtain the correct behaviour of the ionic distribution but such Cauchy boundary conditions do not give a unique solution of the PB equation 69 and thus the correct gradient is sacrificed in the boundary condition approach. Note that true non-periodicity is exclusive to real-space methods because the potential doesn't have to be periodic in the desired direction. In plane-wave calculations, z 0 = z m /2 needs to be chosen to make the potential periodic in all directions and extrapolation of the potential is needed as discussed in Ref. 19.

III. IMPLEMENTATION
The theoretical concepts derived above for the classical continuum electrolyte are implemented as part of the GPAW code 70,71 which uses the projector augmented wave method (PAW) 73 for replacing the rapidly oscillating KS wave functions by smooth pseudo wave functions near the nuclear regions. More specifically, the Simple Continuum Model based on Volumetric Data (SCMVD) 51 already implemented in GPAW is augmented with new features for treating charged (periodic) system using modified Poisson-Boltzmann models and the boundary condition method closely related to the ESM 19 .
A. Poisson-Boltzmann model with SCMVD SCMVD belongs to a class of implicit solvation models where a quantum subsystem (often called the solute) is immersed in a dielectric continuum. As in all implicit solvation model, the choice of the cavity, where the solute is placed, is critical to retain a chemically meaningful framework. A common approach is to use a dielectric constant ✏ which is unity inside the cavity and goes to the bulk value of the solvent in the liquid phase. This approach was taken in the original SCMVD article 51 and used here as well. The spatially changing dielectric is where g(r) is the e↵ective pair distribution function for of the solvent in the presence of the solute. Since our focus is to combine classical continuum and electronic GC-DFT for describing electrochemical interfaces, it is appealing to treat the dielectric as a functional of the pair-density function which is widely used in the description of classical fluids 7 .
We modify the original SCMVD dielectric to account for local screening by the ions in solution and to include changes in the water dipole orientation due to the external potential created DFT region. Together these terms can be accounted by the ionic decrement 50,52,57 ✏[⇢ ± (r)g(r)] = 1 + (✏[⇢ ± (r)] 1)g(r) (57) where ✏[⇢ ± ] replaces ✏ 1 to account for the dielectric treatment due to ion accumulation. To make a clear connection between the cavity g and chemical quantities, the form of the e↵ective pair distribution in the original SCMVD work was motivated by the Kirkwood-Bu↵ equation relating the pair distribution to the partial molar volume v 1 M and isothermal compressibility  as While several functional forms can be chosen for the effective g, it is desirable that the function goes smoothly from 1 to ✏ 1 and can be di↵erentiated several times. From a conceptual point, the cavity should reflect the solvent/ion distribution around the molecule (see also discussion in the SM and section IV). From classical DFT it is known that the pair-distribution function depends on pair-wise interactions v(r) and external potential (r) as g(r) ⇠ exp[ (v(r) + (r))] 7 . It has also been noted that a dielectric depending explicitly of the electron density can lead to numerical problems. Furthermore, a density-independent cavity has recently been shown to have comparable or even better accuracy for solvation energies than a density-dependent one 56,74 . Based on these considerations, the e↵ective potential in the original SCMVD implementation was chosen to have the form where u(r) the e↵ective pair-potential. We stress that ✏[g(r)] does not depend on the electron density which circumvents numerical instabilities and makes electrostatic potential particularly simple. This point was also stressed in the recently developed soft-sphere continuum method 74 and its adaptation 56 . SCMVD and the softsphere cavities are very similar 74 .
Once the form of the dielectric and the cavity have been chosen, the final quantity to define within SCMVD are the non-electrostatic terms consisting of cavitation, repulsion, dispersion and changes in rotational/translation free energies of ,i.e., the chemical interactions presented in Eq. (38). These all relate to the geometry of the cavity and depend on the cavity surface area which in SCMVD is defined as S[g] = R dr|rg(r)|. Following this definition 51,54 , the chemical free energy is The next step towards treating an electrolyte in contact with a surface is to add the ions to the solution. We note that the implicit solvation model of an electrolyte is best suited for situations where specific adsorption of the ions does not occur as has been discussed in the theory section. In this case the electrolyte forms a di↵use double layer, also called the Gouy-Chapmann layer (GC). In the previous section we showed how the GC model can be made more general by including the finite size e↵ecs, the Stern layer and ion solvation terms. The e↵ect of both terms is to prevent unphysical charge build-up near the surface and also account for ion specific interactions. Still they are not capable of treating specific adsorption of ions and to account for this the ionic density of the electrolyte is excluded from the cavity. This is e ciently achieved by multiplying the ion distribution with a function that goes to zero within the cavity 51,54,56 . In the absence of a separate Stern exclusion layer, natural choice is the e↵ective potential g(r). The the Stern layer layer is present, it is simple to make the replacement g(r) ! (g(r) + f (a)). With a further specification of whether or not the Stern layer is accounted for, the ion concentration including the exclusion is If ion solvation e↵ects are neglected, ion distributions have the form of the GC For small potentials, the GC model can be linearized to yield GC and its linearized form are the most widely used and simplest approaches to include the electrolyte in the context implicit solvation models. 22,28 Both are based on the assumption that the ions are point-like which leads overestimation of ionic concentration close to the charge surface at high potentials and/or high (> 1M ) concentrations. The simplest way to account for the finite-size of the ions is the Bikerman-Freise equation 48 for a symmetric binary salt Ion specific e↵ects can also be incorporated using the ion decrement model of Eqs. (40). Inserting the decrement term from Eq. (36) and utilizing Eq. (57) leads a Poisson-Boltzmann equation with both ion decrement and finite size e↵ects which is inserted in (40). At this point all the energetic terms entering the free energy of Eq. (42) within PB-SCMVD have been defined as the dielectric does not depend on the electron density. Since the potential is very similar to the original SCMVD, details for handling this potential within the projector augmented wave method as done in GPAW, are detailed in the original article 51 .
The last step in a fully self-consistent treatment is the derivation of forces acting on the nuclei in the presence of the electrolyte. The forces are found from F A = r A ⌦. The kinetic and exchange-correlation energies are the same in standard DFT and, thus, also the forces due to these terms remain unchanged. The electrostatic energy, however, is altered by the dielectric medium and ions in solution. The force due to the electrostatics is 54,75 : Both force terms were derived in the original SCMVD article 51 . The only di↵erences are that the electrostatic potential includes contributions from the ions in solution and that the dielectric can depend on the ion density through Eq. (57). The first integral can be calculated using integration by parts to yield 54,75 : Computation of the e↵ective pair density distribution di↵erential was derived previously 51 and, for completeness, we present it here as well where v denotes the direction of the di↵erential. Force contribution from chemical interactions,G chem , have been detailed previously 51 . However, it serves to note that in the PAW formulation, only quantities localized within the augmentation spheres (r c ) contribute to the force (see e.g. Section 1.7 of Ref. 76). Therefore, the term r A |r (r)| 2 = 0 for r > r c . Also the choice of the cavity ensures that ✏(r) = 1 and ⇢(r) ± = 0 when r < r c . Therefore, the additional force contribution from the ions and the continuum in PAW is the second term on the last line of Eq. (66) and all other terms are analogous to a gas-phase GPAW-DFT calculation.

B. Solving the generalized Poisson-Boltzmann equation
The non-linear character of the mPB equation places pressure on the Poisson solver. In GPAW, all potentials and wave functions are presented on a real-space grid and the Poisson equation is solved using multigrid techniques 70,71 . As GPAW works in real-space, the simulation cell can be chosen to be either truly periodic or non-periodic. In all the calculations where the mPB models are applied are at least partially non-periodic; slab, wire, and molecule calculations are 2D, 1D and 0D periodic, respectively.
To meet the increased complexity arising from the nonlinearity of the PB di↵erential equation, the standard Poisson solver in GPAW was augmented with the recently developed algorithms of Ref. 27. Specifically, we have implemented the self-consistent PB algorithms with Selfconsistent iterative procedure (SC) (Algorithms 1 and 3 of Ref. 27). We used a mixing parameter of 0.6.
Following the discussion in Section II H, charged system can be modelled either by enforcing charge neutrality using a (chemically motivated) Lagrange multiplier or by modifying the boundary conditions of the underlying Poisson solver. When using the Lagrange method, Dirichlet boundary conditions with (boundary) = 0 are chosen in the non-periodic directions. In the modified boundary conditions method, Dirichlet boundary conditions with (boundary) = 0 is chosen for the active side and the other boundary is determined from the dipole correction 72 . The practical implementation of the boundary condition method is detailed in the Supporting material.

IV. COMPUTATIONAL DETAILS
All calculations are carried out using the GPAW code 70,71,77 and the PBE 78 functional is utilized. The grid-spacing is set to 0.18Å and the geometry of all systems is fully optimized until the force is < 0.05eV /Å. Convergence criterion for the electronic SCF is 10 5 eV /atom. The PB model is considered converged when the maximum di↵erence between to ionic density between successive iterations is below 10 12Å 3 . 4⇥4⇥1 k-point sampling was used in all calculations. The zdirection of the simulation cell is 3 1 D , where  1 D is the Debye screening length. For the ion decrement and finite ion size mPB models were taken from Ref. 50: the ion diameter for Na (F) is 7.16 (7.04)Å while the linear decrement coe cient is 8. (5.)dm 3 /mol. When including the Stern layer, the reported values used the ionic exclusion radius smaller than the ion diameter: 3.5Å. Smaller and larger exclusion radii were also tested but the results were similar and therefore not shown. The chosen Stern radius used for NaF is close to that of NaCl as determined from AIMD in Ref. 26. For all, but the modified boundary condition method, a symmetric simulation cell with dielectric/electrolyte on both sides was used. This is accounted for in the calculation of the capacitance.
The cavity shape parameters to form g(r) and the effective surface tension are taken from Ref. 51. The parameters are fitted for main group elements and might not be suitable for the metallic surfaces considered in this work. The computational focus of this work is not on quantitative description of energetics but on qualitative performance of the mPB models. The main parameter a↵ecting the electrostatic potential and double layer capacitance considered in this work is the van der Waals radius of the elements used for defining the size of the cavity. In the original SCMVD implementation and the recent SJM model the radii by Bondi 79 were used, and the Au radius is 1.66Å. However, the Bondi van der Waals radii are determined for non-metallic compounds or metal-organic molecules and might therefore be unsuitable for the treatment of the metal-water interface. Indeed, crystallographic data on metallic Au systems indicates that the van der Waals radius is larger and around 1.9Å. 80 During this study we found that the double-layer capacitance sensitively depends on the used van der Waals radius. We therefore determined the distance between the metal-water distance from a planar distribution function (PDF) obtained from an ab initio molecular dynamics simulation as detailed in the Supplementary material. From the PDF we can conclude that the Au-H half-peak is found at distance of 2Å. For comparison, all data was computed using the van der Waals radii 1.66 and 2.00Å. A more complete parametrization for metallic surfaces is currently in process.
While all calculations were performed in the canonical ensemble, a simple grand canonical calculator for optimizing the charge required to reach a set potential in successive fixed charge calculations in the spirit of SJM 30 was also implemented. Also, a constant-µ method by Otani 81 was also implemented in the ASE code 77 . Testing e ciency and performance of di↵erent constant potential methods is the subject of on-going work but here we focus on the description of the electrode-electrolyte interface under various computational schemes in the canonical ensemble.

V. STUDIED SYSTEMS
Using the outlined classical continuum grand canonical scheme for electrons and the electrolyte, we study the di↵erential capacitance (C d ) of well-defined single crystal metal electrodes in aqueous solutions. The di↵eren-tial capacitance is a sensitive but indirect measure of the electrode-electrolyte interaction. As shown in SM Fig. (1) the experimentally extracted di↵erential capacitance for Au(210) in di↵erent concentrations of NaF, the shape and magnitude of the capacitance curve depends sensitively on the used concentration and the surface charge state of the electrode. The di↵erential capacitance is defined as where m is the surface charge density of the metal. A particular point where m = 0 is the potential of zero charge, PZC. Now it most be noted that the capacitance is specific to a given electrolyte and as such also PZC depends on the electrolyte and its concentration. Since the mPB model is suitable for treating the di↵use electrical double layer, an electrolyte without specific adsorption should be chosen. Using experimental values for PZC in di↵erent solvents and electrolytes allows one the connect the PB solvent with a real electrolyte by finding S of the electrolyte.
It is important to bear in mind that in linear and GC mPB models all ions are point-like while in the BF-model the ions have a finite size but neither model e↵ects the dielectric and only the pure solvent dielectric constant is utilized. Therefore it is expected that S is equal for all ions in the GC and BF models. On the other hand, in the dielectric decrement mPB dielectric constant does not equal that of the pure solvent but depends linearly on the ion concentration. Thus, S in the dielectric decrement model is in principle ion and concentration dependent. The e↵ect of ion density and decrement of the dielectric constant to the surface potentials are of the order ⇠ ±10mV for 1:1 electrolytes when the concentration is under 6 mol/dm 3 . 82 Therefore, and the reference electrode are equal for all solvent models considered in this work (however, the value may depend on the used vdW radii defining the cavity).
In addition to the capacitances we also study the interfacial electrostatic potential obtained from di↵erent schemes. Correct behaviour of the electrostatic potentials at di↵erent charges and potentials is paramount for describing the electrochemical cell.
The surface capacitances and the interfacial potential was studied for the Au(210) in 0.5 M NaF. Two di↵erent cavity radii were tested: i) 1.66Å Bondi vdW radius as done in the original SCMVD article 51 and in the recent SJM model 30 and ii) 2.00Å obtained in this work from AIMD simulations. Several di↵erent methods for treating the continuum ion distribution models, namely linearized, Gouy-Chapmann, Bikerman-Fraise, and iondecrement mPB were tested. For all the models di↵erent neutralization schemes were also studied and we used the Donnan Lagrangian method of Eq. (50), uniform neutralization in the dielectric, and modified double layer scheme of equation (54) with a smooth step function in order to reach the bulk ion concentration at the cell boundary. Also simplified ion distributions namely the uniform jellium in the dielectric in the spirit of SJM 30 and Gaussian charge sheets as done in Ref. 83 were tested. The modified boundary condition method without neutralizing ions in the spirit of ESM 19 was studied for comparison.

A. Di↵erential capacitance
Here we focus on the most important features of the computed surface capacitance profiles from di↵erent methods. The detailed surface-charge -electrode potential plots for several tested combinations are presented in the supplementary material.
The main finding is that, regardless of the mPB model, only constant capacitances are obtained from the capacitance calculations. This is in contrast to the experimentally obtained curve which has a distinct bell-like shape peaking at the PZC around 65µF/cm 2 and declining to 25µF/cm 2 at around ±0.5V from PZC as shown in the SM Fig. (1). It is unlikely that the di↵erences are caused by the choice the F-anion as almost identical di↵erential capacitances are obtained for Au(210) in NaBF 4 solutions. 84 The constant value for the simple models such as the Gaussian and jellium 30 as well as linearized 28,56 , point-like GC is expected in a linear dielectric model. However, the quantitatively similar linear behaviour for finite-size Bikerman-Freise model has not been shown before. Also the inclusion of the Stern layer does remedy the linear behavior in the used linear dielectric model. It is surprising that the ion-decrement model also predicts constant di↵erential capacitance when coupled with SCMVD and DFT as this model can produce both maxima and minima of the capacitance in simple classical models 50,52 .A plausible reason for poor performance of the dielectric decrement model is the need to enforce charge neutrality as discussed below. Even if the capacitance is not improved by increasing the complexity and including more ion-specificity, the energetics are expected to be improved in the more complicated mPB models 26 .
While the shape of the di↵erential capacitance does not depend on the chosen method, the absolute values of the capacitance are more sensitive. When using the Bondi ion radius of 1.66Å, the capacitance is ⇠ 65 µF/cm 2 for the Gaussian sheet and linearized, GC, and Bikerman-Fraise PB models. In the jellium and ion decrement models the ion distribution is more di↵use and the differential capacitance is ⇠ 55µF/cm 2 . Capacitance from modified boundary condition is much lower: 30µF/cm 2 . When using a larger Au radius (2.00Å) determined from AIMD simulations, smaller capacitances are observed. As before, smaller capacitances are obtained from jellium and ion decrement models but the spread is smaller and the capacitance is within 34 40µF/cm 2 . The boundary condition method is still an outlier giving the capitance of 22µF/cm 2 . Finally, we note that the capacitances obtained with the larger vdW radius are typical range of ⇠ 20 30µF/cm 2 are observed for other mPB models 22,28,56 and explicit ice-like -water structures. 85,86 .
From the study of di↵erent neutralization schemes, one can conclude that the obtained capacitances depend mildly on how the system is neutralized. The modified double-layer and Donnan lagrange neutralization methods give nearly identical results while the uniform dielectric neutralization produces slightly smaller capacitances. There is an important exception, though: in the ion decrement model the neutralization method changes the overall ion distribution and therefore a↵ects the dielectric in an unexpected way: even for the modified double-layer and Donnan method, the neutralization introduces a smooth ion density in the simulation cell and changes in the dielectric due to decrement results to an almost constant decrease instead of pronounced spatial spatial dependence. An example is shown in Figure 3.

B. Electrostatic potentials from di↵erent schemes
The interfacial potential is an important factor in electrochemical systems a↵ecting both reaction thermodynamics and kinetics. To relate the electrostatic potential and electrochemical potential of electrons from the computational models to experimentally utilized reference electrodes, schemes from Section.II G need to be used. Instead of referencing the potential against reference electrodes, computational studies 56,87 of electrochemical systems have used PZC as an internal reference which is afterwards 22,28,30,88 calibrated to the the SHE reference via the experimentel PZC on the SHE scale. The calibration value depends on the used solvent model: for models where dielectric cavity depends on the electron density as ⇠ erfc[log(n(r)/n c )] such as CANDLE 88 , VaspSOL 28 and linear PCM 22 , absolute SHE is calibrated be between 4.4 and 4.7V but for the soft-sphere or SCMVD type of models such a calibration is not available.
The experimental value 89 for PZC of Au(210) in a 0.5 M NaF is 0.144V vs. SHE. In our modified SCMVD model, the absolute potential is simply E F as discussed in Section.II G. By studying the two di↵erent radii for the dielectric/cavity (at which the dielectric and ionic density di↵erence also start di↵ering from zero) we notice that the Fermi-level of the neutral Au(210) depends sensitively on the used radius; for R = 1.66Å the E F ( m = 0) = 4.27V and for R = 2.00Å the E F ( m = 0) = 4.71V . Thus, calibration for R = 1.66Å is 4.41V and for R = 2.00Å 4.84V . Both are close to the values from electron density dependent cavities and experimentally determined absolute SHE value which varies from 4.4 to 4.85V 90 . As discussed, the calibration is closely linked to the used radius which in turn determines the computed di↵erential double layer capacitance. Therefore a judicious choice would to choose the vdW radius for the dielectric/cavity which reproduces the experimentally obtained capacitance and afterwards fix the calibration between the Fermi-level and PZC from comparison experiment and computations.
Given the dependence of the absolute potential on the used vdW radius, we report the electrostatic potential profiles in Fig. 4 with respect to the PZC of a given model. First, the electrostatic potential obtained from the boundary conditions method decreases to zero slower compared to the other methods and the gradient is nonzero, as expected. Similar behaviour was observed also with the similar ESM 19 method. Of the used mPB models, linearized PB, GC, and Bikerman-Fraise give almost identical potential curves and they almost overlap. When the Stern layer is added (shown for Bikerman-Fraise), screening is delayed but otherwise similar to the underlaying mPB. The electrostatic potential from Gaussian charge sheet depends sensitively on the position and width of the sheet. Here, it is seen to screen the charge faster than the mPB methods, whereas the uniform jellium screens the charge slower. As discussed above, the charge distribution from the ion decrement is rather homogeneous and also the resulting potential is also rather similar to the one obtained from the jellium models or when the Stern layer is included. However, the ion decrement model is the only model for which the interfacial electrostatic potentials are not symmetric for positive and negative potentials. At positive potential the ion decrement model predicts slightly stronger screening which reflects the higher dielectric, i.e. smaller dielectric decrement, of the fluoride anions building at the surface at positive potentials.

VII. CONCLUSIONS
We have presented a general density functional theory framework for modeling electrochemical interfaces at fixed electrode and ion potentials. The developed FIG. 4. The interfacial electrostatic potentials from di↵erent schemes. Above the PZC the potential is set to 0.5V and below to 0.5V vs. PZC. The inset highlights the long-range behavior. The top most atom is located at 16.5Å. The Donnan charge neutralization was used and the vdW radius is 2.00Å.
approach provides a rigorous DFT approach in the ion/electron grand canonical ensemble relevant to electrochemical systems. This allows one to perform atomistic simulations of electrochemical interfaces using the experimentally meaningful control variables, namely the electrode and ion chemical potentials, solvent and temperature. Such simulations will advance the understanding and optimization of electrochemical interfaces.
Besides the general framework, a systematic coarsegraining from a fully quantum mechanical nuclei and electrons, to classical nuclei of the electrode, ionic and solvent nuclei interacting with electrons and to continuum electrolyte is presented. The detailed derivation of the coarse-graining on di↵erent levels shows transparently the approximations involved. This allows straightforward assessment of the capability, validity, and applicability of various schemes for modelling the electrochemical interfaces. The performed derivation enables also the removal of the coarse-graining if needed; this freedom can be utilized to study, e.g., nuclear tunnelling and specific adsorption, at a fixed electrode potential. Another specific example of the general applicability of the theory is an extension to treat adiabatic and nonadiabatic electron transfer kinetics as a function of the electron and ion chemical potentials currently under development by us.
On the simplest level of theory, the presented coarsegraining provides a well-defined approach to continuum solvent models and Poisson-Boltzmann models of the ionic double layer. To compare di↵erent continuum models, several mPB models, including the Stern layer, finitesize e↵ects and dielectric decrement, were implemented in the GPAW code and tested. As detailed, other mPB models can be motivated and implemented by using func-tional minimization of the DFT grand free energy functional. We present an extensive discussion and development of simulating charged systems in contact with a continuum dielectric and mPB double-layer is also presented.
The tested mPB models vary from simple homogeneous jellium models 30 and neutralizing Gaussian background charges to common point-like PB and to mPB models with specific ion contributions including the ion size and dielectric decrement e↵ects. All the included mPB approaces are based on the linear dielectric SCMVD 51 dielectric model of the liquid. Calculations for the interfacial electrostatic potential distribution show that the jellium model screens less than other mPB models while the mPB apart from ionic decrement produce almost identical electrostatic potentials. The ionic decrement model screens the charge rather slowly but introduces slight asymmetry between positive and negative potentials. The long-range behavior of the ESM 19 -like boundary condition method di↵ers significantly from the mPB models.
Computed double-layer capacitance profiles for Au(210) in a NaF electrolyte show that the mPB models with a linear dielectric are incapable of capturing the shape of experimentally determined double layer capacitance curve. All the utilized mPB models including linearized, point-like Gouy-Chapmann, finite-ion Bikerman-Freise and the ion decrement model yield almost constant double layer capacitance rather than an experimentally obtained bell-shaped capacitance within the tested potential range in a 0.5 M NaF solution. Furthermore, the absolute capacitance values are very sensitive to the details on the underlying solvent cavity.
Thus, the overall outcome is that a linear dielectric model even when combined with more complex mPB electrolyte description cannot qualitatively reproduce the experimentally measured capacitance; at this level already jellium-like electrolytes result in almost identical double-layer capacitances and interfacial electrostatic potential similar to the more complicated mPB models. On the other hand, non-linear continuum solvents 29 have shown promising results for describing the double-layer capacitance 56 .
Our results show that even addition of ion-specific interactions into simple continuum electrolyte models is unable to provide a fully satisfactory description of electrolyte-electrode interfaces. Instead, approaches beyond linear dielectric continuum models such as including explicit solvent molecules/ions or the dipolar response of the liquid to external potentials, should be investigated. The presented theoretical framework allows a natural way to include important interactions within GC-DFT and to develop better computational models for electrocatalysis. Finally, while we have presented a rigorous and well-defined GC-DFT coarse-graining scheme and applied it to continuum electrolytes, there are serious caveats that need to be recognized and addressed as the use of linear continuum models for modelling the electrochemical interfaces keeps increasing.