Fully harnessing electrochemical interfaces for reactions requires a detailed understanding of solvent effects in the electrochemical double layer. Predicting the significant impact of solvents on entropic and electronic properties of electrochemical interfaces has remained an open challenge of computational electrochemistry. Using molecular dynamics simulations of silver–water and silver–acetonitrile interfaces, we show that switching the solvent changes the signs for both the charge of maximum capacitance (CMC) and charge of maximum entropy (CME). Contrasting the capacitance and CME behavior of these two interfaces, we demonstrate that the preferred orientation of the solvent molecule and the corresponding charge density determine the sign of the CMC and CME and, hence, the qualitatively different charge asymmetry of the electrochemical interface.

The molecular structure of the electrochemical double layer plays a key role in catalysis by determining both the access of the reaction species to the electrode and the subsequent reaction mechanisms.1 Solvent molecules can block reactants, facilitate access to the electrode, or act as reactants in industrially relevant electrochemical reactions, such as the hydrogen evolution reaction, the oxygen reduction reaction, and the hydrogenation of acetonitrile for chemical manufacturing.2 The molecular structure of the double layer changes as the electrode charges, changing the capacitance3,4 and entropy.5 We have recently shown using molecular dynamics (MD) simulations of aqueous interfaces that both capacitance and entropy are asymmetric with electrode charge, with charges of maximum capacitance (CMC) and entropy (CME) both negative, but with distinct magnitudes.6,7 However, the relationship between the molecular structure of the double layer and the charges of maximum capacitance and entropy is yet to be unraveled.

Here, we predict CMC and CME from MD simulations of Ag(100) interfaces with water and acetonitrile to identify the solvent effects on these quantities. We first compare MD simulations with an explicit electrolyte against those with a pure-solvent + continuum electrolyte (i.e., an explicit solvent layer followed by a region with solvent plus a continuum electrolyte) and observe only modest changes in both the CMC and CME, in agreement with small changes of the potential of maximum entropy with electrolyte in experiments.8,9 In contrast, we find significant contributions to both the CMC and CME from the solvent and find that their origins lie in the same molecular phenomenon. Molecular properties of the solvent determine the sign of the CMC and CME, with both flipping from negative for water to positive for acetonitrile.

Molecular dynamics simulations of an aqueous NaF–Ag(100) interface find that the electrochemical capacitance increases for negatively charged interfaces, reaching a maximum at a CMC of −3.7 μC/cm2 [Fig. 1(a)].7 Here, we replace the explicit NaF electrolyte with a continuum electrolyte beyond an inner-layer width parameterized by z2 (discussed below). Figure 1(a) shows that z2 controls the overall magnitude of capacitance, but does not affect the shape of the capacitance curve, with the maximum remaining at a CMC of −3.1 μC/cm2. Continuum electrolyte predictions for a physically reasonable inner-layer width of z2 ≈ 3 Å closely match explicit electrolyte predictions: the CMC is slightly smaller, but the peak shape and the most of the charge asymmetry are captured well. This indicates that the inner layer capacitance is determined by the solvent and is not modified substantially by the ions in a non-adsorbing electrolyte.3 

