A major challenge in the molecular simulation of electric double layer capacitors (EDLCs) is the choice of an appropriate model for the electrode. Typically, in such simulations the electrode surface is modeled using a uniform fixed charge on each of the electrode atoms, which ignores the electrode response to local charge fluctuations in the electrolyte solution. In this work, we evaluate and compare this Fixed Charge Method (FCM) with the more realistic Constant Potential Method (CPM), [S. K. Reed et al, J. Chem. Phys.126, 084704 (2007)], in which the electrode charges fluctuate in order to maintain constant electric potential in each electrode. For this comparison, we utilize a simplified LiClO4-acetonitrile/graphite EDLC. At low potential difference (ΔΨ ⩽ 2 V), the two methods yield essentially identical results for ion and solvent density profiles; however, significant differences appear at higher ΔΨ. At ΔΨ ⩾ 4 V, the CPM ion density profiles show significant enhancement (over FCM) of “inner-sphere adsorbed” Li+ ions very close to the electrode surface. The ability of the CPM electrode to respond to local charge fluctuations in the electrolyte is seen to significantly lower the energy (and barrier) for the approach of Li+ ions to the electrode surface.

Electric double-layer capacitors (EDLCs) are non-Faradaic, high power-density devices that have wide application in energy storage. Together with pseudo-capacitors, EDLCs make up a class of energy storage devices called supercapacitors. The energy storage and release mechanism in EDLCs is rapid and possesses a long cycle life due to the physical nature of the charging/discharging process—in the charging process, ions in the electrolyte solution aggregate at the interface to form an electric double layer, which in turn induces charges on the electrode surfaces. Over the past decade, accompanying the recent bloom of novel electrode (i.e., nanoporous materials) and electrolyte materials (i.e., ionic liquids), EDLCs have been the subject of numerous experimental studies.1,2 From these studies, the performance of EDLCs is affected by a number of different factors including, but not limited to, electrolyte composition,3,4 interface structure,5,6 and surface area.7,8 As a result of these efforts, the performance of EDLCs has been significantly enhanced, extending their applicability.

These experimental studies have been complemented by a number of theoretical/computational studies ranging from atomistic simulation to mesoscale continuum modeling.9 Analytical continuum models have been developed to describe the electrode/electrolyte interface, for example, the Gouy-Chapman-Stern model.10,11 Classical Density Functional Theories (cDFT)12 have also been applied to estimate properties of EDLCs, which are more accurate than continuum models but considerably less computationally demanding than atomistic simulation. Atomistic simulations, however, have an advantage over the continuum models and cDFT because they provide a molecular-level description of the structure, dynamics, and thermodynamics of EDLCs. Molecular simulation tools, such as ab initio molecular dynamics (AIMD) simulation13 and classical molecular dynamics (MD) simulation are widely used to investigate EDLC interfacial phenomena at the molecular level.

In molecular simulations of EDLCs, the modeling of the electrode is a particular challenge because of the difficulty in defining a consistent classical atomistic model for a conductor. In many MD studies of EDLCs, the electrode atoms are assumed to carry a uniform fixed charge. This Fixed Charge Method (FCM), however, neglects charge fluctuations on the electrode induced by local density fluctuations in the electrolyte solution.14–18 To explicitly take into account such fluctuations, the Constant Potential Method (CPM) was developed by Reed et al.19 This method is based on earlier work of Siepmann and Sprik20 in which the constraint of constant electrode potential was enforced on average using an extended Hamiltonian approach (similar to the Nosé method21 for constant temperature); however, in the CPM the constraint is applied instantaneously at every step. Additional corrections to the CPM were added later by Gingrich and Wilson.22 In the CPM, the electric potential Ψi on each electrode atom is constrained at each simulation step to be equal to a preset applied external potential V, which is constant over a given electrode. This constraint leads to the following equation for the charge, qi, on each electrode atom (where i indexes the atoms in the electrode):

