To work with this tutorial, you should be familiar with the following topics:

- Setting up and running simulations in ESPResSo - creating particles,
incorporating interactions.
If you are unfamiliar with this, you can go through the tutorial
in the
`lennard_jones`

folder. - Basic knowledge of classical electrostatics: dipoles, surface and image charges.
- Reduced units and how to relate them to physical quantities, see the ESPResSo user guide.

Ionic liquid (IL) based capacitors have long been established as promising candidates in the area of efficient energy storage devices due to their extraordinary capacitance, which is why they are also termed super-capacitors. Typically, in these setups a fluid consisting of mobile charge carriers is placed between two electrodes and thus gets polarized upon application of an external potential formic an electric double layer at the interfaces [1]. Electric double-layer capacitors (EDLCs) can be constructed from electrodes of various geometries and materials where energy is stored by potential driven adsorption of counterions on the surface of the electrodes and forming the double layer. Thus, conducting, high surface area electrode materials can further maximize the energy per volume.

In this tutorial we are going to investigate the ionic layer formation between
two conducting dielectric walls in presence of an applied voltage using **ESPResSo**'s
"Electrostatic layer correction with image charges" (ELC-IC) method
[2].
We employ a primitive model of an aqueous salt solution, where the solvent is
treated implicitly as a homogeneous dielectric background.
Thus, within the limits of a mean-field treatment and for not too small slit
widths, we can compare our findings with the analytical Gouy–Chapmann
solution for a single charged plane since additivity is assumed to hold if the
potential of both walls is screened sufficiently.
While this mean-field formalism properly describes the behavior of Coulomb
fluids composed of monovalent ions at low concentrations in the vicinity of
weakly charged interfaces, for strongly charged systems, where correlation and
finite size effects begin to dominate, Poisson–Boltzmann theory falls
inadequate.
Our goal in this tutorial is to demonstrate how coarse grained implicit solvent
simulations can corroborate on some of the approximations and hint on
extensions/deviations.

The inclusion of dielectric inhomogeneities appearing at the conducting interfaces demands to take into account image effects that involve the full solution of the Poisson equation. This is dealt with in a computationally cost-effective way using the ELC-IC method to treat the image charge effects in the presence of 2D dielectric bounding interfaces.

Charged surfaces in contact with a liquid containing free charges (ions) attract oppositely charged ions that form a diffusive electric double layer. The competition between electrostatic interactions and entropy of ions in solution determines the exact distribution of mobile ions close to charged membranes. Gouy [3] and Chapman [4] derived in the early 20th century the analytic solution for the case of a single planar wall within the mean-field approximation of the Poisson–Boltzmann (PB) equation. We will use it to describe our two-electrode system, which is justified if the electrodes are so far apart that one surface does not influence the ion distribution in front of the other surface.

In the case of a monovalent electrolyte, double integrating the PB equation and employing the corresponding boundary conditions for the charged plane located at $z=0$ yields the electrostatic potential: $$\phi(z) = -2\ln\left[ \frac{1-\tanh(\phi_\mathrm{s}/4)e^{-\kappa_\mathrm{D} z}} {1+\tanh(\phi_\mathrm{s}/4)e^{-\kappa_\mathrm{D}z}} \right].$$ Here, $\phi_\mathrm{s}=\phi(z=0)=$ const is the surface potential such that $\phi(z\rightarrow \infty)=0$. $\kappa_\mathrm{D} = \lambda_\mathrm{D}^{-1}$ is the inverse Debye screening length given by $$ \lambda_\mathrm{D} = \left(\frac{\varepsilon \, k_{\mathrm B} T}{\sum_{j = 1}^N n_j^0 \, q_j^2}\right)^{1/2}, $$ where $n_j^{(0)}$ and $q_j$ are the equilibrium number density and charge of the $j$-th ion species. For the monovalent salt this can conveniently be expressed in terms of the Bjerrum length $\ell_\mathrm{B}$ and the equilibrium salt concentration $\rho^{(0)}=\sum_j \rho_j^{(0)}$, $$ \lambda_{\mathrm D} = \left(4 \pi \, \ell_\mathrm{B} \, \rho^{(0)}/e\right)^{-1/2} . $$