Briefly, to predict electrochemical capacitance and entropy from pure-solvent MD simulations + continuum electrolyte, we apply an electric field to the solvent using charge densities ±σ on two electrodes in a simulation cell that is periodic in x, y and non-periodic in the z-direction [Fig. 1(b)]. (See the supplementary material for further simulation details.) We then fit the linear electrostatic potential profile in the central bulk region to extract the field-dependent dielectric constant ϵb(σ) from the slope and the intercepts ϕint(+σ) and ϕint(−σ) on both electrodes [Fig. 1(c)], where σ = D, the externally applied electric field throughout the solvent region. From these quantities and the vacuum permittivity ϵ0, we can reconstruct the effective inner layer potential ϕil(σ) = ϕint(σ) − z1σ/ϵ0 + z2σ/(ϵ0ϵb(σ)), where z1 accounts for the distance between the electrode atoms and its charge response plane (Δz in Ref. 6), while z2 is the width of the solvent-only “inner-layer” region of the electrochemical interface as in the Gouy–Chapman–Stern model. Note that we use the bulk permittivity ϵb(σ) only to define the interface potential ϕint(σ) as a way of neatly separating bulk and interface contributions between these two charge dependent functions. The actual permittivity changes dramatically near the interface [where the potential oscillates in Fig. 1(c)] and is captured within ϕint(σ) in our approach. Finally, we calculate the overall potential in the electrochemical interface as ϕ(σ) = ϕil(σ) + ϕd(σ), where ϕd(σ) is the potential across the diffuse layer at a specified ionic concentration from the Debye model (see the supplementary material), effectively constructing the potential of the half-cells shown in Fig. 1(d), and predict the differential capacitance C(σ) = (σ)/. This allows us to predict the overall capacitance of ideal electrochemical interfaces for different ionic concentrations starting from a single set of pure-solvent MD simulations. As discussed above, we find that z2 = 3 Å of solvent, plus continuum electrolyte, reasonably captures explicit 1 mol/liter electrolyte predictions, while the electrode-plane distance z1 = 0.46 Å remains unchanged from explicit electrolyte simulations.6,7

Now that we have shown that the solvent MD + continuum electrolyte closely approximates explicit electrolyte MD simulations with z1 = 0.46 Å and z2 = 3.0 Å, we fix these parameters and the continuum electrolyte concentration of 1 mol/liter and investigate the solvent dependence. Figure 2(a) shows that water and acetonitrile exhibit opposite asymmetry in the capacitance, with charges of maximum capacitance (CMC) of (−3.1 ± 0.3) and (+4.4 ± 0.6) μC/cm2, respectively. Next, to predict asymmetry in the entropy of the interface, we simulate heating of the interfaces, directly mimicking experimental measurements of the charge of maximum entropy (CME).10 Specifically, from the thermodynamic relation ∂S/∂σ|T = −∂V/∂T|σ, where V is the interface potential difference, the maximum of entropy S(σ) occurs where ∂V/∂T|σ crosses 0 with a positive slope. From ∂V/∂T|σ shown in Fig. 2(b), computed from a finite difference of MD simulations at 298 and 318 K,6 we see that water and acetonitrile also exhibit opposite asymmetry in entropy with CMEs of (−5.8 ± 0.3) and (+4.7 ± 0.6) μC/cm2, respectively. The uncertainty in the capacitance [shaded regions in Figs. 2(a) and 2(b)], and hence in the CMC and CME, is calculated by resampling from five independent MD simulations at each charge, as discussed in Ref. 6. Note that the sign of the CMC and CME is identical for each solvent. Their magnitudes differ for water, but agree within the uncertainties for acetonitrile.

To break down contributions to the capacitance and entropy predictions, Figs. 2(c) and 2(d), respectively, show the bulk dielectric constant ϵb from the slope and the intercept ϕint of MD electrostatic potential profiles [defined in Fig. 1(c)]. The bulk dielectric constant decreases symmetrically with the increasing charge magnitude due to dielectric saturation with the applied electric field E = σ/ϵ0. This saturation is captured well by the Booth equation for bulk water and acetonitrile11,12 and is relevant for electrochemical interfaces with ϵb reducing by more than a factor of 2 over the range of relevant electrode charges.13 This bulk ϵb is relevant in determining the diffuse-layer and, hence, the total capacitance, but is distinct from the changes in the interface dielectric constant. These changes are captured within ϕint, which encapsulates all interfacial contributions in the present approach. Since ϵb is symmetric by definition, all the asymmetry in the capacitance and entropy is necessarily from the intercept potential, ϕint. Note that ϕint at zero charge has opposite signs for water and acetonitrile, increases monotonically with charge, and crosses zero at negative and positive signs, respectively. However, the zero-crossing of ϕint is not directly connected to CMC or CME; the capacitance contribution from the interface is related to the second charge derivative 2ϕint/∂σ2, while the entropy contribution is related to the temperature derivative ∂ϕint/∂T.