\begin{equation}V=\Psi _{i}=\frac{\partial U}{\partial q_{i}},\end{equation}
V=Ψi=Uqi,
(1)

where U is the total Coulomb energy of the system.

The structure of the Coulomb energy expression is such that Eq. (1) is a system of linear equations for qi and can be solved with standard linear algebra techniques. To guarantee that the linear system corresponding to Eq. (1) has a solution, the electrode point charges are generally replaced with a narrow Gaussian charge distribution. For a detailed study of the optimal choice for Gaussian width, see Gingrich.23 

Several studies of EDLCs employing the CPM have been reported. Merlet et al.24–27 studied a nanoporous carbon electrode in contact with electrolyte consisting of an ionic liquid or an ionic liquid/acetonitrile mixture. Vatamanu and co-workers28–31 investigated ionic liquid electrolytes with carbon or gold electrodes using the smooth particle mesh Ewald (SMPE)32,33 method to simplify the calculation. The hydration of metal-electrode surfaces was examined by Limmer et al.34,35

In all of these studies, detailed comparisons of the CPM with the fixed charge method have been lacking. In a recent paper, Merlet et al.36 examined the differences between CPM and FCM simulations as measured by the relaxation kinetics in EDLC with nanoporous carbide-derived carbon electrode and the electrolyte structure at interface in EDLC with planar graphite electrode. It was showed that CPM predicted more reasonable relaxation time than FCM. In their study of electrolyte structure, there were some quantitative differences between the results of the two methods, but the qualitative features were unchanged for these ionic-liquid based EDLCs.

In this work, we study an organic electrolyte/salt-based EDLC, namely an LiClO4/acetonitrile electrolyte at a graphite electrode. To compare the results for this system using the CPM and FCM, several structural aspects normally reported in EDLC simulations are studied for comparison, specifically, the particle and charge density profiles near the electrodes and the solvation structure of the cation (Li+) both in bulk and near the surface.

In our simulations, the atoms of the electrolyte solution are placed between two carbon electrodes, each consisting of three graphite layers. The simulation geometry is shown in Fig. 1, which shows a snapshot from a simulation at 298 K with a potential difference, ΔΨ of 2 V. For the production runs, the distance between the two inner-most electrode layers (labeled L1 and R1, respectively, in Fig. 1) is 6.365 nm, which is far beyond the Debye length (0.2 nm). The electrolyte between the electrodes consists of 588 acetonitrile molecules and 32 Li+/

${\rm ClO}_4^-$
ClO 4 pairs, corresponding to a LiClO4 concentration of 1.00 M. The total dimensions of the cell are 2.951, 2.982, and 8.040 nm. The positions of the electrode carbon atoms are held fixed during the simulation.

FIG. 1.

Simulation snapshot at 298 K with ΔΨ = 2V: negative electrode is on left and positive is on right. The color of each electrode atom indicates its charge (refer to color scale bar, unit e). The electrolyte solution is shown between the two electrodes. Orange spheres: Li+; red spheres: O in

${\rm ClO}_4^-$
ClO 4⁠; cyan spheres: Cl in
${\rm ClO}_4^-$
ClO 4
; transparent stick models: acetonitrile. For clarity, in this figure, the distance between L1 and R1 (5.43 nm) is smaller than that used in the production runs.

FIG. 1.

Simulation snapshot at 298 K with ΔΨ = 2V: negative electrode is on left and positive is on right. The color of each electrode atom indicates its charge (refer to color scale bar, unit e). The electrolyte solution is shown between the two electrodes. Orange spheres: Li+; red spheres: O in

${\rm ClO}_4^-$
ClO 4⁠; cyan spheres: Cl in
${\rm ClO}_4^-$
ClO 4
; transparent stick models: acetonitrile. For clarity, in this figure, the distance between L1 and R1 (5.43 nm) is smaller than that used in the production runs.

