Reliable first-principles calculations of electrochemical processes require accurate prediction of the interfacial capacitance, a challenge for current computationally efficient continuum solvation methodologies. We develop a model for the double layer of a metallic electrode that reproduces the features of the experimental capacitance of Ag(100) in a non-adsorbing, aqueous electrolyte, including a broad hump in the capacitance near the potential of zero charge and a dip in the capacitance under conditions of low ionic strength. Using this model, we identify the necessary characteristics of a solvation model suitable for first-principles electrochemistry of metal surfaces in non-adsorbing, aqueous electrolytes: dielectric and ionic nonlinearity, and a dielectric-only region at the interface. The dielectric nonlinearity, caused by the saturation of dipole rotational response in water, creates the capacitance hump, while ionic nonlinearity, caused by the compactness of the diffuse layer, generates the capacitance dip seen at low ionic strength. We show that none of the previously developed solvation models simultaneously meet all these criteria. We design the nonlinear electrochemical soft-sphere solvation model which both captures the capacitance features observed experimentally and serves as a general-purpose continuum solvation model.

The electrochemical double layer plays a central role in the chemistry of a wide variety of processes including electrocatalysis, biocatalysis, and corrosion. Over the last century, the fundamental structure and properties of the electrochemical double layer have been characterized extensively, facilitated by advances in spectroscopy,1,2 microscopy,3 and single-crystal electrochemistry.4 The charge and electric field distributions in the double layer vary with applied potential in complex ways, exhibiting a rich interplay between the surface structure, electrolyte, and adsorbates.5 Measuring, understanding, and quantitatively modeling these distributions is vital because they can strongly affect chemical reactions at the electrochemical interface.

The charge and electric field distributions collectively contribute to the overall change in charge with potential or electrochemical capacitance of the double layer. This capacitance frequently contains signatures of ion adsorption and chemical reactions at the interface. Fundamental studies of the intrinsic response of the double layer therefore employ single-crystalline metal surfaces and a non-adsorbing electrolyte. Figure 1 shows the experimental capacitance curves of Ag(100) in the aqueous KPF6 electrolyte,6 and the archetypal capacitance of a metal surface in the non-adsorbing aqueous electrolyte. These curves exhibit a characteristic wide “hump” shape that is nearly independent of electrolyte concentration, and a narrower “dip” at the potential of zero charge (PZC) that becomes increasingly pronounced at low electrolyte concentrations.7–10 

FIG. 1.

Comparison of capacitance of Ag(100) in KPF6 with experimental measurements6 for (a) linear solvation models including SCCS,26 linearPCM,22 and CANDLE25 at a fixed ionic concentration of 1M (M = mol/l) and (b) nonlinear solvation models including the original parametrization of NonlinearPCM22 and its refit version.30 The linear models show none of the potential or concentration dependence features of the experimental data,30 while refit NonlinearPCM exhibits only a modest capacitance reduction near the PZC at low ionic concentration.

FIG. 1.

Comparison of capacitance of Ag(100) in KPF6 with experimental measurements6 for (a) linear solvation models including SCCS,26 linearPCM,22 and CANDLE25 at a fixed ionic concentration of 1M (M = mol/l) and (b) nonlinear solvation models including the original parametrization of NonlinearPCM22 and its refit version.30 The linear models show none of the potential or concentration dependence features of the experimental data,30 while refit NonlinearPCM exhibits only a modest capacitance reduction near the PZC at low ionic concentration.

Close modal

The concentration-dependent dip in the electrochemical capacitance can be explained by the Gouy-Chapman-Stern (GCS) model,11 one of the simplest models for the interface. This model contains two capacitance contributions in series: a concentration-independent contribution from the dielectric region devoid of ions next to the electrode (Helmholtz capacitance), and the strong concentration and potential-dependent contribution from ions that are further away (diffuse capacitance). The latter ionic contribution is strongly nonlinear with the potential and grows exponentially as the potential deviates from the neutral value, producing the dip feature in the net series capacitance. However, this model does not exhibit the reduction in capacitance further away from the potential of zero charge (the hump).

More generally, the total series capacitance observed experimentally can be decomposed into a diffuse contribution and a concentration-independent but potential-dependent contribution containing the hump.12,13 Theories to explain this hump have historically included dielectric saturation of water, adsorption of ions, and ion packing.9,10,14–17

