# Constant Chemical Potential-Quantum Mechanical-Molecular Dynamics
simulations of the Graphene-electrolyte double layer

Nicodemo Di Pasquale Corresponding author<EMAIL_ADDRESS>Department of Chemical Engineering, Brunel University London, Uxbridge, UB8
3PH, United Kingdom  Aaron R. Finney  Department of Chemical Engineering,
University College London, London, WC1E 7JE, United Kingdom  Joshua Elliott
Department of Chemical Engineering, University of Manchester, Manchester, M13
9PL, United Kingdom  Diamond Light Source, Harwell Science and Innovation
Park, Didcot, Oxfordshire OX11 8UQ, United Kingdom  Paola Carbone  Department
of Chemical Engineering, University of Manchester, Manchester, M13 9PL, United
Kingdom  Matteo Salvalaglio  Department of Chemical Engineering, University
College London, London, WC1E 7JE, United Kingdom

(October 2021)

###### Abstract

We present the coupling of two frameworks—the pseudo-open boundary simulation
method known as constant potential Molecular Dynamics simulations (C$\mu$MD),
combined with QMMD calculations—to describe the properties of graphene
electrodes in contact with electrolytes.

The resulting C$\mu$QMMD model was then applied to three ionic solutions
(LiCl, NaCl and KCl in water) at bulk solution concentrations ranging from 0.5
M up to 6 M in contact with a charged graphene electrode. The new approach we
are describing here provides a simulation protocol to control the
concentration of the electrolyte solutions while including the effects of a
fully polarizable electrode surface. Thanks to this coupling, we are able to
accurately model both the electrode and solution side of the double layer and
provide a thorough analysis of the properties of electrolytes at charged
interfaces, such as the screening ability of the electrolyte and the
electrostatic potential profile. We also report the calculation of the
integral electrochemical double layer capacitance in the whole range of
concentrations analysed for each ionic species, while the QM simulations
provide access to the differential and integral quantum capacitance. We
highlight how subtle features, such as the adsorption of potassium at the
interface or the tendency of the ions to form clusters, emerge from our
simulations, contribute to explaining the ability of graphene to store charge
and suggest implications for desalination.

##  1 Introduction

Interest in graphene-based devices has grown in recent years, thanks of the
versatility and physical characteristics of this new material, in particular
for applications in which it is in contact with an electrolyte solution. Use
of nanoporous graphene as a membrane for water desalination [1, 2] is one
important example. The presence of pores of equal size to the electrolytes
allows the selective passage of water through the membrane. Combined with the
atomic scale thickness of graphene, this can lead to the creation of
desalination membranes with higher performances than common polymer-based ones
[3]. Another promising technologically relevant applications is the use of
graphene electrodes in electrochemical double layer (super)capacitor (EDLC)
devices[4, 5, 6]. In fact, graphene [7, 8, 9, 10], porous activated carbon
[11] and carbon nanotube [12, 13] electrodes potentially have relatively high
charge storage capacity and a favourable specific energy to power ratio, due
to rapid charge-discharge cycling [8] controlled by changes of an applied
potential, together with lifetimes that can reach millions of cycles [11].

Typically, charge storage at carbonaceous electrodes is a non-faradaic
process, where mobile ionic species accumulate at the interface between the
electrode and the liquid phase. An important class of systems of this kind,
which has gained lots of attention recently, is represented by cheap and easy-
to-prepare aqueous-based electrolytes in contact with a graphene electrode
[6]. Carbon-based EDLCs with aqueous-based electrolytes do not generally
suffer from electrochemical degradation, can be non-toxic, and provide an
attractive alternative solution at the problem of energy storage compared with
traditional battery devices. Combined with a longer lifetime and high power
density,[14] these energy storage systems could be increasingly applied to
power small electronic devices and for acceleration and breaking in electrical
vehicles [5].

Several experimental works were undertaken to understand the physicochemical
properties of neutral and charged graphene interfaces in contact with
electrolyte solutions and the nature of these systems charge storage capacity
[15, 16, 17]. However, the delicate balance between hydration-free energy and
surface effects, which regulate the physisorption of ionic species at
surfaces, resulted in conflicting experimental findings (see [18] for a more
detailed account). For instance, there are reports both supporting the
conclusion that the capacitance of graphene films is ion-independent [16], as
well as contrasting observations suggesting that basal capacitance is instead
ion-specific (with, for example, a greater propensity for Na+ and K+
adsorption over Li+ adsorption at negatively charged electrodes in the case of
group I cations)[17]. Atomic-scale defects in the graphitic surface, its
topography, dimensionality and chemical modifications are difficult to control
and have non-negligible effects in experimental measurements. As an example,
mechanical cutting produces structural defects known as “dangling bonds” which
modifies the measured capacitance of the sample [15, 19]. In this respect, a
model of the graphene interface and its interactions with an electrolyte
solution can exclude all the spurious effects coming from uncontrolled defects
and chemical modification of the surface. Molecular modelling and simulations
can help to improve understanding of the mechanisms involved in such complex
systems and guide the interpretation of experimental results.

Many key features of supercapacitive devices are underpinned by the properties
of the electrochemical double layer, and their responses to the charging of
the electrode. Gouy-Chapman theory [20, 21] describes the double layer as a
diffuse charged layer in the solution that compensates an applied surface
charge on the electrode. Modifications to this model include the adsorption of
counter-ions at the surface in the so-called Stern layer [22]. The development
of a mean-field theory based on the Poisson-Boltzmann lattice-gas model [23]
has shown that features not present in the Gouy-Chapman theory, such as steric
effects, ion correlations, and preferential adsorption [24, 25, 26] need to be
accounted for in order to correctly describe the interactions between the ions
and the electrode. Mechanistic insight for these kinds of effects and how they
control charge storage can be gained by atomistic simulations of the
graphene/electrolyte interface; these also enable the evaluation of ensemble
properties, such as the free energy of adsorption of the ions at the interface
[27]. Furthermore, simulations can establish the effect of solution
concentration on ion accumulation at the electrode, their interfacial
structure, and dynamical properties.

In order to compare simulations with a macroscopic system, this adsorption
should ideally be modelled in the presence of bulk electroneutral solution
with fixed composition to ensure a constant driving force for the adsorption
at a charged surface. This can be obtained for example as shown in Finney et
al. [28], where the authors performed MD simulations using constant chemical
potential MD simulations, C$\mu$MD [29], which mimics open-boundary
conditions. With C$\mu$MD, the authors simulated NaCl(aq) with concentrations
spanning $\sim 0.1-10$ M at graphite surfaces. Their results indicate that the
interface charge screening behaviour is a function of bulk solution
concentration, with a transition (at $\sim 1$M) from diffuse charge screening,
qualitatively consistent with the picture from simple mean field models, to a
complex multi-layered structuring that systematically either over or under
screens the surface potential. The multiple charged layers result from ion
finite-size effects, over-compensation of the surface charge by oppositely
charged ions closest to the surface, and non-idealities in solution, i.e.,
when the hypothesis of non-interaction between oppositely charged ions breaks
down for large ions concentrations [30]. This last effect also has
consequences on the conductance of the ions, which deviates from the
prediction of the Nerst-Einstein equations [31].

Together with a constant driving force for ion adsorption from the bulk,
another important effect to consider in the description of such systems is the
the polarisation of the electrode exerted by the adsorbing electrolytes [18].
Classical simulations typically model the non-bonded interactions between
atoms within the electrolyte and atoms belonging to the interface using
additive pairwise potentials such as the Lennard-Jones potential and Coulomb
interactions between fixed point atom charges. Polarisation can be introduced
using e.g., oscillating charge models, or by fitting short-range potentials to
binding energies obtained from ab initio methods [27, 32, 33]. However, these
models may not accurately capture the complex many-body effect associated with
charge polarization at the electrode-solution interface. Another way to
include polarisation in classical MD simulations is through the constant
potential method developed in [34]. This constant potential method has been
successfully deployed to describe the properties of the electrochemical double
layer of aqueous electrolytes and ionic liquids in contact with metal
electrodes such as Au and Cu. Despite its successes, one of the key
approximations of the constant potential method is that the electrode is fully
metallic and can perfectly screen charges, which is not the case for
(semimetallic) graphene [18].