Close modal

To model the molecular interactions we employ a variety of literature force fields. For acetonitrile, we use the united atom model of Edwards et al.37 For the ions, we use the forcefield of Eilmes and Kubisiak,38 excluding the polarizability terms. The interaction parameters for the graphite electrode carbon atoms are taken from Ref. 39. Lorentz-Berthelot mixing rules are used to construct all cross interactions.

In this slab geometry, we define the z-axis as the direction normal to the electrodes and apply periodic boundary conditions only in the x-y plane (parallel to the graphite layers). Unlike the original CPM, which used 2d-periodic Ewald sums,32,40 3d-periodic Ewald sums with shape corrections41 were used in this work to improve the calculation speed, with a volume factor set to 3. The correction term to the usual 3d-periodic energy expression is given by (see Eq. (A.21) in Ref. 23)

\begin{equation}\frac{2\pi }{V}\left(\sum _{i}q_{i}z_{i}\right)^{2} .\end{equation}
2πViqizi2.
(2)

In studies using FCM, the uniform charge on each electrode atom either is arbitrarily chosen15,16,18 or is estimated using a time-consuming trial and error procedure to yield the specified electric potential difference,17 which is calculated by numerically integrating the Poisson equation using the charge density profile. The CPM is able to predict the explicit average charge on each electrode atom at a given potential difference; therefore, for consistency, in our FCM simulations we set the charge on each electrode atom to the average charge per atom obtained by the CPM calculations at the same potential difference. Otherwise, all other force-field parameters in the FCM simulations are identical to those used in the CPM simulations.

All simulations were performed using the molecular-dynamics simulation code LAMMPS,42 modified to implement the CPM, using a time step of 1 fs. Constant NVT conditions are enforced using a Nosé-Hoover thermostat with a relaxation time of 100 fs and a temperature of 298 K. The cutoffs for all non-bonded interactions are 1.4 nm. For the Ewald sums, an accuracy (relative RMS error in per-atom forces) is set to 10−8. The parameter of the Gaussian electrode charge distribution (see Eq. (S2) in the supplementary material43) is set to 19.79 nm−1, which is the same as in Ref. 19. All results reported here are statistical averages taken from runs of 25–30 ns in length, each preceded by 2 ns of equilibration. Further details as to the CPM method and our implementation in LAMMPS can be found in the supplementary material.43 

Unlike the FCM, the CPM allows the individual atom charges on the electrode to fluctuate in response to local charge rearrangements in the electrolyte. The total charges on each layer for several potential differences are shown in Table I. The net electrode charge for each simulation does not show values statistically significant from zero, as expected. In a perfect conductor, the charges on the electrode atoms would be concentrated entirely on the electrode surface layers (L1 and R1 in Fig. 1); however, in the CPM the charge distribution on the electrode is approximated by discrete point charges centered on the electrode itself, so a small amount of charge is found on the second layer and to a much lesser extent the third.

Table I.

Average total charge on each electrode layer.a

Layer0 V2 V4 V5 V
L3 0.00(2) 0.03(2) 0.06(2) 0.06(2) 
L2 0.02(5) 0.33(5) 0.64(4) 0.75(5) 
L1 −0.1(4) −3.4(4) −6.6(4) −8.4(5) 
R1 0.0(4) 3.3(4) 6.5(4) 8.4(5) 
R2 0.01(5) −0.29(4) −0.58(5) −0.75(6) 
R3 0.00(2) −0.02(2) −0.03(2) −0.06(2) 
Total 0.01(6) 0.02(6) 0.04(6) 0.01(6) 
Layer0 V2 V4 V5 V
L3 0.00(2) 0.03(2) 0.06(2) 0.06(2) 
L2 0.02(5) 0.33(5) 0.64(4) 0.75(5) 
L1 −0.1(4) −3.4(4) −6.6(4) −8.4(5) 
R1 0.0(4) 3.3(4) 6.5(4) 8.4(5) 
R2 0.01(5) −0.29(4) −0.58(5) −0.75(6) 
R3 0.00(2) −0.02(2) −0.03(2) −0.06(2) 
Total 0.01(6) 0.02(6) 0.04(6) 0.01(6) 
a

