Improving accuracy of electrochemical capacitance and solvation energetics in first-principles calculations

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,³ and single-crystal electrochemistry.⁴ 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.⁵ 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 KPF₆ electrolyte,⁶ 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 

The concentration-dependent dip in the electrochemical capacitance can be explained by the Gouy-Chapman-Stern (GCS) model,¹¹ 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,²¹ 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 KPF₆, 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 x₂ in GCS theory, or they can share the same cavity (x₂ = 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 charges²⁵ and for capturing capacitance of metallic surfaces.³⁰ The recent soft-sphere solvation model²¹ 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 (x₂ ≠ 0) and nonlinear responses of the dielectric and ions, resulting in the NESS model.

We start by comparing the electrochemical capacitance of Ag(100) in the non-adsorbing KPF₆ electrolyte against the capacitance predicted by previously developed solvation models.⁶ For each solvation model, we perform a series of grand-canonical DFT calculations³³ in the JDFTx plane-wave DFT software³⁴ 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 E_h (see the supplementary material for convergence checks), the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional,³⁵ and Garrity-Bennett-Rabe-Vanderbilt (GBRV) ultrasoft pseudopotentials³⁶ with a 20 E_h (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 interactions³⁷ to isolate periodic images along the third direction.²⁹ 

First, Fig. 1(a) compares the experimental capacitance of Ag(100) with solvation models that solve the linearized Poisson-Boltzmann equation. Of these, the SCCS²⁶ and LinearPCM²² 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.³⁰ (Note that the commonly used VASPsol model²⁴ re-implements and hence produces identical results to LinearPCM.²²) The CANDLE²⁵ 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.²² 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 capacitance³⁰ 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.²³ 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 x₁, and the ions are excluded from a solvent region of width x₂. Assuming the metal is at potential ϕ₀ and extends until x = 0, the potential ϕ(x) for x > 0 satisfies the nonlinear Poisson-Boltzmann equation

with

and with boundary conditions ϕ(0) = ϕ₀ and ϕ(∞) = 0. Here, ϵ(E) is the field-dependent dielectric constant of liquid water, N_ion is the concentration of a 1:1 electrolyte, κ is the inverse Debye screening length, T is the temperature, e is the elementary charge, k_B 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 D₀ = 36 μC/cm², which accurately reproduces Extended Simple Point Charge (SPC/E) molecular dynamics predictions for ϵ(E).³⁸ We solve differential equation (1) using a 1D finite-element method for a series of ϕ₀, calculate the metal surface charge density $σ=−4π∂ϕ∂xx=0$ for each, and extract the differential capacitance ∂σ/∂ϕ₀ 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,¹⁷ where the effective κ² is instead given by

For the ion packing fraction $η≡4π3(Rcat3+Ran3)Nion$⁠, we use an upper-bound cation radius of R_cat = 2.79 Å for K⁺ including its hydration shell,³⁹ and an anion radius of R_an = 2.41 Å for the $PF6−$ ion.⁴⁰ 

Figure 2(b) shows the capacitance predicted by this toy model for the case of the linear dielectric, ϵ(0). With x₁ = 0, the model reduces to GCS theory and the predicted capacitance is much higher than experiment for a typical value of x₂ = 2.6 Å. A small vacuum region, x₁ = 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 x₂ = 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 x₂ = 2.6 Å in Fig. 2(d). When x₂ = 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 κ² 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,¹⁷ 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 model²² already includes ion packing effects at a level comparable to mPB theory¹⁷ and yet it did not produce an observable hump in Fig. 1.

Figure 2(e) shows that the nonlinear toy model with non-zero x₁ and x₂ 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 (x₂ ≠ 0). (Continuum solvation models already include x₁ ≠ 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 x₂ ≠ 0. Additionally, our previous work^25,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.²¹ In this model, the cavity shape function is determined as

which goes smoothly from zero in the vicinity of atoms located at positions $r→i$ to one in the solvent region more than the radius R_i 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)⁴⁴ 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 NonlinearPCM²² and including an additional ionic cavity so that x₂ 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 n_c 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 x₂. All remaining details of the model are identical to NonlinearPCM and share most of its implementation in JDFTx,³⁴ 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⁻⁵, 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.²² It still falls short of CANDLE,²⁵ 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.

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 x₂ 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.

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 x₁ 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.

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.⁴⁴ 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 parametrizations^{22,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.