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, as described in the
**ESPResSo**user guide

This tutorial introduces some basic concepts for simulating charges close to an
electrode interface using **ESPResSo**.
In this first part, we focus on the interaction of a single ion pair confined in
a narrow metallic slit pore using the ICC$^\star$-algorithm
[1] for the computation of the surface polarization.
Here, we verify the strong deviation from a Coulomb-like interaction:
In metallic confinement, the ion pair interaction energy is screened
exponentially due to the presence of induced charges on the slit walls.
[2]

The normal component of electric field across a surface dividing two dielectric media in presence of a surface charge density $\sigma$ is discontinuous and follows as $(\varepsilon_1\vec{E}_1 - \varepsilon_2\vec{E}_2).\hat{n}=-\sigma(\vec{r})$. This expression describes the jump in the electric field across the material interface going from a dielectric medium $\varepsilon_1$ to another one, $\varepsilon_2$.

While in the case of non-polarizable materials ($\varepsilon_1 = \varepsilon_2 = 1$), this jump is only related to surface charges and the potential is continuous across the interface, for polarizable materials, the polarization field $\vec{P}$ will also give a contribution. In order to solve for the electric field in presence of a jump of the dielectric constant across an interface, one must know the electric fields on both sides.

Another approach is to replace this two-domain problem by an equivalent one
without the explicit presence of the dielectric jump.
This is achieved by introducing an additional fictitious charge, i.e. an induced
charge density $\sigma_{\mathrm{ind}}$ on the surface.
With this well known "method of image charges", it is sufficient to know the
electric field on one side of the interface.
**ESPResSo** provides the "Induced Charge Calculation with fast Coulomb Solvers"
(ICC$^\star$) algorithm [1] which employs a numerical scheme for solving the boundary integrals and the induced charge.

*Note*: Apart from ICC$^\star$ that solves for image charges spatially resolved, **ESPResSo** offers the "Electrostatic layer
correction with image charges" (ELC-IC) method [3], for
planar 2D+h partially periodic systems with dielectric interfaces that accounts for laterally averaged surface charge.
The tutorial on *Basic simulation of electrodes in ESPResSo part II*
addresses this in detail.

The Green's function for two point charges in a dielectric slab can be solved analytically [2]. In the metallic limit ($\varepsilon_2 \to\infty$) the dielectric contrast is $$ \Delta = \frac{\varepsilon_1 - \varepsilon_2} {\varepsilon_1 + \varepsilon_2} = -1 .$$ If the ions are placed in the center of a slab of width $w$ and a distance $r$ away from each other, the Green's function accounting for all image charges simplifies to [4] $$ 4 \pi \varepsilon_0 \varepsilon_r w \mathcal{G}(x) = \sum_{n=-\infty}^\infty \frac{-1^n}{\sqrt{x^2+n^2}} ,$$ where we have introduced the scaled separation $x = r/w$.

For $x\to 0$ the term for $n = 0$ dominates and one recovers $$ \lim_{x\to 0} 4 \pi \varepsilon_0 \varepsilon_r w \mathcal{G}(x) = \frac{1}{x},$$ which is the classical Coulomb law. Contrary, for large distances $x \to \infty$ one finds $$ \lim_{x\to \infty} 4 \pi \varepsilon_0 \varepsilon_r w \mathcal{G}(x) = \sqrt{\frac{8}{x}} e^{-\pi x},$$ i.e. the interaction decays exponentially. Such screened electrostatic interactions might explain unusual features concerning the nano-confined ionic liquids employed for supercapacitors referred to 'super-ionic states'.

We will explore this interaction numerically using ICC$^\star$ in the following.