The numbers in parentheses represent 95% confidence intervals in the last significant figure.

Figure 2 is a log-linear plot of the probability distribution of individual electrode atom charges, p(q), on the inner electrode layers for various potential differences. For the inner layer of the positively charged electrode (R1), the charge distribution is well described by a Gaussian distribution over the entire range of potential differences studied (ΔΨ = 0 V to 5 V). For the negatively charged electrode (L1), this is true at the lower potential differences (up to ΔΨ = 2 V), but the charge distribution takes on a bimodal structure at higher potential differences (ΔΨ = 4 V and 5 V). A similar non-Gaussian charge distribution has also been seen in simulations at the negative electrode of a H2O/Pt system;34 however, that system differs somewhat from our simulated system in that no ions are present. The bimodal charge distribution at the higher potential differences in this system highlights an important difference between the CPM and the FCM. Unlike the FCM, the electrode charges in the CPM can adjust to respond to local fluctuations in the electrolyte/ion charge density. As we will see in Sec. III B, this second peak in p(q) seen at high potential differences is due to the presence of Li+ ions near the electrode surface, which induce higher than average local charges on the adjacent electrode.

FIG. 2.

Distribution of electrode atom charges for the inner electrode layers at various potential differences: L1 (solid lines) and R1 (dashed lines).

FIG. 2.

Distribution of electrode atom charges for the inner electrode layers at various potential differences: L1 (solid lines) and R1 (dashed lines).

Close modal

To better understand the origin of the bimodal charge distributions in the CPM at high potential difference, we examine the ion density profiles at the interface and compare the results obtained using CPM and FCM (Fig. 3).

FIG. 3.

Ion densities near each of the two electrodes (L1 and R1) for various applied potential differences (a) L1; ΔΨ = 2 V, (b) R1; ΔΨ = 2 V, (c) L1; ΔΨ = 4 V (d) R1; ΔΨ = 4 V (e) L1; ΔΨ = 5V, (f) R1; ΔΨ = 5 V. The negatively and positively charged electrodes (L1 and R1) are at z = 0 and 6.365 nm, respectively.

FIG. 3.

Ion densities near each of the two electrodes (L1 and R1) for various applied potential differences (a) L1; ΔΨ = 2 V, (b) R1; ΔΨ = 2 V, (c) L1; ΔΨ = 4 V (d) R1; ΔΨ = 4 V (e) L1; ΔΨ = 5V, (f) R1; ΔΨ = 5 V. The negatively and positively charged electrodes (L1 and R1) are at z = 0 and 6.365 nm, respectively.

Close modal

At the lowest non-zero potential difference (ΔΨ = 2 V), the ion density profiles for the CPM and FCM methods show relatively minor quantitative differences in peak heights, but otherwise the density peak positions and overall structure are identical. At higher potential differences (4 V and 5 V), however, qualitative differences emerge. At ΔΨ = 4 V, a peak in the Li+ density profile near the negatively charged electrode (L1) at 0.22 nm appears in the CPM calculation, but is absent in the FCM results. This peak represents inner-sphere adsorbed Li+ ions that no longer possess a full solvation shell of acetonitrile, but are also “partially solvated” by the electrode. The existence of these inner-sphere adsorbed ions is made possible because the CPM allows for more physically correct electrode charge distribution - one that can respond to the presence of the nearby Li+ ion with locally larger-than-average negative electrode charges. The partial desolvation was reported to be highly related with the confinement of ions, which fundamentally affects the capacitance.44 

