Density-Potential Functional Theory of Electrochemical Double Layers: Calibration on the Ag(111)-KPF6 System and Parametric Analysis

The density-potential functional theory (DPFT) of electrochemical double layer (EDL) is upgraded by adopting (generalized) gradient approximations for kinetic, exchange, and correlation functionals of metal electrons. A new numerical scheme that is more stable and converges faster is proposed to solve the DPFT model. The DPFT model is calibrated with existing differential double-layer capacitance (Cdl) data of the EDL at Ag(111)-KPF6 aqueous interface at five concentrations at room temperature. Metal electronic effects are essential to explain why the two peaks of the camel-shaped Cdl curves are almost symmetric in spite of the size difference of the hydrated cations and anions. A systematic parametric analysis is then conducted in terms of key EDL properties, including the potential of zero charge and the differential capacitance. The parametric analysis, on the one hand, elucidates how quantum mechanical behaviors of metal electrons as well as interactions between metal electrons and the electrolyte solution impact the EDL properties and, on the other hand, identifies key parameters of the DPFT model, which should be calibrated using first-principles calculations and/or advanced experiments in the future.


■ INTRODUCTION
Over the past few years, we have been developing a densitypotential functional theory (DPFT) for electrochemical double layers (EDLs), so named because the grand potential of EDLs is a hybrid functional in terms of the electron density and the electric potential. 1−3 Broadly speaking, DPFT belongs to semiclassical theories of EDLs, along with joint density functional theory (JDFT) 4−9 and DFT-reference interaction site model (DFT-RISM). 10−13 The key feature distinguishing DPFT from other semiclassical models is that the kinetic energy of electrons is expressed as an explicit functional of electron density rather than calculated from wave functions of an auxiliary noninteracting electronic system in the Kohn−Sham DFT. Such an idea can be traced back to the Thomas−Fermi theory of an electron gas in the 1920s, which is widely termed orbital-free DFT (OFDFT) nowadays. 14,15 OFDFT has many developments and applications, notably including the proof of stability of matter using the Thomas− Fermi−von Weizsacker (TFvW) theory of electronic gas, 16 theory of work function and surface energy of metal surfaces, 17−19 computational modeling of materials, 14,15,20 and hydrodynamic theory of quantum plasmonics. 21−24 Inspired by earlier works on metal surfaces, 17−19 the jellium model of EDL developed in the 1980s also employed the OFDFT to describe metal electronic effects. 25−29 As the kinetic energy functional can be, at best, approximated, OFDFT can by no means surpass Kohn−Sham DFT in terms of accuracy in computing material properties, though encouraging progress of improving the OFDFT accuracy has been reported. 30−33 That said, OFDFT has unique merits that make it attractive to the problem of modeling EDLs. First and foremost, combined with a classical DFT of the electrolyte solution, OFDFT allows a grand canonical, namely, constantpotential description of EDLs. Second, the low computational cost of OFDFT allows us to model realistic EDLs at nanoparticles that will be beyond the reach of Kohn−Sham DFT in years to come. In short, OFDFT offers a viable approach to tackle two grand challenges of modeling EDLs, namely, treatment of the electrode potential and simulation of EDLs of realistic scales. EDL, in turn, provides an exciting playground for OFDFT to realize its full potential.
In addition to the electronic structure, equally important to EDL models is the electrolyte solution, which can be treated on varying levels of complexity, including the Debye−Huckel theory, the Poisson−Boltzmann theory, and its variants such as the modified Poisson−Boltzmann theory, and integral equation theories such as the reference interaction site model, see recent review articles. 34−36 In DPFT, the free energy functional of the electrolyte solution was derived rigorously from statistical field theory and is exactly on the mean-field level. 2 Compared to the Poisson−Boltzmann theory that is most frequently used in joint DFTs, the obtained free energy functional further considers ion size effects, orientational polarization of solvent molecules, and other short-range interactions using a reference system.
DPFT is different from earlier jellium models of EDL in several aspects. First, previous jellium models were developed and solved under constant-charge conditions, namely, conditions where the number of metal electrons is preset. 25−29 Therefore, the double-layer capacitance calculated from these jellium models was often plotted as a function of the electrode charge. On the contrary, DPFT is a constant-potential model which is able to simulate EDLs as a function of electrode potential. Second, previous jellium models inherited Smith's treatment of the metal surface. 17 Specifically, the metal electrons were described using the Thomas−Fermi−Dirac−Wagner theory without gradient terms, and the electron density was solved approximately using trail functions. In comparison, DPFT upgrades the description of metal electrons by including gradient terms, and solves for the electron density and the conjugate electric potential self-consistently using a numerical scheme. Note in passing that the trail functions are limited to planar, one-dimensional (1D) EDLs, while the numerical scheme is, in principle, applicable to EDLs of arbitrary structures. Third, as for the electrolyte solution, previous jellium models used either the Poisson−Boltzmann theory or the integral equation theory of liquid, while DPFT introduces a general statistical field theory framework into which many effects and interactions that were neglected before can be incorporated. This paper is the fourth one in the series. The first paper introduced the first version of DPFT that used the Thomas− Fermi−Weizsacker theory for the kinetic energy functional, the Dirac−Wigner (DW) theory for the exchange−correlation energy functional, and a modified Poisson−Boltzmann theory for the electrolyte solution. 1 Constant-potential simulations of metal−solution interfaces were conducted, from which the potential of zero charge (pzc) was obtained and analyzed in terms of electrolyte effects. In the second paper, we introduced the discrete structure of metal cationic cores, considered specific adsorption of electrolyte ions using a model Hamiltonian approach, and compared the model with experimental data measured on Ag(111)-aqueous KPF 6 solution interfaces. 3 In the third paper, we focused on improving the treatment of the electrolyte solution and derived a comprehensive mean-field free energy functional by following the statistical field theory procedure. 2 The purpose of this paper is 4-fold. First, we modify the theory by upgrading the exchange−correlation energy functional from the local density approximation level to the generalized gradient approximation (GGA) level and using Morse potentials to describe metal−solution interactions. Second, we improve the numerical stability and efficiency of solving DPFT by proposing a new numerical scheme. Specifically, we solve for the cubic root of the electron density rather than the electron density itself, which avoids numerical problems encountered in calculating terms of the electron density with fractional exponents. In addition, we implement the model in COMSOL Multiphysics, rather than Matlab previously, because it is more convenient to model EDLs of various structures in COMSOL Multiphysics. Third, we calibrate the model with high-quality experimental data of the EDL at the Ag(111)-KPF 6 aqueous interface at five concentrations. 37 Comparison between this model and classical EDL models reveals the essential roles of metal electronic effects in understanding detailed features of the experimental data. Fourth, using the calibrated model as the baseline, we conduct a systematic parametric analysis of the model, guiding where further improvements of DPFT should be targeted.