To connect CMC and CME to the molecular behavior of the solvent at the interface, Fig. 3 compares the solvent charge density distribution and its temperature derivative at the interface for various electrode charges. Note that the solvent response is expected to dominate the overall charge response in the inner-layer region shown here, an idealized interface of a silver electrode and weakly interacting, non-adsorbing electrolyte. The charge density of water closest to the electrode is positive at zero electrode charge [gray line in Fig. 3(a)], indicating an interfacial water orientation with the hydrogen slightly closer to the electrode on average. This slight charge imbalance, with a small positive (negative) charge density near (away from) the electrode, is consistent with experimental predictions of water molecules in a flat orientation at neutral silver–water interfaces.14,15 This charge density increases in magnitude for negative electrode charges, but for positive electrode charges decreases in magnitude first, then crosses zero (near σ = 8 μC/cm2), and finally increases in magnitude with oxygen now facing the electrode. Additionally, the location of the charge density is much closer to the electrode when the hydrogen faces the negative electrode, compared to when the oxygen is next to the positive electrode. This molecular behavior leads to a stronger solvent response to negative charges, leading to a negative CMC, as discussed previously.6 

The above behavior is reversed overall from water to acetonitrile. The solvent charge nearest the neutral electrode [gray line in Fig. 3(c)] is negative, indicating the CN end of the molecule facing the electrode. For a large enough magnitude of σ < 0, this charge density flips positive, with CH3 facing the electrode. This behavior is in agreement with Sum Frequency Generation (SFG) measurements on Pt electrodes, where acetonitrile binds strongly,16 and reversibly flips between CH3 and CN toward the metal at low and high potentials, respectively.17 For acetonitrile, the negative (CN) solvent charge approaches positively charged electrodes more closely than the positive (CH3) solvent charge approaches negative electrodes, leading to a stronger charge response for σ > 0 and, hence, a positive CMC. Thus, the opposite sign of CMC for water and acetonitrile stems from the closer electrode approach of the positive H end of water and the negative CN end of acetonitrile, compared to negative O and positive CH3 ends, respectively.

To understand the entropy response and CME, we need to compare the behavior of ρ and /dT for both solvents. To relate the charge of maximum entropy to the molecular interaction at the interface, we first consider the neutral interface. For neutral interfaces of both water and acetonitrile, the peaks of ρ closest to the electrode in Figs. 3(a) and 3(c) line up with the zero crossings of /dT in Figs. 3(b) and 3(d) (connected by gray vertical lines). This signature, along with the direction of zero crossing, corresponds to a charge density that moves closer to the electrode with the increasing temperature. This is expected from a thermodynamic standpoint: as the temperature increases, more thermal energy is available to overcome the repulsive potential at the interface. In contrast, for highly charged electrodes with σ = ±15 μC/cm2, the peaks of ρ line up with opposite peaks of /dT (green and magenta vertical lines), which indicates that the ρ peak primarily broadens with the increasing temperature. The commonality of these behaviors leads to the sign of ρ and /dT nearest the interface being linked. This leads to the CMC and CME having the same sign, in general, as they stem from the same molecular orientation preferences, but they could differ in magnitude as that depends on further details of the overall complex charge distributions and temperature derivatives in Fig. 3.

Finally, to aid direct experimental comparisons, Fig. 4 plots the simulated capacitance as a function of electrode potential (instead of electrode charge in previous figures). Additionally, we present results for several ionic concentrations in the diffuse capacitance contribution, approximated using the Gouy–Chapman model (essentially the nonlinear Poisson–Boltzmann model). For both water and acetonitrile, the capacitance drop is centered at the potential of zero charge when the ionic concentration is sufficiently low, irrespective of the asymmetry of the solvent response, which affects the shape of the capacitance curves at larger magnitudes of potential. However, capacitance curves for both water18 and acetonitrile appear to be narrower than those observed experimentally.