Partially periodic ionic systems with dielectric interfaces are very often encountered in the study of energy materials or bio-macromolecular and membrane studies. These systems usually exhibit a confinement along one ($z$) direction, where the confining boundary or interface imposes a dielectric discontinuity, while the other $x$-$y$ directions are treated periodic. To investigate such a 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) [3]. The latter corrects for the contributions from the periodic images in the constrained 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. The method relies on an empty vacuum slab in the simulation box in the $z$-direction perpendicular to slab. While in theory this can become relatively small (typically 10% of the box length), its choice in practice strongly affects the performance due to the tuning of the P3M parameters to obtain the desired accuracy.

We furthermore employ ICC$^\star$ to solve the Poisson equation for an inhomogeneous dieletric: $$ \nabla (\varepsilon \nabla \phi)=-4\pi \rho$$

The image charge formalism can be derived as follows:

- Integrate the latter expression at the boundary over an infinitesimally wide pillbox, which will give the induced surface charge in this infinitesimal segment as (Gauss law): $$q_{\mathrm{ind}} = \frac{1}{4\pi} \oint\, dA\, \cdot \varepsilon\nabla \phi = \frac{A}{4\pi}(\varepsilon_1\vec{E}_1 \cdot \hat{n}-\varepsilon_2\vec{E}_2 \cdot\hat{n})$$
- The electric field in region 1 at the closest proximity of the interface, $\vec{E}_{1}$, can be written as a sum of electric field contributions from the surface charge $\sigma$ and the external electric field $\vec{E}$: $$ \vec{E}_{1} =\vec{E} + 2\pi/\varepsilon_1\sigma\hat{n} $$
- Combining this with the previous expression, the induced charge can be written in terms of the dielectric mismatch $\Delta$ and the electric field as: $$\sigma = \frac{\varepsilon_1}{2\pi} \frac{\varepsilon_1-\varepsilon_2}{\varepsilon_1+\varepsilon_2}\vec{E} \cdot \hat{n} =: \frac{\varepsilon_1}{2\pi} \Delta \, \vec{E} \cdot \hat{n}$$

The basic idea of the ICC$^\star$ formalism now is to employ a discretization of the surface by means of spatially fixed ICC particles. The charge of each ICC particle is not fixed but rather iterated using the expressions for $\vec{E}_{1}$ and $\sigma$ above until a self-consistent solution is found.

We first import all ESPResSo features and external modules.

In [1]:

```
import matplotlib.pyplot as plt
import tqdm
import numpy as np
import espressomd
import espressomd.electrostatics
import espressomd.electrostatic_extensions
espressomd.assert_features(['ELECTROSTATICS'])
plt.rcParams.update({'font.size': 18})
```

We need to define the system dimensions and some physical parameters related to length, time and energy scales of our system. All physical parameters are defined in reduced units of length ($\sigma=1$; particle size), mass ($m=1$; particle mass), arbitrary time (we do not consider particle dynamics) and elementary charge ($e=1$).

Another important length scale is the Bjerrum Length, which is the length at which the electrostatic energy between two elementary charges is comparable to the thermal energy $k_\mathrm{B}T$. It is defined as $$\ell_\mathrm{B}=\frac{1}{4\pi\varepsilon_0\varepsilon_r k_\mathrm{B}T}.$$

In our case, if we choose the ion size ($\sigma$) in simulation units to represent a typical value for mono-atomic salt, 0.3 nm in real units, then the Bjerrum length of water at room temperature, $\ell_\mathrm{B}=0.71 \,\mathrm{nm}$ is $\ell_\mathrm{B}\sim 2$ in simulations units.

In [2]:

```
# Box dimensions
# To construct a narrow slit Lz << (Lx , Ly)
box_l_x = 100.
box_l_y = 100.
box_l_z = 5.
# Additional space for ELC
ELC_GAP = 6*box_l_z
system = espressomd.System(box_l=[box_l_x, box_l_y, box_l_z + ELC_GAP])
system.time_step = 0.01
system.cell_system.skin = 0.4
# Elementary charge
q = 1.0
# Interaction parameters for P3M with ELC
BJERRUM_LENGTH = 2.0 # Electrostatic prefactor passed to P3M ; prefactor=lB KBT/e2
ACCURACY = 1e-7 # P3M force accuracy
MAX_PW_ERROR = 1e-7 # maximum pairwise error in ELC
ICC_EPSILON_DOMAIN = 1. # epsilon inside the slit
ICC_EPSILON_WALLS = 1e5 # epsilon outside the slit. Very large to mimic metal
ICC_CONVERGENCE = 1e-3 # ICC numeric/performance parameters
ICC_RELAXATION = 0.95
ICC_MAX_ITERATIONS = 1000
# Lennard-Jones parameters
LJ_SIGMA = 1.0
LJ_EPSILON = 1.0
# Particle parameters
TYPES = {"Cation": 0, "Anion": 1 ,"Electrodes": 2}
charges = {"Cation": q, "Anion": -q }
# Test particles to measure forces
p1 = system.part.add(pos=[box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges["Cation"], fix=3*[True])
p2 = system.part.add(pos=[3.0*box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges["Anion"], fix=3*[True])
```

First, we define our 3D electrostatics solver (P3M)

In [3]:

```
p3m = espressomd.electrostatics.P3M(
prefactor=BJERRUM_LENGTH,
accuracy=ACCURACY,
check_neutrality=False,
mesh=[100, 100, 150],
cao=5
)
```

In [4]:

```
elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)
```

Next, we set up the ICC particles on both electrodes

In [5]:

```
ICC_PARTCL_NUMBER_DENSITY = 1.
icc_partcls_bottom = []
icc_partcls_top = []
```

- Using the (area) density of ICC particles defined in the cell above, calculate the x/y positions of the particles for a uniform, quadratic grid.
- Add fixed particles on the electrodes. Make sure to use the correct
`type`

. Give the top (bottom) plate a total charge of $+1$ ($-1$). - Store the created particles in lists
`icc_partcls_bottom`

,`icc_partcls_top`

.

In [6]:

```
line_density = np.sqrt(ICC_PARTCL_NUMBER_DENSITY)
xs = np.linspace(0, system.box_l[0], num=int(round(system.box_l[0] * line_density)), endpoint=False)
ys = np.linspace(0, system.box_l[1], num=int(round(system.box_l[1] * line_density)), endpoint=False)
n_partcls_each_electrode = len(xs) * len(ys)
# Bottom electrode
for x in xs:
for y in ys:
icc_partcls_bottom.append(system.part.add(pos=[x, y, 0.], q=-1. / n_partcls_each_electrode,
type=TYPES["Electrodes"], fix=3*[True]))
# Top electrode
for x in xs:
for y in ys:
icc_partcls_top.append(system.part.add(pos=[x, y, box_l_z], q=1. / n_partcls_each_electrode,
type=TYPES["Electrodes"], fix=3*[True]))
```

- Set
`elc`

as`system.electrostatics.solver`

- Create an ICC object and set it as
`system.electrostatics.extension`

- ICC variables are defined in the second code cell from the top.
- Make sure to not mark our test particles
`p1`

and`p2`

(with ids 0 and 1) as ICC particles.

In [7]:

```
system.electrostatics.solver = elc
n_icc_partcls = len(icc_partcls_top) + len(icc_partcls_bottom)
# Common area, sigma and metallic epsilon
icc_areas = 1. / ICC_PARTCL_NUMBER_DENSITY * np.ones(n_icc_partcls)
icc_sigmas = np.zeros(n_icc_partcls)
icc_epsilons = ICC_EPSILON_WALLS * np.ones(n_icc_partcls)
icc_normals = np.repeat([[0, 0, 1], [0, 0, -1]], repeats=n_icc_partcls//2, axis=0)
icc = espressomd.electrostatic_extensions.ICC(
first_id=min(system.part.select(type=TYPES["Electrodes"]).id),
n_icc=n_icc_partcls,
convergence=ICC_CONVERGENCE,
relaxation=ICC_RELAXATION,
ext_field=[0, 0, 0],
max_iterations=ICC_MAX_ITERATIONS,
eps_out=ICC_EPSILON_DOMAIN,
normals=icc_normals,
areas=icc_areas,
sigmas=icc_sigmas,
epsilons=icc_epsilons
)
system.electrostatics.extension = icc
```