Capturing complex capacitance behavior in calculations that include electronic structure effects of the surface and the reaction species is essential for understanding double layer effects on adsorbates and chemical reactions.18–20 However, accurately capturing capacitance features with strictly ab initio methods would require ab initio molecular dynamics (AIMD) calculations that are unfeasibly large because of the large number of solvent molecules and electrolyte ions necessary for meaningful statistics, and the long equilibration time of the electrolyte distribution. A tractable alternative is to couple density-functional theory (DFT) calculations of the electronic-structure sensitive reactants and surface with a continuum representation of the electrolyte. A continuum solvation model suitable for this purpose must then capture the effect of non-adsorbing electrolytes, while the electrode and specific adsorbates are included explicitly in the DFT calculation.

In this paper, we demonstrate the shortcomings of current state-of-the-art solvation models for metallic electrodes in non-adsorbing, aqueous electrolytes and find that they do not reproduce the capacitance features that are experimentally observed. We then use a simple extension of the Gouy-Chapman-Stern model to illustrate the need for a dielectric-only region at the interface and nonlinear dielectric effects in the bulk electrolyte. We formulate the Nonlinear Electrochemical Soft-Sphere (NESS) solvation model, an extension of the recently developed soft-sphere (atomic-sphere-based) solvation model,21 that includes the necessary nonlinear effects and dielectric-only region. We show that NESS, used in conjunction with DFT calculations of the metal surface, can capture the magnitude, shape, and ionic concentration dependence of the measured electrochemical capacitance of Ag(100) immersed in aqueous KPF6, while simultaneously remaining accurate for solvation energies of ions and neutral molecules.

Continuum solvation models approximate the effect of a liquid environment by placing the solute (treated quantum-mechanically using DFT) in a cavity within a dielectric continuum. Ions in an electrolyte additionally contribute a Debye screening response in the continuum. Beyond this ansatz, solvation models differ in the construction and parametrization of the cavity, as well as the treatment of the continuum dielectric and ionic response. Broadly, determining the cavity from the solute electron density (iso-density) and/or electrostatic potential involves a few easily determined parameters, whereas determining the cavity using atom-centered spheres involves many parameters and more complexity but offers greater flexibility. Additionally, the approximations used to model the dielectric and ionic continuum regions can be either linear or nonlinear in response to the applied potential. Nonlinearity of the dielectric response arises from saturation of solvent dipoles aligning to the local electric field, while nonlinearity of the ionic response arises from the Boltzmann occupation factor of the ions at the local potential. The latter exponential increase in local ion concentration results in increased localization of ions near the electrode and increased diffuse contribution to the capacitance, giving rise to the capacitance dip as discussed for GCS theory above. Finally, the ionic and dielectric cavities can be separated by a distance typically denoted by x2 in GCS theory, or they can share the same cavity (x2 = 0).

Table I summarizes the above choices for the predominant solvation models used for electrochemistry, most of which adopt linear dielectric and ionic responses, precluding treatment of the capacitance dip and hump, as we show below. Additionally, simultaneous accuracy for properties of molecules, ions, and surfaces is challenging in iso-density models because of different ideal cavity sizes for capturing solvation energies of positive and negative solute charges25 and for capturing capacitance of metallic surfaces.30 The recent soft-sphere solvation model21 brings the added flexibility of atomic parametrization, traditionally employed only in finite molecular calculations,31,32 to calculations for periodic systems including solid-liquid interfaces. This model has so far not included ionic response or nonlinearity and has not yet been applied to electrochemical systems. We take advantage of the flexibility of the atomic parametrization of the soft-sphere model and add the key ingredients required for accurate electrochemical predictions that we show below: a separate ionic cavity (x2 ≠ 0) and nonlinear responses of the dielectric and ions, resulting in the NESS model.

TABLE I.

Categorization of solvation models for electrochemistry by linearity of fluid dielectric (ϵ) and ionic (κ) response, cavity parametrization type, and presence of a separate ionic cavity (x2 ≠ 0). Here L = linear, NL = nonlinear, n = solute electron density (iso-density models), ϕ = solute electrostatic potential and atom = atomic spheres.

ResponseCavities
ModelϵκTypex2 ≠ 0
LinearPCM22,23 = VASPsol24  n No 
CANDLE25  n, ϕ No 
SCCS26,27 n No 
Anderson and Jinnouchi28  NL n Yes 
Dabo et al.29  NL n Yes 
NonlinearPCM22,30 NL NL n No 
Soft-sphere21  None Atom No 
NESS (this work) NL NL Atom Yes 
ResponseCavities
ModelϵκTypex2 ≠ 0
LinearPCM22,23 = VASPsol24  n No 
CANDLE25  n, ϕ No 
SCCS26,27 n No 
Anderson and Jinnouchi28  NL n Yes 
Dabo et al.29  NL n Yes 
NonlinearPCM22,30 NL NL n No 
Soft-sphere21  None Atom No 
NESS (this work) NL NL Atom Yes 