The cationic and anionic density profiles then follow from the Boltzmann distribution as: $$n_{\pm}(z)=n_\pm^{(0)}\left(\frac{1\pm\gamma e^{-\kappa_\mathrm{D}z}} {1\mp\gamma e^{-\kappa_\mathrm{D}z}} \right)^2$$ Here, $\gamma$ is associated with the surface potential as $\phi_\mathrm{s}=-4\tanh^{-1}(\gamma)$. At large z, where the potential decays to zero, the ionic profiles tend to their bulk (reservoir) densities, $n_\pm(z\to\infty) = n_\pm^{(0)}$

The relation between the surface potential and the surface charge induced on the electrode is given by the Grahame Equation [5] : $$ \sigma = \sinh(\phi_\mathrm{s}/2) \sqrt{\frac{2 n_\mathrm{b}}{\pi \ell_\mathrm{B}}} $$ The latter expression thus yields the differential capacitance $C=\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}\phi_\mathrm{s}}$ within the mean-field solution for non-overlapping double layers.

In this tutorial we employ a parallel plate capacitor setup with constant potential boundary conditions which needs to be treated appropriately by the electrostatics solver. To simulate a two-dimensional partially periodic system, we combine the efficient scaling behavior of three-dimensional mesh-based solvers (typically $\mathcal{O}(N \log N)$ for P3M) with the Electrostatic Layer Correction (ELC) [1]. The latter removes the contributions from the periodic images in the non-periodic direction and its numerical cost grows linear with the number of point charges $N$, hence the performance overall depends on the underlying 3D Coulomb solver. Furthermore, ELC can be extended straightforwardly to metallic boundary conditions (or any other dielectric contrast) by using the method of image charges, which is referred to as the “Electrostatic Layer Correction with Image Charges” (ELC-IC) approach used in this tutorial.

A voltage difference can be applied between the electrodes by following considerations: The total potential drop $\Delta \phi$ across the simulated system is readily obtained from the ion distribution and integrating twice over the one-dimensional Poisson equation: $$-\varepsilon_{0}\varepsilon_{r}\phi_\mathrm{ind}(z)=\iint_{0}^{z}\rho(z^{\prime})(z-z^{\prime})dz^{\prime}$$ Here, the subscript 'ind' indicates that this is the potential due to the induced inhomogeneous charge distribution. In order to set up a constant potential difference $\Delta \phi$, a homogeneous electric field is superimposed such that $$ \Delta \phi = \Delta \phi_\mathrm{ind} + \Delta \phi_\mathrm{bat},$$ where $\Delta \phi_\mathrm{bat}$ corresponds to the potential of a (virtually) applied battery. In practice, the linear electric field in $E_z^\mathrm{(bat)}=-\phi_\mathrm{bat}/L_z$ in the $z$-direction normal to the surface that one needs to apply can be calculated straightforwardly, as the corresponding contribution from the induced charges is known: $$ E_z^\mathrm{(ind)} = \frac{1}{\varepsilon_0 \varepsilon_r L^2 L_z} P_z$$ Here, $L$ denotes the lateral system size, $L_z$ the distance between the plates and $P_z = \sum_i q_i z_i$ is the total dipole moment in $z$-direction due to the charges $q_i$. Then, to maintain $\Delta \phi$, a force $F_z^\mathrm{bat} = q E_z^\mathrm{(bat)}$ is applied on all charges. Since ELC already calculates $P_z$, the constant potential correction requires no additional computational effort.

*Note*: Apart from ELC-IC, **ESPResSo** also provides the ICC$^\star$ method
[2] which employs an iterative numerical scheme with
discretized surface particles to solve the boundary integrals at the dieletcric
interface.
The tutorial on *Basic simulation of electrodes in ESPResSo part I* addresses this
in detail.

First we import all ESPResSo features and external modules

In [1]:

```
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as constants
import scipy.stats
import tqdm
import espressomd
import espressomd.electrostatics
import espressomd.electrostatic_extensions
import espressomd.observables
import espressomd.accumulators
import espressomd.shapes
espressomd.assert_features(['WCA', 'ELECTROSTATICS'])
rng = np.random.default_rng(42)
plt.rcParams.update({'font.size': 18})
```

We need to define system dimensions and some physical parameters related to length, time and energy scales of our system. As discussed in previous tutorials, all physical parameters are defined in terms of a length $\sigma$, mass $m$ and time $t$ and unit of charge $q$. Since we are not explicitly interested in the dynamics of the system, we set the mass to $m=1$ (particle mass). For convenience, we choose the elementary charge as fundamental unit ($q=1e$) and $\sigma = 1 \,\mathrm{nm}$.

With this we can now define the fundamental parameters of our system:

In [2]:

```
# water at room temperature
EPSILON_R = 78.4 # Relative dielectric constant of water
TEMPERATURE = 300.0 # Temperature in Kelvin
BJERRUM_LENGTH = constants.elementary_charge**2 / constants.nano / \
(4 * np.pi * constants.epsilon_0 * EPSILON_R * constants.Boltzmann * TEMPERATURE)
# BERRUM_LENGTH of water at room temperature is 0.71 nm; electrostatic prefactor passed to P3M KBT/e2
# Lennard-Jones parameters
LJ_SIGMA = 0.3 # Particle size in nanometers
LJ_EPSILON = 1.0
CONCENTRATION = 1e-2 # desired concentration in mol/l
DISTANCE = 10 # 10 Debye lengths
N_IONPAIRS = 500
POTENTIAL_DIFF = 5.0
# Elementary charge
q = 1.
types = {"Cation": 0, "Anion": 1, "Electrodes": 2}
charges = {"Cation": q, "Anion": -q}
```

We want to make use of the optimal performance of **ESPResSo** in this tutorial,
which is roughly at 1000 particles/core.
Thus, we fixed above the number of ion pairs to `N_IONPAIRS = 500`

.

To be able to employ the analytical solution for the single plate also for the
double layer capacitor setup, the two electrodes need to be sufficiently far
away such that the additivity of the two surface potentials holds. In practice,
a separation of $d=10\lambda_\mathrm{D}$ is a good choice, represented by
`DISTANCE = 10`

.

Our choice of $c=10\,\mathrm{mmol}$ is a compromise between a sufficiently low
concentration for the PB theory to hold and not too large distances $d$ such
that the equilibration/diffusion of the ions is sufficiently fast
(`CONCENTRATION = 1e-2`

).

Note that in order to obtain results that we can interpret easily, we explicitly set a unit system using nanometers as length-scale above. The corresponding ion size of about 0.3 nm is a typical value for a simple salt; this, however, is in sharp contrast to the mean-field assumption of point-like ions. The latter are not easily studied within Molecular Dynamics simulations due to the required small time steps and are better suited for Monte Carlo type simulations. We instead focus here on analyzing deviations from PB theory due to the finite ion size.

- write a function
`get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS)`

that returns the lateral and normal box lengths`box_l_xy`

and`box_l_z`

(in nanometers) for the given parameters.

**Hint:** To account for the finite ion size and the wall interaction it is
useful to define the effective separation $d^\prime = d-2\sigma$, such that the
concentration is $\rho = N/(A \cdot d^\prime)$.

In [3]:

```
def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):
"""
For a given number of particles, determine the lateral area of the box
to match the desired concentration.
"""
# concentration is in mol/l, convert to 1/sigma**3
rho = concentration * (constants.Avogadro / constants.liter) * constants.nano**3
debye_length = (4. * np.pi * BJERRUM_LENGTH * rho * 2.)**(-1. / 2.) # desired Debye length in nm
l_z = distance * debye_length
box_volume = n_ionpairs / rho
area = box_volume / (l_z - 2. * LJ_SIGMA) # account for finite ion size in density calculation
l_xy = np.sqrt(area)
return l_xy, l_z
```