This underestimated width of capacitance curves may be caused by missing chemisorption interactions in our simulation. Specifically, a classical force field may not be able to adequately represent the intrinsically quantum-mechanical strong attraction of the water to the metal interface, beyond the average metal–molecule interactions captured by parameterizing the wall interaction to DFT (see the supplementary material). AIMD simulations of the aqueous double layer of Pt(111)19 suggest that chemisorption of water onto the electrode surface leads to the experimentally observed bell-shaped capacitance feature.20 The chemisorption of water to silver interfaces is believed to be weaker than that of platinum, but the electronic interaction is still strong enough to significantly shift the potential of zero charge (PZC) away from its expected value based on the work function,21 and so it may still play an important role in determining the capacitance. Note, however, that ab initio molecular dynamics simulations are typically coupled with a capacitance model for the bulk solvent that is constant with potential19 and thus may not adequately describe solvent contributions to the capacitance that are captured by longer and larger classical MD simulations presented here.

For acetonitrile, we note that qualitatively the Ag–acetonitrile interface for both NaF and LiClO4 electrolytes has increasing capacitance with potential across the potential of zero charge.22,23 A more quantitative comparison of the capacitance curve width is difficult for acetonitrile because of experimental challenges with the chemical stability of the interface, including oxidation of the surface and contamination by solvent impurities.22,23

In conclusion, we related the CMC and CME to the molecular behavior of the solvent and compared the opposite-charged behavior of the silver interface with acetonitrile and water. Continuum treatment of ions does not significantly change the CMC or CME compared to explicit MD simulations of ions, allowing us to systematically investigate solvent effects in water and acetonitrile and setting the stage for expanding this analysis to several more solvents in future work. The natural orientation of solvent molecules at the neutral electrode leads to a closer charge approach to electrodes of one charge sign compared to the other, which directly leads to the common sign of the CMC and CME. The overall shape of capacitance from these classical MD simulations does not match experiment perfectly likely due to missing solvent chemisorption. Future work is needed to combine such effects captured naturally in AIMD with the longer sampling times and larger bulk regions captured in our classical MD simulations. This includes extending the present approach to a pure-solvent AIMD + continuum electrolyte as a strategy to reduce the required time scales.24 Combinations of electronic DFT with classical MD,25 AIMD simulations and the detailed solvent response models developed here hold promise for a more complete picture of electrified electrochemical interfaces.

See the supplementary material for details of the molecular dynamics simulations in the LAMMPS,26 including force field details for water (SPC/E)27 and acetonitrile;28 parameterization of silver–solvent interactions from electronic density functional theory calculations in JDFTx;29 and input files and analysis scripts for all results presented above.

Note: Certain software is identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology nor does it imply that the software identified is necessarily the best available for the purpose.

R.S. acknowledges the support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0022247. Calculations were carried out at the Center for Computational Innovations at Rensselaer Polytechnic Institute, at NIST, and at the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award ERCAP0020105.

The authors have no conflicts to disclose.

Ravishankar Sundararaman: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Kathleen Schwarz: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal).

The data that support the findings of this study are available within the supplementary material.

