Basic simulation of electrodes in ESPResSo part I: ion-pair in a narrow metallic slit-like confinement using ICC$^\star$

Prerequisites

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

Introduction

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]

Theoretical Background

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.

Green's function for charges in a dielectric slab

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.

2D+h periodic systems, dielectric interfaces and Induced Charge Computation with ICC$^\star$

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.

1. System setup

We first import all ESPResSo features and external modules.

import matplotlib.pyplot as plt
import tqdm
import numpy as np

import espressomd
import espressomd.electrostatics
import espressomd.electrostatic_extensions

espressomd.assert_features(['ELECTROSTATICS', 'EXTERNAL_FORCES'])
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.

# 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 = 100

# 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])

Setup of electrostatic interactions

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

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

Task

• Set up ELC with p3m as its actor.

# SOLUTION CELL
elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR, check_neutrality=False)

Next, we set up the ICC particles on both electrodes

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

Task

• 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.

# SOLUTION CELL
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]))

Task

• Set elc as system.electrostatics.solver
• Create an ICC object and set it as system.electrostatics.extension

Hints

• 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.

# SOLUTION CELL
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

2. Calculation of the forces

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 [00:43<00:00,  4.32s/it]

3. Analysis and Interpretation of the data

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.

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

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()

References

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