We start by comparing the electrochemical capacitance of Ag(100) in the non-adsorbing KPF6 electrolyte against the capacitance predicted by previously developed solvation models.6 For each solvation model, we perform a series of grand-canonical DFT calculations33 in the JDFTx plane-wave DFT software34 for 25 uniformly spaced potentials in a range extending 0.8 V on either side of the neutral potential (PZC) and calculate the differential capacitance as a finite-difference derivative. We use a seven-layer Ag(100) slab with 12 × 12 × 1 k-point sampling, cold smearing with width σ = 0.01 Eh (see the supplementary material for convergence checks), the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional,35 and Garrity-Bennett-Rabe-Vanderbilt (GBRV) ultrasoft pseudopotentials36 with a 20 Eh (Hartree) plane-wave cutoff. We set the vacuum/electrolyte spacing to at least three times the Debye screening length for each ionic concentration, additionally employing truncated Coulomb interactions37 to isolate periodic images along the third direction.29 

First, Fig. 1(a) compares the experimental capacitance of Ag(100) with solvation models that solve the linearized Poisson-Boltzmann equation. Of these, the SCCS26 and LinearPCM22 models lead to a capacitance well below the experimental values, with a slight linear increase with potential, capturing none of the potential or electrolyte-concentration dependence, as expected and discussed previously.30 (Note that the commonly used VASPsol model24 re-implements and hence produces identical results to LinearPCM.22) The CANDLE25 model exhibits higher capacitance for negative potentials, along with an unphysical spike when transitioning from positive to negative potentials, because it adjusts the cavity size based on solute electric field to account for asymmetry in solvation of cations and anions. (See the supplementary material for a more detailed discussion on the capacitance of CANDLE.)

Next, Fig. 1(b) shows the effect of including dielectric saturation from the rotational response of water molecules as well as the response of electrolyte ions based on the nonlinear Poisson-Boltzmann equation, as in the NonlinearPCM model.22 The overall capacitance of this model is too low because of a large cavity that results in too low a Helmholtz capacitance, which also dominates the series capacitance and masks nonlinear effects in the diffuse capacitance. Refitting this model’s cavity to match the typical magnitude of the experimental capacitance30 exhibits some reduction in capacitance near the PZC for lower ionic concentrations. However, the capacitance still varies almost linearly with potential, with a larger positive slope. This asymmetry with respect to the potential of zero charge results from the electron-density parametrization of the solvent cavity: at positive potentials, the positively charged surface has a lower electron density tail extending from the surface, which pulls the cavity closer and increases the capacitance.23 This is a particularly large effect for the refit version of NonlinearPCM because the cavity is determined by a higher electron density threshold, which occurs relatively closer to the metal atoms.

To understand the limitations of existing DFT-based solvation approaches and identify the key pieces necessary to correctly capture electrochemical capacitance, we next study a model system illustrated in Fig. 2(a) that minimally extends GCS theory. The model is composed of a continuum plane of charge representing the metal, a continuum solvent, and a continuum electrolyte. The metal and solvent are separated by a small gap of width x1, and the ions are excluded from a solvent region of width x2. Assuming the metal is at potential ϕ0 and extends until x = 0, the potential ϕ(x) for x > 0 satisfies the nonlinear Poisson-Boltzmann equation

ϵ(x,|ϕ|)ϕ+κ2(x,ϕ)ϕ=0,
(1)

with

ϵ(x,E)1+Θ(xx1)(ϵ(E)1),
(2)
κ2(x,ϕ)Θ(xx1x2)8πNione2sinhϕkBTϕ
(3)

and with boundary conditions ϕ(0) = ϕ0 and ϕ() = 0. Here, ϵ(E) is the field-dependent dielectric constant of liquid water, Nion is the concentration of a 1:1 electrolyte, κ is the inverse Debye screening length, T is the temperature, e is the elementary charge, kB is the Boltzmann constant, and Θ(x) is the Heaviside step function. We use the implicit form ϵ(E)=ϵb/(1+4(D/D0)2+3(D/D0)4) with Dϵ(E)E, bulk dielectric constant ϵb = 78.4, and fit parameter D0 = 36 μC/cm2, which accurately reproduces Extended Simple Point Charge (SPC/E) molecular dynamics predictions for ϵ(E).38 We solve differential equation (1) using a 1D finite-element method for a series of ϕ0, calculate the metal surface charge density σ=4πϕxx=0 for each, and extract the differential capacitance ∂σ/∂ϕ0 from a finite difference derivative. In order to consider if ion packing is responsible for the hump, we also solve the modified Poisson-Boltzmann (mPB) equation,17 where the effective κ2 is instead given by