As the potential difference is increased beyond 4 V, the inner-sphere adsorbed ion peak grows substantially. For the FCM, we see this peak at 5 V, but it has an amplitude that is much smaller than that predicted by the CPM. In addition, Fig. 3 also shows that the height of the Li+ density peak at 0.69 nm, representing the first fully solvated outer-sphere adsorbed Li+ layer in the EDLC, is overestimated at ΔΨ = 4 V and 5 V in the FCM. This is due to the migration of Li+ ion density from the inner-sphere adsorbed peak at 0.21 nm, compared with the CPM.

In contrast to Li+, the

${\rm ClO}_4^-$
ClO 4 density profiles are nearly identical for both methods at all studied potential differences. The charge density in
${\rm ClO}_4^-$
ClO 4
is spread out over a larger radius and does not create as large a local charge concentration in the electrolyte solution to perturb the electrode charge distribution enough to affect the density profiles.

In Fig. 4 we plot the solvent (acetonitrile, center of mass) and charge density profiles. Unlike the ion density profiles, the density profiles for acetonitrile are identical for both methods up to a potential difference of 4V. However, at ΔΨ = 5V the acetonitrile peak closest to the negative electrode (L1) splits into two peaks in the CPM, a feature not seen in the FCM calculation. This is because the large amount of electrode solvated Li+ seen in the CPM pulls the acetonitrile molecules in its solvation shell closer to the surface.

FIG. 4.

Acetonitirle density and charge profiles at various applied potentials (a) ΔΨ = 2 V, (b) ΔΨ = 4 V, (c) ΔΨ = 5 V. Vertical lines show the position of the electrode inner layers L1 (left) and R1 (right).

FIG. 4.

Acetonitirle density and charge profiles at various applied potentials (a) ΔΨ = 2 V, (b) ΔΨ = 4 V, (c) ΔΨ = 5 V. Vertical lines show the position of the electrode inner layers L1 (left) and R1 (right).

Close modal

The structure of the electric double layer in this system is best illustrated by the charge density profile (see Figure 4). For this quantity the results for both methods (CPM and FCM) are very similar with the only difference coming at high potential difference, where small charge peaks corresponding to the electrode-solvated Li+ are present in the CPM result, but not in the FCM.

In this work, the principal difference between the results derived from the CPM and FCM is the appearance in the CPM of Li+ ions very close to the negatively charged electrode at high potentials. This inner-sphere adsorbed Li+ also is seen in the FCM, but only at very low concentration at the highest potential difference studied. To examine this feature more closely, we analyze the solvation structure of Li+ by plotting in Fig. 5 (for ΔΨ = 4 V) the CPM radial distribution function (RDF) between the Li+ ion and the acetonitrile N for three Li+ distances from the electrode surface: 0.21, 0.69, and 4.50 nm. These correspond to the inner-sphere adsorbed Li+, the first layer of outer-sphere adsorbed Li+ and the bulk conditions, respectively. Also plotted in Fig. 5 are the number of solvent molecules (NS) within a given radius—for the bulk ions, the first plateau in this quantity corresponds to the equilibrium coordination number. The results of RDF and NS for the FCM or other potential differences are identical to that shown for the CPM at 4 V (Fig. 5)—except for the fact that no Li+ ions at 0.21 nm for the CPM at ΔΨ = 2 V and for the FCM at 2 V and 4 V are detected and therefore no RDF (nor NS) data was obtained.

FIG. 5.

Radial distribution functions (RDF, solid lines) and corresponding number of solvent molecules (NS, dashed lines) between Li+ and N atom in acetonitrile for ΔΨ = 4 V. In the legend, the number indicates the distance (nm) from the negative electrode surface to the center of bin (with width ±0.02 nm) containing the Li+ ion. The insets show typical snapshots of the Li+ solvation structure for each distance studied.

FIG. 5.