■ MODIFICATIONS TO DPFT
We exposit the DPFT framework in the Supporting Information of this paper and mention modifications made in this paper herein. In the quantum part, the exchange−correlation functional is upgraded from the Dirac−Wigner (DW) functional, a local density approximation, to the Perdew−Burke−Ernzerhof (PBE) functional, 38 a generalized gradient approximation, see eqs S7−S10. The kinetic energy is described using the Thomas− Fermi−von Weizsacker (TFvW) functional as in the previous work. 2,3 In Figure 1, we compare two sets of exchange− correlation functionals, including the PBE functional and the DW functional used in previous work, 3   where e au = 27.2 eV is the energy in atomic units, t TF is the kinetic energy functional, see eq S6, u X 0 is the exchange energy functional of a homogeneous electron gas, eq S8, and u C 0 is the correlation energy functional of a homogeneous electron gas, eq S9. n ̅ e is the dimensionless electron density. In this paper, an overbar denotes a dimensionless number density normalized to the reference number density (a 0 ) −3 , where a 0 is the reference length, the Bohr radius. It is shown that the chemical energy of electrons μ e increases nearly linearly after a brief decrease as the metal cationic charge density n ̅ cc 0 increases. The difference between these two functionals is within 0.3 eV in the examined range of n ̅ cc 0 . 17 The difference transitions from negative to positive as n ̅ cc 0 increases. Following Shandilya, Schwarz, and Sundararaman, 39 we use Morse potentials, expressed in eq S12, instead of truncated power potentials in our previous work, 2,3 to describe metal− solution interactions. A Morse potential contains three parameters, namely, a well depth D l , a coefficient controlling the well width β l , and an equilibrium bond distance B l , where the subscript l indexes different solution components (s for solvent, c for cation, a for anion). We are about to show in the following section how these parameters can be determined using firstprinciples calculations.
Previously, we solved differential equations in terms of n ̅ e . 2,3 It was found that careful treatments are needed for a stable and convergent numerical solution. Specifically, zero or negative n ̅ e must be avoided as n ̅ e occurs in the denominator, and fractional exponents of n ̅ e are required in many places. A positive-defined n ̅ e , n ̅ e mod = |n ̅ e | + n ̅ e lb with n ̅ e lb being a positive lower bound was used for this end. In this work, we propose a new numerical scheme in which we solve for ψ = (n ̅ e ) 1/3 rather than n ̅ e . This way, we can get rid of issues related to fractional exponents of n ̅ e . The new controlling equation in terms of ψ is given in eq S39.
The modified model is solved in COMSOL Multiphysics. The two controlling equations, eqs S39 and S40, are implemented as a coefficient form partial differential equation. A step-by-step tutorial of constructing the model from scratch is provided in the Supporting Information of this article.