κ2(x,ϕ)Θ(xx1x2)8πNione2sinhϕkBTϕ1+4ηsinh2ϕ2kBT.
(4)

For the ion packing fraction η4π3(Rcat3+Ran3)Nion, we use an upper-bound cation radius of Rcat = 2.79 Å for K+ including its hydration shell,39 and an anion radius of Ran = 2.41 Å for the PF6 ion.40 

FIG. 2.

Experimental capacitance of Ag(100)6 compared against the toy model illustrated in (a) containing vacuum and solvent gaps of width x1 and x2, respectively, between the metal electrode and electrolyte. (b) Effect of x1 on model capacitance for a linear dielectric ϵ(0), where x1 = 0 is GCS theory. [(c) and (d)] Effect of nonlinearity (ϵ(E)) and ion packing as described by the modified Poisson-Boltzmann (mPB) theory on model capacitance for (c) x2 = 0 and (d) x2 = 2.6 Å; both nonzero solvent gap and dielectric saturation are needed to match experimental shape and mPB theory does not produce a hump to the extent observed in experiment. (e) Electrolyte concentration dependence of the capacitance: nonlinear toy model with x1 and x2 reproduces hump and dip features well.

FIG. 2.

Experimental capacitance of Ag(100)6 compared against the toy model illustrated in (a) containing vacuum and solvent gaps of width x1 and x2, respectively, between the metal electrode and electrolyte. (b) Effect of x1 on model capacitance for a linear dielectric ϵ(0), where x1 = 0 is GCS theory. [(c) and (d)] Effect of nonlinearity (ϵ(E)) and ion packing as described by the modified Poisson-Boltzmann (mPB) theory on model capacitance for (c) x2 = 0 and (d) x2 = 2.6 Å; both nonzero solvent gap and dielectric saturation are needed to match experimental shape and mPB theory does not produce a hump to the extent observed in experiment. (e) Electrolyte concentration dependence of the capacitance: nonlinear toy model with x1 and x2 reproduces hump and dip features well.

Close modal

Figure 2(b) shows the capacitance predicted by this toy model for the case of the linear dielectric, ϵ(0). With x1 = 0, the model reduces to GCS theory and the predicted capacitance is much higher than experiment for a typical value of x2 = 2.6 Å. A small vacuum region, x1 = 0.1 Å, is necessary to make the magnitude of the capacitance reasonable; note that this small vacuum gap has the same capacitance as a solvent gap x2 = 7.8 Å due to the bulk dielectric constant of water, ϵb = 78.4. With a linear dielectric, GCS theory and this new model both produce a constant capacitance far from PZC. Figures 2(c) and 2(d) demonstrate that introducing ϵ(E), the nonlinear dielectric, which includes dielectric saturation, produces the hump feature (the reduction of capacitance far from PZC). Also note that the shape of the hump agrees best with experimental measurements for a theory which includes a solvent-only region, characterized by non-zero x2 = 2.6 Å in Fig. 2(d). When x2 = 0, as in Fig. 2(c), the width and the height of the capacitance hump are both too large compared to the experiment. Figures 2(c) and 2(d) also show the effects of ion packing, as described by mPB theory with a modified κ2 as discussed above. Inclusion of ion packing produces a hump-like feature, but with a substantially smaller magnitude than observed in experiment. While this effect may be dominant for ionic liquids,17 dielectric saturation is the more important effect for the non-adsorbing aqueous electrolyte because of water’s high dielectric constant, the small size of the ions, and the fact that water molecules (rather than ions) are closest to the surface. Additionally, the NonlinearPCM solvation model22 already includes ion packing effects at a level comparable to mPB theory17 and yet it did not produce an observable hump in Fig. 1.