On the other hand, a full Quantum Mechanical (QM) treatment of the
interactions between the electrolyte and the substrate is still unfeasible,
due to the length (tens of nm) and time (hundreds of ns) scales required for
modeling the effect of the aqueous electrolytes. However, while the full QM
model of the electrode/electrolyte system is out of reach, QM calculations can
be used to compute a set of atomic partial charges on the electrode in the
presence of the electrostatic potential arising from the position of the
electrolyte atoms. This is exactly the spirit of our QMMD scheme, where QM
calculations are coupled to MD simulations at fixed intervals of time
integration. As such, the surface atom partial charges within the classical
force field are updated on the fly. In a more recent development Machine
Learning models have recently proven to be a viable option in tuning the
surface polarization if the scope of the system becomes too large for QM
simulations. This is achieved by replacing the QM calculations with a Neural
Network (NN) model trained to reproduce results from a wide range of QM
calculations with varying distributions of electrolytes in solution. The NN
acts as a polarizable-like force field, combining fast classical MD
simulations with more accurate QM calculations of the interface polarization
[35].

This present work leverages the QMMD framework introduced in [18] and the
C$\mu$MD introduced in [29, 28]. The approach simultaneously captures surface
polarization and concentration effects that can modify the structure and
composition of the electrochemical double layer. We use the resulting
C$\mu$QMMD protocol to examine interfaces between aqueous alkali chloride
solutions at different concentrations with a graphene electrode surface,
elucidating complex interfacial structure, dynamics, and electrochemical
properties.

This paper is organized as follows: we first provide a brief overview of the
QMMD and C$\mu$MD protocols, pointing to the relevant literature for the
interested reader; we present the systems to which we apply the C$\mu$QMMD
framework: a charged graphene electrode in contact with three different
electrolyte solution, NaCl(aq), LiCl(aq), KCl(aq) at different concentrations.
We derive the electrical properties of the interface in terms of the screening
factor and electrical potential and calculate the total integral capacitance
of this system by deriving the quantum and electrical double layer
capacitance. Finally, we discuss the effects of complex solute speciation on
the performance of graphene-electrolyte devices and draw some conclusions
regarding this new proposed simulation scheme.

##  2 Computational Models

In order to capture the dynamic polarization of a charged graphene surface in
response to the evolving configuration of an electrolyte at a prescribed
concentration, we coupled the classical C$\mu$MD simulation to the electronic
structure theory calculations at regular time intervals. We will give a more
detailed account of both models (C$\mu$MD and QMMD) in the following sections,
while here we will only discuss their coupling.

A sketch of the sequence of the operations involved is given in Figure 1. All
the operations shown in Figure 1 are obtained through an in-house python
wrapper. During the MD time integration obtained with GROMACS 2018.4 MD
package [36], ion positions are passed to the Plumed software (v. 2.7) [37]
patched with GROMACS, to compute the C$\mu$MD forces (see section 2.1 for more
details). After the evolution of the atom positions, the final configuration
of the electrolyte is extracted to compute the electrostatic potential. In
turn, this latter quantity is used as input for the QM calculations obtained
with the Dftb+ software package [38]. From the QM results, the distribution of
the charges on the graphene is extracted (see section 2.2 for more details)
and used as input for the new iteration of the loop.

Figure 1: The computational workflow adopted in this work highlighting the two
“black boxes” (the MD software and the QM software) in the blue squares and
the operations included in the python wrapper (red squares).

###  2.1 C$\mu$MD Model

Figure 2: Example configuration from a C$\mu$QMMD simulation of KCl(aq) in
contact with graphene in this work projected onto simulation $x,z$ dimensions.
K+, Cl-, O of water and C of graphene are shown by the pink, cyan, red and
grey spheres. The blue lines highlight the C$\mu$QMMD control and reservoir
regions, which also indicate the simulation cell boundaries. An extended
vacuum region, around 8 nm in $z$, is truncated in the image.

The graphene electrode we considered is located at $z=0$ and is in contact
with an electrolyte slab of thickness 8 nm. A further 8 nm of vacuum separates
the system from its periodically repeating images. The electrolyte phase is
divided into three regions: the first region starts at the graphene electrode
up to a distance of 4 nm. The second one is the control region, which is used
to control the concentrations. The third region is the reservoir region which
provides the reservoir of ions to adjust the concentration of the electrolytes
in the other regions. Figure 2 provides an example of the set-up adopted in
this work, where we highlighted the different C$\mu$MD simulation cell
regions.

The control of the concentration of the ions in solution is obtained by
applying a force at the edge of the reservoir region according to a continuous
function of the form,

$F_{i}^{\mu}(z)=k_{i}(n_{i}^{\mathrm{CR}}-n_{i}^{0})\left[\frac{1}{4\omega}\left(1+\mathrm{cosh}\left(\frac{z-z_{F}}{\omega}\right)\right)^{-1}\right].$
(1)

Here, $\omega$ was set to 0.2 nm, and represents the width of the force region
(between the control and reservoir regions highlighted by the blue lines in
Figure 2) while $k$ was $2\times 10^{4}$ kJ mol-1 nm-1, giving the correct
densities in the bulk (see [28] for a discussion on these parameters). $n^{0}$
is the target ion number density, while $n^{\mathrm{CR}}$ is the density
calculated instantaneously during time integration in the control region.
Finally, $z_{F}$ is the position in $z$ where the C$\mu$MD forces are applied.
In our simulations this is set to 5.5 nm beyond the graphene surface. Using
this approach, the densities of cations and anions are constrained in the
control region to maintain target concentrations of 0.5, 2.0, 3.0, 4.0, 4.4
and 6 M. At each MD time-step ion positions are passed to Plumed in order to
compute the C$\mu$MD forces only acting on those ions in the region of
$z_{F}$. No external forces are applied to the ions outside of this region,
and any local change in the ion density at the interface results from the
physical interactions between graphene and the solution.

###  2.2 QMMD Model

The generality of electronic structure theory and its ability to reproduce the
electronic charge density distribution in semiconductors, metals, and
semimetals implies that that the QMMD approach can describe both long- and
short-ranged redistribution of the surface charge induced by the presence of
the electrolyte. Within each iteration (see Figure 1) of our scheme, the fully
classical system is taken as input for a quantum mechanical calculation. The
simulation box is partitioned into surface atoms whose electronic structure is
explicitly treated, and electrolyte atoms that are converted into a set of
point charges. The point charges take the values of the partial charges
contained in the classical force field and form the background electrostatic
potential during the computation of the electron structure. Upon derivation of
the electronic structure, partitioning of the charge density via Mulliken
population analysis yields the set of surface atom partial charges, which are
then passed to the classical force field. Finally, a short MD trajectory on
the order of several picoseconds can then be carried out (in the presence of
the quantum mechanically polarized surface) to generate the electrolyte
configuration for the following iteration. In our simulations we employ a
coupling between QM and MD calculations of 5 ps. We previously found for this
class of systems that 5 ps represents a good compromise in terms of
computational accuracy of the computed charges (0.004 $e$) vs computing time
when compared with a QMMD simulation where the charges were updated at every
MD time step [18].