In [4]:

```
box_l_xy, box_l_z = get_box_dimension(CONCENTRATION, DISTANCE, N_IONPAIRS)
# useful quantities for the following calculations
DEBYE_LENGTH = box_l_z / DISTANCE # in units of nm
rho = N_IONPAIRS / (box_l_xy * box_l_xy * box_l_z) # in units of 1/nm^3
```

We now can create the **ESPResSo** system.

Note that for ELC to work properly, we need to add a gap of `ELC_GAP`

in the
non-periodic direction.
The precise value highly affects the performance due to the tuning of the P3M
electrostatic solver.
For $d=10\lambda$ a gap value of $6 L_z$ is a good choice.

We also set the time-step $dt = 0.01 \tau$, which is limited by the choice of $\sigma$ and $\tau$ in the repulsive WCA interaction.

In [5]:

```
ELC_GAP = 6. * box_l_z
system = espressomd.System(box_l=[box_l_xy, box_l_xy, box_l_z + ELC_GAP])
system.time_step = 0.01
```

- add two wall constraints at $z=0$ and $z=L_z$ to stop particles from
crossing the boundaries and model the electrodes.
Refer to
espressomd.constraints.ShapeBasedConstraint
and its
wall constraint
in the documentation to set up constraints and the
`types`

dictionary for the particle type.

In [6]:

```
# Bottom wall, normal pointing in the +z direction
floor = espressomd.shapes.Wall(normal=[0, 0, 1])
c1 = system.constraints.add(
particle_type=types["Electrodes"], penetrable=False, shape=floor)
# Top wall, normal pointing in the -z direction
ceiling = espressomd.shapes.Wall(normal=[0, 0, -1], dist=-box_l_z)
c2 = system.constraints.add(
particle_type=types["Electrodes"], penetrable=False, shape=ceiling)
```

In [7]:

```
offset = LJ_SIGMA # avoid unfavorable overlap at close distance to the walls
init_part_btw_z1 = offset
init_part_btw_z2 = box_l_z - offset
ion_pos = np.empty((3), dtype=float)
for i in range (N_IONPAIRS):
ion_pos[0] = rng.random(1) * system.box_l[0]
ion_pos[1] = rng.random(1) * system.box_l[1]
ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1
system.part.add(pos=ion_pos, type=types["Cation"], q=charges["Cation"])
for i in range (N_IONPAIRS):
ion_pos[0] = rng.random(1) * system.box_l[0]
ion_pos[1] = rng.random(1) * system.box_l[1]
ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1
system.part.add(pos=ion_pos, type=types["Anion"], q=charges["Anion"])
```

- For excluded volume interactions, add a WCA potential.

Refer to the documentation to set up the WCA interaction under Non-bonded section.

In [8]:

```
for t1 in types.values():
for t2 in types.values():
system.non_bonded_inter[t1, t2].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)
```

For the (2D+h) electrostatic with dielectrics we choose the ELC-IC with P3M.

Refer the documentation to set up ELCIC with P3M under the electrostatics section.

As later we will study different potential drops between the electrodes, write a
function that sets up the electrostatic solver for a given value
`POTENTIAL_DIFF.`

This function will take care of tuning the P3M and ELC parameters.
For our purposes, an accuracy of $10^{-3}$ is sufficient.

- Write a function
`setup_electrostatic_solver(potential_diff)`

that returns the ELC instance.

In [9]:

```
def setup_electrostatic_solver(potential_diff):
delta_mid_top = -1. # (Fully metallic case both -1)
delta_mid_bot = -1.
p3m_accuracy = 1e-3
elc_accuracy = 1e-3
p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,
accuracy=p3m_accuracy,
mesh=[12, 12, 48], # pinned for tuning reproducibility
cao=4, # pinned for tuning reproducibility
tune=True,
verbose=False)
elc = espressomd.electrostatics.ELC(actor=p3m,
gap_size=ELC_GAP,
const_pot=True,
pot_diff=potential_diff,
maxPWerror=elc_accuracy,
delta_mid_bot=delta_mid_bot,
delta_mid_top=delta_mid_top)
return elc
```