Figure 2(e) shows that the nonlinear toy model with non-zero x1 and x2 has both the expected capacitance maximum and dips depending on the ionic concentration. However, relatively minor disagreements remain with the experimental curve, including the missing plateaus and slight asymmetry far from PZC, speculated to arise from specific adsorption of water or ions.10,41,42 Understanding the origin of these secondary features seen experimentally and capturing them in a solvation model will be the subject of future work. For instance, future studies could consider the interplay between dielectric saturation and electrolyte packing using a method such as classical density-functional theory with molecular orientations for the solvent.38,43

Here, we focus on obtaining a continuum solvation model for DFT that captures the hump and dip features qualitatively, and Fig. 2 indicates that such a solvation model must be nonlinear in both the dielectric and ionic responses, and it must have a solvent region that excludes ions (x2 ≠ 0). (Continuum solvation models already include x1 ≠ 0 due to the gap between the electron density of the solute and the solvent cavity.) Table I illustrates that none of the previously developed solvation models simultaneously include both nonlinearities and x2 ≠ 0. Additionally, our previous work25,30 found that electron density-based parametrization of solvation model cavities results in limited simultaneous accuracy for properties of molecules, ions, and metal surfaces.

We therefore construct a solvation model with cavities based on atom-centered spheres, which is also suitable for calculations of extended systems in the plane-wave basis, starting from the recently proposed soft-sphere solvation model.21 In this model, the cavity shape function is determined as

s(r)=i12erfcRi|rri|σ,
(5)

which goes smoothly from zero in the vicinity of atoms located at positions ri to one in the solvent region more than the radius Ri away, with a transition width σ set to 0.5 bohr. The sphere radii Ri=fRi0, where f is a scale parameter fit to solvation energies, and Ri0 are a set of standard atomic radii, chosen to be the universal force field (UFF)44 van der Waals radii in Ref. 21.

The soft-sphere solvation model includes only linear dielectric response and neglects Debye screening from the electrolyte and hence cannot correctly capture electrochemical capacitance. Therefore, we construct the Nonlinear Electrochemical Soft-Sphere (NESS) solvation model by using the soft-sphere cavity definition in NonlinearPCM22 and including an additional ionic cavity so that x2 is non-zero. Specifically we set the spatial modulation s(r) of the nonlinear dielectric response in NonlinearPCM, exactly as given by (5), (instead of a threshold nc on the electron density n(r) as in Ref. 22). Further, we set the spatial modulation of the nonlinear ionic response to a separate cavity function sion(r), which is still given by (5) except with modified radii Riion=fRi0+x2, which introduces an electrolyte-excluded solvent region of width x2. All remaining details of the model are identical to NonlinearPCM and share most of its implementation in JDFTx,34 including the nonlinear response functions, the term with effective surface tension t to capture non-electrostatic free energy contributions and the algorithms to converge the solvation free energy and to handle charged solutes in periodic boundary conditions. (See Refs. 22 and 33 for further details.)

We perform solvation free energy calculations of a large set of neutral solutes, cations, and anions, including ionic screening (1M electrolyte) for the ion calculations. We fit the empirical parameters f (cavity scale factor) and t (effective tension = α + γ in Ref. 21’s notation) to the solvation free energies. When fit to neutral solutes, we find optimum parameters f = 1.00 and t = 1.02 × 10−5, which we use for subsequent calculations here. Table II compares the accuracy of this parametrization for solvation energies of neutral solutes, cations, and anions with that of several previous solvation models. (See the supplementary material for further details on the parametrization, including fits to cation or anion solvation energies alone.) Our new NESS model performs comparably to the linear soft-sphere model and, in its parametrization to neutral solutes, is substantially more accurate than the original NonlinearPCM.22 It still falls short of CANDLE,25 which is able to deliver the best overall accuracy for solvation free energies within a single parametrization because it accounts for charge asymmetry. However, that feature of CANDLE results in unphysical capacitance predictions, as shown in Fig. 1 and discussed above and in the supplementary material.

TABLE II.

Comparison of mean absolute errors (MAEs) in aqueous solvation energies of the same set of 240 neutral molecules, 51 cations, and 55 anions between several prominent grid/plane-wave based solvation models. CANDLE yields the lowest MAE in a single parametrization (in boldface below) due to charge-asymmetry correction but that leads to unphysical features in the capacitance (Fig. 1). The NESS model parametrized to neutral solvation energies is reasonably accurate for solvation energies (highlighted below), while simultaneously predicting electrochemical capacitance in good agreement with experiment (Fig. 3). (See the supplementary material for additional comparison of ionic parametrizations of the SCCS, soft-sphere and NESS models.)

