To work with this tutorial, you should be familiar with the following topics:
lennard_jones
folder.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:
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.
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.
# 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)
p3m = espressomd.electrostatics.P3M(
prefactor=BJERRUM_LENGTH,
accuracy=ACCURACY,
check_neutrality=False,
mesh=[100, 100, 150],
cao=5
)
Task
p3m
as its actor.elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)
Next, we set up the ICC particles on both electrodes
ICC_PARTCL_NUMBER_DENSITY = 1.
icc_partcls_bottom = []
icc_partcls_top = []
Task
type
. Give the top (bottom) plate a total charge of $+1$ ($-1$). icc_partcls_bottom
, icc_partcls_top
.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
elc
as system.electrostatics.solver
system.electrostatics.extension
Hints
p1
and p2
(with ids 0 and 1) as ICC particles.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
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:34<00:00, 3.49s/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.
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()
[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.