In [8]:

```
N_AXIAL_POINTS = 10
r = np.logspace(0., box_l_z / 4., N_AXIAL_POINTS)
elc_forces_axial = np.empty((N_AXIAL_POINTS, 2))
n_icc_per_electrode = len(icc_partcls_top)
p1.pos = [0., box_l_y / 2., box_l_z / 2.]
for i in tqdm.trange(N_AXIAL_POINTS):
p2.pos = [r[i], box_l_y / 2., box_l_z / 2.]
system.integrator.run(0, recalc_forces=True)
elc_forces_axial[i, 0] = p1.f[0]
elc_forces_axial[i, 1] = p2.f[0]
# reset ICC charges to ensure charge neutrality
for part in icc_partcls_top:
part.q = 1. / n_icc_per_electrode
for part in icc_partcls_bottom:
part.q = -1. / n_icc_per_electrode
```

100%|██████████| 10/10 [01:00<00:00, 6.08s/it]

With this we can now compare the force between the two ions to the analytical prediction. To evaluate the infinite series we truncate at $n=1000$, which already is well converged.

In [9]:

```
def analytic_force_centered(r,w):
def summand(x, n):
return (-1)**n * x / (x**2 + n**2)**(3. / 2.)
def do_sum(x):
limit = 1000
accumulator = 0.
for n in range(-limit + 1, limit + 1):
accumulator += summand(x, n)
return accumulator
x = r / w
prefactor = BJERRUM_LENGTH / w**2
F = do_sum(x) * prefactor
return F
def coulomb_force(x):
prefactor = BJERRUM_LENGTH
E = prefactor / x**2
return E
def exponential_force(r,w):
x = r / w
prefactor = BJERRUM_LENGTH
E = prefactor * np.sqrt(2.) * (1. / x)**(3. / 2.) * np.exp(-np.pi * x)
return E
```

In [10]:

```
fig = plt.figure(figsize=(10, 6))
x = np.logspace(-0.25, 1.45, 100)
plt.plot(x / BJERRUM_LENGTH, analytic_force_centered(x, box_l_z), color='b', label="analytic", marker='')
plt.plot(x / BJERRUM_LENGTH, coulomb_force(x), color='g', ls='--', label='Coulomb')
plt.plot(x / BJERRUM_LENGTH, exponential_force(x, box_l_z), color='r', ls='--', label='Exponential')
plt.plot(r / BJERRUM_LENGTH, -elc_forces_axial[:,1], color='r', label="sim (p2)", marker='o', ls='')
plt.plot(r / BJERRUM_LENGTH, +elc_forces_axial[:,0], color='k', label="sim (p1)", marker='x', ls='')
plt.xlabel(r'$r \, [\ell_\mathrm{B}]$')
plt.ylabel(r'$F \, [k_\mathrm{B}T / \sigma$]')
plt.loglog()
plt.legend()
plt.show()
```

[1] 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.

[2] Kondrat, S.; Feng, G.; Bresme, F.; Urbakh, M.; Kornyshev, A. A. Theory and Simulations of Ionic Liquids in Nanoconfinement. Chem. Rev. 2023, 123 (10), 6668–6715. https://doi.org/10.1021/acs.chemrev.2c00728.

[3] 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.

[4] Loche, P.; Ayaz, C.; Wolde-Kidan, A.; Schlaich, A.; Netz, R. R. Universal and Nonuniversal Aspects of Electrostatics in Aqueous Nanoconfinement. J. Phys. Chem. B 2020, 124 (21), 4365–4371. https://doi.org/10.1021/acs.jpcb.0c01967.