MAE (kcal/mol)
ModelNeutralCationsAnionsAll
Linear dielectric 
LinearPCM22 = VASPsol24  1.27 2.10 15.1 3.59 
SCCS26  1.20 2.55 17.4 3.97 
CANDLE25  1.27 2.62 3.46 1.81 
Soft-sphere21  1.25 6.02 10.4 3.41 
Nonlinear dielectric 
Orig. nonlinearPCM22  1.28 16.1 27.0 7.55 
Refit nonlinearPCM30  1.44 14.4 27.6 7.51 
NESS (this work) 1.29 3.67 9.58 2.96 
MAE (kcal/mol)
ModelNeutralCationsAnionsAll
Linear dielectric 
LinearPCM22 = VASPsol24  1.27 2.10 15.1 3.59 
SCCS26  1.20 2.55 17.4 3.97 
CANDLE25  1.27 2.62 3.46 1.81 
Soft-sphere21  1.25 6.02 10.4 3.41 
Nonlinear dielectric 
Orig. nonlinearPCM22  1.28 16.1 27.0 7.55 
Refit nonlinearPCM30  1.44 14.4 27.6 7.51 
NESS (this work) 1.29 3.67 9.58 2.96 

Using the new NESS model parametrization, we perform capacitance calculations of the Ag(100) surface. We note that for the capacitance calculations, the silver atom radius in the original soft-sphere model is small enough so that the fluid leaks into the interior of the silver slab. For these calculations, we have modified our model to exclude the solvent from the interior of the slab. Figure 3 shows that the NESS model predicts the electrochemical capacitance of Ag(100) in much better agreement with experiment than any previous solvation model, including both the hump and the dip features. However, the capacitance maximum is not centered at the PZC as was the case for the toy model in Fig. 2(e). In this solvation model, the cavity is independent of the electron density, while the electrons extend further from the metal atoms with increasing negative charge, reducing the gap, and increasing the capacitance. As x2 increases, the contribution of this gap capacitance to the total series capacitance becomes less important, reducing the asymmetry of the hump, as seen in Fig. 3. In reality, the solvent likely moves further away with increasing negative charge, due to excess electrons (as captured by the iso-density cavities), resulting in a more symmetric hump.

FIG. 3.

NESS solvation model predicted capacitance for x2 = 0.5 and 2.6 Å, exhibiting both the hump and dip features of the experimental capacitance.

FIG. 3.

NESS solvation model predicted capacitance for x2 = 0.5 and 2.6 Å, exhibiting both the hump and dip features of the experimental capacitance.

Close modal

To more conclusively identify the source of the asymmetry in the capacitance curves from the NESS solvation model, we use the electron density results from the NESS calculations as inputs for the x1 parameter in our toy model. Figure 4(a) depicts the induced charge at the Ag(100) surface as a function of applied potential, assigning the distance between the electron density peak and the cavity onset as the x1eff parameter, plotted in Fig. 4(b) as a function of the electrode surface charge density σ. Figure 4(c) depicts the results of the toy model, using the σ-dependent x1eff results from Fig. 4(b). The toy model capacitance in Fig. 4(c) and the NESS capacitance of Fig. 3 are in very good agreement, confirming that the asymmetry in the NESS model is caused by the dielectric cavity remaining spatially fixed while the induced charge density shifts with potential.

FIG. 4.

(a) Induced charge density at the Ag(100) surface as a function of applied potential, from bottom to top: +0.8 V to −0.8 V relative to the PZC, as in Fig. 3. Dotted line indicates the location of dielectric cavity turn-on, and the solid line indicates the maximum of the induced charge density. (b) The separation between these locations in part (a) is the effective vacuum distance, x1eff between the metal and the solvent as defined in the toy model above, but which now depends on the electrode surface charge density σ. (c) Capacitance curves predicted by the toy model with this σ-dependent x1eff.

FIG. 4.

(a) Induced charge density at the Ag(100) surface as a function of applied potential, from bottom to top: +0.8 V to −0.8 V relative to the PZC, as in Fig. 3. Dotted line indicates the location of dielectric cavity turn-on, and the solid line indicates the maximum of the induced charge density. (b) The separation between these locations in part (a) is the effective vacuum distance, x1eff between the metal and the solvent as defined in the toy model above, but which now depends on the electrode surface charge density σ. (c) Capacitance curves predicted by the toy model with this σ-dependent x1eff.

Close modal