In practice, in order to describe the electronic structure of solid-
electrolyte interfaces on the length scales required, we leverage the self-
consistent charge Density Functional Tight-Binding (SCC-DFTB) [39] approach,
which is an approximation to Kohn-Sham Density Functional Theory. The
empirical description in our Dftb+ calculations of the interactions between
the C atoms in the surface are described by the mio-1-1 parameter set. The SCC
charge threshold and Fermi temperature have been set to $1\times 10^{-2}$
Hartree and 300 K, respectively. Whereas, on first inspection, these criteria
can be considered loose and should not be adopted for the calculation of the
total electronic energy, rigorous testing in our previous works [18, 40] found
that they provide a sufficiently accurate description of the surface charge
distribution with respect to fully converged simulations, at a fraction of the
computational cost. Finally, to compute the partial charges passed to the
graphene force field at each MD step, we perform a Mulliken population
analysis [41], which gives reasonable results for this class of systems[18,
40].

###  2.3 Simulations Details

In our simulations, we consider a graphene electrode composed of 336 carbon
atoms in contact with aqueous electrolyte solutions. We investigated three
electrolyte systems, NaCl, KCl and LiCl at concentrations ranging from 0.5 M
to 6 M. However, due to the solubility limits of the KCl(aq) [42, 43], we
limit the investigated concentrations to 4.4 M for the KCl system.

Our simulations are carried out at constant surface charge, which makes it
difficult to draw comparisons across different electrodes since the potential
applied is not necessarily constant. As such, when we compute the capacitance,
we use the potential drop of the neutral electrode as a reference. This
approach has been applied previously to compare the properties of the
electrochemical double layer for different electrolytes [44]. Each operating
condition was therefore repeated for two different total charges of the
electrode: a charged graphene layer with a constant charge on the surface [45]
$\sigma$ of -0.449 e nm-2 (-0.0719 C/$m^{2}$) and a neutral one ($\sigma=0$).

Structural analyses of the solutions are carried out using PLUMED by post-
processing the simulation trajectories. The first-shell coordination numbers
for cations with anions ($N_{\mathrm{X-Cl}}$) and cations with water oxygen
atoms ($N_{\mathrm{X-Ow}}$) were computed using a continuous switching
function:

$N=\frac{1}{M}\sum_{i}^{M}\mathbf{e}^{\left(\frac{-(r-d_{0})^{2}}{2r_{0}^{2}}\right)}$
(2)

where $r$ are distances between atoms, $r_{0}$ was 0.01 nm, and $d_{0}$ was
chosen such that the function goes smoothly from one to zero at the position
of the first minimum in radial distribution functions for the cations with
anions and water oxygen atoms. This ensured that a conservative definition of
first-shell coordination was adopted in the analyses. Coordination numbers
were evaluated in 1.3 nm regions in $z$ closest to the graphene surface and
3.5 nm from the surface, representing the double layer and bulk solution
regions, respectively. The first coordination sphere distributions for ions
were used to construct a graph of ion-ion contacts using the NetworkX Python
library [46]. This allowed us to identify and compute the size of the ion
clusters formed. Ion clusters at the interface and within the bulk were
identified by sampling the regions defined for computing the coordination
numbers.

Molecular dynamics calculations in the NVT ensemble are carried out using
GROMACS [47, 48], version 2018.4. The leapfrog algorithm with a timestep of 1
fs was used to integrate the equations of motion at a constant temperature of
298.15 K, controlled with the Nosé-Hoover thermostat, with a relaxation time
of 0.1 ps. Long-range electrostatic interactions were treated using the
particle-mesh Ewald approach, with a cut-off of 1.4 nm. Non-bonded
interactions were computed using a Lennard-Jones 12-6 potential, truncated
smoothly at 1.0 nm using a switch function starting at a distance of 0.99 nm.
In all simulations, graphene carbon atoms were frozen, and water was modelled
using the SPC/E model [49] with the SETTLE algorithm used to maintain rigid
molecule geometries [50]. This choice is compatible with the Werder water-
graphene parameters that reproduce the experimentally measured graphene/water
contact angle [7, 51]. Ion force field parameters (for K+, Li+, Na+, Cl+),
also compatible with the SPC/E model, are taken from the work of Joung and
Cheatham [52]. In order to prevent water molecules and ions from escaping the
solution into the vacuum space, we added a fixed wall above the reservoir,
interacting with water molecules and ions only through a short-range Lennard-
Jones potential.

We equilibrated each system for 20 ns followed by 130 ns production runs to
collect data for subsequent analyses of the steady-state structure of the
interface. In all analyses discussed below, mean values and standard
deviations (error bars) are obtained via averaging performed using 5 ns
windows.

##  3 Results and Discussion

Thanks to the simulation protocol implemented, electroneutral solutions with
fixed ion concentrations can be maintained in the Control Region in Figure 2,
representing bulk solutions in equilibrium with the electrode-solution
interfaces. This allows us to compare the behaviour of different electrolytes
while controlling the electrolyte background concentration.

###  3.1 Density Profiles

We start this section by reporting in Figure 3 the concentration of the
different ionic species in solution as a function of the $z$-coordinate,
corresponding to the simulation cell direction orthogonal to the surface of
negatively charged graphene electrodes. As expected, these profiles show
preferential adsorption of cations at the electrode surfaces. For Na+ and Li+,
a sharp density peak is observed at a distance of $0.5$ nm from graphene,
followed by a second, less pronounced peak at $0.75$ nm. At the highest
concentrations, a third cation peak emerges around $1.15$ nm, which is more
pronounced for Li+. In contrast, in the case of K+, a small peak at $0.3$ nm
is followed by a much larger and relatively diffuse density peak at $0.6$ nm.
This is due to specific adsorption of the larger cation at the carbon surface,
a small number of which partially dehydrate to directly coordinate to carbon.

The difference in the $z$-density profiles for the different systems is less
notable when considering Cl- with respect to cations. At the lowest bulk
concentrations, there is a monotonically increasing density which reaches bulk
values around 1.5 nm from the graphene interface. As the concentration rises,
further density peaks are observed close to the carbon substrate, determined
by the emergence of a multi-layered electrical double-layer structure,
consistent with previously reported results [28, 18]. In such double-layer
configurations, adjacent solution layers, rich in cations or anions, arise at
the interface due to ion crowding (as in the case of the cations that are
attracted towards the negatively charged surface of the electrode) and ion
correlation (the localized positive excess charge in the closest layers to the
electrode, in turn, attracts the anions).

The results reported in Figure 3 are consistent with those of [40] with
NaCl(aq) and LiCl(aq) systems displaying, qualitatively, the same solution
side double layer structure. The case of KCl(aq) differs somewhat, with the
same position of the two first two peaks for Figure 3(e) in both cases, but a
different intensity compared with [40]. In turn, this intensity difference can
be due to the different classical force fields used for water, carbon and ions
as well as the use of scaled ionic charges not considered here. Besides these
rather minor differences, the results presented here seem to be robust with
respect to the chosen classical model. However, other results in the
literature (see [53]) show clear qualitative differences (in particular for
the KCl(aq) system where no adsorption is observed), most likely due to the
lack of dynamic polarization considered for the graphene electrodes.

(a) Na+

(b) K+

(c) Li+

(d) Cl- (NaCl)

(e) Cl- (KCl)

(f) Cl- (LiCl)

Figure 3:  Molar (M) density of the cations (top row) and the corresponding
anions (bottom row) for the three systems considered in this work. Black,
green, magenta, red, orange, blue curves correspond to bulk solution
concentrations 0.5, 2, 3, 4, 4.4 and 6 M, respectively.

###  3.2 Electrical Double Layer Properties

In this section we will derive and analyze the electrical properties of the
electrode-electrolyte systems considered in this work.

#### Electrode Charge Screeninig

We begin by considering the screening factor [28] $f$ defined as:

$f(z)=-\int_{0}^{z}\frac{\rho_{ions}(z^{\prime})}{\sigma}\mbox{d}z^{\prime}$
(3)