■ MODEL PARAMETERIZATION AND CALIBRATION
A basic set of model parameters are obtained from calibrating the model using experimental data measured on the Ag(111)-KPF 6 aqueous interface by Valette. 37 This interface is selected for model calibration for two reasons. First, PF 6 − is found to be a very weakly specifically adsorbing anion on Ag(111), rendering that this interface can well mimic an ideally polarizable interface in a wide potential range that is treated in this model. Specific adsorption and chemisorption of ions were treated using a Newns−Anderson Hamiltonian with several additional parameters in our previous work. 3 These Hamiltonian parameters can be determined from Kohn−Sham DFT calculations. Calibration of the DPFT model with chemisorption will be reported in a future part of this series. Second, Valette reported C dl curves at a series of electrolyte concentrations, allowing us to calibrate the model using a set of C dl curves rather than a single one.
The model describes a one-dimensional EDL. The metal cationic cores are represented by a uniform background with a positive charge density n ̅ cc 0 , namely, representing the number of electrons of a silver atom, and a Ag = 4.08 Å is the length of the cubic closed-packed cell of Ag, which contains four silver atoms, and a 0 = 0.529 Å is the Bohr radius. n ̅ cc 0 is normalized to the reference number density (a 0 ) −3 . x ̅ = x/a 0 is the dimensionless coordinate. This model considers all electrons and does not require pseudopotentials for the metal cationic cores anymore. In an all-electrode model, we should use the vacuum permittivity ϵ 0 in the Poisson equation for the metal phase.
Within a 1D approximation, we do not consider the discreteness of metal cationic cores, which can be best treated in a three-dimensional (3D) model. We designate the coordinate origin, x ̅ = 0, in the metal bulk phase, and x ̅ M as the dimensionless thickness of the metal region. x ̅ M should be large enough, say 20 (∼1 nm), such that the metal bulk is reached at x ̅ = 0. In the metal bulk, the gradients of electron density and electric potential, denoted ϕ, are zero, Here, an overbar denotes variables and operators in the  because metal electrons cannot go beyond several Å from the metal surface, and the electric potential is designated as the reference = 0 (5) In this model, the electric potential of the metal, ϕ M , is not used explicitly in boundary conditions. Instead, constant-potential simulations of the present model are realized by actually varying μẽ, which is tantamount to varying ϕ M up to a constant of e 0 , Note in passing that it is more fundamental to control μẽ than ϕ M because a voltmeter measures the difference in μẽ, not ϕ M between two electrodes. 40 We will use −μẽ, which is equivalent to ϕ M up to some constants, to represent electrode potential in the figures of this article. In the subsequent discussion, we will use electrode potential and −μ̃e interchangeably. After specifying the setup, boundary conditions, and operating conditions of the model, we now explain the acquisition of model parameters. The metal parameters, n ̅ cc 0 and ϵ op,M , are already determined for the EDL, and there is no freedom to tune them. We use the PBE functional with recommended gradient coefficients θ X = 0.1235 and θ C = 0.046, and the only free parameter in electronic functionals is the gradient coefficient in the kinetic energy, θ T , which is to be determined from calibration with experimental data.
Parameters of the electrolyte solution are adopted from the literature. Molecular dynamics (MD) simulations have determined the radii of hydrated K + and PF 6 − to be within 4− 5 41 and 3−4 Å, 42 respectively. We use the upper limits, 5 and 4 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Å, herein. The length of the cubic lattice accommodating a water molecule is determined to be 3.11 Å, ensuring a bulk concentration of c s b = 55M. The relative permittivity of the bulk solution is ϵ b = 78.5. Therefore, knowing the optical dielectric constant ϵ op,S of the electrolyte solution, we can calculate the effective dipole moment of water molecules to be where k B is the Boltzmann constant, T is the temperature, and N A is the Avogadro constant. The only free parameter of the electrolyte solution is thus ϵ op,S , which is to be determined from calibration with experimental data. The continuous transition of the dielectric constant from ϵ op,M = 1 in the metal phase to ϵ op,S in the solution phase is empirically described as with a coefficient β op controlling the width of the transition region, and erf (x ̅ ) is the error function. Parameters of the Morse potentials describing metal− solution interactions can be acquired from Kohn−Sham DFT calculations. Water adsorption on Ag(111) has been studied by Le et al. using ab initio molecular dynamics (AIMD). 43 The binding energy of water is 0.25 eV with a bond length of 2.82 Å. So, the Morse potential of water reads w s ( The coefficient β s can be determined by fitting the whole potential profile, which is not available in the AIMD work of Le et al., and we use β s = 1 as an initial guess. Since hydrated K + (subscript c) and PF 6 − (subscript a) are very weakly adsorbing ions, D a and D c should be much smaller than D s , while B a and B c are larger than B s .
To summarize, unknown model parameters in the current model are θ T , ϵ op,S , β l (l = s,c,a,op), D l (l = c,a), and B l (l = c,a). All of these unknown model parameters except θ T can be, in principle, determined from Kohn−Sham DFT calculations. Herein, these unknown model parameters are estimated by comparing model-based and experimental C dl curves at five concentrations. The impact of these parameters on the model results will be gauged in a parametric analysis in the next section.
The model calculates C dl by differentiating the surface free charge σ free with respect to electrode potential where the second equal sign transforms the electrode potential to the electrochemical potential of electrons as expressed in eq 6, the third equal sign is due to the definition of 2 results from the balance of dimensionality, and the fourth equal sign corresponds to another way of calculating σ free . These two ways of calculating σ free are equivalent because the total system is electroneutral, ∫ dx ̅ (n ̅ e + n ̅ a − n ̅ c − n ̅ cc ) = 0.
At low concentrations, the C dl curve has a Gouy−Chapman minimum that is located at the potential of zero charge (pzc). The specific value of μẽ corresponding to the pzc, denoted μẽ ,pzc , is decomposed into two parts, where ϕ M,pzc is the inner potential of the bulk metal relative to the bulk solution under this condition, and μ e is the chemical potential of electrons in the bulk metal, as expressed in eq 1. Figure 2 shows a comparison between model-based and experimental results of C dl of Ag(111)-KPF 6 aqueous interface at five concentrations (100, 40, 20, 10, 5 mM) at 25 ± 2°C. 37 The model-based results at five concentrations are calculated using the same set of tunable model parameters θ T = 2.08, ϵ op,S = 4, β l (l = s,c,a,op) = 1, D l (l = c,a) = D s /6, and B l (l = c,a) = 7.56. The value of θ T is greater than 5/3 given by the original von Weiszacker theory, which is not surprising as the latter theory is exact only for a single-electron system. 44 ϵ op,S falls into the usual range between 3 and 6 for the dielectric constant of the inner layer of the EDL. 45 When converted to Å −1 , β l corresponds to 0.53 Å −1 , which is close to values between 0.4 and 1.1 Å −1 for a similar system. 39 B l represents a reasonable bond length of 4 Å. The model can well reproduce the well-known camel shape of the C dl curve and the lifting trend of the Gouy−Chapman minimum with increasing concentration. Model-based and experimental data are plotted in the same graph in Figure S1, where noticeable differences between them are observed. Experimental values are overall larger than model-based values, which could be attributed to an underestimated roughness factor. The two peaks that signify overcrowding of counterions when the EDL is highly charged are slightly broadened in the model. In addition, model-based values are much lower than experimental values when −μ̃e < 3.8 eV. An improved agreement between the model and experiments can be expected by global optimization of model parameters, which are tuned by hand in this work.
Interestingly enough, the two peaks almost have the same height, though hydrated K + (the left peak) is bigger than PF 6 − (the right peak). This implies that the left peak will be higher than the right one if cations and anions are of the same size, as will be shown later. This asymmetry is ascribed to metal electronic effects in this model. By contrast, without considering metal electronic effects, the classical GCS model with symmetric size would lead to two symmetrical peaks, e.g., see Figure 7 in ref  Figure 3 displays the model results for the Ag(111)-0.1M KPF 6 aqueous interface at five electrode potentials spanning from negative to positive σ free . In these plots, the metal edge is relocated at x = 0. Figure 3a shows the overall distribution of the electron density normalized to its bulk value, and (b) gives a close-up of the spillover region just outside the interface. The electron tail stretches out more at lower electrode potentials, viz., less negative μ̃e when electrons have a higher potential energy. An overshoot in the electron density is observed near x = 0. Figure 3c exhibits the distribution of electric potential ϕ, and a close-up just outside the metal is shown in Figure 3d. A change in the sign of ϕ signifies a change in the sign of σ free and a change in the identity of counterions when varying μẽ. At more negative μẽ, namely, more positive electrode potential, ϕ is more positive, attracting anions into the near-surface region and repelling cations out of the region, as shown in Figure 3e,f for the normalized density of anions and cations, respectively.
All of the above phenomena are already known from classical GCS models 46,48−50 and in our previous works. 2,3 The use of a Morse potential that is deeper and closer to the metal surface for water results in a dense packing of water molecules on the metal surface, Figure 3g. Hydrated counterions reside outside of this water layer; see a schematic diagram in Figure 3i. The steric effects of counterions diminish water density beyond the first water layer. The total permittivity is shown in Figure 3h, which follows the distribution of water density. We see a transition from ϵ op,M = 1 in the metal phase to ϵ op,S = 4 around x = 2 Å. Even though the densities of water and ions are zero around x = 2 Å, ϵ op,S is larger than 1 to account for the screening effects of the electronic cloud of the first water layer.

■ PARAMETRIC ANALYSIS
We now use the calibrated model to explore the influence of selected model parameters. The parametric analysis will be conducted in three groups. The first group contains parameters related to the electronic structure, including gradient coefficients in electronic functionals, θ T and θ XC , and the metal cationic charge density n ̅ cc 0 . The second group focuses on the influence of solvent, including the distribution of optical permittivity near the metal surface, described by ϵ op (x ̅ ) in eq 7, the dipole moment of solvent, and the Morse potential of metal−solvent interactions. The last group consists of parameters related to ions, including ion size, ion concentration, and parameters in Morse potentials of metal−ion interactions.
In addition to C dl defined in eq 8, the moment of the electron density distribution will be used in the subsequent analysis, defined as )xd ed e cc (10) where n ̅ cc is given by eq 2. A more positive M ed represents a more extended distribution of electron density, equivalently, a larger potential drop across the metal surface if the permittivity distribution is fixed. Electronic Structure Parameters. Figure 4a,c shows the C dl curve as a function of −μẽ at three values of θ T and θ XC , respectively. Other parameters except the one under evaluation have their base values in this one-factor-at-a-time study. According to eq 6, varying −μẽ is equivalent to varying the inner potential of the electrode ϕ M . We do not analyze θ X and θ C independently because the gradient terms in the exchange and correlation functionals are combined into one term with one composite gradient coefficient θ XC = θ X − π 2 θ C /3.
The C dl curves are shifted along the −μ̃e axis with varying θ T and θ XC . Specifically, the C dl curve is shifted to the right side at more positive θ T and more negative θ XC . Accordingly, the pzc is more positive at more positive θ T and more negative θ XC . As the chemical potential of electrons μ e is independent of θ T and θ XC , as seen from eq 1, variation in the pzc is solely due to variation in ϕ M,pzc , according to eq 9, which is determined by the distribution of the electron density at the metal surface. Figure 4b,d shows the dimensionless moment of the electron density distribution, M ed , defined in eq 10. The general trend is that M ed decreases with more positive electrode potential. That is to say, the electron tail is shrunk at more positive electrode potential when the electrochemical potential of electrons is lower. In addition, M ed increases at more positive θ T and more negative θ XC . The opposite influence of θ T and θ XC is the consequence of opposite signs of the kinetic and exchange− correlation energy of a homogeneous electron gas. The electron density is distributed in such a profile that the total energy is lowest. The gradient terms in the total energy grow at more positive θ T and more negative θ XC if the electron density distribution is kept the same. To counteract this increasing trend, the electron density distribution should become more even, namely, M ed increases.
Different metals and different facets of a metal have different values of n ̅ cc 0 . Figure 5 shows that the C dl curve is elevated at larger n ̅ cc 0 because there are more electrons to screen the electric field in the double layer. A nonmonotonic dependence of μẽ ,pzc on n ̅ cc 0 is observed in Figure 5. According to eq 9, μẽ ,pzc is the difference between the chemical part μ e and the electrostatic part e 0 ϕ M,pzc . The chemical potential μ e grows with n ̅ cc 0 in the examined range, as shown in Figure 1. ϕ M,pzc also grows with n ̅ cc 0 . Therefore, the nonmonotonic dependence of μẽ ,pzc results from the balance of two monotonically increasing trends. It is important to understand that we have used the same θ T . At this level of theory, θ T is an empirical parameter that should be calibrated for each metal. In addition, it is about to show that μẽ ,pzc is markedly impacted by ϵ op (x ̅ ), which depends on the solvent properties.
Effect of Solvent Properties. Figure 6 exhibits the influence of solvent parameters, including the solvent dipole moment (p s ) in Figure 6a,b, and the metal−solvent bond length in the Morse potential B s in Figure 6c,d. The C dl curve is elevated up with increasing p s in the unit of Debye (D), as expected. In addition, the pzc shifts from −4.1 eV at p s = 3D to −4.06 eV at p s = 4D and further to −4.02 eV at p s = 5D. At larger p s , Figure 6b shows that the electron density distribution is shrunk when μ̃e is more negative than −3.8 eV but expanded when μ̃e is more positive than −3.8 eV. Overall, M ed drops faster with increasing electrode potential when p s is higher.
At high electrode potentials, ϕ is positive near the metal surface, see Figure 3d. Solvent with larger p s can screen the electric field more efficiently, leading to a less positive ϕ in the electron tail region. This means that chemical potential μ e is higher at a fixed μẽ when p s is larger. According to Figure 1a, μ e decreases with electron density in the electron tail with n e < 0.003. Therefore, when p s is larger, n e should decrease to compensate for the decreased ϕ in the electron tail region at a fixed μẽ.
C dl grows when the metal−solvent bond length in the Morse potential B s is shortened, as shown in Figure 6c. This phenomenon can be understood readily using Helmholtz's model of EDL. The capacitance of a planar capacitor is larger when the gap between two end plates constituting the planar capacitor is narrower. When B s is increased to 7a 0 , namely, when the bond length is similar for solvent and ions, we find a spike of the cation density distribution around 0.35 nm before the metal surface at μẽ = −3.8 eV, Figure 6d. This spike indicates that cations as the counterions can break into the first solvent layer when the local electric potential is negative enough, as schematically shown in Figure 6e. Such a competition between solvent molecules and counterions was also revealed in our previous work. 3 Shatla et al. recently studied the pzc of the Au(111)nonaqueous solution interfaces. 51 Their study shows that the pzc is 0.31 V Ag/Ag + in DMSO with p s = 1.97D, −0.01 V Ag/Ag + in DMSO with p s = 3.96D, and −0.1 V Ag/Ag + in PC with p s = 4.9D. This trend that the pzc shifts to negative values when p s increases is consistent with the model results in Figure 6a. However, the pzc of Ag(111)-ACN aqueous interface with p s = 3.93D is 0.24 V Ag/Ag + , which breaks the above trend.
It is important to note that the model results in Figure 6a are obtained by only changing p s while keeping all other parameters unchanged. In experiments, when one changes the solvent, many parameters change at the same time. Therefore, the discrepancy between the model and experiments shall not be surprising. For example, we show in Figure 7 that the pzc can be markedly influenced by other parameters related to the solvent. We change ϵ op,S in eq 7 and keep all other parameters unchanged. Figure 7a shows that the C dl curve is elevated, and the pzc decreases when ϵ op,S is increased from 3.6 to 4.4. The total permittivity is shown in Figure 7d. When ϵ op,S is lower, the potential drop across the metal surface at the pzc is more positive, as shown in Figure 7b. Since the potential in the solution bulk is taken as the reference, the inner potential in the metal bulk is higher when ϵ op,S is lower, Figure 7b. A higher inner potential of the metal leads to a more negative μ̃e ,pzc , namely, a more positive pzc, since the chemical potential μ e in the metal bulk does not change with ϵ op,S . As μ̃e ,pzc is less negative at larger ϵ op,S , electron distribution is more extended outside the metal, Figure 7c.
Effect of Ion Properties. The ion concentration is one of the easiest parameters tunable in experiments. In Figure 8, we examine how the C dl curve depends on the ion concentration while all other parameters have their base values. Compared to Figure 2, we include a more concentrated electrolyte solution of 1M. The camel shape with the minima at the pzc is transitioned to a volcano shape when the ion concentration increases. The camel−volcano transition has been clearly exposited by Kornyshev,48 and many others; see a recent review. 49 However, these classical EDL models would predict the volcano peak is also at the pzc where the camel minima are located. By contrast, the present model shows a left shift of the volcano peak relative to the pzc, which is a manifestation of metal electronic effects.
The model uses different sizes for cations and anions. Specifically, the radius of hydrated K + and PF 6 − are 5 and 4 Å, respectively. According to classical EDL models, 46,48−50 the size asymmetry will lead to asymmetric camel shapes. For the present case, one would expect that the peak at lower electrode potential that signifies crowding of larger cations would be smaller than that at more positive electrode potential. By contrast, experimental data and model results show a nearly symmetric camel shape. In Figure 8b, we show the model results when the anion radius varies between 3, 4, and 5 Å. The peak at more positive electrode potential grows when the anion radius is smaller. When cations and anions have the same radius of 5 Å, the peak at lower electrode potential is higher, which would not be expected from classical EDL models. This again reflects the essential role of metal electrons. At lower electrode potentials, electron density is extended more outside the metal, see Figure  3b, enhancing the capability of screening the electric field emanated from the metal, namely, increasing the C dl .
The influence of the well depth of Morse potential on C dl is examined in Figure 9a,b for anions and cations, respectively. The anion effect is seen on the right peak at more positive electrode potentials when the metal surface is positively charged, while the cation effect is seen on the left side at more negative electrode potentials when the metal surface is negatively charged. Around the pzc, C dl grows when the well depth is larger for both cations and anions because it effectively enriches the local concentration. This has been recently discussed in detail by Doblhoff− Dier and Koper in the context of the EDL at Pt-aqueous solution interfaces. 52 A shift in the pzc is observed for both cases when the well depth is changed. This reminds us of the fact that the pzc is an interfacial property that is determined by both parties constituting the interface, namely, the metal and the nearby electrolyte solution. A change in the local reaction environment in the EDL would change the pzc, even though the cations and anions are not specifically adsorbed. Note in passing that Valette used the shift in the pzc at different ion concentrations as a descriptor of ion-specific adsorption. 53−55 The model results presented in Figure 9 reveal that it might be problematic to use the concentration-dependent shift of the pzc as a descriptor of specific adsorption.