Radial distribution functions (RDF, solid lines) and corresponding number of solvent molecules (NS, dashed lines) between Li+ and N atom in acetonitrile for ΔΨ = 4 V. In the legend, the number indicates the distance (nm) from the negative electrode surface to the center of bin (with width ±0.02 nm) containing the Li+ ion. The insets show typical snapshots of the Li+ solvation structure for each distance studied.

Close modal

At z = 0.69 nm and 4.50 nm, the Li+ ion is fully solvated by the acetonitrile solvent and the first solvation shell peaks for both distances are nearly identical with a peak distance of 0.22 nm. In addition, the coordination number (given by the plateau value of NS) of each Li+ ion is equal to about 5.0 for the “bulk” case (4.50 nm), just slightly smaller than that for the first fully solvated peak with 5.4, presumably due to the enhanced acetonitrile density near the electrode region (see Fig. 4). The value for the “bulk” is consistent with that seen in other simulation studies of bulk acetonitrile-lithium salt solutions.45,46

For Li+ ions at 0.21 nm from the electrode, however, the solvent coordination is significantly reduced from the bulk value to 3.1 due to partial solvation by the electrode atoms themselves. The bottom inset in Fig. 5 shows a representative snapshot from the CPM simulation (ΔΨ = 4 V) of an electrode-solvated ion and its nearby environment. In the CPM, because the charges on the electrode can fluctuate individually in response to the local electrolyte charge distribution, the presence of the Li+ ion very close to the electrode induces larger than average negative charges on the nearby electrode ions, as seen in the inset. The ability of the electrode to respond to local fluctuations is necessary for the stabilization of the electrode-solvated Li+ below ΔΨ ≲ 4 V, as no such ions are seen in the FCM simulations in this range. At 5 V, such electrode-solvated ions are seen in the FCM, but at significantly lower concentration than in the CPM, indicating a much higher energy for such configurations in the FCM, relative to the CPM. We have confirmed that sampling in the simulations is sufficient, as we see multiple crossings and recrossings of Li+ ions between the positions at 0.69 nm and 0.21 nm.

In this work, two methods, the CPM and the FCM, are compared in a simulation of a model for a LiClO4-acetonitrile/graphite electric double-layer capacitor. The major difference between the two methods is that in the CPM the charges on the individual electrode atoms can fluctuate in response to local fluctuations in electrolyte charge density, whereas those charges for the FCM are static. For this system, there are no measurable differences between the results of these two methods at low potential differences (ΔΨ ⩽ 2 V); however, at larger potential differences significant qualitative differences emerge in the EDLC ion spatial distribution.

At a potential difference of 4 V, a new peak in the Li+ density profile appears in the CPM calculation at 0.21 nm away from the electrode. This peak — absent in the FCM calculation — corresponds to a inner-sphere adsorbed Li+ ion that is close enough to the electrode to have lost some of its acetonitrile solvation shell and is partially solvated by the electrode atoms. For this electrode-solvated Li+ ion the acetonitrile coordination number is reduced from about 5.0 for a fully solvated ion to 3.1. This partial solvation is possible because in the CPM the electrode atom charges can respond to local fluctuations in the electrolyte charge density—in this case negative charges build up in the CPM electrode near the lithium ion. This ability lowers the energy of the “electrode-solvated” Li+ ion relative to that of a fixed-charge electrode (as in the FCM). This close approach of Li+ to the electrode is possible in the FCM, but only at higher electrode potential differences and a considerably lower concentrations than when CPM is used. The energetics of the approach of a Li+ ion to the electrode is important in many pseudo-capacitor applications, and, as our calculations show, the CPM would be preferable to the FCM in the calculation of the barrier energy of this process because it more accurately represents the fluctuating charges on the electrode.

The authors thank Dr. David Limmer for helpful discussions on the constant potential method. This work is supported as part of the Molecularly Engineered Energy Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001342.