where $\sigma$ is the superficial charge of the electrode interface and
$\rho_{ions}(z)$ is the density charge of ions only, which is considered a
function of just the $z$ coordinate, i.e., it is averaged over the $x$ and $y$
coordinates.

The screening factor represents the extent to which the electrolyte phase
electrically screens the charged interface. When $f$ converges to a value of
one, the charge on the electrode is entirely shielded by the electrolyte. By
considering only the ions in the calculations of $f$, we can compare their
screening potential to predictions of simple mean field models. The
integration shown in Equation 3 and Equation 5 is performed numerically. Data
are first smoothed by applying the Savitzky-Golay [54] finite impulse response
smoothing filter of order 3 with a window width of 5 points, implemented in
Matlab. The smoothed curves obtained are then integrated using the trapezoidal
rule. Error bars are computed by error propagation through the integration
procedure.

The screening factors for all systems are reported in Figure 4. When the
concentration of the ions is below 1 M, an under-screening near the interface
can be observed. $f$ increases smoothly to a value of one at around $z=2$ nm.
This is qualitatively consistent with the predictions of Gouy-Chapman’s
theory. For higher concentrations, however, $f$ transitions to over-screening
at relatively small values of the $z$ coordinate. The over-screening,
highlighted by the first peak at $z\approx 0.6$ nm reported in Figures 4(a),
4(b) and 4(c), depends both on the particular ion and the bulk concentration.
In particular, the LiCl system has the strongest over-screening effect on the
electrode across the entire concentration range considered. Over-screening is
a well-known effect for ionic liquids [25] and is usually not considered
important in the electrolyte solutions, as this was only apparent at
relatively high concentrations [18, 28, 53]. The fact that over-screening
appears for higher concentration of the solute, in turn, can be linked
directly to the structuring of the ions near the interface observed in Figure
3. With the increase in concentration, the density of the cations closest to
the electrode increases with respect to their value in the solution bulk (see
Figure 3). The excess charge associated with this ion accumulation is balanced
in adjacent solution layers until the average bulk density is reached [14].
This description is consistent with our observations, where lithium and sodium
show a high degree of structuring near the interface relative to potassium
(i.e., multiple ion density peaks are observed, accompanied by a significant
over-screening effect). In contrast, potassium, with the lowest degree of
structuring near the interface, shows the smallest over-screening among the
three ion solutions considered. Moreover, for potassium, we observe a
variation in the slope of the screening factor when $z\approx 0.5$ nm, which
increases (becoming more pronounced) as a function of concentration. This
additional feature in the screening factor, absent in NaCl and LiCl, can be
explained by the direct coordination of the K+ (i.e., through the first
coordination sphere) to carbon atoms (as also observed in [40]), as opposed to
the behaviour of the cations in LiCl(aq) and NaCl(aq) systems (see the first
peak at $\approx$0.35 nm in Figure 3(b) with respect to the first peak at
$\approx$0.5 nm in Figures 3(c) and 3(a)).

(a) KCl

(b) LiCl

(c) NaCl

Figure 4:  Screening factor as defined in Equation 3 for the three systems
considered using ion solution charge densities only. We included only a subset
of the concentrations for clarity, and the results for all the concentrations
are reported in the SI (see Fig.S3 of the SI).

#### Electrode Polarisation

The coordination of the K+ with the carbon atoms on the graphene electrode is
shown qualitatively in Figure 3(b) for the lowest (0.5 M) and the highest
concentration (4.4 M) considered here. The plots in Figure 5 represents a
single snapshot in the 150 ns long simulation with the highest number of
potassium cations in direct contact with the interface (i.e., at a distance of
$0.26$ nm from the interface). As expected, the number of K+ in direct contact
with the interface increases as the bulk concentration of the cations
increases, consistent with the observation in Figure 4 for the short-distance
(from the electrode) behaviour of the screening factor, which increases with
concentration.

The accumulation of K+ in the nearby region to the negative electrode (see
Figure 3(b)) results in an increased non-uniformity of the partial charge
distribution on the electrode, with higher negative charges located on the
carbons closer to the coordinated K+. This, in turn, demonstrates how
polarisation effects are important to be considered in systems where direct
coordination of electrolytes to the electrode may occur.

(a) 0.5 M

(b) 4.4 M

Figure 5:  Representative plot of the computed Mulliken charges on the
graphene sheet charged with 4 and in contact with kCl solutions at different
concentrations. Circled X’s mark the coordinates of K ions directly adsorbed
on the surface.

(a) KCl

(b) LiCl

(c) NaCl

Figure 6:  Electrostatic potential as defined in Equation 5 for the three
systems considered. We included only a subset of the concentrations for
clarity. We included the results for all the concentrations in the SI (see
Fig.S3 of the SI).

#### Electrical Potential in the Double Layer

We calculated the electrical field $E(z)$ and the electrical potential,
$\psi(z)$ in the direction orthogonal to the interface using the Poisson
equation:

$-\frac{\mbox{d}^{2}\psi(z)}{\mbox{d}z^{2}}=\frac{\mbox{d}E(z)}{\mbox{d}z}=\frac{\rho(z)}{\epsilon(z)}$
(4)

where $\rho(z)$ is the charge density calculated for all atoms on the
perpendicular axis and we defined $\epsilon(z)=\epsilon_{r}(z)\epsilon_{0}$,
the product of the permittivity in vacuum $\epsilon_{0}$ and relative
permittivity $\epsilon_{r}$. It was reported that this latter quantity could
be a function of the distance from the electrode [55], a function of the
concentration of the electrolyte [16], or possibly both. Given such
uncertainties, we consider a constant relative permittivity equal to one in
this work. The electrical potential, $\psi$(z), is obtained from Equation 4 by
integrating twice with respect to the $z$-coordinate:

$\psi(z)=-\int_{0}^{z}\int_{0}^{z^{\prime}}\frac{\rho(\zeta)}{\epsilon(\zeta)}\mbox{d}\zeta\mbox{d}z^{\prime}$
(5)

The two integration constants in Equation 5 are chosen to set the
electrostatic field and potential equal to zero in the bulk, which amounts to
considering the bulk as the reference for the calculation of the electrostatic
potential.

The results of Equation 5 are reported in Figure 6 for a selection of
concentrations (see Figure S.3 of the SM for the entire range of
concentrations). In stark contrast to the exponential behaviour predicted by
models based on the Gouy-Chapman double layer theory which treats the solvent
medium as a continuum with known dielectric, atom/molecule finite-size effects
give rise to an undulating $\psi(z)$ function in the interfacial region at all
concentrations and in all systems. When calculating the charge distribution,
we include all solution species, including water molecules partial charges.
Hence, it is unsurprising that the structuring of ions and water molecules at
the interface gives rise to a significant departure from the predictions of
simple mean field models. Indeed, these finite size effects are a well-
reported feature of electrode-electrolyte systems [56, 57].

From a relatively large negative value of the potential at the electrode, the
(partial) charges of ions and water give rise to fluctuations that attenuate
at larger values of $z$, where the bulk solution behaviour is recovered.
Generally, increasing the bulk solution concentration increase the amplitude
of $\psi(z)$ fluctuations. Furthermore, it is evident from Figure 6(b) and
Figure 6(c) that the crowding of ions in the double-layer increases with
concentration as the positions of peaks and minima in $z$ shift to lower
values, a feature also observed by Finney et al. [28] with graphite and which
was related to changes in the screening factor. This concentration dependence
is less apparent in the case of KCl(aq), where the value of $\psi(z)$ at the
first maximum is less susceptible to changes in the concentration as opposed
to NaCl(aq) and LiCl(aq).

#### Electrical Double Layer Capacitance

The total capacitance $C_{TOT}$ in these kinds of systems is usually
considered as composed of three independent components combined in series: the
Electrochemical Double-Layer Capacitance (EDLC), $C_{EDL}$, and the quantum
capacitance (or the space charge capacitance)($C_{Q}$), depending on the
spatial distribution of the charges on the graphene [40]. The total
capacitance is then given by