In [10]:

```
# Set up steepest descent integration
system.integrator.set_steepest_descent(f_max=0., gamma=30., max_displacement=0.01)
system.integrator.run(100)
particle_forces = np.linalg.norm(system.part.all().f, axis=1)
assert np.max(particle_forces) < 1.
```

We now switch to a Velocity Verlet integrator and set up a Langevin thermostat. Note, that we only analyze static properties, thus the damping and temperature chosen here only determine the simulation time towards the equilibrium distribution.

In [11]:

```
# Equilibration parameters
system.time_step = 0.003 # this time step limits particle movement to ~5% sigma per integration step
system.integrator.set_vv()
system.thermostat.set_langevin(kT=1., gamma=2., seed=42)
system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)
times = []
e_kin = []
for i in tqdm.trange(100):
energy = system.analysis.energy()
e_kin.append(energy['kinetic'])
times.append(system.time)
system.integrator.run(10)
```

100%|██████████| 100/100 [00:18<00:00, 5.51it/s]

In [12]:

```
# Plot the convergence of the total energy
plt.figure(figsize=(10, 6))
plt.plot(times, len(e_kin) * [N_IONPAIRS * 2 * 3. / 2.], label='heat bath')
plt.plot(times, e_kin, label='kinetic energy')
plt.xlabel('Simulation time')
plt.ylabel('Energy')
plt.legend()
plt.show()
```

Convergence after $t \sim 5$ time units.

Now we let ions diffuse to the electrodes.

In [13]:

```
N_SAMPLES_EQUIL = 20
STEPS_PER_EQUIL = 200
system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)
for tm in tqdm.trange(N_SAMPLES_EQUIL):
system.integrator.run(STEPS_PER_EQUIL)
```

100%|██████████| 20/20 [00:59<00:00, 2.95s/it]

We now need to set up an espressomd.observables.DensityProfile observable to calculate the anion and cation density profiles.

The time average is obtained through a espressomd.accumulators.MeanVarianceCalculator.

- Write a function
`setup_densityprofile_accumulators(bin_width)`

that returns the`bin_centers`

and the accumulators for both ion species in the $z$-range $[0,d]$. Since we are not estimating errors in this tutorial, the choice of`delta_N`

is rather arbitrary and does not affect the results. In practice, a typical value is`delta_N=20`

.

In [14]:

```
def setup_densityprofile_accumulators(bin_width):
cations = system.part.select(type=types["Cation"])
anions = system.part.select(type=types["Anion"])
n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))
density_profile_cation = espressomd.observables.DensityProfile(
ids=cations.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,
max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)
density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(
obs=density_profile_cation, delta_N=20)
density_profile_anion = espressomd.observables.DensityProfile(
ids=anions.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,
max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)
density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(
obs=density_profile_anion, delta_N=20)
zs = density_profile_anion.bin_centers()[0, 0, :, 2]
return zs, density_accumulator_cation, density_accumulator_anion
```

In [15]:

```
zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(
bin_width=DEBYE_LENGTH / 10.)
```

Now we take some measurement sampling the density profiles.

In [16]:

```
N_SAMPLES_PROD = 20
STEPS_PER_SAMPLE = 200
# Add the accumulators
system.auto_update_accumulators.add(density_accumulator_cation)
system.auto_update_accumulators.add(density_accumulator_anion)
for tm in tqdm.trange(N_SAMPLES_PROD):
system.integrator.run(STEPS_PER_SAMPLE)
cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]
anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]
```

100%|██████████| 20/20 [00:57<00:00, 2.89s/it]

Since we assume additivity, the total ion density follows from $$ \rho (z) = \rho_+(z) - \rho_+ (d-z) + \rho_-(z) - \rho_-(d-z) .$$

In [17]:

```
def gouy_chapman_potential(x, debye_length, phi_0):
kappa = 1. / debye_length
return 2. * np.log((1. + np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))) \
/ (1. - np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))))
def gouy_chapman_density(x, c0, debye_length, phi_0):
phi = gouy_chapman_potential(x, debye_length, phi_0)
return c0 / 2. * np.exp(-phi)
```

In [18]:

```
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(16, 4), nrows=1, ncols=3, sharey=True)
fig.subplots_adjust(wspace=0)
x = np.linspace(LJ_SIGMA, box_l_z-LJ_SIGMA, 100)
anion_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \
+ gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2.
cation_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) \
+ gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.))/2.
ax1.set_title('Cation')
ax2.set_title('Anion')
ax3.set_title('Total')
ax1.plot(x, cation_profile_analytic, label='analytic')
ax2.plot(x, anion_profile_analytic, label='analytic')
ax3.plot(x, cation_profile_analytic + anion_profile_analytic, label='analytic')
ax1.plot(zs[1:-1], cation_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax2.plot(zs[1:-1], anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax3.plot(zs[1:-1], cation_profile_mean[1:-1] + anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax1.legend(loc='upper center')
ax2.legend(loc='upper center')
ax3.legend(loc='upper center')
ax2.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax1.set_ylabel(r'$\rho(z)\,\mathrm{[mol/l]}$')
plt.show()
```

In [19]:

```
sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:len(zs)//2]) * (zs[1] - zs[0])
sigma_right = np.sum((cation_profile_mean-anion_profile_mean)[len(zs)//2:]) * (zs[1] - zs[0])
def grahame_sigma(phi):
return np.sinh(phi / 4.) * np.sqrt(2. * rho / (np.pi * BJERRUM_LENGTH))
sigma_grahame = grahame_sigma(POTENTIAL_DIFF)
print(f'simulation: {sigma_right:.3f} e/nm^2')
print(f'grahame: {sigma_grahame:.3f} e/nm^2')
print(f'relative deviation: {abs(1. - sigma_right/sigma_grahame) * 100.:.0f}%')
```

simulation: 0.076 e/nm^2 grahame: 0.117 e/nm^2 relative deviation: 35%

In [20]:

```
# plot the electric field
fig, ax = plt.subplots(figsize=(10, 6))
dz_SI = (zs[1] - zs[0]) * constants.nano
chargedensity = (cation_profile_mean - anion_profile_mean) * constants.elementary_charge / constants.nano**3
E_SI = 1. / (EPSILON_R * constants.epsilon_0) * np.cumsum(chargedensity * dz_SI)
# integration constant: zero field in the center
E_SI -= E_SI.min()
E = E_SI / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE) / constants.nano)
ax2 = plt.twinx()
ax.plot(zs, E_SI)
ax2.plot(zs, E)
ax.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax.set_ylabel(r'$E_\mathrm{ind}\,\mathrm{[V/m]}$')
ax2.set_ylabel(r'$E_\mathrm{ind}\,\mathrm{[(k_\mathrm{B}T/e)/nm]}$')
plt.show()
```

The electric potential can be calculated from $\phi(z) = \int_0^z -E(z^\prime)\,\mathrm{d}z^\prime$.

In [21]:

```
# plot the elecrostatic potential
fig, ax = plt.subplots(figsize=(10, 6))
ax2 = ax.twinx()
phi_SI = -np.cumsum(E_SI * dz_SI)
phi = phi_SI * (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))
ax.plot(zs, phi_SI)
ax.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax.set_ylabel(r'$\phi\,[V]$')
ax2.set_ylabel(r'$\phi\,[k_\mathrm{B}T/e]$')
ax2.axhline(-5, ls='--', color='k')
ax2.axhline(0, ls='--', color='k')
ax.set_xlim(0, 10. * DEBYE_LENGTH)
plt.show()
```

In [22]:

```
measured_potential_difference = -(phi[-1] - phi[0])
print(f'applied voltage: {POTENTIAL_DIFF:.2f} V')
print(f'measured voltage: {measured_potential_difference:.2f} V')
print(f'relative deviation: {abs(1. - measured_potential_difference / POTENTIAL_DIFF) * 100:.0f}%')
```