The NESS solvation model brings the computational electrochemistry community closer to developing an all-purpose solvation model for ions, molecules, and electrified surfaces, but there are still many opportunities for improvement and deeper understanding. For instance, improvements to the dielectric cavity definition and to the dielectric response at the interface have potential to increase accuracy. For the dielectric cavity, the original soft-sphere model and NESS employ an atomic radius for each element, currently set based on vdW radii from the UFF force field.44 These radii have been carefully evaluated for the few atoms present in organic solutes for which extensive solvation free energy data are readily available. This choice fortuitously gave good results for Ag(100) capacitance, but that might not be the case for all metal surfaces, and for all properties. Additionally, the solvent cavity presently remains fixed with potential and charge on the metal surface; understanding and incorporating the variation of water-metal distance with surface charge could result in further improvements in accuracy for the capacitance. We note that previous solvation models with isodensity cavity parametrizations22,25–28,30,45 mitigate some of this asymmetry but have not yet simultaneously captured solvation energetics and electrochemical capacitance within a single parametrization. Therefore, future work may need to employ a hybrid approach that combines the flexibility of atomic-sphere based cavities with the response to electron density inherent to isodensity cavities. Additionally, for the dielectric response at the interface and near ions, we have assumed that the response of the water remains unchanged from bulk; future improvements could include modulation of the dielectric response adjacent to the solute.

In summary, we find that a solvation model must contain nonlinearities in both the ionic and dielectric responses, along with a region of space containing only the dielectric, in order to capture the features seen in the experimental capacitance curves of metals such as Ag(100) immersed in non-adsorbing electrolytes. By including these characteristics, the NESS solvation model substantially improves on previous solvation models for electrochemical systems, qualitatively capturing key features in the experimental electrochemical capacitance.

See supplementary material for the convergence of the capacitance with the number of layers of Ag, number of k-points, and electron density smearing width. Additionally, the supplementary material contains a discussion of the capacitance behavior of the CANDLE solvation model and parametrization details of the NESS solvation model developed in this paper.

R.S. acknowledges start-up funding from the Department of Materials Science and Engineering at Rensselaer Polytechnic Institute and computational resources provided by the Center for Computational Innovations (CCI) at Rensselaer Polytechnic Institute. K.L.W.’s use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357.

1.
S. D.
Fried
and
S. G.
Boxer
,
Acc. Chem. Res.
48
,
998
(
2015
).
2.
M.
Fleischmann
,
P. J.
Hendra
, and
A. J.
McQuillan
,
J. Chem. Soc., Chem. Commun.
0
(
3
),
80
(
1973
).
3.
J.
Wiechers
,
T.
Twomey
,
D.
Kolb
, and
R.
Behm
,
J. Electroanal. Chem. Interfacial Electrochem.
248
,
451
(
1988
).
4.
J.
Clavilier
,
J. Electroanal. Chem.
107
,
205
(
1980
).
5.
A.
Baskin
and
D.
Prendergast
,
J. Electrochem. Soc.
164
,
E3438
(
2017
).
6.
G.
Valette
,
J. Electroanal. Chem. Interfacial Electrochem.
138
,
37
(
1982
).
7.
A.
Hamelin
, in
Trends in Interfacial Electrochemistry
, NATO ASI Series, edited by
A.
Silva
(
Springer Netherlands
,
1986
), Vol. 179, Chap. 5, pp.
83
102
.
8.
9.
S.
Trasatti
,
Mod. Aspects Electrochem.
13
,
81
(
1979
).
10.
B.
Damaskin
and
A.
Frumkin
,
Electrochim. Acta
19
,
173
(
1974
).
11.
O.
Stern
,
Z. Elektrochem. Angew. Phys. Chem.
30
,
508
(
1924
).
13.
R.
Parsons
and
F.
Zobel
,
J. Electroanal. Chem.
9
,
333
(
1965
).
14.
W.
Schmickler
and
D.
Henderson
,
J. Chem. Phys.
85
,
1650
(
1986
).
15.
R.
Saradha
and
M.
Sangaranarayanan
,
J. Phys. Chem. B
102
,
5468
(
1998
).
16.
S.
Amokrane
and
J.
Badiali
, in
Modern Aspects of Electrochemistry
, edited by
J.
Bockris
(
Springer
,
Boston, MA
,
1992
), Vol. 22.
17.
T.
Markovich
,
D.
Andelman
, and
R.
Podgornik
, in
Handbook of Lipid Membranes: Molecular, Functional, and Materials Aspects
, edited by
C.
Safinya
and
J.
Radler
(
Taylor & Francis
,
2014
), Chap. 9.
18.
H.
Xiao
,
T.
Cheng
,
W. A.
Goddard
 III