1.
M.
Yan
,
Y.
Kawamata
, and
P. S.
Baran
,
Chem. Rev.
117
,
13230
(
2017
).
2.
R.
Xia
,
D.
Tian
,
S.
Kattel
,
B.
Hasa
,
H.
Shin
,
X.
Ma
,
J. G.
Chen
, and
F.
Jiao
,
Nat. Commun.
,
12
1949
(
2021
).
3.
S.
Park
and
J. G.
McDaniel
,
J. Phys. Chem. C
126
,
16461
(
2022
).
4.
M. H.
Motevaselian
and
N. R.
Aluru
,
ACS Nano
14
,
12761
(
2020
).
5.
J. A.
Harrison
,
J. E. B.
Randles
, and
D. J.
Schiffrin
,
J. Electroanal. Chem. Interfacial Electrochem.
48
,
359
(
1973
).
6.
A.
Shandilya
,
K.
Schwarz
, and
R.
Sundararaman
,
J. Chem. Phys.
156
,
014705
(
2022
).
7.
A.
Shandilya
,
K.
Schwarz
, and
R.
Sundararaman
,
J. Chem. Phys.
156
,
129901
(
2022
).
8.
R. W.
Haid
,
X.
Ding
,
T. K.
Sarpey
,
A. S.
Bandarenka
, and
B.
Garlyyev
,
Curr. Opin. Electrochem.
32
,
100882
(
2022
).
9.
X.
Ding
,
D.
Scieszka
,
S.
Watzele
,
S.
Xue
,
B.
Garlyyev
,
R. W.
Haid
, and
A. S.
Bandarenka
,
ChemElectroChem
9
,
e202101088
(
2022
).
10.
V.
Climent
,
B. A.
Coles
, and
R. G.
Compton
,
J. Phys. Chem. B
106
,
5258
(
2002
).
11.
I. N.
Daniels
,
Z.
Wang
, and
B. B.
Laird
,
J. Phys. Chem. C
121
,
1025
1031
(
2017
).
12.
F.
Booth
,
J. Chem. Phys.
19
,
391
(
1951
).
13.
K.
Schwarz
and
R.
Sundararaman
,
Surf. Sci. Rep.
75
,
100492
(
2020
).
14.
Z. D.
Schultz
,
S. K.
Shaw
, and
A. A.
Gewirth
,
J. Am. Chem. Soc.
127
,
15916
(
2005
).
15.
K.-i.
Ataka
,
T.
Yotsuyanagi
, and
M.
Osawa
,
J. Phys. Chem.
100
,
10664
10672
(
1996
).
16.
O. A.
Petrii
and
I. G.
Khomchenko
,
J. Electroanal. Chem. Interfacial Electrochem.
106
,
277
(
1980
).
17.
S.
Baldelli
,
G.
Mailhot
,
P.
Ross
,
Y.-R.
Shen
, and
G. A.
Somorjai
,
J. Phys. Chem. B
105
,
654
(
2001
).
18.
G.
Valette
,
J. Electroanal. Chem. Interfacial Electrochem.
138
,
37
(
1982
).
19.
J.-B.
Le
,
Q.-Y.
Fan
,
J.-Q.
Li
, and
J.
Cheng
,
Sci. Adv.
6
,
eabb1219
(
2020
).
20.
K.
Ojha
,
K.
Doblhoff-Dier
, and
M. T. M.
Koper
,
Proc. Natl. Acad. Sci. U. S. A.
119
,
e21160161
(
2021
).
21.
J.
Le
,
M.
Iannuzzi
,
A.
Cuesta
, and
J.
Cheng
,
Phys. Rev. Lett.
119
,
016801
(
2017
).
22.
V. A.
Safonov
,
M. A.
Choba
, and
O. A.
Petrii
,
J. Electroanal. Chem.
808
,
278
(
2018
).
23.
V. A.
Safonov
,
A. G.
Krivenko
, and
M. A.
Choba
,
Electrochim. Acta
53
,
4859
4866
(
2008
).
24.
R.
Sundararaman
,
D.
Vigil-Fowler
, and
K.
Schwarz
,
Chem. Rev.
122
,
10651
(
2022
).
25.
S.-J.
Shin
,
D. H.
Kim
,
G.
Bae
,
S.
Ringe
,
H.
Choi
,
H.-K.
Lim
,
C. H.
Choi
, and
H.
Kim
,
Nat Commun.
13
,
174
(
2022
).
26.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
27.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
28.
V. A.
Koverga
,
O. M.
Korsun
,
O. N.
Kalugin
,
B. A.
Marekha
, and
A.
Idrissi
,
J. Mol. Liq.
233
,
251
(
2017
).
29.
R.
Sundararaman
,
K.
Letchworth-Weaver
,
K. A.
Schwarz
,
D.
Gunceler
,
Y.
Ozhabes
, and
T. A.
Arias
,
SoftwareX
6
,
278
(
2017
).

Supplementary Material