diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt
index 508bda4caf..479f3491e5 100644
--- a/doc/tutorials/CMakeLists.txt
+++ b/doc/tutorials/CMakeLists.txt
@@ -114,6 +114,7 @@ add_subdirectory(visualization)
add_subdirectory(ferrofluid)
add_subdirectory(constant_pH)
add_subdirectory(widom_insertion)
+add_subdirectory(electrodes)
configure_file(Readme.md ${CMAKE_CURRENT_BINARY_DIR} COPYONLY)
configure_file(convert.py ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/doc/tutorials/electrodes/CMakeLists.txt b/doc/tutorials/electrodes/CMakeLists.txt
new file mode 100644
index 0000000000..4ba4352845
--- /dev/null
+++ b/doc/tutorials/electrodes/CMakeLists.txt
@@ -0,0 +1,26 @@
+#
+# Copyright (C) 2020-2022 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+configure_tutorial_target(TARGET tutorial_electrodes DEPENDS
+ electrodes_part1.ipynb electrodes_part2.ipynb)
+
+nb_export(TARGET tutorial_electrodes SUFFIX "1" FILE "electrodes_part1.ipynb"
+ HTML_RUN)
+nb_export(TARGET tutorial_electrodes SUFFIX "2" FILE "electrodes_part2.ipynb"
+ HTML_RUN)
diff --git a/doc/tutorials/electrodes/NotesForTutor.md b/doc/tutorials/electrodes/NotesForTutor.md
new file mode 100644
index 0000000000..47497b59da
--- /dev/null
+++ b/doc/tutorials/electrodes/NotesForTutor.md
@@ -0,0 +1,32 @@
+# Simulations of electrodes in ESPResSO
+
+## Physics learning goals
+
+### Part 1
+
+* Give a short recap about image charges, dielectric media, ...
+* Nano-confinement can exhibit a broad variety of interesting effects that can
+ be studied with computer simulations!
+* Electrostatics in nano-confinement: concept of a Green's function
+* Discrete image charges: ICC\*
+
+### Part 2
+
+* Nano-confined charged liquids as super-capacitors
+* Advanced Poisson-Boltzmann theory: Gouy-Chapman, Graham equation
+* Limits of PB: finite ion size, correlations, ...
+* Coarse-grained view: Surface charge density via ELC-IC
+* How to apply a potential difference in the simulation.
+
+After the tutorial, students should be able to:
+
+* Explain how ESPResSo implements 2D periodic electrostatics.
+* What are the limitations of the mean-field PB description.
+* How to evaluate the differential capacitance from simulations.
+* The basic idea of super-ionic states.
+
+## ESPResSo learning goals
+
+* Setting up ELC, ICC and ELC-IC simulations
+* Why is choosing the ELC-gap a crucial parameter?
+* Setting up observables and accumulators for the density profiles.
diff --git a/doc/tutorials/electrodes/electrodes_part1.ipynb b/doc/tutorials/electrodes/electrodes_part1.ipynb
new file mode 100644
index 0000000000..acda4631df
--- /dev/null
+++ b/doc/tutorials/electrodes/electrodes_part1.ipynb
@@ -0,0 +1,572 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Basic simulation of electrodes in ESPResSo part I: ion-pair in a narrow metallic slit-like confinement using ICC$^\\star$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Prerequisites\n",
+ "\n",
+ "To work with this tutorial, you should be familiar with the following topics:\n",
+ "\n",
+ "- Setting up and running simulations in ESPResSo - creating particles,\n",
+ " incorporating interactions.\n",
+ " If you are unfamiliar with this, you can go through the tutorial\n",
+ " in the `lennard_jones` folder.\n",
+ "- Basic knowledge of classical electrostatics:\n",
+ " dipoles, surface and image charges\n",
+ "- Reduced units, as described in the **ESPResSo** [user guide](https://espressomd.github.io/doc/introduction.html#on-units)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Introduction\n",
+ "\n",
+ "This tutorial introduces some basic concepts for simulating charges close to an\n",
+ "electrode interface using **ESPResSo**.\n",
+ "In this first part, we focus on the interaction of a single ion pair confined in\n",
+ "a narrow metallic slit pore using the ICC$^\\star$-algorithm\n",
+ "[1] for the computation of the surface polarization.\n",
+ "Here, we verify the strong deviation from a Coulomb-like interaction:\n",
+ "In metallic confinement, the ion pair interaction energy is screened\n",
+ "exponentially due to the presence of induced charges on the slit walls.\n",
+ "[2] "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Theoretical Background \n",
+ "\n",
+ "The normal component of electric field across a surface dividing two dielectric\n",
+ "media in presence of a surface charge density $\\sigma$ is discontinuous and follows as\n",
+ "$(\\varepsilon_1\\vec{E}_1 - \\varepsilon_2\\vec{E}_2).\\hat{n}=-\\sigma(\\vec{r})$.\n",
+ "This expression describes the jump in the electric field across the material\n",
+ "interface going from a dielectric medium $\\varepsilon_1$ to another one,\n",
+ "$\\varepsilon_2$.\n",
+ "\n",
+ "While in the case of non-polarizable materials ($\\varepsilon_1 = \\varepsilon_2 = 1$),\n",
+ "this jump is only related to surface charges and the\n",
+ "potential is continuous across the interface, for polarizable materials, the\n",
+ "polarization field $\\vec{P}$ will also give a contribution. \n",
+ "In order to solve for the electric field in presence of a jump of the dielectric constant\n",
+ "across an interface, one must know the electric fields on both sides. \n",
+ "\n",
+ "Another approach is to replace this two-domain problem by an equivalent one\n",
+ "without the explicit presence of the dielectric jump.\n",
+ "This is achieved by introducing an additional fictitious charge, i.e. an induced\n",
+ "charge density $\\sigma_{\\mathrm{ind}}$ on the surface. \n",
+ "With this well known \"method of image charges\", it is sufficient to know the\n",
+ "electric field on one side of the interface. \n",
+ "**ESPResSo** provides the \"Induced Charge Calculation with fast Coulomb Solvers\"\n",
+ "(ICC$^\\star$) algorithm [1] which employs a numerical scheme for solving the boundary integrals and the induced charge. \n",
+ "\n",
+ "*Note*: Apart from ICC$^\\star$ that solves for image charges spatially resolved, **ESPResSo** offers the \"Electrostatic layer\n",
+ "correction with image charges\" (ELC-IC) method [3], for\n",
+ "planar 2D+h partially periodic systems with dielectric interfaces that accounts for laterally averaged surface charge.\n",
+ "The tutorial on *Basic simulation of electrodes in ESPResSo part II*\n",
+ "addresses this in detail."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Green's function for charges in a dielectric slab\n",
+ "\n",
+ "The Green's function for two point charges in a dielectric slab can be solved\n",
+ "analytically [2].\n",
+ "In the metallic limit ($\\varepsilon_2 \\to\\infty$) the dielectric contrast is\n",
+ "$$ \\Delta = \\frac{\\varepsilon_1 - \\varepsilon_2} {\\varepsilon_1 + \\varepsilon_2} = -1 .$$\n",
+ "If the ions are placed in the center of a slab of width $w$ and a distance $r$\n",
+ "away from each other, the Green's function accounting for all image charges\n",
+ "simplifies to [4] \n",
+ "$$ 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\sum_{n=-\\infty}^\\infty \\frac{-1^n}{\\sqrt{x^2+n^2}} ,$$\n",
+ "where we have introduced the scaled separation $x = r/w$.\n",
+ "\n",
+ "For $x\\to 0$ the term for $n = 0$ dominates and one recovers\n",
+ "$$ \\lim_{x\\to 0} 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\frac{1}{x},$$\n",
+ "which is the classical Coulomb law.\n",
+ "Contrary, for large distances $x \\to \\infty$ one finds\n",
+ "$$ \\lim_{x\\to \\infty} 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\sqrt{\\frac{8}{x}} e^{-\\pi x},$$\n",
+ "i.e. the interaction decays exponentially.\n",
+ "Such screened electrostatic interactions might explain unusual features\n",
+ "concerning the nano-confined ionic liquids employed for supercapacitors referred\n",
+ "to 'super-ionic states'.\n",
+ "\n",
+ "We will explore this interaction numerically using ICC$^\\star$ in the following."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2D+h periodic systems, dielectric interfaces and Induced Charge Computation with ICC$^\\star$\n",
+ "\n",
+ "Partially periodic ionic systems with dielectric interfaces are very often\n",
+ "encountered in the study of energy materials or bio-macromolecular and membrane\n",
+ "studies. \n",
+ "These systems usually exhibit a confinement along one ($z$) direction, where the\n",
+ "confining boundary or interface imposes a dielectric discontinuity, while the\n",
+ "other $x$-$y$ directions are treated periodic. \n",
+ "To investigate such a partially periodic system, we combine the efficient scaling behavior\n",
+ "of three-dimensional mesh-based solvers (typically\n",
+ "$\\mathcal{O}(N \\log N)$ for P3M) with the Electrostatic Layer Correction (ELC)\n",
+ "[3].\n",
+ "The latter corrects for the contributions from the periodic images in the\n",
+ "constrained direction and its numerical cost grows linear with the number of\n",
+ "point charges $N$, hence the performance overall depends on the underlying 3D\n",
+ "Coulomb solver.\n",
+ "The method relies on an empty vacuum slab in the simulation box in the\n",
+ "$z$-direction perpendicular to slab.\n",
+ "While in theory this can become relatively small (typically 10% of the box\n",
+ "length), its choice in practice strongly affects the performance due to the\n",
+ "tuning of the P3M parameters to obtain the desired accuracy.\n",
+ "\n",
+ "We furthermore employ ICC$^\\star$ to solve the Poisson equation for an\n",
+ "inhomogeneous dieletric:\n",
+ "$$ \\nabla (\\varepsilon \\nabla \\phi)=-4\\pi \\rho$$\n",
+ "\n",
+ "The image charge formalism can be derived as follows:\n",
+ "- Integrate the latter expression at the boundary over an infinitesimally wide\n",
+ "pillbox, which will give the induced surface charge in this infinitesimal\n",
+ "segment as (Gauss law):\n",
+ "$$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})$$\n",
+ "- The electric field in region 1 at the closest proximity of the interface, $\\vec{E}_{1}$,\n",
+ "can be written as a sum of electric field contributions from the surface charge\n",
+ "$\\sigma$ and the external electric field $\\vec{E}$:\n",
+ "$$ \\vec{E}_{1} =\\vec{E} + 2\\pi/\\varepsilon_1\\sigma\\hat{n} $$\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:\n",
+ "$$\\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}$$\n",
+ "\n",
+ "The basic idea of the ICC$^\\star$ formalism now is to employ a discretization\n",
+ "of the surface by means of spatially fixed ICC particles.\n",
+ "The charge of each ICC particle is not fixed but rather iterated using the\n",
+ "expressions for $\\vec{E}_{1}$ and $\\sigma$ above until a self-consistent\n",
+ "solution is found."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1. System setup \n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We first import all ESPResSo features and external modules."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import tqdm\n",
+ "import numpy as np\n",
+ "\n",
+ "import espressomd\n",
+ "import espressomd.electrostatics\n",
+ "import espressomd.electrostatic_extensions\n",
+ "\n",
+ "espressomd.assert_features(['ELECTROSTATICS'])\n",
+ "plt.rcParams.update({'font.size': 18})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We need to define the system dimensions and some physical parameters related to\n",
+ "length, time and energy scales of our system.\n",
+ "All physical parameters are defined in reduced units of length ($\\sigma=1$;\n",
+ "particle size), mass ($m=1$; particle mass), arbitrary time (we do not consider particle dynamics) and\n",
+ "elementary charge ($e=1$).\n",
+ "\n",
+ "Another important length scale is the Bjerrum Length, which is the length at\n",
+ "which the electrostatic energy between two elementary charges is comparable to\n",
+ "the thermal energy $k_\\mathrm{B}T$.\n",
+ "It is defined as\n",
+ "$$\\ell_\\mathrm{B}=\\frac{1}{4\\pi\\varepsilon_0\\varepsilon_r k_\\mathrm{B}T}.$$\n",
+ "\n",
+ "In our case, if we choose the ion size ($\\sigma$) in simulation units to represent a\n",
+ "typical value for mono-atomic salt, 0.3 nm in real units, then the\n",
+ "Bjerrum length of water at room temperature, $\\ell_\\mathrm{B}=0.71 \\,\\mathrm{nm}$ is\n",
+ "$\\ell_\\mathrm{B}\\sim 2$ in simulations units."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Box dimensions\n",
+ "# To construct a narrow slit Lz << (Lx , Ly)\n",
+ "box_l_x = 100.\n",
+ "box_l_y = 100.\n",
+ "box_l_z = 5.\n",
+ "\n",
+ "# Additional space for ELC\n",
+ "ELC_GAP = 6*box_l_z\n",
+ "\n",
+ "system = espressomd.System(box_l=[box_l_x, box_l_y, box_l_z + ELC_GAP])\n",
+ "\n",
+ "system.time_step = 0.01\n",
+ "system.cell_system.skin = 0.4\n",
+ "\n",
+ "# Elementary charge \n",
+ "q = 1.0 \n",
+ "\n",
+ "# Interaction parameters for P3M with ELC\n",
+ "BJERRUM_LENGTH = 2.0 # Electrostatic prefactor passed to P3M ; prefactor=lB KBT/e2 \n",
+ "ACCURACY = 1e-7 # P3M force accuracy \n",
+ "MAX_PW_ERROR = 1e-7 # maximum pairwise error in ELC\n",
+ "ICC_EPSILON_DOMAIN = 1. # epsilon inside the slit\n",
+ "ICC_EPSILON_WALLS = 1e5 # epsilon outside the slit. Very large to mimic metal\n",
+ "ICC_CONVERGENCE = 1e-3 # ICC numeric/performance parameters\n",
+ "ICC_RELAXATION = 0.95\n",
+ "ICC_MAX_ITERATIONS = 1000\n",
+ "\n",
+ "# Lennard-Jones parameters\n",
+ "LJ_SIGMA = 1.0\n",
+ "LJ_EPSILON = 1.0 \n",
+ "\n",
+ "# Particle parameters\n",
+ "TYPES = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n",
+ "charges = {\"Cation\": q, \"Anion\": -q }\n",
+ "\n",
+ "# Test particles to measure forces\n",
+ "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])\n",
+ "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])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Setup of electrostatic interactions\n",
+ "First, we define our 3D electrostatics solver (P3M)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p3m = espressomd.electrostatics.P3M(\n",
+ " prefactor=BJERRUM_LENGTH,\n",
+ " accuracy=ACCURACY,\n",
+ " check_neutrality=False,\n",
+ " mesh=[100, 100, 150],\n",
+ " cao=5\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "### Task\n",
+ "\n",
+ "* Set up [ELC](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.ELC) with ``p3m`` as its actor."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, we set up the ICC particles on both electrodes"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ICC_PARTCL_NUMBER_DENSITY = 1.\n",
+ "icc_partcls_bottom = []\n",
+ "icc_partcls_top = []"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "### TASK\n",
+ "\n",
+ "* 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. \n",
+ "* Add fixed particles on the electrodes. Make sure to use the correct ``type``. Give the top (bottom) plate a total charge of $+1$ ($-1$). \n",
+ "* Store the created particles in lists ``icc_partcls_bottom``, ``icc_partcls_top``."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "line_density = np.sqrt(ICC_PARTCL_NUMBER_DENSITY)\n",
+ "xs = np.linspace(0, system.box_l[0], num=int(round(system.box_l[0] * line_density)), endpoint=False)\n",
+ "ys = np.linspace(0, system.box_l[1], num=int(round(system.box_l[1] * line_density)), endpoint=False)\n",
+ "n_partcls_each_electrode = len(xs) * len(ys)\n",
+ "\n",
+ "# Bottom electrode\n",
+ "for x in xs:\n",
+ " for y in ys:\n",
+ " icc_partcls_bottom.append(system.part.add(pos=[x, y, 0.], q=-1. / n_partcls_each_electrode,\n",
+ " type=TYPES[\"Electrodes\"], fix=3*[True]))\n",
+ "# Top electrode\n",
+ "for x in xs:\n",
+ " for y in ys:\n",
+ " icc_partcls_top.append(system.part.add(pos=[x, y, box_l_z], q=1. / n_partcls_each_electrode,\n",
+ " type=TYPES[\"Electrodes\"], fix=3*[True]))\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "### Task\n",
+ "\n",
+ "* Set ``elc`` as ``system.electrostatics.solver``\n",
+ "* Create an [ICC object]((https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatic_extensions.ICC) and set it as ``system.electrostatics.extension``\n",
+ "\n",
+ "### Hints\n",
+ "\n",
+ "* ICC variables are defined in the second code cell from the top.\n",
+ "* Make sure to not mark our test particles ``p1`` and ``p2`` (with ids 0 and 1) as ICC particles."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "system.electrostatics.solver = elc\n",
+ "\n",
+ "n_icc_partcls = len(icc_partcls_top) + len(icc_partcls_bottom)\n",
+ "\n",
+ "# Common area, sigma and metallic epsilon\n",
+ "icc_areas = 1. / ICC_PARTCL_NUMBER_DENSITY * np.ones(n_icc_partcls)\n",
+ "icc_sigmas = np.zeros(n_icc_partcls)\n",
+ "icc_epsilons = ICC_EPSILON_WALLS * np.ones(n_icc_partcls)\n",
+ "icc_normals = np.repeat([[0, 0, 1], [0, 0, -1]], repeats=n_icc_partcls//2, axis=0)\n",
+ "\n",
+ "icc = espressomd.electrostatic_extensions.ICC(\n",
+ " first_id=min(system.part.select(type=TYPES[\"Electrodes\"]).id),\n",
+ " n_icc=n_icc_partcls,\n",
+ " convergence=ICC_CONVERGENCE,\n",
+ " relaxation=ICC_RELAXATION,\n",
+ " ext_field=[0, 0, 0],\n",
+ " max_iterations=ICC_MAX_ITERATIONS,\n",
+ " eps_out=ICC_EPSILON_DOMAIN,\n",
+ " normals=icc_normals,\n",
+ " areas=icc_areas,\n",
+ " sigmas=icc_sigmas,\n",
+ " epsilons=icc_epsilons\n",
+ ")\n",
+ "system.electrostatics.extension = icc\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2. Calculation of the forces"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N_AXIAL_POINTS = 10\n",
+ "r = np.logspace(0., box_l_z / 4., N_AXIAL_POINTS) \n",
+ "elc_forces_axial = np.empty((N_AXIAL_POINTS, 2))\n",
+ "n_icc_per_electrode = len(icc_partcls_top)\n",
+ "\n",
+ "p1.pos = [0., box_l_y / 2., box_l_z / 2.]\n",
+ "\n",
+ "for i in tqdm.trange(N_AXIAL_POINTS):\n",
+ " p2.pos = [r[i], box_l_y / 2., box_l_z / 2.]\n",
+ "\n",
+ " system.integrator.run(0, recalc_forces=True)\n",
+ " elc_forces_axial[i, 0] = p1.f[0]\n",
+ " elc_forces_axial[i, 1] = p2.f[0]\n",
+ " \n",
+ " # reset ICC charges to ensure charge neutrality \n",
+ " for part in icc_partcls_top:\n",
+ " part.q = 1. / n_icc_per_electrode\n",
+ " for part in icc_partcls_bottom:\n",
+ " part.q = -1. / n_icc_per_electrode"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 3. Analysis and Interpretation of the data\n",
+ "\n",
+ "With this we can now compare the force between the two ions to the analytical prediction.\n",
+ "To evaluate the infinite series we truncate at $n=1000$, which already is well converged."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def analytic_force_centered(r,w):\n",
+ " def summand(x, n):\n",
+ " return (-1)**n * x / (x**2 + n**2)**(3. / 2.)\n",
+ " \n",
+ " def do_sum(x):\n",
+ " limit = 1000\n",
+ " accumulator = 0.\n",
+ " for n in range(-limit + 1, limit + 1):\n",
+ " accumulator += summand(x, n)\n",
+ " return accumulator\n",
+ "\n",
+ " x = r / w\n",
+ " prefactor = BJERRUM_LENGTH / w**2\n",
+ " F = do_sum(x) * prefactor\n",
+ " return F\n",
+ "\n",
+ "def coulomb_force(x):\n",
+ " prefactor = BJERRUM_LENGTH\n",
+ " E = prefactor / x**2\n",
+ " return E\n",
+ "\n",
+ "def exponential_force(r,w):\n",
+ " x = r / w\n",
+ " prefactor = BJERRUM_LENGTH\n",
+ " E = prefactor * np.sqrt(2.) * (1. / x)**(3. / 2.) * np.exp(-np.pi * x)\n",
+ " return E"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig = plt.figure(figsize=(10, 6))\n",
+ "\n",
+ "x = np.logspace(-0.25, 1.45, 100)\n",
+ "plt.plot(x / BJERRUM_LENGTH, analytic_force_centered(x, box_l_z), color='b', label=\"analytic\", marker='')\n",
+ "plt.plot(x / BJERRUM_LENGTH, coulomb_force(x), color='g', ls='--', label='Coulomb')\n",
+ "plt.plot(x / BJERRUM_LENGTH, exponential_force(x, box_l_z), color='r', ls='--', label='Exponential')\n",
+ "\n",
+ "plt.plot(r / BJERRUM_LENGTH, -elc_forces_axial[:,1], color='r', label=\"sim (p2)\", marker='o', ls='')\n",
+ "plt.plot(r / BJERRUM_LENGTH, +elc_forces_axial[:,0], color='k', label=\"sim (p1)\", marker='x', ls='')\n",
+ "\n",
+ "plt.xlabel(r'$r \\, [\\ell_\\mathrm{B}]$')\n",
+ "plt.ylabel(r'$F \\, [k_\\mathrm{B}T / \\sigma$]')\n",
+ "plt.loglog()\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## References\n",
+ "\n",
+ "[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.\n",
+ " \n",
+ "[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.\n",
+ "\n",
+ "[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.\n",
+ "\n",
+ "[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."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb
new file mode 100644
index 0000000000..35f7723156
--- /dev/null
+++ b/doc/tutorials/electrodes/electrodes_part2.ipynb
@@ -0,0 +1,1088 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Basic simulation of electrodes in ESPResSo part II: Electrolyte capacitor and Poisson–Boltzmann theory"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Prerequisites\n",
+ "\n",
+ "To work with this tutorial, you should be familiar with the following topics:\n",
+ "\n",
+ "1. Setting up and running simulations in ESPResSo - creating particles,\n",
+ " incorporating interactions.\n",
+ " If you are unfamiliar with this, you can go through the tutorial\n",
+ " in the `lennard_jones` folder.\n",
+ "2. Basic knowledge of classical electrostatics:\n",
+ " dipoles, surface and image charges.\n",
+ "3. Reduced units and how to relate them to physical quantities, see the ESPResSo\n",
+ " [user guide](https://espressomd.github.io/doc/introduction.html#on-units)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Introduction\n",
+ "\n",
+ "Ionic liquid (IL) based capacitors have long been established as promising\n",
+ "candidates in the area of efficient energy storage devices due to their\n",
+ "extraordinary capacitance, which is why they are also termed super-capacitors.\n",
+ "Typically, in these setups a fluid consisting of mobile charge carriers is\n",
+ "placed between two electrodes and thus gets polarized upon application of an\n",
+ "external potential formic an electric double layer at the interfaces\n",
+ "[1].\n",
+ "Electric double-layer capacitors (EDLCs) can be constructed from electrodes of\n",
+ "various geometries and materials where energy is stored by potential driven\n",
+ "adsorption of counterions on the surface of the electrodes and forming the\n",
+ "double layer. \n",
+ "Thus, conducting, high surface area electrode materials can further maximize the\n",
+ "energy per volume.\n",
+ "\n",
+ "In this tutorial we are going to investigate the ionic layer formation between\n",
+ "two conducting dielectric walls in presence of an applied voltage using **ESPResSo**'s\n",
+ "\"Electrostatic layer correction with image charges\" (ELC-IC) method\n",
+ "[2].\n",
+ "We employ a primitive model of an aqueous salt solution, where the solvent is\n",
+ "treated implicitly as a homogeneous dielectric background.\n",
+ "Thus, within the limits of a mean-field treatment and for not too small slit\n",
+ "widths, we can compare our findings with the analytical Gouy–Chapmann\n",
+ "solution for a single charged plane since additivity is assumed to hold if the\n",
+ "potential of both walls is screened sufficiently.\n",
+ "While this mean-field formalism properly describes the behavior of Coulomb\n",
+ "fluids composed of monovalent ions at low concentrations in the vicinity of\n",
+ "weakly charged interfaces, for strongly charged systems, where correlation and\n",
+ "finite size effects begin to dominate, Poisson–Boltzmann theory falls\n",
+ "inadequate.\n",
+ "Our goal in this tutorial is to demonstrate how coarse grained implicit solvent\n",
+ "simulations can corroborate on some of the approximations and hint on\n",
+ "extensions/deviations.\n",
+ "\n",
+ "The inclusion of dielectric inhomogeneities appearing at the conducting\n",
+ "interfaces demands to take into account image effects that involve the full\n",
+ "solution of the Poisson equation.\n",
+ "This is dealt with in a computationally cost-effective way using the ELC-IC method to\n",
+ "treat the image charge effects in the presence of 2D dielectric bounding\n",
+ "interfaces. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Theoretical Background \n",
+ "\n",
+ "### Poisson–Boltzmann Theory\n",
+ "\n",
+ "Charged surfaces in contact with a liquid containing free charges (ions) attract\n",
+ "oppositely charged ions that form a diffusive electric double layer. \n",
+ "The competition between electrostatic interactions and entropy of ions in\n",
+ "solution determines the exact distribution of mobile ions close to charged\n",
+ "membranes.\n",
+ "Gouy [3] and Chapman [4] derived in the\n",
+ "early 20th century the analytic solution for the case of a single planar wall\n",
+ "within the mean-field approximation of the Poisson–Boltzmann (PB) equation. \n",
+ "We will use it to describe our two-electrode system, which is justified if the electrodes are \n",
+ "so far apart that one surface does not influence the ion distribution in front of the other surface.\n",
+ "\n",
+ "In the case of a monovalent electrolyte, double integrating the PB equation and\n",
+ "employing the corresponding boundary conditions for the charged plane located at\n",
+ "$z=0$ yields the electrostatic potential:\n",
+ "$$\\phi(z) = -2\\ln\\left[\n",
+ " \\frac{1-\\tanh(\\phi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D} z}}\n",
+ " {1+\\tanh(\\phi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D}z}} \\right].$$ \n",
+ "Here, $\\phi_\\mathrm{s}=\\phi(z=0)=$ const is the surface potential such\n",
+ "that $\\phi(z\\rightarrow \\infty)=0$.\n",
+ "$\\kappa_\\mathrm{D} = \\lambda_\\mathrm{D}^{-1}$ is the inverse Debye screening length given by\n",
+ "$$ \\lambda_\\mathrm{D} = \\left(\\frac{\\varepsilon \\, k_{\\mathrm B} T}{\\sum_{j = 1}^N n_j^0 \\, q_j^2}\\right)^{1/2}, $$\n",
+ "where $n_j^{(0)}$ and $q_j$ are the equilibrium number density and charge of the\n",
+ "$j$-th ion species.\n",
+ "For the monovalent salt this can conveniently be expressed in terms of the\n",
+ "Bjerrum length $\\ell_\\mathrm{B}$ and the equilibrium salt concentration\n",
+ "$\\rho^{(0)}=\\sum_j \\rho_j^{(0)}$,\n",
+ "$$ \\lambda_{\\mathrm D} = \\left(4 \\pi \\, \\ell_\\mathrm{B} \\, \\rho^{(0)}/e\\right)^{-1/2} . $$\n",
+ "\n",
+ "The cationic and anionic density profiles then follow from the Boltzmann\n",
+ "distribution as:\n",
+ "$$n_{\\pm}(z)=n_\\pm^{(0)}\\left(\\frac{1\\pm\\gamma e^{-\\kappa_\\mathrm{D}z}}\n",
+ " {1\\mp\\gamma e^{-\\kappa_\\mathrm{D}z}} \\right)^2$$\n",
+ "Here, $\\gamma$ is associated with the surface potential as\n",
+ "$\\phi_\\mathrm{s}=-4\\tanh^{-1}(\\gamma)$.\n",
+ "At large z, where the potential decays to zero, the ionic profiles tend to their\n",
+ "bulk (reservoir) densities, $n_\\pm(z\\to\\infty) = n_\\pm^{(0)}$\n",
+ "\n",
+ "The relation between the surface potential and the surface charge induced on the\n",
+ "electrode is given by the Grahame Equation [5] :\n",
+ "$$ \\sigma = \\sinh(\\phi_\\mathrm{s}/2) \\sqrt{\\frac{2 n_\\mathrm{b}}{\\pi \\ell_\\mathrm{B}}} $$\n",
+ "The latter expression thus yields the differential capacitance\n",
+ "$C=\\displaystyle\\frac{\\mathrm{d}\\sigma}{\\mathrm{d}\\phi_\\mathrm{s}}$ within the mean-field\n",
+ "solution for non-overlapping double layers."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## ELC-IC for 2D+h periodic systems with dielectric interfaces\n",
+ "\n",
+ "In this tutorial we employ a parallel plate capacitor setup with constant\n",
+ "potential boundary conditions which needs to be treated appropriately by the\n",
+ "electrostatics solver.\n",
+ "To simulate a two-dimensional partially periodic system, we combine the efficient\n",
+ "scaling behavior of three-dimensional mesh-based solvers (typically\n",
+ "$\\mathcal{O}(N \\log N)$ for P3M) with the Electrostatic Layer Correction (ELC)\n",
+ "[1].\n",
+ "The latter removes the contributions from the periodic images in the\n",
+ "non-periodic direction and its numerical cost grows linear with the number of\n",
+ "point charges $N$, hence the performance overall depends on the underlying 3D\n",
+ "Coulomb solver.\n",
+ "Furthermore, ELC can be extended straightforwardly to metallic boundary\n",
+ "conditions (or any other dielectric contrast) by using the method of image charges,\n",
+ "which is referred to as the “Electrostatic Layer Correction with Image Charges”\n",
+ "(ELC-IC) approach used in this tutorial.\n",
+ "\n",
+ "A voltage difference can be applied between the electrodes by following\n",
+ "considerations:\n",
+ "The total potential drop $\\Delta \\phi$ across the simulated system is readily\n",
+ "obtained from the ion distribution and integrating twice over the one-dimensional Poisson equation:\n",
+ "$$-\\varepsilon_{0}\\varepsilon_{r}\\phi_\\mathrm{ind}(z)=\\iint_{0}^{z}\\rho(z^{\\prime})(z-z^{\\prime})dz^{\\prime}$$\n",
+ "Here, the subscript 'ind' indicates that this is the potential due to the\n",
+ "induced inhomogeneous charge distribution.\n",
+ "In order to set up a constant potential difference $\\Delta \\phi$, a homogeneous\n",
+ "electric field is superimposed such that\n",
+ "$$ \\Delta \\phi = \\Delta \\phi_\\mathrm{ind} + \\Delta \\phi_\\mathrm{bat},$$\n",
+ "where $\\Delta \\phi_\\mathrm{bat}$ corresponds to the potential of a (virtually)\n",
+ "applied battery.\n",
+ "In practice, the linear electric field in $E_z^\\mathrm{(bat)}=-\\phi_\\mathrm{bat}/L_z$\n",
+ "in the $z$-direction normal to the surface that one needs to apply can be\n",
+ "calculated straightforwardly, as the corresponding contribution from the induced\n",
+ "charges is known:\n",
+ "$$ E_z^\\mathrm{(ind)} = \\frac{1}{\\varepsilon_0 \\varepsilon_r L^2 L_z} P_z$$\n",
+ "Here, $L$ denotes the lateral system size, $L_z$ the distance between the plates\n",
+ "and $P_z = \\sum_i q_i z_i$ is the total dipole moment in $z$-direction due to\n",
+ "the charges $q_i$.\n",
+ "Then, to maintain $\\Delta \\phi$, a force $F_z^\\mathrm{bat} = q E_z^\\mathrm{(bat)}$\n",
+ "is applied on all charges.\n",
+ "Since ELC already calculates $P_z$, the constant potential correction requires no\n",
+ "additional computational effort.\n",
+ "\n",
+ "*Note*: Apart from ELC-IC, **ESPResSo** also provides the ICC$^\\star$ method\n",
+ "[2] which employs an iterative numerical scheme with\n",
+ "discretized surface particles to solve the boundary integrals at the dieletcric\n",
+ "interface.\n",
+ "The tutorial on *Basic simulation of electrodes in ESPResSo part I* addresses this\n",
+ "in detail."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 1. System setup \n",
+ "\n",
+ "First we import all ESPResSo features and external modules"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import scipy.constants as constants\n",
+ "import scipy.stats\n",
+ "import tqdm\n",
+ "\n",
+ "import espressomd\n",
+ "import espressomd.electrostatics\n",
+ "import espressomd.electrostatic_extensions\n",
+ "import espressomd.observables\n",
+ "import espressomd.accumulators\n",
+ "import espressomd.shapes\n",
+ "\n",
+ "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n",
+ "rng = np.random.default_rng(42)\n",
+ "plt.rcParams.update({'font.size': 18})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We need to define system dimensions and some physical parameters related to\n",
+ "length, time and energy scales of our system.\n",
+ "As discussed in previous tutorials, all physical parameters are defined in terms\n",
+ "of a length $\\sigma$, mass $m$ and time $t$ and unit of charge $q$.\n",
+ "Since we are not explicitly interested in the dynamics of the system, we set the\n",
+ "mass to $m=1$ (particle mass).\n",
+ "For convenience, we choose the elementary charge as fundamental unit ($q=1e$)\n",
+ "and $\\sigma = 1 \\,\\mathrm{nm}$.\n",
+ "\n",
+ "With this we can now define the fundamental parameters of our system:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# water at room temperature\n",
+ "EPSILON_R = 78.4 # Relative dielectric constant of water\n",
+ "TEMPERATURE = 300.0 # Temperature in Kelvin\n",
+ "BJERRUM_LENGTH = constants.elementary_charge**2 / constants.nano / \\\n",
+ " (4 * np.pi * constants.epsilon_0 * EPSILON_R * constants.Boltzmann * TEMPERATURE)\n",
+ "# BERRUM_LENGTH of water at room temperature is 0.71 nm; electrostatic prefactor passed to P3M KBT/e2 \n",
+ "\n",
+ "# Lennard-Jones parameters\n",
+ "LJ_SIGMA = 0.3 # Particle size in nanometers\n",
+ "LJ_EPSILON = 1.0\n",
+ "\n",
+ "CONCENTRATION = 1e-2 # desired concentration in mol/l\n",
+ "DISTANCE = 10 # 10 Debye lengths\n",
+ "N_IONPAIRS = 500\n",
+ "\n",
+ "POTENTIAL_DIFF = 5.0\n",
+ "\n",
+ "# Elementary charge \n",
+ "q = 1.\n",
+ "types = {\"Cation\": 0, \"Anion\": 1, \"Electrodes\": 2}\n",
+ "charges = {\"Cation\": q, \"Anion\": -q}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 1.1 Setting up the box dimensions and create system\n",
+ "\n",
+ "We want to make use of the optimal performance of **ESPResSo** in this tutorial,\n",
+ "which is roughly at 1000 particles/core.\n",
+ "Thus, we fixed above the number of ion pairs to `N_IONPAIRS = 500`.\n",
+ "\n",
+ "To be able to employ the analytical solution for the single plate also for the\n",
+ "double layer capacitor setup, the two electrodes need to be sufficiently far\n",
+ "away such that the additivity of the two surface potentials holds. In practice,\n",
+ "a separation of $d=10\\lambda_\\mathrm{D}$ is a good choice, represented by \n",
+ "`DISTANCE = 10`.\n",
+ "\n",
+ "Our choice of $c=10\\,\\mathrm{mmol}$ is a compromise between a sufficiently low\n",
+ "concentration for the PB theory to hold and not too large distances $d$ such\n",
+ "that the equilibration/diffusion of the ions is sufficiently fast\n",
+ "(`CONCENTRATION = 1e-2`).\n",
+ "\n",
+ "Note that in order to obtain results that we can interpret easily, we explicitly\n",
+ "set a unit system using nanometers as length-scale above.\n",
+ "The corresponding ion size of about 0.3 nm is a\n",
+ "typical value for a simple salt; this, however, is in sharp contrast to the\n",
+ "mean-field assumption of point-like ions.\n",
+ "The latter are not easily studied within Molecular Dynamics simulations due to\n",
+ "the required small time steps and are better suited for Monte Carlo type\n",
+ "simulations.\n",
+ "We instead focus here on analyzing deviations from PB theory due to the finite\n",
+ "ion size."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "### Task\n",
+ "\n",
+ "* write a function \n",
+ "`get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS)`\n",
+ "that returns the lateral and normal box lengths `box_l_xy` and `box_l_z` (in\n",
+ "nanometers) for the given parameters.\n",
+ "\n",
+ "**Hint:** To account for the finite ion size and the wall interaction it is\n",
+ "useful to define the effective separation $d^\\prime = d-2\\sigma$, such that the\n",
+ "concentration is $\\rho = N/(A \\cdot d^\\prime)$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):\n",
+ " \"\"\"\n",
+ " For a given number of particles, determine the lateral area of the box\n",
+ " to match the desired concentration.\n",
+ " \"\"\"\n",
+ "\n",
+ " # concentration is in mol/l, convert to 1/sigma**3\n",
+ " rho = concentration * (constants.Avogadro / constants.liter) * constants.nano**3\n",
+ " debye_length = (4. * np.pi * BJERRUM_LENGTH * rho * 2.)**(-1. / 2.) # desired Debye length in nm\n",
+ " l_z = distance * debye_length\n",
+ " \n",
+ " box_volume = n_ionpairs / rho\n",
+ " area = box_volume / (l_z - 2. * LJ_SIGMA) # account for finite ion size in density calculation\n",
+ " l_xy = np.sqrt(area)\n",
+ "\n",
+ " return l_xy, l_z\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "box_l_xy, box_l_z = get_box_dimension(CONCENTRATION, DISTANCE, N_IONPAIRS)\n",
+ "\n",
+ "# useful quantities for the following calculations\n",
+ "DEBYE_LENGTH = box_l_z / DISTANCE # in units of nm\n",
+ "rho = N_IONPAIRS / (box_l_xy * box_l_xy * box_l_z) # in units of 1/nm^3"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now can create the **ESPResSo** system.\n",
+ "\n",
+ "Note that for ELC to work properly, we need to add a gap of `ELC_GAP` in the\n",
+ "non-periodic direction.\n",
+ "The precise value highly affects the performance due to the tuning of the P3M\n",
+ "electrostatic solver.\n",
+ "For $d=10\\lambda$ a gap value of $6 L_z$ is a good choice.\n",
+ "\n",
+ "We also set the time-step $dt = 0.01 \\tau$, which is limited by the choice of\n",
+ "$\\sigma$ and $\\tau$ in the repulsive WCA interaction."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ELC_GAP = 6. * box_l_z\n",
+ "system = espressomd.System(box_l=[box_l_xy, box_l_xy, box_l_z + ELC_GAP])\n",
+ "system.time_step = 0.01"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 1.2 Set up the double-layer capacitor\n",
+ "\n",
+ "We now set up an electrolyte solution made of monovalent cations and anions\n",
+ "between two metallic electrodes at constant potential. \n",
+ "\n",
+ "#### 1.2.1 Electrode walls "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "### Task\n",
+ "* add two wall constraints at $z=0$ and $z=L_z$ to stop particles from\n",
+ "crossing the boundaries and model the electrodes.\n",
+ "Refer to \n",
+ "[espressomd.constraints.ShapeBasedConstraint](https://espressomd.github.io/doc/espressomd.html#espressomd.constraints.ShapeBasedConstraint)\n",
+ "and its\n",
+ "[wall constraint](https://espressomd.github.io/doc/constraints.html?highlight=constraint#wall)\n",
+ "in the documentation to set up constraints and the `types` dictionary for the\n",
+ "particle type."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "# Bottom wall, normal pointing in the +z direction \n",
+ "floor = espressomd.shapes.Wall(normal=[0, 0, 1])\n",
+ "c1 = system.constraints.add(\n",
+ " particle_type=types[\"Electrodes\"], penetrable=False, shape=floor)\n",
+ "\n",
+ "# Top wall, normal pointing in the -z direction\n",
+ "ceiling = espressomd.shapes.Wall(normal=[0, 0, -1], dist=-box_l_z) \n",
+ "c2 = system.constraints.add(\n",
+ " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceiling)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "#### 1.2.2 Add particles for the ions\n",
+ "\n",
+ "### Task\n",
+ "\n",
+ "* place ion pairs at random positions between the electrodes.\n",
+ "\n",
+ "Note, that unfavorable overlap can be avoided by placing the particles in the\n",
+ "interval $[\\sigma, d-\\sigma]$ in the $z$-direction only."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "offset = LJ_SIGMA # avoid unfavorable overlap at close distance to the walls\n",
+ "init_part_btw_z1 = offset \n",
+ "init_part_btw_z2 = box_l_z - offset\n",
+ "ion_pos = np.empty((3), dtype=float)\n",
+ "\n",
+ "for i in range (N_IONPAIRS):\n",
+ " ion_pos[0] = rng.random(1) * system.box_l[0]\n",
+ " ion_pos[1] = rng.random(1) * system.box_l[1]\n",
+ " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n",
+ " system.part.add(pos=ion_pos, type=types[\"Cation\"], q=charges[\"Cation\"])\n",
+ " \n",
+ "for i in range (N_IONPAIRS):\n",
+ " ion_pos[0] = rng.random(1) * system.box_l[0]\n",
+ " ion_pos[1] = rng.random(1) * system.box_l[1]\n",
+ " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n",
+ " system.part.add(pos=ion_pos, type=types[\"Anion\"], q=charges[\"Anion\"])\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "#### 1.2.3 Add interactions:\n",
+ "\n",
+ "### Task\n",
+ "\n",
+ "* For excluded volume interactions, add a WCA potential. \n",
+ "\n",
+ "Refer to the documentation to set up the\n",
+ "[WCA interaction](https://espressomd.github.io/doc/espressomd.html#espressomd.interactions.WCAInteraction) \n",
+ "under [Non-bonded](https://espressomd.github.io/doc/inter_non-bonded.html)\n",
+ "section."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "for t1 in types.values():\n",
+ " for t2 in types.values():\n",
+ " system.non_bonded_inter[t1, t2].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "For the (2D+h) electrostatic with dielectrics we choose the ELC-IC with P3M.\n",
+ "\n",
+ "Refer the documentation to set up\n",
+ "[ELCIC with P3M](https://espressomd.github.io/doc/electrostatics.html#electrostatic-layer-correction-elc)\n",
+ "under the [electrostatics](https://espressomd.github.io/doc/electrostatics.html)\n",
+ "section. \n",
+ "\n",
+ "As later we will study different potential drops between the electrodes, write a\n",
+ "function that sets up the electrostatic solver for a given value\n",
+ "`POTENTIAL_DIFF.`\n",
+ "This function will take care of tuning the P3M and ELC parameters.\n",
+ "For our purposes, an accuracy of $10^{-3}$ is sufficient.\n",
+ "\n",
+ "### Task\n",
+ "\n",
+ "* Write a function `setup_electrostatic_solver(potential_diff)` that\n",
+ "returns the ELC instance."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "def setup_electrostatic_solver(potential_diff):\n",
+ " delta_mid_top = -1. # (Fully metallic case both -1) \n",
+ " delta_mid_bot = -1.\n",
+ " accuracy = 1e-3\n",
+ " elc_accuracy = 1e-3\n",
+ " p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,\n",
+ " accuracy=accuracy,\n",
+ " tune=True,\n",
+ " verbose=False)\n",
+ " elc = espressomd.electrostatics.ELC(actor=p3m,\n",
+ " gap_size=ELC_GAP,\n",
+ " const_pot=True,\n",
+ " pot_diff=potential_diff,\n",
+ " maxPWerror=elc_accuracy,\n",
+ " delta_mid_bot=delta_mid_bot,\n",
+ " delta_mid_top=delta_mid_top)\n",
+ " return elc\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now add the solver to the system:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2. Equilibration\n",
+ "\n",
+ "### 2.1 Steepest descent\n",
+ "\n",
+ "Before we can start the simulation, we need to remove the overlap between particles.\n",
+ "For this, we use the steepest descent integrator.\n",
+ "Afterwards, we switch to a Velocity Verlet integrator and set up a Langevin thermostat.\n",
+ "Note, that we only analyze static properties, thus the damping and temperature chosen\n",
+ "here only determine the simulation time towards the equilibrium distribution."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "system.integrator.set_steepest_descent(f_max=10, gamma=50.0,\n",
+ " max_displacement=0.02)\n",
+ "system.integrator.run(250)\n",
+ "system.integrator.set_vv()\n",
+ "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Equilibrate the ion distribution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Equlibration parameters\n",
+ "STEPS_PER_SAMPLE = 200\n",
+ "N_SAMPLES_EQUIL = 50\n",
+ "N_PART = 2 * N_IONPAIRS\n",
+ "\n",
+ "times = np.zeros(N_SAMPLES_EQUIL)\n",
+ "e_total = np.zeros_like(times)\n",
+ "e_kin = np.zeros_like(times)\n",
+ "\n",
+ "for i in tqdm.trange(N_SAMPLES_EQUIL):\n",
+ " times[i] = system.time\n",
+ " energy = system.analysis.energy()\n",
+ " e_total[i] = energy['total']\n",
+ " e_kin[i] = energy['kinetic']\n",
+ " system.integrator.run(STEPS_PER_SAMPLE)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot the convergence of the total energy\n",
+ "plt.figure(figsize=(10, 6))\n",
+ "plt.plot(times, e_total, label='total')\n",
+ "plt.plot(times, e_kin, label='kinetic')\n",
+ "plt.xlabel('Simulation time')\n",
+ "plt.ylabel('Energy')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Convergence after $t\\sim50$ time units."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
+ "source": [
+ "## 3. Calculate and analyze ion profile\n",
+ "\n",
+ "### 3.1 Set up the density accumulators\n",
+ "\n",
+ "We now need to set up an \n",
+ "[espressomd.observables.DensityProfile](https://espressomd.github.io/doc/espressomd.html#espressomd.observables.DensityProfile)\n",
+ "observable to calculate the anion and cation density profiles.\n",
+ "\n",
+ "The time average is obtained through a\n",
+ "[espressomd.accumulators.MeanVarianceCalculator](espressomd.accumulators.MeanVarianceCalculator).\n",
+ "\n",
+ "### Task\n",
+ "\n",
+ "* Write a function `setup_densityprofile_accumulators(bin_width)` that returns the\n",
+ "`bin_centers` and the accumulators for both ion species in the $z$-range $[0,d]$.\n",
+ "Since we are not estimating errors in this tutorial, the choice of `delta_N` is\n",
+ "rather arbitrary and does not affect the results. In practice, a typical value is\n",
+ "`delta_N=20`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "def setup_densityprofile_accumulators(bin_width):\n",
+ " cations = system.part.select(type=types[\"Cation\"]) \n",
+ " anions = system.part.select(type=types[\"Anion\"])\n",
+ " n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))\n",
+ " density_profile_cation = espressomd.observables.DensityProfile(\n",
+ " 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,\n",
+ " max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)\n",
+ " density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(\n",
+ " obs=density_profile_cation, delta_N=20)\n",
+ " density_profile_anion = espressomd.observables.DensityProfile(\n",
+ " 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,\n",
+ " max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)\n",
+ " density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(\n",
+ " obs=density_profile_anion, delta_N=20)\n",
+ " zs = density_profile_anion.bin_centers()[0, 0, :, 2]\n",
+ " return zs, density_accumulator_cation, density_accumulator_anion\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(\n",
+ " bin_width=DEBYE_LENGTH / 10.)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 3.2 Run the simulation\n",
+ "\n",
+ "Now we take some measurement sampling the density profiles."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N_SAMPLES_PROD = 10\n",
+ "\n",
+ "# Add the accumulators\n",
+ "system.auto_update_accumulators.clear()\n",
+ "system.auto_update_accumulators.add(density_accumulator_cation)\n",
+ "system.auto_update_accumulators.add(density_accumulator_anion)\n",
+ " \n",
+ "times = []\n",
+ "e_total = []\n",
+ "for tm in tqdm.trange(N_SAMPLES_PROD):\n",
+ " system.integrator.run(STEPS_PER_SAMPLE)\n",
+ " times.append(system.time)\n",
+ " energy = system.analysis.energy()\n",
+ " e_total.append(energy['total'])\n",
+ "\n",
+ "cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n",
+ "anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Compare to analytical prediction\n",
+ "\n",
+ "Since we assume additivity, the total ion density follows from\n",
+ "$$ \\rho (z) = \\rho_+(z) - \\rho_+ (d-z) + \\rho_-(z) - \\rho_-(d-z) .$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def gouy_chapman_potential(x, debye_length, phi_0):\n",
+ " kappa = 1. / debye_length\n",
+ " return 2. * np.log((1. + np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))) \\\n",
+ " / (1. - np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))))\n",
+ "\n",
+ "def gouy_chapman_density(x, c0, debye_length, phi_0):\n",
+ " phi = gouy_chapman_potential(x, debye_length, phi_0)\n",
+ " return c0 / 2. * np.exp(-phi)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, (ax1, ax2, ax3) = plt.subplots(figsize=(16, 4), nrows=1, ncols=3, sharey=True)\n",
+ "fig.subplots_adjust(wspace=0)\n",
+ "\n",
+ "x = np.linspace(LJ_SIGMA, box_l_z-LJ_SIGMA, 100)\n",
+ "anion_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n",
+ " + gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2.\n",
+ "cation_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) \\\n",
+ " + gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.))/2.\n",
+ "\n",
+ "ax1.set_title('Cation')\n",
+ "ax2.set_title('Anion')\n",
+ "ax3.set_title('Total')\n",
+ "\n",
+ "ax1.plot(x, cation_profile_analytic, label='analytic')\n",
+ "ax2.plot(x, anion_profile_analytic, label='analytic')\n",
+ "ax3.plot(x, cation_profile_analytic + anion_profile_analytic, label='analytic')\n",
+ "\n",
+ "ax1.plot(zs[1:-1], cation_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n",
+ "ax2.plot(zs[1:-1], anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n",
+ "ax3.plot(zs[1:-1], cation_profile_mean[1:-1] + anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n",
+ "\n",
+ "ax1.legend(loc='upper center')\n",
+ "ax2.legend(loc='upper center')\n",
+ "ax3.legend(loc='upper center')\n",
+ "\n",
+ "ax2.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n",
+ "ax1.set_ylabel(r'$\\rho(z)\\,\\mathrm{[mol/l]}$')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We see good agreement between our simulation and the meanfield solution of Guy and Chapman. Low density and reasonably low potential make the assumptions of the analytical approach justified."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now check how well the surface charge agrees with Grahame's equation.\n",
+ "To this end we calculate \n",
+ "$$\\sigma = \\int_0^{d/2} \\rho(z) \\,\\mathrm{d}z .$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:len(zs)//2]) * (zs[1] - zs[0])\n",
+ "sigma_right = np.sum((cation_profile_mean-anion_profile_mean)[len(zs)//2:]) * (zs[1] - zs[0])\n",
+ "\n",
+ "def grahame_sigma(phi):\n",
+ " return np.sinh(phi / 4.) * np.sqrt(2. * rho / (np.pi * BJERRUM_LENGTH))\n",
+ "\n",
+ "sigma_grahame = grahame_sigma(POTENTIAL_DIFF)\n",
+ "print(f'simulation: {sigma_right:.3f} e/nm^2')\n",
+ "print(f'grahame: {sigma_grahame:.3f} e/nm^2')\n",
+ "print(f'relative deviation: {abs(1. - sigma_right/sigma_grahame) * 100.:.0f}%')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The electric field is readily obtained from the integral \n",
+ "$$E(z) = \\int_0^{z} \\frac{1}{\\varepsilon_0 \\varepsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot the electric field\n",
+ "fig, ax = plt.subplots(figsize=(10, 6))\n",
+ "\n",
+ "dz_SI = (zs[1] - zs[0]) * constants.nano\n",
+ "chargedensity = (cation_profile_mean - anion_profile_mean) * constants.elementary_charge / constants.nano**3 \n",
+ "E_SI = 1. / (EPSILON_R * constants.epsilon_0) * np.cumsum(chargedensity * dz_SI)\n",
+ "# integration constant: zero field in the center\n",
+ "E_SI -= E_SI.min()\n",
+ "E = E_SI / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE) / constants.nano)\n",
+ "ax2 = plt.twinx()\n",
+ "\n",
+ "ax.plot(zs, E_SI)\n",
+ "ax2.plot(zs, E)\n",
+ "ax.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n",
+ "ax.set_ylabel(r'$E_\\mathrm{ind}\\,\\mathrm{[V/m]}$')\n",
+ "ax2.set_ylabel(r'$E_\\mathrm{ind}\\,\\mathrm{[(k_\\mathrm{B}T/e)/nm]}$')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We see that the electric field reduces to 0 in the middle of the channel, justifying the assumption that the two electrodes are far enough apart to not influence each other."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The electric potential can be calculated from $\\phi(z) = \\int_0^z -E(z^\\prime)\\,\\mathrm{d}z^\\prime$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot the elecrostatic potential\n",
+ "fig, ax = plt.subplots(figsize=(10, 6))\n",
+ "ax2 = ax.twinx()\n",
+ "phi_SI = -np.cumsum(E_SI * dz_SI)\n",
+ "phi = phi_SI * (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))\n",
+ "ax.plot(zs, phi_SI)\n",
+ "ax.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n",
+ "ax.set_ylabel(r'$\\phi\\,[V]$')\n",
+ "ax2.set_ylabel(r'$\\phi\\,[k_\\mathrm{B}T/e]$')\n",
+ "ax2.axhline(-5, ls='--', color='k')\n",
+ "ax2.axhline(0, ls='--', color='k')\n",
+ "ax.set_xlim(0, 10. * DEBYE_LENGTH)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "measured_potential_difference = -(phi[-1] - phi[0])\n",
+ "print(f'applied voltage: {POTENTIAL_DIFF:.2f} V')\n",
+ "print(f'measured voltage: {measured_potential_difference:.2f} V')\n",
+ "print(f'relative deviation: {abs(1. - measured_potential_difference / POTENTIAL_DIFF) * 100:.0f}%')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 4. Differential capacitance\n",
+ "\n",
+ "With the above knowledge, we can now assess the \n",
+ "differential capacitance of the system, by changing the applied voltage\n",
+ "difference and determining the corresponding surface charge density."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sigma_vs_phi = []\n",
+ "MIN_PHI = 0.5\n",
+ "MAX_PHI = 10\n",
+ "N_PHI = 7\n",
+ "N_SAMPLES_EQUIL_CAP = 5\n",
+ "N_SAMPLES_CAP = 5\n",
+ "\n",
+ "# sample from high to low potential to improve sampling\n",
+ "for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n",
+ " system.electrostatics.solver = setup_electrostatic_solver(potential_diff)\n",
+ " system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)\n",
+ " sigmas = []\n",
+ "\n",
+ " for tm in range(N_SAMPLES_CAP):\n",
+ " zs, density_accumulator_cation, density_accumulator_anion = \\\n",
+ " setup_densityprofile_accumulators(bin_width=DEBYE_LENGTH / 10.)\n",
+ "\n",
+ " system.auto_update_accumulators.clear()\n",
+ " system.auto_update_accumulators.add(density_accumulator_cation)\n",
+ " system.auto_update_accumulators.add(density_accumulator_anion)\n",
+ " system.integrator.run(STEPS_PER_SAMPLE)\n",
+ "\n",
+ " cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n",
+ " anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]\n",
+ "\n",
+ " sigmas.append(np.sum((cation_profile_mean - anion_profile_mean)[:int(len(zs) / 2.)]) * (zs[1] - zs[0]))\n",
+ "\n",
+ " sigma_vs_phi.append([potential_diff, np.mean(sigmas), scipy.stats.sem(sigmas)]) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fig, ax = plt.subplots(figsize=(10, 6))\n",
+ "sigma_vs_phi = np.array(sigma_vs_phi)\n",
+ "x = np.linspace(0, sigma_vs_phi[:,0].max())\n",
+ "phi_SI = sigma_vs_phi[:,0] / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))\n",
+ "plt.errorbar(-sigma_vs_phi[:,1] * constants.elementary_charge / constants.nano**2,\n",
+ " phi_SI, xerr=sigma_vs_phi[:,2] * scipy.constants.elementary_charge / scipy.constants.nano**2,\n",
+ " fmt='o',label='Simulation')\n",
+ "plt.plot(grahame_sigma(x) * constants.elementary_charge / constants.nano**2,\n",
+ " x / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE)), label='Grahame')\n",
+ "x = np.linspace(0, ax.get_ylim()[1])\n",
+ "plt.plot(EPSILON_R * constants.epsilon_0 * x / 2. / (DEBYE_LENGTH * constants.nano), x, label='linear PB',\n",
+ " ls='--')\n",
+ "plt.xlabel(r'$\\sigma\\,\\mathrm{[C/m^2]}$')\n",
+ "plt.ylabel(r'$\\phi_\\mathrm{s}\\,\\mathrm{[V]}$')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For small potential drops, one observes the expected Poisson–Boltzmann behavior. It also agrees with the linearized solution $\\sigma(\\phi_\\mathrm{s}) = \\varepsilon_r\\varepsilon_0 \\frac{\\phi_\\mathrm{s}}{2 \\lambda_\\mathrm{D}}$.\n",
+ "However, we observe already for potentials $\\sim 0.1\\,\\mathrm{V} = 4\\,k_\\mathrm{B}T / e$ a significant deviation which can be attributed to the fact that our ions are of finite size and thus the surface charge saturates."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## References\n",
+ "\n",
+ "[1] Conway, B. E. Electrochemical Supercapacitors; Springer US: Boston, MA, 1999. https://doi.org/10.1007/978-1-4757-3058-6.\n",
+ "\n",
+ "[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.\n",
+ "\n",
+ "[3] Gouy, G. Constitution of the Electric Charge at the Surface of an Electrolyte. J. Phys. 1910, 9 (4), 457–467.\n",
+ "\n",
+ "[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.\n",
+ "\n",
+ "[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.\n",
+ "\n",
+ "[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."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt
index e93755f95a..b79ce7c963 100644
--- a/testsuite/scripts/tutorials/CMakeLists.txt
+++ b/testsuite/scripts/tutorials/CMakeLists.txt
@@ -61,6 +61,8 @@ tutorial_test(FILE test_ferrofluid_1.py)
tutorial_test(FILE test_ferrofluid_2.py)
tutorial_test(FILE test_ferrofluid_3.py)
tutorial_test(FILE test_constant_pH__ideal.py)
+tutorial_test(FILE test_electrodes_1.py)
+tutorial_test(FILE test_electrodes_2.py)
tutorial_test(FILE test_constant_pH__interactions.py)
tutorial_test(FILE test_widom_insertion.py)
diff --git a/testsuite/scripts/tutorials/test_electrodes_1.py b/testsuite/scripts/tutorials/test_electrodes_1.py
new file mode 100644
index 0000000000..eaa490e7a9
--- /dev/null
+++ b/testsuite/scripts/tutorials/test_electrodes_1.py
@@ -0,0 +1,44 @@
+#
+# Copyright (C) 2023 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+import unittest as ut
+import importlib_wrapper
+import numpy as np
+
+tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
+ "@TUTORIALS_DIR@/electrodes/electrodes_part1.py", N_AXIAL_POINTS=6,
+ script_suffix="@TEST_SUFFIX@")
+
+
+@skipIfMissingFeatures
+class Tutorial(ut.TestCase):
+
+ def test_force(self):
+ ref_force = tutorial.analytic_force_centered(
+ tutorial.r, tutorial.box_l_z)
+ msg = "The force for particle 1 doesn't agree with analytical expression"
+ np.testing.assert_allclose(np.log(tutorial.elc_forces_axial[:, 0]),
+ np.log(ref_force), rtol=0.05, err_msg=msg)
+ msg = "The force for particle 2 doesn't agree with analytical expression"
+ np.testing.assert_allclose(np.log(-tutorial.elc_forces_axial[:, 1]),
+ np.log(ref_force), rtol=0.05, err_msg=msg)
+
+
+if __name__ == "__main__":
+ ut.main()
diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py
new file mode 100644
index 0000000000..5654db1ed8
--- /dev/null
+++ b/testsuite/scripts/tutorials/test_electrodes_2.py
@@ -0,0 +1,72 @@
+#
+# Copyright (C) 2023 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+
+import unittest as ut
+import importlib_wrapper
+import numpy as np
+from scipy import constants
+
+params = {'N_SAMPLES_EQUIL': 25, 'N_SAMPLES_PROD': 5,
+ 'N_SAMPLES_EQUIL_CAP': 5, 'N_SAMPLES_CAP': 1,
+ 'MIN_PHI': 1, 'MAX_PHI': 2.5, 'N_PHI': 4}
+
+tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
+ "@TUTORIALS_DIR@/electrodes/electrodes_part2.py",
+ script_suffix="@TEST_SUFFIX@", **params)
+
+
+@skipIfMissingFeatures
+class Tutorial(ut.TestCase):
+
+ def test_potential_difference(self):
+ # Test that the applied potential difference equals the one from
+ # integrating Poisson equation
+ msg = 'The potential difference is not equal to the one from integrating Poisson equation.'
+ self.assertAlmostEqual(
+ tutorial.measured_potential_difference / tutorial.POTENTIAL_DIFF,
+ 1., delta=0.25, msg=msg)
+
+ def test_charge_profile(self):
+ # Roughly test the profile, deviations are expected!!
+ charge_profile = (
+ tutorial.cation_profile_mean +
+ tutorial.anion_profile_mean)
+ analytic = (
+ tutorial.cation_profile_analytic +
+ tutorial.anion_profile_analytic)
+ msg = 'The density profile is not sufficiently equal to PB theory.'
+ np.testing.assert_allclose(
+ charge_profile[1:-1],
+ analytic[1:-1],
+ rtol=5e-2,
+ atol=5e-2,
+ err_msg=msg)
+
+ def test_capacitance(self):
+ # For low potentials the capacitance should be in line with Grahame/DH
+ # equilibration performance limiting
+ grahame = -tutorial.sigma_vs_phi[:, 0] / (
+ constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE))
+ msg = 'The capacitance at low potentials should be in line with Grahame/DH.'
+ np.testing.assert_allclose(
+ grahame, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg)
+
+
+if __name__ == "__main__":
+ ut.main()