$\frac{1}{C_{TOT}}=\frac{1}{C_{EDL}}+\frac{1}{C_{Q}}$ (6)

From Figure 6 we can easily derive the potential drop, $\Delta\psi$, across
the interface as*** A more precise notation for the potential drop across the
interface would have been $\Delta\Delta\psi=\Delta\psi^{-}-\Delta\psi_{ref}$.
[40] $\Delta\psi=\Delta\psi^{-}-\Delta\psi_{ref}$ where $\Delta\psi^{-}$ and
$\Delta\psi_{ref}$ represent the potential drop at the interface with respect
to the bulk for the charged and neutral electrodes, respectively. As a
reference for the calculation of the potential drop, we use the potential at
the interface in a neutral electrode with all other conditions unchanged. We
report the calculation of the potential across the system for a neutral
electrode in the SI (see Figure S.4 of the SI) along with the potential drop
at the interface ($\Delta\psi_{ref}$) (see Table S.2 of the SI). With this
definition of the potential drop, the EDLC can be obtained, as

$C_{EDL}=\frac{\sigma}{\Delta\psi}$ (7)

The quantum capacitance instead is obtained by calculating the differential
quantum capacitance $C_{Q}^{diff}$ according to [40]:

$C_{Q}^{diff}(\psi)=\frac{e^{2}}{4k_{B}T}\int_{-\infty}^{\infty}\left[D(E)\mbox{sech}^{2}\left(E+\psi\right)\right]\mbox{d}E$
(8)

Where $e$ is the electron charge, $E$ is the energy relative to the Fermi
level, $D(E)$ is the density of states at a given energy, $k_{B}$ is the
Boltzmann constant, and $T$ is the temperature. By integrating the
differential quantum capacitance with respect to the potential $\psi$ up to
the potential drop $\Delta\psi$ calculated for each system, we obtain the
integral quantum capacitance $C_{Q}$:

$C_{Q}=\frac{1}{\Delta\psi}\int_{0}^{\Delta\psi}C_{Q}^{diff}(\psi)\mbox{d}\psi$
(9)

For more detailed information about the calculation of the quantum capacitance
we refer the reader to our previous work [40].

concentration | $\Delta\psi$ | $C_{EDL}$ | $C_{Q}$ | $C_{TOT}$
---|---|---|---|---
| KCl
0.5 | -1.03 | 6.95 | 10.56 | 4.19
2.0 | -1.01 | 7.10 | 10.31 | 4.20
3.0 | -1.00 | 7.16 | 10.23 | 4.21
4.0 | -0.984 | 7.28 | 9.81 | 4.18
4.4 | -0.988 | 7.25 | 10.07 | 4.22
| LiCl
0.5 | -1.05 | 6.82 | 10.78 | 4.18
2.0 | -1.01 | 7.09 | 10.31 | 4.20
3.0 | -1.00 | 7.16 | 10.21 | 4.21
4.0 | -1.00 | 7.16 | 10.21 | 4.21
4.4 | -1.09 | 6.57 | 11.25 | 4.15
6.0 | -1.02 | 7.02 | 10.44 | 4.20
| NaCl
0.5 | -1.05 | 6.82 | 10.78 | 4.18
2.0 | -1.02 | 7.02 | 10.44 | 4.20
3.0 | -0.997 | 7.18 | 10.17 | 4.21
4.0 | -0.996 | 7.19 | 10.15 | 4.21
4.4 | -0.986 | 7.26 | 10.05 | 4.22
6.0 | -1.00 | 7.16 | 10.21 | 4.21
Table 1: Electrostatic potential drop ($\Delta\psi$) across the interface (in
V), Electrochemical double layer capacitance $C_{EDL}$, Quantum Capacitance
$C_{Q}$, and total capacitance $C_{TOT}$ (in $\mu$F cm-2) for each
concentration considered (in M).

The results for $C_{Q}$, $C_{EDL}$, and $C_{TOT}$ for all of the systems
considered are reported in Table 1. The data show that the total capacitance
is practically constant across all the concentration range and for all
solution types. The largest variation in $C_{TOT}$ we obtained among all the
systems is $\approx$2% (between the LiCl(aq) and KCl(aq) at 4.4 M). This
result contrasts with the different behaviour of the three cations in solution
and near the electrode interfaces, as highlighted in the discussion of the
number density of ionic species at the interface (see Figure 3) their
screening effect on the charge of the electrode (Figure 4), and as further
discussed in the following section in relation to their clustering properties.

An important point we want to highlight here is that such differences in the
behaviour of the cation in solution can be correctly captured through the use
of a simulation protocol that combines the pseudo-open boundary condition,
i.e., C$\mu$MD to maintain constant composition electroneutral bulk solutions
beyond the double layer, and the quantum mechanical description for the
distribution of partial charges of the electrode. However, while the
capacitance is a critical parameter for this kind of system’s applications as
supercapacitors, we showed here that the physics of the interfaces between
graphene electrodes and electrolytes is much richer than the one captured by
such quantity.

###  3.3 Ion Association

An often overlooked effect in systems in alkali chloride solutions is the
tendency for ions to associate, forming clusters. In particular, even simple
salt solutions exhibit significant non-ideal behaviour at high concentrations.
Recent experiments [58] and simulations [59] have shown that extended liquid-
like clusters exist in bulk NaCl(aq) at high concentrations and the extent of
these ionic networks is promoted in the double layer at carbon surfaces [28].
Since the effectiveness of the graphene-electrolyte devices often depends on
the ability to ‘build up the double layer’ (i.e., accumulate ions from the
bulk solution in the interfacial region), the structure and mobility of ion
species can be essential to this.

#### Ion Clusters

To identify and characterise ion associates in the simulations in this work,
pairwise RDFs were computed (see Fig. S.1 of the Supplementary Material,
(SM)), and the first minima in these informed truncation distances ($r_{c}$)
for first-sphere ion-ion coordination. $r_{c}=0.29$, 0.34 and 0.39 nm for Li–
Na– and K–Cl, respectively, reflecting the different sizes of the cations.
Clusters were identified as fully connected networks in the graph of adjacent
ion-ion connections according to this geometric criteria, regardless of their
total charge or lifetime.

Figure 7 provides the average first-sphere coordination number between cations
and O of water (see Figure 7(a)) as well as cations and anions for all
systems, calculated using Equation 2.

(a) Cation/Water

(b) Cation/Anion

Figure 7:  Coordination number for the different systems at the different
concentrations.

The results shown in Figure 7 indicate no significant surface effect on the
coordination of cations with water or chloride when ions in the interface
($0<z<2.5$ nm) and bulk ($2.5<z<4.5$ nm) regions were investigated. There is a
slight increase in the mean cation-anion coordination, and a concomitant
decrease in cation-water coordination, at the interface compared to the bulk;
however, this difference is within the margin of error. Generally, the effect
of increasing concentration is to increase the number of cation-anion
contacts, particularly for KCl(aq), where the coordination number is more than
double that of the other systems for all concentrations (and with Li–Cl
coordination being negligible even at 6 M). From the largest to smallest
variation in the coordination number we can write
$\mbox{K}^{+}\rightarrow\mbox{Na}^{+}\rightarrow\mbox{Li}^{+}$. This trend
follows the decrease of the ion radius and it is likely due to the stronger
binding of water in the solvation spheres of smaller cations. Furthermore, the
average cation-water coordination number is unchanging with a concentration
within the margin of error.

In simulations of NaCl(aq) in contact with graphite, [28] the substrate was
found to increase cation-anion correlations in the double layer with respect
to the bulk, particularly beyond $5$ M. It is important to note that different
models (due to the different system) were used and also that system size
likely plays a role in the extent that clusters can grow (both in e.g., the
system-size dependence of the availability of ions to form associates and the
extent to which finite-size and percolating clusters may form in effectively
confined canonical systems.)