■ CONCLUSIONS
We have modified the density-potential functional theory (DPFT) of electrochemical interfaces using a generalized gradient approximation for exchange−correlation energy and Morse potentials for metal−solvent/ions interactions. We have calibrated an EDL model based on DPFT using experimental differential double-layer capacitance (C dl ) curves measured on Ag(111)-KPF 6 aqueous solution by Valette. The calibrated model has then been subject to a parametric study in which parameters of metal electrons, solvent, and ions are examined. The difference in chemical potential of electrons between the PBE functional and the Dirac−Wagner functional is within 0.3 eV when the charge density of metal cationic cores (n ̅ cc 0 ) varies from 10 −4 to 10 −1 . The difference in work function and potential of zero charge is thus mainly due to the interfacial potential drop and the electron density distribution across the metal surface. In this regard, the gradient coefficients in kinetic and exchange− correlation functionals, namely, θ T and θ XC , as well as the permittivity near the metal surface, characterized by ϵ op,S , are important. A larger θ T , a more negative θ XC , and a lower ϵ op,S lead to a more extended metal electron tail, resulting in a larger interfacial potential drop and a more positive potential of zero charge (pzc). C dl at the pzc is larger when the charge density of metal cationic cores (n ̅ cc 0 ) is larger. However, the pzc varies nonmonotonically with n ̅ cc 0 . The C dl curve is elevated, while the pzc decreases when the solvent dipole moment (p s ) is greater. Changing the well depth of Morse potential affects different peaks of the C dl curve for cations and anions and also shifts the pzc. The parametric analysis indicates that θ T , the permittivity within 1 nm from the metal surface, and metal−solution interactions are key parameters of the present model, which  should be calibrated carefully with first-principles methods and experiments.
This study also reveals several aspects of the difference between the present model and classical EDL models. In terms of physics, the present model contains quantum mechanical effects of metal electrons, as well as short-range interactions between the metal and solvent/ions; both are not considered in many classical EDL models. This allows this model to calculate the pzc, which is usually taken as a potential reference in classical EDL models. In terms of phenomenology, the present model and classical EDL models differ in describing the dependence of C dl on the ion size and concentration. Both models give a camel−volcano transition of the shape of C dl with increasing the ion concentration. Additionally, the present model shows that the peak of the volcano peak at high ion concentration does not coincide with the minimum of the camel-shaped C dl . Classical EDL models would lead to a camel-shaped C dl curve with two symmetric peaks when cations and anions have the same size. However, the present model gives a camel-shaped C dl curve with a higher left peak for the case of symmetric size. Both discrepancies between this model and EDL models are attributed to metal electronic effects.