applied voltage: 5.00 V measured voltage: 4.79 V relative deviation: 4%

With the above knowledge, we can now assess the differential capacitance of the system, by changing the applied voltage difference and determining the corresponding surface charge density.

In [23]:

```
sigma_vs_phi = []
MIN_PHI = 0.5
MAX_PHI = 10
N_PHI = 7
N_SAMPLES_EQUIL_CAP = 15
N_SAMPLES_CAP = 5
# sample from high to low potential to improve sampling
for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):
system.auto_update_accumulators.clear()
system.electrostatics.solver = setup_electrostatic_solver(potential_diff)
system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)
sigmas = []
for tm in range(N_SAMPLES_CAP):
zs, density_accumulator_cation, density_accumulator_anion = \
setup_densityprofile_accumulators(bin_width=DEBYE_LENGTH / 10.)
system.auto_update_accumulators.clear()
system.auto_update_accumulators.add(density_accumulator_cation)
system.auto_update_accumulators.add(density_accumulator_anion)
system.integrator.run(STEPS_PER_SAMPLE)
cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]
anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]
sigmas.append(np.sum((cation_profile_mean - anion_profile_mean)[:int(len(zs) / 2.)]) * (zs[1] - zs[0]))
sigma_vs_phi.append([potential_diff, np.mean(sigmas), scipy.stats.sem(sigmas)])
```

100%|██████████| 7/7 [06:36<00:00, 56.66s/it]

In [24]:

```
fig, ax = plt.subplots(figsize=(10, 6))
sigma_vs_phi = np.array(sigma_vs_phi)
x = np.linspace(0, sigma_vs_phi[:,0].max())
phi_SI = sigma_vs_phi[:,0] / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))
plt.errorbar(-sigma_vs_phi[:,1] * constants.elementary_charge / constants.nano**2,
phi_SI, xerr=sigma_vs_phi[:,2] * scipy.constants.elementary_charge / scipy.constants.nano**2,
fmt='o',label='Simulation')
plt.plot(grahame_sigma(x) * constants.elementary_charge / constants.nano**2,
x / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE)), label='Grahame')
x = np.linspace(0, ax.get_ylim()[1])
plt.plot(EPSILON_R * constants.epsilon_0 * x / 2. / (DEBYE_LENGTH * constants.nano), x, label='linear PB',
ls='--')
plt.xlabel(r'$\sigma\,\mathrm{[C/m^2]}$')
plt.ylabel(r'$\phi_\mathrm{s}\,\mathrm{[V]}$')
plt.legend()
plt.show()
```

[1] Conway, B. E. Electrochemical Supercapacitors; Springer US: Boston, MA, 1999. https://doi.org/10.1007/978-1-4757-3058-6.

[2] Tyagi, S.; Arnold, A.; Holm, C. Electrostatic Layer Correction with Image Charges: A Linear Scaling Method to Treat Slab 2D+h Systems with Dielectric Interfaces. J. Chem. Phys. 2008, 129 (20), 204102. https://doi.org/10.1063/1.3021064.

[3] Gouy, G. Constitution of the Electric Charge at the Surface of an Electrolyte. J. Phys. 1910, 9 (4), 457–467.

[4] Chapman, D. L. A Contribution to the Theory of Electrocapillarity. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1913, 25 (148), 475. https://doi.org/10.1080/14786440408634187.

[5] Grahame, D. C. The Electrical Double Layer and the Theory of Electrocapillarity. Chem. Rev. 1947, 41 (3), 441–501. https://doi.org/10.1021/cr60130a002.

[6] Tyagi, S.; Süzen, M.; Sega, M.; Barbosa, M.; Kantorovich, S. S.; Holm, C. An Iterative, Fast, Linear-Scaling Method for Computing Induced Charges on Arbitrary Dielectric Boundaries. J. Chem. Phys. 2010, 132 (15), 154112. https://doi.org/10.1063/1.3376011.