The change in coordination for different salts is reflected in the cluster
size probability distributions presented in Figure 8 for the case of 4.4 M (we
report the results for the entire range of concentrations in Figure S.5 of the
SM). There is a clear difference in the extent to which clusters can grow,
with lithium forming clusters containing at most four ions and potassium
forming much larger networks containing as many as 35 ions. Even at the
highest concentrations, the majority of the Li+ are dispersed in solution,
fully solvated in their first shell. A snapshot of a configuration obtained
during the simulation of KCl at 4.4 M is shown in Figure 8. Although the most
probable clusters contain only a few ions (for clusters composed of five ion
units, we obtained a relative frequency of 0.01), larger species do contribute
to the charge storage capacity and must be considered. What we observe is a
stronger tendency of the potassium to associate into large aggregates—albeit
ones which are highly dynamic on the timescales of the simulations—compared to
sodium or lithium.

Since the KCl(aq) system shows the formation of large aggregates of ions, it
is interesting to study the relative frequency of the charge of these
aggregates. In Figure 9 we plot the 2-dimensional histogram showing the
relative frequencies of the charge vs. the cluster size for the KCl(aq)
system. The histogram is skewed towards positive charges, with the appearance
of clusters containing an excess of positive charge as large as +7e, although
the majority of the clusters are neutral.

Figure 8:  On the left: Histogram of the relative frequency of the cluster of
different sizes for the concentration of 4.4 M. In the inset, the same
quantity for the 0.5 M case. On the right: an example of a cluster composed of
26 ions for the KCl system at 4.4 M.  Figure 9:  2-dimensional histogram
(charge VS dimension of the clusters) for the KCl(aq) system at the largest
concentration considered (4.4 M).

#### Ion Mobilities

As well as a high capacity to store charge, an optimal charge storage device
must also be a good electrical conductor. In this regard, we might expect the
conductivity of solutions to decrease when clusters are present. Indeed, this
can be perceived as a relative decrease in the activity of charge carriers due
to increasingly non-ideal solutions. To test this, we calculated the
conductivity of bulk NaCl(aq) solutions with concentrations ranging from 1–10
M from the ion diffusion coefficients calculated by Finney and Salvalaglio in
finite size systems and in the dilute limit [59, 60]. To determine
conductivity we use the Nernst-Einstein equation:

$\sigma_{NE}=\frac{e^{2}}{Vk_{\mathrm{B}}T}(N_{+}z_{+}^{2}D_{+}+N_{-}z_{-}^{2}D_{-})$
(10)

where $e$, $V$, $k_{\mathrm{B}}$ and $T$ are the elementary charge, simulation
cell volume, Boltzmann’s constant and temperature, respectively. $N$ and $D$
are the total number of ions and diffusion coefficients for ions with charge
indicated by the subscript. Furthermore, we assume that, given the highly
dynamic nature of the clusters observed in solution, the valency of ionic
species, $z$, is equal to one.

Figure 10: Solution conductivities, $\sigma_{NE}$, of bulk NaCl(aq) solutions
calculated for a range of concentrations. To this aim, the Nernst-Einstein
equation was adopted where ion diffusion coefficients were determined from
simulations at finite concentration, $D_{ion}$ (blue), or from a single
simulation at the dilute limit, $D_{ion}^{0}$ (red). Dashed lines are a guide
for the eye, while error bars indicate uncertainties in the conductivities
associated with the calculated $D$ value from Refs. [59] and [60].

Figure 10 provides the solution conductivities for NaCl(aq) where either the
diffusion of ions in finite concentration simulations was used ($D_{ion}$) or
the diffusion of ions in the dilute limit ($D_{ion}^{0}$) was considered. For
the latter, ions are assumed to be completely, dispersed as association, even
beyond the second solvation sphere, did not occur in simulations at the dilute
limit. For the estimate of $D_{ion}^{0}$, Finney and Salvalaglio [60]
performed extended simulations of a single cation and anion in 4,000 water
molecules; here, $D_{+}^{0}=1.223\pm 0.005\times 10^{-5}$ cm2 s-1 and
$D_{-}^{0}=1.282\pm 0.008\times 10^{-5}$ cm2 s-1. In all cases, diffusion
coefficients were corrected to account for simulation finite size effects
[61]. Unsurprisingly, a linear correlation in $\sigma_{NE}$ as a function of
concentration is found when a constant $D_{ion}^{0}$ is used for the diffusion
of ions, independent of concentration. This is inaccurate at relatively high
concentrations, given the simulation and experimental observations of ion-ion
correlations[59, 58].

When accounting for the non-idealities in the solution and the formation of
clusters explicitly in the diffusion of ions, we find that the solution
conductivity reaches an upper limit between 4 and 5 M. At the lowest
concentrations (1–2 M), the conductivity from finite concentration and dilute
simulation data agree, and the simulation predictions match well with
experimental measurements [62]. A crossover in the conductivity behaviour from
the ‘pseudo-ideal’ to non-ideal regime occurs between 2 and 3 M. Therefore,
over a wide concentration range up to the salt solubility, non-idealities will
affect the performance of electrical devices; depending upon the chosen
application, electrolytes should be chosen to minimize these effects.

##  4 Conclusions

In this work, we presented an extended set of simulations describing the
interface between three different electrolyte solutions - (KCl(aq), LiCl(aq),
and NaCl(aq)) - in contact with the surface of a negatively charged graphene
electrode. To investigate these systems, we combined QM/MD and C$\mu$MD
methodologies into a new simulation framework. QM/MD models of the graphene
electrode in contact with an electrolyte enabled the explicit coupling of the
electrode polarizability with the instantaneous configuration of the
electrolyte. The latter was maintained in equilibrium with a liquid phase at
constant bulk concentration thanks to the C$\mu$MD model, which mimics open-
boundary conditions.

We performed a thorough analysis of the interaction of the ions with the
electrode by showing the different behaviour of the three cations in the
double layer, focusing on K+, which, according to our results, is able to
directly adsorb at the electrode surface at shorter distances compared to Li+
and Na+, modifying the screening effect of the solution.

Calculations of the integral capacitance indicated no concentration dependence
or specific ion effects, with a total capacitance of around 4.2 µF cm-2 across
all systems. However, the lack of variation in capacitance hides the rich
electrolyte solution behaviour, particularly for the ions close to the
electrode. We showed, for example, that large KCl clusters emerge in solution,
which might be important when considering properties associated with ion
mobility and charge transfer.

Our results indicate that accurate models of the interface - able to account
for the position-dependent non-ideality of electrolyte solutions - better
capture the configurational and dynamical details underpinning the
electrochemical behavior of interfaces at the atomistic level, and that is
often overshadowed by the calculation of aggregated quantities such as the
integral capacitance. We plan to extend our calculations to include a range of
charged electrodes, both positive and negative, and further investigate ion
dynamics in solution.

##  5 Acknowledgements

We acknowledge the support provided by the IT Services use of the
Computational Shared Facility (CSF) and at the University of Manchester. NDP,
JDE and PC thank the European Union’s Horizon 2020 research and innovation
programme project VIMMP under Grant Agreement No. 760907. ARF and MS
acknowledge funding from an EPSRC Programme Grant (Grant EP/R018820/1), which
funds the Crystallisation in the Real World consortium.