, and
R.
Sundararaman
,
J. Am. Chem. Soc.
138
,
483
(
2016
).
19.
A. T. B.
Jason
,
D.
Goodpaster
, and
M.
Head-Gordon
,
J. Phys. Chem. Lett.
7
,
1471
(
2016
).
20.
A. J.
Garza
,
A. T.
Bell
, and
M.
Head-Gordon
,
J. Phys. Chem. Lett.
9
,
601
(
2018
).
21.
G.
Fisicaro
,
L.
Genovese
,
O.
Andreussi
,
S.
Mandal
,
N. N.
Nair
,
N.
Marzari
, and
S.
Goedecker
,
J. Chem. Theory Comput.
13
,
3829
(
2017
).
22.
D.
Gunceler
,
K.
Letchworth-Weaver
,
R.
Sundararaman
,
K. A.
Schwarz
, and
T. A.
Arias
,
Modell. Simul. Mater. Sci. Eng.
21
,
074005
(
2013
).
23.
K.
Letchworth-Weaver
and
T. A.
Arias
,
Phys. Rev. B
86
,
075140
(
2012
).
24.
K.
Mathew
,
R.
Sundararaman
,
K.
Letchworth-Weaver
,
T. A.
Arias
, and
R. G.
Hennig
,
J. Chem. Phys.
140
,
084106
(
2014
).
25.
R.
Sundararaman
and
W.
Goddard
,
J. Chem. Phys.
142
,
064107
(
2015
).
26.
O.
Andreussi
,
I.
Dabo
, and
N.
Marzari
,
J. Chem. Phys.
136
,
064102
(
2012
).
27.
C.
Dupont
,
O.
Andreussi
, and
N.
Marzari
,
J. Chem. Phys.
139
,
214110
(
2013
).
28.
R.
Jinnouchi
and
A.
Anderson
,
Phys. Rev. B
77
,
245417
(
2008
).
29.
I.
Dabo
,
Y.
Li
,
N.
Bonnet
, and
N.
Marzari
, in
Fuel Cell Science: Theory, Fundamentals, and Biocatalysis
, edited by
A.
Wieckowski
and
J.
Norskov
(
Wiley
,
2010
), Chap. 13.
30.
R.
Sundararaman
and
K.
Schwarz
,
J. Chem. Phys.
146
,
084111
(
2017
).
31.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
,
Chem. Rev.
105
,
2999
(
2005
).
32.
A. V.
Marenich
,
C. J.
Cramer
, and
D. G.
Truhlar
,
J. Phys. Chem. B
113
,
6378
(
2009
).
33.
R.
Sundararaman
,
W. A.
Goddard
 III
, and
T. A.
Arias
,
J. Chem. Phys.
146
,
114104
(
2017
).
34.
R.
Sundararaman
,
K.
Letchworth-Weaver
,
K.
Schwarz
,
D.
Gunceler
,
Y.
Ozhabes
, and
T. A.
Arias
,
SoftwareX
6
,
278
(
2017
).
35.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
36.
K. F.
Garrity
,
J. W.
Bennett
,
K.
Rabe
, and
D.
Vanderbilt
,
Comput. Mater. Sci.
81
,
446
(
2014
).
37.
R.
Sundararaman
and
T.
Arias
,
Phys. Rev. B
87
,
165122
(
2013
).
38.
R.
Sundararaman
,
K.
Letchworth-Weaver
, and
T. A.
Arias
,
J. Chem. Phys.
140
,
144504
(
2014
).
39.
40.
H. K.
Roobottom
,
H. D. B.
Jenkins
,
J.
Passmore
, and
L.
Glasser
,
J. Chem. Educ.
76
,
1570
(
1999
).
41.
S.
Trasatti
,
J. Electroanal. Chem.
150
,
1
(
1983
).
42.
G.
Valette
,
J. Electroanal. Chem.
122
,
285
(
1981
).
43.
R.
Sundararaman
and
T.
Arias
,
Comput. Phys. Commun.
185
,
818
(
2014
).
44.
A. K.
Rappe
,
C. J.
Casewit
,
K. S.
Colwell
,
W. A.
Goddard
, and
W. M.
Skiff
,
J. Am. Chem. Soc.
114
,
10024
(
1992
).
45.
J. L.
Fattebert
and
F.
Gygi
,
Int. J. Quantum Chem.
93
,
139
(
2003
).

Supplementary Material