1.
N. S.
Choi
,
Z.
Chen
,
S. A.
Freunberger
,
X.
Ji
,
Y.-K.
Sun
,
K.
Amine
,
G.
Yushin
,
L. F.
Nazar
,
J.
Cho
, and
P. G.
Bruce
,
Angew. Chem.
51
,
9994
(
2012
).
2.
V.
Aravindan
,
J.
Gnanaraj
,
Y. S.
Lee
, and
S.
Madhavi
, “
Insertion-type electrodes for nonaqueous Li-ion capacitors
,”
Chem. Rev.
(published online).
3.
T.
Sato
,
G.
Masuda
, and
K.
Takagi
,
Electrochim. Acta
49
,
3603
(
2004
).
4.
K.
Yuyama
,
G.
Masuda
,
H.
Yoshida
, and
T.
Sato
,
J. Power Sources
162
,
1401
(
2006
).
5.
D.
Wei
,
M. R.
Astley
,
N.
Harris
,
R.
White
,
T.
Ryhänen
, and
J.
Kivioja
,
Nanoscale
6
,
9536
(
2014
).
6.
H.
Ji
,
X.
Zhao
,
Z.
Qiao
,
J.
Jung
,
Y.
Zhu
,
Y.
Lu
,
L. L.
Zhang
,
A. H.
MacDonald
, and
R. S.
Ruoff
,
Nat. Commun.
5
,
3317
(
2014
).
7.
C.
Largeot
,
C.
Portet
,
J.
Chmiola
,
P. L.
Taberna
,
Y.
Gogotsi
, and
P.
Simon
,
J. Am. Chem. Soc.
130
,
2730
(
2008
).
8.
H.
Itoi
,
H.
Nishihara
,
T.
Kogure
, and
T.
Kyotani
,
J. Am. Chem. Soc.
133
,
1165
(
2011
).
9.
M. V.
Fedorov
and
A. A.
Kornyshev
,
Chem. Rev.
114
,
2978
(
2014
).
10.
H.
Wang
and
L.
Pilon
,
Electrochim. Acta
64
,
130
(
2012
).
11.
H.
Wang
,
A.
Thiele
, and
L.
Pilon
,
J. Phys. Chem C
117
,
18286
(
2013
).
12.
D. E.
Jiang
and
J.
Wu
,
J. Phys. Chem. Lett.
4
,
1260
(
2013
).
13.
Y.
Ando
,
Y.
Gohda
, and
S.
Tsuneyuki
,
Chem. Phys. Lett.
556
,
9
(
2013
).
14.
L.
Yang
,
B. H.
Fishbine
,
A.
Migliori
, and
L. R.
Pratt
,
J. Am. Chem. Soc.
131
,
12373
(
2009
).
15.
L.
Yang
,
B. H.
Fishbine
,
A.
Migliori
, and
L. R.
Pratt
,
J. Chem. Phys.
132
,
044701
(
2010
).
16.
Y.
Shim
,
H. J.
Kim
, and
Y.
Jung
,
Faraday Discuss.
154
,
249
(
2012
).
17.
G.
Feng
and
P. T.
Cummings
,
J. Phys. Chem. Lett.
2
,
2859
(
2011
).
18.
G.
Feng
,
S.
Li
,
J. S.
Atchison
,
V.
Presser
, and
P. T.
Cummings
,
J. Phys. Chem. C
117
,
9178
(
2013
).
19.
S. K.
Reed
,
O. J.
Lanning
, and
P. A.
Madden
,
J. Chem. Phys.
126
,
084704
(
2007
).
20.
J. I.
Siepmann
and
M.
Sprik
,
J. Chem. Phys.
102
,
511
(
1995
).
21.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
22.
T. R.
Gingrich
and
M.
Wilson
,
Chem. Phys. Lett.
500
,
178
(
2010
).
23.
T. R.
Gingrich
, “
Simulating surface charge effects in carbon nanotube templated ionic crystal growth
,” M.S. thesis (
University of Oxford
,
2010
).
24.
C.
Merlet
,
M.
Salanne
,
B.
Rotenberg
, and
P. A.
Madden
,
J. Phys. Chem. C
115
,
16613
(
2011
).
25.
C.
Merlet
,
B.
Rotenberg
,
P. A.
Madden
,
P. L.
Taberna
,
P.
Simon
,
Y.
Gogotsi
, and
M.
Salanne
,
Nat. Mater.
11
,
306
(
2012
).
26.
C.
Merlet
,
M.
Salanne
, and
B.
Rotenberg
,
J. Phys. Chem. C
116
,
7687
(
2012
).
27.
C.
Merlet
,
M.
Salanne
,
B.
Rotenberg
, and
P. A.
Madden
,
Electrochim. Acta
101
,
262
(
2013
).
28.
J.
Vatamanu
,
O.
Borodin
, and
G. D.
Smith
,
J. Am. Chem. Soc.
132
,
14825
(
2010
).
29.
J.
Vatamanu
,
O.
Borodin
,
D.
Bedrov
, and
G. D.
Smith
,
J. Phys. Chem. C
116
,
7940
(
2012
).
30.
J.
Vatamanu
,
O.
Borodin
, and
G. D.
Smith
,
J. Phys. Chem C
116
,
1114
(
2012
).
31.
J.
Vatamanu
,
Z.
Hu
,
D.
Bedrov
,
C.
Perez
, and
Y.
Gogotsi
,
J. Phys. Chem. Lett.
4
,
2829
(
2013
).
32.
M.
Kawata
and
M.
Mikami
,
Chem. Phys. Lett.
340
,
157
(
2001
).
33.
M.
Kawata
,
M.
Mikami
, and
U.
Nagashima
,
J. Chem. Phys.
116
,
3430
(
2002
).
34.
D. T.
Limmer
,
C.
Merlet
,
M.
Salanne
,
D.
Chandler
,
P. A.
Madden
,
R.
van Roij
, and
B.
Rotenberg
,
Phys. Rev. Lett.
111
,
106102
(
2013
).
35.
D. T.
Limmer
,
A. P.
Willard
,
P.
Madden
, and
D.
Chandler
,
Proc. Natl. Acad. Sci. U.S.A.
110
,
4200
(
2013
).
36.
C.
Merlet
,
C.
Péan
,
B.
Rotenberg
,
P. A.
Madden
,
P.
Simon
, and
M.
Salanne
,
J. Phys. Chem. Lett.
4
,
264
(
2013
).
37.
D. M. F.
Edwards
,
P. A.
Madden
, and
I. R.
Mcdonald
,
Mol. Phys.
51
,
1141
(
1984
).
38.
A.
Eilmes
and
P.
Kubisiak
,
J. Phys. Chem. B
115
,
14938
(
2011
).
39.
M. W.
Cole
and
J. R.
Klein
,
Surf. Sci.
124
,
547
(
1983
).
40.
S. W.
De Leeuw
and
J. W.
Perram
,
Mol. Phys.
37
,
1313
(
1979
).
41.
I. C.
Yeh
and
M. L.
Berkowitz
,
J. Chem. Phys.
111
,
3155
(
1999
).
42.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
43.
See supplementary material at http://dx.doi.org/10.1063/1.4899176 for details of constant potential method and its implementation in LAMMPS.
44.
C.
Merlet
,
C.
Péan
,
B.
Rotenberg
,
P. A.
Madden
,
B.
Daffos
,
P. L.
Taberna
,
P.
Simon
, and
M.
Salanne
,
Nat. Commun.
4
,
2701
(
2013
).
45.
S.
Oliveira Costa
and
J.
López Cascales
,
J. Electroanal. Chem.
644
,
13
(
2010
).
46.
D.
Spångberg
and
K.
Hermansson
,
Chem. Phys.
300
,
165
(
2004
).

Supplementary Material