## References

  * Cohen-Tanugi and Grossman [2012] David Cohen-Tanugi and Jeffrey C Grossman.  Water desalination across nanoporous graphene.  _Nano letters_ , 12(7):3602–3608, 2012.
  * Heiranian et al. [2021] Mohammad Heiranian, Yechan Noh, and Narayana R Aluru.  Dynamic and weak electric double layers in ultrathin nanopores.  _The Journal of Chemical Physics_ , 154(13):134703, 2021.
  * Surwade et al. [2015] Sumedh P Surwade, Sergei N Smirnov, Ivan V Vlassiouk, Raymond R Unocic, Gabriel M Veith, Sheng Dai, and Shannon M Mahurin.  Water desalination using nanoporous single-layer graphene.  _Nature nanotechnology_ , 10(5):459–464, 2015\.
  * Simon and Gogotsi [2008] Patrice Simon and Yury Gogotsi.  Materials for electrochemical capacitors.  In _Nanoscience and technology: a collection of reviews from Nature journals_ , pages 320–329. World Scientific, 2008.
  * Wang et al. [2021] Yifan Wang, Lin Zhang, Haoqing Hou, Wenhui Xu, Gaigai Duan, Shuijian He, Kunming Liu, and Shaohua Jiang.  Recent progress in carbon-based materials for supercapacitor electrodes: a review.  _Journal of Materials Science_ , 56(1):173–200, 2021.
  * Elliott et al. [2022a] J. Elliott, A. A. Papaderakis, R. Dryfe, and P. Carbone.  The electrochemical double layer at the graphene/aqueous electrolyte interface: what we can learn from simulations, experiments, and theory.  _Journal of Materials Chemistry C_ , 2022a.
  * Wang et al. [2009] Yan Wang, Zhiqiang Shi, Yi Huang, Yanfeng Ma, Chengyang Wang, Mingming Chen, and Yongsheng Chen.  Supercapacitor devices based on graphene materials.  _Journal of Physical Chemistry C_ , 113(30):13103–13107, 2009.
  * Liu et al. [2010] Chenguang Liu, Zhenning Yu, David Neff, Aruna Zhamu, and Bor Z Jang.  Graphene-based supercapacitor with an ultrahigh energy density.  _Nano letters_ , 10(12):4863–4868, 2010.
  * Yu et al. [2010] Aiping Yu, Isaac Roes, Aaron Davies, and Zhongwei Chen.  Ultrathin, transparent, and flexible graphene films for supercapacitor application.  _Applied physics letters_ , 96(25):253105, 2010\.
  * Zhang et al. [2010] Li Li Zhang, Rui Zhou, and XS Zhao.  Graphene-based materials as supercapacitor electrodes.  _Journal of Materials Chemistry_ , 20(29):5983–5992, 2010.
  * Zhu et al. [2011] Y. Zhu, S. Murali, M. D. Stoller, K. J. Ganesh, W. Cai, P. J. Ferreira, A. Pirkle, R. M. Wallace, K. A. Cychosz, M. Thommes, D. Su, E. A. Stach, and R. S. Ruoff.  Carbon-based supercapacitors produced by activation of graphene.  _science_ , 332(6037):1537–1541, 2011.
  * An et al. [2001] Kay Hyeok An, Won Seok Kim, Young Soo Park, Young Chul Choi, Seung Mi Lee, Dong Chul Chung, Dong Jae Bae, Seong Chu Lim, and Young Hee Lee.  Supercapacitors using single-walled carbon nanotube electrodes.  _Advanced Materials_ , 13(7):497–500, 2001.
  * Yang et al. [2019] Zhoufei Yang, Jiarui Tian, Zefang Yin, Chaojie Cui, Weizhong Qian, and Fei Wei.  Carbon nanotube-and graphene-based nanomaterials and applications in high-voltage supercapacitor: A review.  _Carbon_ , 141:467–480, 2019.
  * Merlet et al. [2013] C. Merlet, B. Rotenberg, P. A. Madden, and M. Salanne.  Computer simulations of ionic liquids at electrochemical interfaces.  _Physical Chemistry Chemical Physics_ , 15(38):15781–15792, 2013.
  * Iamprasertkun et al. [2019] P. Iamprasertkun, W. Hirunpinyopas, A. Keerthi, B. Wang, B. Radha, M. A. Bissett, and R. A. W. Dryfe.  Capacitance of basal plane and edge-oriented highly ordered pyrolytic graphite: specific ion effects.  _The Journal of Physical Chemistry Letters_ , 10(3):617–623, 2019.
  * Yang et al. [2017] H. Yang, J. Yang, Z. Bo, X. Chen, X. Shuai, J. Kong, J. Yan, and K. Cen.  Kinetic-dominated charging mechanism within representative aqueous electrolyte-based electric double-layer capacitors.  _The Journal of Physical Chemistry Letters_ , 8(15):3703–3710, 2017.
  * Qu et al. [2008] QT Qu, B Wang, LC Yang, Y Shi, S Tian, and YP Wu.  Study on electrochemical performance of activated carbon in aqueous li2so4, na2so4 and k2so4 electrolytes.  _Electrochemistry Communications_ , 10(10):1652–1655, 2008.
  * Elliott et al. [2020] J. D. Elliott, A. Troisi, and P. Carbone.  A QM/MD coupling method to model the ion-induced polarization of graphene.  _J. Chem. Theory Comput._ , 16(8):5253–5263, 2020\.
  * Wang et al. [2008] Xinran Wang, Scott M Tabakman, and Hongjie Dai.  Atomic layer deposition of metal oxides on pristine and functionalized graphene.  _Journal of the American Chemical Society_ , 130(26):8152–8153, 2008.
  * Gouy [1910] MJJPTA Gouy.  Sur la constitution de la charge électrique à la surface d’un électrolyte.  _J. Phys. Theor. Appl._ , 9(1):457–468, 1910\.
  * Chapman [1913] D. L. Chapman.  LI. A contribution to the theory of electrocapillarity.  _The London, Edinburgh, and Dublin philosophical magazine and journal of science_ , 25(148):475–481, 1913.
  * Stern [1924] Otto Stern.  Zur theorie der elektrolytischen doppelschicht.  _Zeitschrift für Elektrochemie und angewandte physikalische Chemie_ , 30(21-22):508–516, 1924.
  * Popović and Šiber [2013] Marko Popović and Antonio Šiber.  Lattice-gas poisson-boltzmann approach for sterically asymmetric electrolytes.  _Physical Review E_ , 88(2):022302, 2013.
  * Borukhov et al. [1997] Itamar Borukhov, David Andelman, and Henri Orland.  Steric effects in electrolytes: A modified poisson-boltzmann equation.  _Physical review letters_ , 79(3):435, 1997.
  * Fedorov and Kornyshev [2008] M. V. Fedorov and A. A. Kornyshev.  Ionic liquid near a charged wall: Structure and capacitance of electrical double layer.  _The Journal of Physical Chemistry B_ , 112(38):11868–11872, aug 2008.  doi: 10.1021/jp803440q.
  * Howard et al. [2010] J. J. Howard, J. S. Perkyns, and B. M. Pettitt.  The behavior of ions near a charged wall-dependence on ion size, concentration, and surface charge.  _The journal of physical chemistry B_ , 114(18):6074–6083, 2010.
  * Misra and Blankschtein [2021a] Rahul Prasanna Misra and Daniel Blankschtein.  Ion adsorption at solid/water interfaces: Establishing the coupled nature of ion–solid and water–solid interactions.  _The Journal of Physical Chemistry C_ , 125(4):2666–2679, 2021a.
  * Finney et al. [2021] Aaron R Finney, Ian J McPherson, Patrick R Unwin, and Matteo Salvalaglio.  Electrochemistry, ion adsorption and dynamics in the double layer: a study of nacl (aq) on graphite.  _Chemical science_ , 12(33):11166–11180, 2021\.
  * Perego et al. [2015] C. Perego, M. Salvalaglio, and M. Parrinello.  Molecular dynamics simulations of solutions at constant chemical potential.  _The Journal of Chemical Physics_ , 142(14):144113, apr 2015.  doi: 10.1063/1.4917200.
  * Doblhoff-Dier and Koper [2021] K. Doblhoff-Dier and M. T. M. Koper.  Modeling the Gouy–Chapman diffuse capacitance with attractive ion–surface interaction.  _The Journal of Physical Chemistry C_ , 125(30):16664–16673, 2021.
  * France-Lanord and Grossman [2019] Arthur France-Lanord and Jeffrey C Grossman.  Correlations from ion pairing and the nernst-einstein equation.  _Physical review letters_ , 122(13):136001, 2019\.
  * Misra and Blankschtein [2021b] Rahul Prasanna Misra and Daniel Blankschtein.  Uncovering a Universal Molecular Mechanism of Salt Ion Adsorption at Solid/Water Interfaces.  _Langmuir_ , 37(2):722–733, January 2021b.
  * Williams et al. [2017] C. D. Williams, J. Dix, A. Troisi, and P. Carbone.  Effective Polarization in Pairwise Potentials at the Graphene–Electrolyte Interface.  _J. Phys. Chem. Lett._ , 8(3):703, February 2017\.
  * Siepmann and Sprik [1995] J Ilja Siepmann and Michiel Sprik.  Influence of surface topology and electrostatic potential on water/electrode systems.  _The Journal of chemical physics_ , 102(1):511, 1995.
  * Di Pasquale and Davidchack [2020] N. Di Pasquale and R. L. Davidchack.  Shuttleworth Equation: A Molecular Simulations Perspective.  _Journal of Chemical Physics_ , 153:154705, 2020.
  * Abraham et al. [2015] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl.  GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.  _SoftwareX_ , 1:19–25, 2015.
  * Tribello et al. [2014] G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi.  Plumed 2: New feathers for an old bird.  _Computer physics communications_ , 185(2):604–613, 2014.
  * Hourahine et al. [2020] B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W.-z. Yu, and T. Frauenheim.  DFTB$+$, a software package for efficient approximate density functional theory based atomistic simulations.  _The Journal of Chemical Physics_ , 152(12):124101, 2020.  doi: 10.1063/1.5143190.
  * Elstner et al. [1998] Marcus Elstner, Dirk Porezag, G Jungnickel, J Elsner, M Haugk, Th Frauenheim, Sandor Suhai, and Gotthard Seifert.  Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties.  _Physical Review B_ , 58(11):7260, 1998.
  * Elliott et al. [2022b] J. D. Elliott, M. Chiricotto, A. Troisi, and P. Carbone.  Do specific ion effects influence the physical chemistry of aqueous graphene-based supercapacitors? perspectives from multiscale qmmd simulations.  _arXiv_ , 2022b.  doi: 10.48550/ARXIV.2203.02469.
  * Mulliken [1955] Robert S Mulliken.  Electronic population analysis on lcao–mo molecular wave functions. i.  _The Journal of Chemical Physics_ , 23(10):1833–1840, 1955.
  * Haynes et al. [2016] William M Haynes, David R Lide, and Thomas J Bruno.  _CRC handbook of chemistry and physics_.  CRC press, 2016.
  * Zeron et al. [2019] I. M. Zeron, J. L. F. Abascal, and C. Vega.  A force field of Li+, Na+, K$+$, Mg2+, Ca2+, Cl-, and SO${}_{4}^{2-}$ in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions.  _The Journal of chemical physics_ , 151(13):134504, 2019.
  * Ho and Striolo [2013] Tuan A Ho and Alberto Striolo.  Capacitance enhancement via electrode patterning.  _The Journal of Chemical Physics_ , 139(20):204708, 2013.
  * Xu et al. [2020] Kui Xu, Hui Shao, Zifeng Lin, Céline Merlet, Guang Feng, Jixin Zhu, and Patrice Simon.  Computational insights into charge storage mechanisms of supercapacitors.  _Energy & environmental materials_, 3(3):235–246, 2020.
  * Hagberg et al. [2008] Aric Hagberg, Pieter Swart, and Daniel S Chult.  Exploring network structure, dynamics, and function using networkx.  Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.
  * Berendsen et al. [1995] H. J. C. Berendsen, D. van der Spoel, and R. van Drunen.  Gromacs: a message-passing parallel molecular dynamics implementation.  _Computer physics communications_ , 91(1-3):43–56, 1995.
  * van der Spoel et al. [2005] D. van der Spoel, E. Lindhal, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen.  Gromacs: Fast, flexible and free.  _J. Comput. Chem._ , 26(16):1701–1718, 2005.
  * Berendsen et al. [1987] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma.  The missing term in effective pair potentials.  _Journal of Physical Chemistry_ , 91(24):6269–6271, 1987.
  * Miyamoto and Kollman [1992] S. Miyamoto and P. A. Kollman.  SETTLE: An analytical version of the SHAKE and RATTLE algorithm for rigid water molecules.  _J. Comput. Chem_ , 13:952–962, 1992.
  * Huang et al. [2012] Yi Huang, Jiajie Liang, and Yongsheng Chen.  An overview of the applications of graphene-based materials in supercapacitors.  _small_ , 8(12):1805–1834, 2012.
  * Joung and Cheatham [2008] I. S. Joung and T. E. III Cheatham.  Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations.  _The journal of physical chemistry B_ , 112(30):9020–9041, 2008.
  * Dočkal et al. [2022] J. Dočkal, M. Lísal, and F. Moučka.  Molecular dynamics of the interfacial solution structure of alkali-halide electrolytes at graphene electrodes.  _Journal of Molecular Liquids_ , 353:118776, 2022.
  * Savitzky and Golay [1964] Abraham Savitzky and Marcel JE Golay.  Smoothing and differentiation of data by simplified least squares procedures.  _Analytical chemistry_ , 36(8):1627–1639, 1964\.
  * Olivieri et al. [2021] Jean-François Olivieri, James T Hynes, and Damien Laage.  Confined water’s dielectric constant reduction is due to the surrounding low dielectric media and not to interfacial molecular ordering.  _The Journal of Physical Chemistry Letters_ , 12(17):4319–4326, 2021.
  * Kornyshev [2007] Alexei A Kornyshev.  Double-layer in ionic liquids: paradigm change?  _The Journal of Physical Chemistry B_ , 111(20):5545–5557, 2007.
  * Snook and van Megen [1981] Ian Snook and William van Megen.  Finite ion size effects in the electrical double layer—a monte carlo study.  _The Journal of Chemical Physics_ , 75(8):4104–4106, 1981.
  * Hwang et al. [2021] Hyerim Hwang, Yong Chan Cho, Sooheyong Lee, Yun-Hee Lee, Seongheun Kim, Yongjae Kim, Wonhyuk Jo, Patrick Duchstein, Dirk Zahn, and Geun Woo Lee.  Hydration breaking and chemical ordering in a levitated NaCl solution droplet beyond the metastable zone width limit: evidence for the early stage of two-step nucleation.  _Chemical Science_ , 12(1):179–187, 2021.  doi: 10.1039/d0sc04817h.  URL https://doi.org/10.1039%2Fd0sc04817h.
  * Finney and Salvalaglio [2022] Aaron R Finney and Matteo Salvalaglio.  Multiple pathways in nacl homogeneous crystal nucleation.  _Faraday Discussions_ , 2022.
  * Finney and Salvalaglio [2021] A. R. Finney and M. Salvalaglio.  Bridging the gap between mesoscopic and molecular models of solid/liquid interfaces out-of-equilibrium.  _arXiv preprint arXiv:2109.00568_ , 2021.
  * Yeh and Hummer [2004] I. Yeh and G. Hummer.  System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions.  _The Journal of Physical Chemistry B_ , 108(40):15873–15879, October 2004.
  * Widodo et al. [2018] C. S. Widodo, H. Sela, and D. R. Santosa.  The effect of nacl concentration on the ionic nacl solutions electrical impedance value using electrochemical impedance spectroscopy methods.  In _AIP Conference Proceedings_ , volume 2021, page 050003. AIP Publishing LLC, 2018.