From 6cc0fdcb0da4a0aaa6edaaba430019b0b8898873 Mon Sep 17 00:00:00 2001 From: Alexander Schlaich Date: Mon, 2 Oct 2023 21:48:51 +0200 Subject: [PATCH 1/3] Add new tutorial on ICC/ELC/ELC-IC Co-authored-by: keerthirk --- doc/tutorials/CMakeLists.txt | 1 + doc/tutorials/electrodes/CMakeLists.txt | 26 + doc/tutorials/electrodes/NotesForTutor.md | 31 + .../electrodes/electrodes_part1.ipynb | 500 +++++++ .../electrodes/electrodes_part2.ipynb | 1240 +++++++++++++++++ testsuite/scripts/tutorials/CMakeLists.txt | 2 + .../scripts/tutorials/test_electrodes_1.py | 40 + .../scripts/tutorials/test_electrodes_2.py | 68 + 8 files changed, 1908 insertions(+) create mode 100644 doc/tutorials/electrodes/CMakeLists.txt create mode 100644 doc/tutorials/electrodes/NotesForTutor.md create mode 100644 doc/tutorials/electrodes/electrodes_part1.ipynb create mode 100644 doc/tutorials/electrodes/electrodes_part2.ipynb create mode 100644 testsuite/scripts/tutorials/test_electrodes_1.py create mode 100644 testsuite/scripts/tutorials/test_electrodes_2.py diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt index 508bda4caf7..479f3491e53 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 00000000000..4ba4352845b --- /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 00000000000..3c87bb760b8 --- /dev/null +++ b/doc/tutorials/electrodes/NotesForTutor.md @@ -0,0 +1,31 @@ +# 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 00000000000..06146bdf413 --- /dev/null +++ b/doc/tutorials/electrodes/electrodes_part1.ipynb @@ -0,0 +1,500 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Basic simulation of electrodes in ESPResSo part I:\n", + "# 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 respective 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 yields a discontinuity, which can be expressed in terms of a finite\n", + "surface charge density \n", + "$(\\epsilon_1\\vec{E}_1 - \\epsilon_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 $\\epsilon_1$ to another one,\n", + "$\\epsilon_2$.\n", + "\n", + "While in the case of non-polarizable materials ($\\epsilon_1 = \\epsilon_2 = 1$),\n", + "this jump is only related to surface charges and dielectric contrast and the\n", + "potential is continuous across the interface, for polarizable materials also the\n", + "polarization field $\\vec{P}$ will give a contribution. \n", + "In order to solve the problem 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_{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\n", + "*scheme for solving the boundary integrals and the induced charge. \n", + "\n", + "*Note*: Apart from ICC $\\star$, **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.\n", + "The tutorial on *Basic simulation of electrodes in ESPResSo part II*\n", + "addresses this in detail.\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 the dielectric contrast is\n", + "$$ \\Delta = \\frac{\\epsilon_1 - \\epsilon_2} {\\epsilon_1 + \\epsilon_2} = -1 .$$\n", + "If the ions are placed in the center of the 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 \\epsilon_0 \\epsilon_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 \\epsilon_0 \\epsilon_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 \\epsilon_0 \\epsilon_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 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 $3-D$\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 (\\epsilon \\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_{ind} = \\frac{1}{4\\pi} \\oint\\, dA\\, \\cdot \\epsilon\\nabla \\phi = \\frac{A}{4\\pi}(\\epsilon_1\\vec{E}_1 \\cdot \\hat{n}-\\epsilon_2\\vec{E}_2 \\cdot\\hat{n})$$\n", + "- The electric field at the closest proximity of the interface, $\\vec{E}_{1/2}$,\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/2} =\\vec{E} + 2\\pi/\\epsilon_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{\\epsilon_1}{2\\pi} \\frac{\\epsilon_1-\\epsilon_2}{\\epsilon_1+\\epsilon_2}\\vec{E} \\cdot \\hat{n} =: \\frac{\\epsilon_1}{2\\pi} \\Delta \\, \\vec{E} \\cdot \\hat{n}$$\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/2}$ 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 espressomd\n", + "import numpy as np\n", + "import espressomd.electrostatics\n", + "import espressomd.electrostatic_extensions\n", + "from espressomd.interactions import *\n", + "\n", + "espressomd.assert_features(['ELECTROSTATICS'])\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from scipy.special import *\n", + "from tqdm import tqdm" + ] + }, + { + "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), time ($t=0.01 \\tau$) 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$ and thus defines the energy scale in our\n", + "system.\n", + "It is defined as\n", + "$$\\ell_\\mathrm{B}=\\frac{1}{4\\pi\\epsilon_0\\epsilon_r k_\\mathrm{B}T}.$$ \n", + "In our case, if we choose the ion size ($\\sigma$) in simulations to represent a\n", + "typical value for mono-atomar 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": [ + "#***************************************************\n", + "# System Setup\n", + "#***************************************************\n", + "\n", + "# 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", + "# ICC* with ELC: 2D electrostatics\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\n", + "system.time_step = 0.01\n", + "system.cell_system.skin = 0.4\n", + "\n", + "# Elementary charge \n", + "q = np.array([1.0]) \n", + "\n", + "# Interaction Parameters for P3M with ELC\n", + "\n", + "BJERRUM_LENGTH = 2.0 # Electrostatic prefactor passed to P3M ; prefactor=lB KBT/e2 \n", + "ACCURACY = 1e-7 # P3M force accuracy \n", + "CHECK_ACCURACY = 1e-7 # maximim pairwise error in ELC\n", + "\n", + "#Lennard-Jones Parameters\n", + "\n", + "LJ_SIGMA = 1.0\n", + "LJ_EPSILON = 1.0 \n", + "\n", + "#Particle parameters\n", + "\n", + "types = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", + "charges = {\"Cation\": q[0], \"Anion\": -q[0] }\n", + "\n", + "p1=system.part.add(pos=[box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges[\"Cation\"])\n", + "print(f\"Cation placed at position: {p1.pos}\")\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\"])\n", + "print(f\"Anion placed at position: {p2.pos}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task:\n", + "Set up electrostatics interactions between the ion pairs using P3M and ICC $\\star$. \n", + "#### Hints:\n", + "- First instantiate a P3M object (using\n", + " [*espressomd.electrostatics.P3M(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.P3M))\n", + " as the 3D Coulomb solver for the system.\n", + "- Since, we are working on a 2D+h partially periodic system, we use ELC in\n", + " combination with P3M. For this purpose, instantiate ELC (using\n", + " [*espressomd.electrostatics.ELC(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.ELC))\n", + " with p3m as actor, which calls necessary initialization routines. \n", + "- Now we account for dielectric interface by adding ICC $^\\star$.\n", + " To this end, first create an ICC object using\n", + " [*espressomd.electrostatic_extensions.ICC(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatic_extensions.ICC).\n", + "- Setting up ICC $^\\star$ requires specifying a fixed number of ICC particles\n", + " and their initial charges.\n", + " Since the ICC particles are normal ESPResSo particles, they need to be added\n", + " at the interfaces using\n", + " [*system.part.add(..)*](https://espressomd.github.io/doc/espressomd.html?#module-espressomd.particle_data).\n", + "\n", + "*Note* - Refer to section\n", + "[**Dielectric interfaces with the ICC algorithm**](https://espressomd.github.io/doc/electrostatics.html#dielectric-interfaces-with-the-icc-star-algorithm)\n", + "in the **ESPResSo** documentation for the basics of an ICC $^\\star$ setup.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "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", + " )\n", + "\n", + "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=CHECK_ACCURACY)\n", + "\n", + "# Set the ICC line density and calculate the number of\n", + "# ICC particles according to the box size\n", + "nicc = 100 # linear density\n", + "nicc_per_electrode = nicc**2 # surface density\n", + "nicc_tot = 2 * nicc_per_electrode\n", + "iccArea = box_l_x * box_l_y / nicc_per_electrode\n", + "lx = box_l_x / nicc\n", + "ly = box_l_y / nicc\n", + "\n", + "# Lists to collect required parameters\n", + "iccNormals = []\n", + "iccAreas = []\n", + "iccSigmas = []\n", + "iccEpsilons = []\n", + "\n", + "# Add the fixed ICC particles:\n", + "\n", + "# Bottom electrode (normal [0, 0, 1])\n", + "for xi in range(nicc):\n", + " for yi in range(nicc):\n", + " system.part.add(pos=[lx * xi, ly * yi, 0.], q=-0.0001,\n", + " type=types[\"Electrodes\"], fix=[True, True, True])\n", + "iccNormals.extend([0, 0, 1] * nicc_per_electrode)\n", + "\n", + "# Top electrode (normal [0, 0, -1])\n", + "for xi in range(nicc):\n", + " for yi in range(nicc):\n", + " system.part.add(pos=[lx * xi, ly * yi, box_l_z], q=0.0001,\n", + " type=types[\"Electrodes\"], fix=[True, True, True])\n", + "iccNormals.extend([0, 0, -1] * nicc_per_electrode)\n", + "\n", + "# Common area, sigma and metallic epsilon\n", + "iccAreas.extend([iccArea] * nicc_tot)\n", + "iccSigmas.extend([0] * nicc_tot)\n", + "iccEpsilons.extend([100000] * nicc_tot)\n", + "\n", + "iccNormals = np.array(iccNormals, dtype=float).reshape(nicc_tot, 3)\n", + "\n", + "icc = espressomd.electrostatic_extensions.ICC(\n", + " first_id=2,\n", + " n_icc=nicc_tot,\n", + " convergence=1e-3,\n", + " relaxation=0.95,\n", + " ext_field=[0, 0, 0],\n", + " max_iterations=1000,\n", + " eps_out=1.0,\n", + " normals=iccNormals,\n", + " areas=np.array(iccAreas, dtype=float),\n", + " sigmas=np.array(iccSigmas, dtype=float),\n", + " epsilons=np.array(iccEpsilons, dtype=float)\n", + " )\n", + "\n", + "system.electrostatics.solver = elc\n", + "system.electrostatics.extension = icc" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Calculation of the forces" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r = np.logspace(0,box_l_z/4.,10) \n", + "elc_forces_axial = np.empty((len(r), 2))\n", + "\n", + "for i, x in enumerate(tqdm(r)):\n", + " p1.pos = [0, box_l_y/2.0, box_l_z/2.0]\n", + " p2.pos = [x, box_l_y/2.0, box_l_z/2.0]\n", + "\n", + " system.integrator.run(0)\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 check passes\n", + " system.part.by_ids(range(2,2+nicc_per_electrode)).q = np.array([-0.0001]*nicc_per_electrode)\n", + " system.part.by_ids(range(2+nicc_per_electrode,2+2*nicc_per_electrode)).q = np.array([0.0001]*nicc_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", + " x = r/w\n", + " prefactor = BJERRUM_LENGTH / w**2\n", + "\n", + " def summand(x, n):\n", + " return (-1)**n * x/(x**2+n**2)**(3./2)\n", + " \n", + " def do_sum(x):\n", + " max = int(1e3)\n", + " sum = 0\n", + " for n in range(-max+1,max+1):\n", + " sum += summand(x,n)\n", + " return sum\n", + "\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", + "plt.plot(r/BJERRUM_LENGTH, -elc_forces_axial[:,1], color='red' , 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", + "x = np.logspace(-.25,1.45,100)\n", + "\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='green', ls='--', label='Coulomb')\n", + "plt.plot(x/BJERRUM_LENGTH, exponential_force(x,box_l_z), color='red', ls='--', label='Exponential')\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()" + ] + }, + { + "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 (ipykernel)", + "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.11.5" + } + }, + "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 00000000000..a22fadbbc40 --- /dev/null +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -0,0 +1,1240 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Basic simulation of electrodes in ESPResSo part II:\n", + "# Electrolyte capacitor and Poisson-Boltzmann theory" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "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 respective 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": { + "deletable": true, + "editable": true + }, + "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 of 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 **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 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 against 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 on the fly.\n", + "This is dealt in a computational cost-effective way using the ELC-IC method to\n", + "treat the image charge effect in the presence of 2D dielectric bounding\n", + "interfaces. " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "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", + "withing the mean-field approximation of the Poisson Boltzmann (PB) equation. \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", + "$$\\psi(z) = -2\\ln\\left[\n", + " \\frac{1-\\tanh(\\psi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D} z}}\n", + " {1+\\tanh(\\psi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D}z}} \\right].$$ \n", + "Here, $\\psi_\\mathrm{s}=\\psi(z=0)=$ const is the surface potential considering\n", + "that $\\psi(z\\rightarrow \\infty)=0$ and\n", + "$\\kappa_\\mathrm{D} = \\lambda_\\mathrm{D}^{-1}$ is the inverse Debye screening length given by\n", + "$$ \\lambda_\\mathrm{D} = \\left(\\frac{\\epsilon \\, 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", + "$\\psi_\\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(\\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 capactitance\n", + "$C=\\frac{\\mathrm{d}\\sigma}{\\mathrm{d}\\phi_\\mathrm{s}}$ within the mean-field\n", + "solution for non-overlapping double layers." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "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 two-dimensional a 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 accounts 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 $3-D$\n", + "Coulomb solver.\n", + "Furthermore, ELC can be extended straightforward to implement 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\n", + "dimensional Poisson's equation:\n", + "$$-\\epsilon_{0}\\epsilon_{r}\\phi_\\mathrm{ind}(z)=\\iint_{0}^{z}\\rho(z^{'})(z-z^{'})dz^{'} $$\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 straightforward as the corresponding contribution from the induced\n", + "charges is known:\n", + "$$ E_z^\\mathrm{(ind)} = \\frac{1}{\\epsilon_0 \\epsilon_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 $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 ab 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": { + "deletable": true, + "editable": true + }, + "source": [ + "## 1. System setup \n", + "\n", + "First we import all ESPResSo features and external modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "import espressomd\n", + "import numpy as np\n", + "rng = np.random.default_rng(42)\n", + "\n", + "import espressomd.electrostatics\n", + "import espressomd.electrostatic_extensions\n", + "from espressomd.interactions import *\n", + "import espressomd.observables\n", + "import espressomd.accumulators\n", + "from espressomd import shapes\n", + "\n", + "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from scipy import constants as const\n", + "from tqdm import tqdm\n", + "from scipy.stats import sem" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "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) and time $t=0.01 \\tau$.\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": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "# water at room temperature\n", + "EPSILON_R = 78.4 # Relative dielectric constant of the water\n", + "TEMPERATURE = 300.0 # Temperature in Kelvin\n", + "BJERRUM_LENGTH = const.elementary_charge**2 / (4*np.pi*const.epsilon_0*EPSILON_R*const.Boltzmann*TEMPERATURE) / const.nano\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 nanometers, not point-like\n", + "LJ_EPSILON = 1.0\n", + "HS_ION_SIZE = 2**(1/6) * LJ_SIGMA\n", + "\n", + "CONCENTRATION = 1e-2 # desired concentration 10 mmol/l\n", + "DISTANCE = 10 # 10 Debye lengths\n", + "N_IONPAIRS = 500\n", + "\n", + "POTENTIAL_DIFF = 5.0\n", + "\n", + "# Elementary charge \n", + "q = np.array([1.0]) \n", + "types = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", + "charges = {\"Cation\": q[0], \"Anion\": -q[0] }" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden", + "solution2_first": true + }, + "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 (`HS_ION_SIZE`) of about 0.33 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.\n", + "\n", + "The first task now is to 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*d^\\prime)$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):\n", + " \"\"\" For a given number of particles, determine the lateral area of the box\n", + " to match the desired concentration \"\"\"\n", + "\n", + " # concentration is in mol/l, convert to 1/sigma**3\n", + " rho = concentration * (const.Avogadro / const.liter) * const.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*HS_ION_SIZE) # account for finite ion size in density calculation\n", + " l_xy = np.sqrt(area)\n", + "\n", + " return l_xy, l_z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "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": { + "deletable": true, + "editable": true + }, + "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 value of `ELC_GAP = 6*box_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": { + "deletable": true, + "editable": true + }, + "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": { + "deletable": true, + "editable": true, + "solution2": "hidden", + "solution2_first": true + }, + "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 \n", + "\n", + "First, we add two wall constraints at $z=0$ and $z=L_z$ to stop particles from\n", + "crossing the boundaries and model the electrodes.\n", + "\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": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "# 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", + "ceil = espressomd.shapes.Wall(normal=[0, 0, -1],\n", + " dist=-box_l_z) \n", + "c2 = system.constraints.add(\n", + " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceil)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "#### 1.2.2 Add particles for the ions\n", + "\n", + "Now, place the 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": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "offset=HS_ION_SIZE # avoid unfavorable overlap at close to the walls\n", + "Init_part_btw_z1=0+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\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "#### 1.2.3 Add interactions:\n", + "For excluded volume interactions, we 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.\n", + "\n", + "Add the corresponding interaction parameters to the system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "for key, val in types.items():\n", + " for key1, val1 in types.items():\n", + " system.non_bonded_inter[val, val1].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true, + "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", + "Write a function `setup_electrostatic_solver(potential_diff)` that\n", + "returns the ELC instance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "def setup_electrostatic_solver(potential_diff):\n", + "\n", + " delta_mid_top = -1.0 #(Fully metallic case both -1) \n", + " delta_mid_bot = -1.0\n", + "\n", + " accuracy = 1e-3\n", + " check_accuracy = 1e-3\n", + " p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,\n", + " accuracy=accuracy, \n", + " tune=True,\n", + " )\n", + " \n", + " elc = espressomd.electrostatics.ELC(actor=p3m,\n", + " gap_size=ELC_GAP,\n", + " const_pot=True,\n", + " pot_diff=potential_diff,\n", + " maxPWerror=check_accuracy,\n", + " delta_mid_bot=delta_mid_bot,\n", + " delta_mid_top=delta_mid_top)\n", + " \n", + " return elc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "Now add the solver to the system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## 2. Equilibration\n", + "\n", + "### 2.1 Steepest descent\n", + "\n", + "Before we can start the simulation, we need to remove the overlap between\n", + "particles to avoid large forces which would crash the simulation.\n", + "\n", + "For this, we use the steepest descent integrator with a relative convergence\n", + "criterion for forces and energies.\n", + "\n", + "After steepest descent, we switch to a Velocity Verlet integrator and set up a\n", + "Langevin thermostat.\n", + "Note, that we only analyze static properties, thus the damping and temperature\n", + "chosen here only determine the relaxation speed towards the equilibrium\n", + "distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "# Relax the overlaps with steepest descent\n", + "\n", + "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() # Switch bach to Velocity Verlet \n", + "\n", + "# Add thermostat \n", + "thermostat_seed = np.random.randint(np.random.randint(1000000))\n", + "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=thermostat_seed)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## Equilibrate the ion distribution\n", + "\n", + "Convergence after $t\\sim25$ time units, possible to run up to $t=100$ here...\n", + "This is a total of 10.000 time steps (~1 minute)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "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(range(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": { + "deletable": true, + "editable": true + }, + "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('t')\n", + "plt.ylabel('E')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true, + "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", + "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": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true, + "solution2": "hidden" + }, + "outputs": [], + "source": [ + "def setup_densityprofile_accumulators(bin_width):\n", + "\n", + " Ion_id=[]\n", + " Cations = system.part.select(type=types[\"Cation\"])\n", + " Cations_id=[]\n", + " for i in Cations:\n", + " Cations_id.append(i.id)\n", + " Ion_id.append(i.id)\n", + " \n", + " Anions = system.part.select(type=types[\"Anion\"])\n", + " Anions_id=[]\n", + " for i in Anions:\n", + " Anions_id.append(i.id)\n", + " Ion_id.append(i.id)\n", + " \n", + " n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))\n", + " \n", + " # Accumulator 1 : observable::Density_Profile\n", + " density_profile_cation = espressomd.observables.DensityProfile(ids=Cations_id,\n", + " n_x_bins=1,\n", + " n_y_bins=1,\n", + " n_z_bins=n_z_bins,\n", + " min_x=0,\n", + " min_y=0,\n", + " min_z=0,\n", + " max_x=system.box_l[0],\n", + " max_y=system.box_l[1],\n", + " max_z=system.box_l[2] - ELC_GAP)\n", + " \n", + " density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(obs=density_profile_cation, delta_N=20)\n", + " \n", + " \n", + " density_profile_anion = espressomd.observables.DensityProfile(ids=Anions_id,\n", + " n_x_bins=1,\n", + " n_y_bins=1,\n", + " n_z_bins=n_z_bins,\n", + " min_x=0,\n", + " min_y=0,\n", + " min_z=0,\n", + " max_x=system.box_l[0],\n", + " max_y=system.box_l[1],\n", + " max_z=system.box_l[2] - ELC_GAP)\n", + " \n", + " density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(obs=density_profile_anion, delta_N=20)\n", + "\n", + " zs = density_profile_anion.bin_centers()[0, 0, :, 2]\n", + " \n", + " return zs, density_accumulator_cation, density_accumulator_anion" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(bin_width = DEBYE_LENGTH/10.)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "### 3.2 Run the simulation\n", + "\n", + "Now we take some measurement sampling the density profiles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "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(range(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": { + "deletable": true, + "editable": true + }, + "source": [ + "### Compare to analytical prediction\n", + "\n", + "Since we assume pair-wise 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": { + "deletable": true, + "editable": true + }, + "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": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 6))\n", + "plt.plot(zs, cation_profile_mean, color='b', label='cation')\n", + "plt.plot(zs, anion_profile_mean, color='r', label='anion')\n", + "plt.plot(zs, cation_profile_mean + anion_profile_mean, color='k', label='total')\n", + "\n", + "x = np.linspace(HS_ION_SIZE, box_l_z-HS_ION_SIZE, 100)\n", + "plt.plot(x, (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2., color='r', ls='--')\n", + "plt.plot(x, (gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='b', ls='--')\n", + "plt.plot(x, (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2 \\\n", + " + (gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='k', ls='--', lw=2)\n", + "\n", + "plt.legend()\n", + "plt.xlabel(r'$z\\,\\mathrm{[nm]}$')\n", + "plt.ylabel(r'$\\rho(z)\\,\\mathrm{[mol/l]}$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "We now check how well the surface charge agrees with Graham 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": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "#Test sigma with Graham equation\n", + "sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:int(len(zs)/2.)]) * (zs[1]-zs[0])\n", + "sigma_right = np.sum((+cation_profile_mean-anion_profile_mean)[int(len(zs)/2.):]) * (zs[1]-zs[0])\n", + "\n", + "def graham_sigma(phi):\n", + " return np.sinh(phi/4.) * np.sqrt(2*rho/(np.pi*BJERRUM_LENGTH))\n", + "sigma_graham = graham_sigma(POTENTIAL_DIFF)\n", + "\n", + "# sigma in e/nm^2\n", + "print('simulation:', sigma_right, 'graham:', sigma_graham, 'relative:', sigma_right/sigma_graham)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "The electric field is readily obtained from the integral \n", + "$$E(z) = \\int_0^{z} \\frac{1}{\\epsilon_0 \\epsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "# plot the electric field\n", + "f, ax = plt.subplots(figsize=(10, 6))\n", + "\n", + "dz_SI = (zs[1]-zs[0])*const.nano\n", + "chargedensity = (cation_profile_mean-anion_profile_mean)*const.elementary_charge/const.nano**3 \n", + "E_SI = 1/(EPSILON_R*const.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 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE) / const.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": { + "deletable": true, + "editable": true + }, + "source": [ + "And the electric potential from $\\phi(z) = \\int_0^z -E(z^\\prime)\\,\\mathrm{d}z^\\prime$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "# plot the elecrostatic potential\n", + "f, ax = plt.subplots(figsize=(10, 6))\n", + "ax2 = ax.twinx()\n", + "phi_SI = -np.cumsum(E_SI*dz_SI)\n", + "phi = phi_SI * (const.elementary_charge / (const.Boltzmann * TEMPERATURE))\n", + "ax.plot(zs, phi_SI)\n", + "ax2.plot(zs, phi)\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", + "ax.axhline(-5 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)))\n", + "ax2.axhline(0, ls='--', color='k')\n", + "ax.axhline(0 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)))\n", + "ax.set_xlim(0, 10*DEBYE_LENGTH)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "measured_potential_difference = -(phi[-1]+phi[0])\n", + "print('applied voltage', POTENTIAL_DIFF, 'measured voltage', measured_potential_difference,\n", + " 'relative deviation:', 1-measured_potential_difference/POTENTIAL_DIFF)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## 4. Differential capacitance\n", + "\n", + "With the above knowledge, we can now easily assess the \n", + "differential capacitance of the system, i.e. change the applied voltage\n", + "difference and determine the corresponding surface charge density." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "sigma_vs_phi = []\n", + "MIN_PHI = 0.5\n", + "MAX_PHI = 14\n", + "N_PHI = 10\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(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n", + "\n", + " system.electrostatics.solver = setup_electrostatic_solver(potential_diff)\n", + "\n", + " times=[]\n", + " e_total=[]\n", + " sigmas = []\n", + " for tm in range(N_SAMPLES_EQUIL_CAP):\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", + " for tm in range(N_SAMPLES_CAP):\n", + "\n", + " zs, density_accumulator_cation, density_accumulator_anion = 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", + "\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, :]\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), sem(sigmas)]) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "f, 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()*1.3)\n", + "phi_SI = sigma_vs_phi[:,0] / (const.elementary_charge / (const.Boltzmann * TEMPERATURE))\n", + "plt.errorbar(-sigma_vs_phi[:,1]*const.elementary_charge/const.nano**2, phi_SI, xerr=sigma_vs_phi[:,2]*const.elementary_charge/const.nano**2, fmt='o',label='Sim')\n", + "plt.plot(graham_sigma(x)*const.elementary_charge/const.nano**2,\n", + " x / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)), label='Graham')\n", + "x = np.linspace(0,ax.get_ylim()[1])\n", + "plt.plot(EPSILON_R*const.epsilon_0*x/2/(DEBYE_LENGTH*const.nano),x, label='linear PB', ls='--')\n", + "plt.xlabel(r'$\\sigma\\,\\mathrm{[C/m^2]}$')\n", + "plt.ylabel(r'$\\Psi_0\\,\\mathrm{[V]}$')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, + "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 (ipykernel)", + "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.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index e93755f95ab..b79ce7c9638 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 00000000000..6724cf8afe6 --- /dev/null +++ b/testsuite/scripts/tutorials/test_electrodes_1.py @@ -0,0 +1,40 @@ +# Copyright (C) 2019-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 . + +import unittest as ut +import importlib_wrapper +import numpy as np + +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/electrodes/electrodes_part1.py", + script_suffix="@TEST_SUFFIX@") + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + + def test_force(self): + msg = 'The force for particle 1 should agree with the analytical expression.' + np.testing.assert_allclose(tutorial.elc_forces_axial[:, 0], + tutorial.analytic_force_centered(tutorial.r / tutorial.BJERRUM_LENGTH, tutorial.box_l_z), rtol=1, err_msg=msg) + msg = 'The force for particle 2 should agree with the analytical expression.' + np.testing.assert_allclose(-tutorial.elc_forces_axial[:, 1], + tutorial.analytic_force_centered(tutorial.r / tutorial.BJERRUM_LENGTH, tutorial.box_l_z), rtol=1, 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 00000000000..355f69e4b09 --- /dev/null +++ b/testsuite/scripts/tutorials/test_electrodes_2.py @@ -0,0 +1,68 @@ +# Copyright (C) 2019-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 . + +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 (within 30 %) + 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.1, 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.gouy_chapman_density(tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, -tutorial.POTENTIAL_DIFF / 2.) + + tutorial.gouy_chapman_density(tutorial.box_l_z - tutorial.HS_ION_SIZE - tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, tutorial.POTENTIAL_DIFF / 2.)) / 2. + msg = 'The density profile is not sufficiently equal to PB theory.' + np.testing.assert_allclose( + charge_profile, + analytic, + rtol=5e-2, + atol=5e-2, + err_msg=msg) + + def test_capacitance(self): + # For low potentials the capacitance should be in line with Graham/DH + # equilibration limiting (2.5 minutes total) + graham = -tutorial.sigma_vs_phi[:, 0] / ( + constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE)) + msg = 'The capacitance at low potentials should be in line with Graham/DH.' + np.testing.assert_allclose( + graham, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg) + + +if __name__ == "__main__": + ut.main() From 25f7267f07c4d74c5fc73cec1754790c08118366 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Mon, 2 Oct 2023 21:50:36 +0200 Subject: [PATCH 2/3] Improve electrode tutorial --- .../electrodes/electrodes_part1.ipynb | 338 ++++++----- .../electrodes/electrodes_part2.ipynb | 542 +++++++----------- .../scripts/tutorials/test_electrodes_1.py | 9 +- .../scripts/tutorials/test_electrodes_2.py | 8 +- 4 files changed, 424 insertions(+), 473 deletions(-) diff --git a/doc/tutorials/electrodes/electrodes_part1.ipynb b/doc/tutorials/electrodes/electrodes_part1.ipynb index 06146bdf413..21e2f8861a3 100644 --- a/doc/tutorials/electrodes/electrodes_part1.ipynb +++ b/doc/tutorials/electrodes/electrodes_part1.ipynb @@ -49,18 +49,17 @@ "## Theoretical Background \n", "\n", "The normal component of electric field across a surface dividing two dielectric\n", - "media yields a discontinuity, which can be expressed in terms of a finite\n", - "surface charge density \n", + "media in presence of a surface charge density $\\sigma$ is discontinuous and follows as\n", "$(\\epsilon_1\\vec{E}_1 - \\epsilon_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 $\\epsilon_1$ to another one,\n", "$\\epsilon_2$.\n", "\n", "While in the case of non-polarizable materials ($\\epsilon_1 = \\epsilon_2 = 1$),\n", - "this jump is only related to surface charges and dielectric contrast and the\n", + "this jump is only related to surface charges and the\n", "potential is continuous across the interface, for polarizable materials also the\n", "polarization field $\\vec{P}$ will give a contribution. \n", - "In order to solve the problem in presence of a jump of the dielectric constant\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", @@ -70,18 +69,19 @@ "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\n", - "*scheme for solving the boundary integrals and the induced charge. \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$, **ESPResSo** offers the \"Electrostatic layer\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.\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.\n", "\n", + "### 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 the dielectric contrast is\n", + "In the metallic limit ($\\epsilon_2 \\to\\infty$) the dielectric contrast is\n", "$$ \\Delta = \\frac{\\epsilon_1 - \\epsilon_2} {\\epsilon_1 + \\epsilon_2} = -1 .$$\n", "If the ions are placed in the center of the slab of width $w$ and a distance $r$\n", "away from each other, the Green's function accounting for all image charges\n", @@ -114,7 +114,7 @@ "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 such a partially periodic system, we combine the efficient scaling behavior\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", @@ -137,10 +137,10 @@ "pillbox, which will give the induced surface charge in this infinitesimal\n", "segment as (Gauss law):\n", "$$q_{ind} = \\frac{1}{4\\pi} \\oint\\, dA\\, \\cdot \\epsilon\\nabla \\phi = \\frac{A}{4\\pi}(\\epsilon_1\\vec{E}_1 \\cdot \\hat{n}-\\epsilon_2\\vec{E}_2 \\cdot\\hat{n})$$\n", - "- The electric field at the closest proximity of the interface, $\\vec{E}_{1/2}$,\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/2} =\\vec{E} + 2\\pi/\\epsilon_1\\sigma\\hat{n} $$\n", + "$$ \\vec{E}_{1} =\\vec{E} + 2\\pi/\\epsilon_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{\\epsilon_1}{2\\pi} \\frac{\\epsilon_1-\\epsilon_2}{\\epsilon_1+\\epsilon_2}\\vec{E} \\cdot \\hat{n} =: \\frac{\\epsilon_1}{2\\pi} \\Delta \\, \\vec{E} \\cdot \\hat{n}$$\n", "\n", @@ -148,7 +148,7 @@ "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/2}$ and $\\sigma$ above until a self-consistent\n", + "expressions for $\\vec{E}_{1}$ and $\\sigma$ above until a self-consistent\n", "solution is found." ] }, @@ -172,17 +172,15 @@ "metadata": {}, "outputs": [], "source": [ - "import espressomd\n", + "import matplotlib.pyplot as plt\n", + "import tqdm\n", "import numpy as np\n", - "import espressomd.electrostatics\n", - "import espressomd.electrostatic_extensions\n", - "from espressomd.interactions import *\n", "\n", + "import espressomd\n", "espressomd.assert_features(['ELECTROSTATICS'])\n", "\n", - "import matplotlib.pyplot as plt\n", - "from scipy.special import *\n", - "from tqdm import tqdm" + "import espressomd.electrostatics\n", + "import espressomd.electrostatic_extensions" ] }, { @@ -192,13 +190,12 @@ "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), time ($t=0.01 \\tau$) and\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$ and thus defines the energy scale in our\n", - "system.\n", + "the thermal energy $k_\\mathrm{B}T$.\n", "It is defined as\n", "$$\\ell_\\mathrm{B}=\\frac{1}{4\\pi\\epsilon_0\\epsilon_r k_\\mathrm{B}T}.$$ \n", "In our case, if we choose the ion size ($\\sigma$) in simulations to represent a\n", @@ -213,48 +210,67 @@ "metadata": {}, "outputs": [], "source": [ - "#***************************************************\n", - "# System Setup\n", - "#***************************************************\n", - "\n", "# 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", - "# ICC* with ELC: 2D electrostatics\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\n", "system.time_step = 0.01\n", "system.cell_system.skin = 0.4\n", "\n", "# Elementary charge \n", - "q = np.array([1.0]) \n", - "\n", - "# Interaction Parameters for P3M with ELC\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", - "CHECK_ACCURACY = 1e-7 # maximim pairwise error in ELC\n", - "\n", - "#Lennard-Jones Parameters\n", - "\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", + "# Particle parameters\n", + "TYPES = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", + "charges = {\"Cation\": q, \"Anion\": -q }\n", "\n", - "types = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", - "charges = {\"Cation\": q[0], \"Anion\": -q[0] }\n", - "\n", - "p1=system.part.add(pos=[box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges[\"Cation\"])\n", - "print(f\"Cation placed at position: {p1.pos}\")\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\"])\n", - "print(f\"Anion placed at position: {p2.pos}\")\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", + " )" ] }, { @@ -264,105 +280,153 @@ "solution2_first": true }, "source": [ - "### Task:\n", - "Set up electrostatics interactions between the ion pairs using P3M and ICC $\\star$. \n", - "#### Hints:\n", - "- First instantiate a P3M object (using\n", - " [*espressomd.electrostatics.P3M(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.P3M))\n", - " as the 3D Coulomb solver for the system.\n", - "- Since, we are working on a 2D+h partially periodic system, we use ELC in\n", - " combination with P3M. For this purpose, instantiate ELC (using\n", - " [*espressomd.electrostatics.ELC(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.ELC))\n", - " with p3m as actor, which calls necessary initialization routines. \n", - "- Now we account for dielectric interface by adding ICC $^\\star$.\n", - " To this end, first create an ICC object using\n", - " [*espressomd.electrostatic_extensions.ICC(...)*](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatic_extensions.ICC).\n", - "- Setting up ICC $^\\star$ requires specifying a fixed number of ICC particles\n", - " and their initial charges.\n", - " Since the ICC particles are normal ESPResSo particles, they need to be added\n", - " at the interfaces using\n", - " [*system.part.add(..)*](https://espressomd.github.io/doc/espressomd.html?#module-espressomd.particle_data).\n", - "\n", - "*Note* - Refer to section\n", - "[**Dielectric interfaces with the ICC algorithm**](https://espressomd.github.io/doc/electrostatics.html#dielectric-interfaces-with-the-icc-star-algorithm)\n", - "in the **ESPResSo** documentation for the basics of an ICC $^\\star$ setup.\n" + "### 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": [ - "p3m = espressomd.electrostatics.P3M(\n", - " prefactor=BJERRUM_LENGTH,\n", - " accuracy=ACCURACY,\n", - " check_neutrality = False,\n", - " mesh = [100,100,150],\n", - " cao = 5\n", - " )\n", - "\n", - "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=CHECK_ACCURACY)\n", - "\n", - "# Set the ICC line density and calculate the number of\n", - "# ICC particles according to the box size\n", - "nicc = 100 # linear density\n", - "nicc_per_electrode = nicc**2 # surface density\n", - "nicc_tot = 2 * nicc_per_electrode\n", - "iccArea = box_l_x * box_l_y / nicc_per_electrode\n", - "lx = box_l_x / nicc\n", - "ly = box_l_y / nicc\n", - "\n", - "# Lists to collect required parameters\n", - "iccNormals = []\n", - "iccAreas = []\n", - "iccSigmas = []\n", - "iccEpsilons = []\n", - "\n", - "# Add the fixed ICC particles:\n", - "\n", - "# Bottom electrode (normal [0, 0, 1])\n", - "for xi in range(nicc):\n", - " for yi in range(nicc):\n", - " system.part.add(pos=[lx * xi, ly * yi, 0.], q=-0.0001,\n", - " type=types[\"Electrodes\"], fix=[True, True, True])\n", - "iccNormals.extend([0, 0, 1] * nicc_per_electrode)\n", - "\n", - "# Top electrode (normal [0, 0, -1])\n", - "for xi in range(nicc):\n", - " for yi in range(nicc):\n", - " system.part.add(pos=[lx * xi, ly * yi, box_l_z], q=0.0001,\n", - " type=types[\"Electrodes\"], fix=[True, True, True])\n", - "iccNormals.extend([0, 0, -1] * nicc_per_electrode)\n", + "### 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", - "iccAreas.extend([iccArea] * nicc_tot)\n", - "iccSigmas.extend([0] * nicc_tot)\n", - "iccEpsilons.extend([100000] * nicc_tot)\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", "\n", - "iccNormals = np.array(iccNormals, dtype=float).reshape(nicc_tot, 3)\n", + "icc_normals = np.array([[0,0,1] for _ in range(n_icc_partcls//2)] + [[0,0,-1] for _ in range(n_icc_partcls//2)])\n", "\n", "icc = espressomd.electrostatic_extensions.ICC(\n", - " first_id=2,\n", - " n_icc=nicc_tot,\n", - " convergence=1e-3,\n", - " relaxation=0.95,\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=1000,\n", - " eps_out=1.0,\n", - " normals=iccNormals,\n", - " areas=np.array(iccAreas, dtype=float),\n", - " sigmas=np.array(iccSigmas, dtype=float),\n", - " epsilons=np.array(iccEpsilons, dtype=float)\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", - "\n", - "system.electrostatics.solver = elc\n", - "system.electrostatics.extension = icc" + "system.electrostatics.extension = icc\n", + "```" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -378,18 +442,22 @@ "source": [ "r = np.logspace(0,box_l_z/4.,10) \n", "elc_forces_axial = np.empty((len(r), 2))\n", + "n_icc_per_electrode = len(icc_partcls_top)\n", + "\n", + "p1.pos = [0, box_l_y/2.0, box_l_z/2.0]\n", "\n", - "for i, x in enumerate(tqdm(r)):\n", - " p1.pos = [0, box_l_y/2.0, box_l_z/2.0]\n", + "for i, x in enumerate(tqdm.tqdm(r)): \n", " p2.pos = [x, box_l_y/2.0, box_l_z/2.0]\n", "\n", " system.integrator.run(0)\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 check passes\n", - " system.part.by_ids(range(2,2+nicc_per_electrode)).q = np.array([-0.0001]*nicc_per_electrode)\n", - " system.part.by_ids(range(2+nicc_per_electrode,2+2*nicc_per_electrode)).q = np.array([0.0001]*nicc_per_electrode) " + " # 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" ] }, { @@ -416,11 +484,11 @@ " return (-1)**n * x/(x**2+n**2)**(3./2)\n", " \n", " def do_sum(x):\n", - " max = int(1e3)\n", - " sum = 0\n", - " for n in range(-max+1,max+1):\n", - " sum += summand(x,n)\n", - " return sum\n", + " limit = int(1e3)\n", + " sum_ = 0\n", + " for n in range(-limit+1,limit+1):\n", + " sum_ += summand(x,n)\n", + " return sum_\n", "\n", " F = do_sum(x) * prefactor\n", " return F\n", @@ -478,7 +546,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -492,7 +560,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb index a22fadbbc40..61a0e18b4c5 100644 --- a/doc/tutorials/electrodes/electrodes_part2.ipynb +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -2,10 +2,7 @@ "cells": [ { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "# Basic simulation of electrodes in ESPResSo part II:\n", "# Electrolyte capacitor and Poisson-Boltzmann theory" @@ -13,10 +10,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## Prerequisites\n", "\n", @@ -34,10 +28,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## Introduction\n", "\n", @@ -48,7 +39,7 @@ "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 of electrodes of\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", @@ -56,13 +47,13 @@ "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 **ESPResSo**'s\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 homogeneous dielectric background.\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 against with the analytical Gouy-Chapmann\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", @@ -76,18 +67,15 @@ "\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 on the fly.\n", - "This is dealt in a computational cost-effective way using the ELC-IC method to\n", - "treat the image charge effect in the presence of 2D dielectric bounding\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": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## Theoretical Background \n", "\n", @@ -101,14 +89,17 @@ "Gouy [3] and Chapman [4] derived in the\n", "early 20th century the analytic solution for the case of a single planar wall\n", "withing 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", - "$$\\psi(z) = -2\\ln\\left[\n", - " \\frac{1-\\tanh(\\psi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D} z}}\n", - " {1+\\tanh(\\psi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D}z}} \\right].$$ \n", - "Here, $\\psi_\\mathrm{s}=\\psi(z=0)=$ const is the surface potential considering\n", - "that $\\psi(z\\rightarrow \\infty)=0$ and\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{\\epsilon \\, 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", @@ -123,9 +114,9 @@ "$$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", - "$\\psi_\\mathrm{s}=-4\\tanh^{-1}(\\gamma)$.\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(\\infty) = n_\\pm^{(0)}$\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", @@ -137,10 +128,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## ELC-IC for 2D+h periodic systems with dielectric interfaces\n", "\n", @@ -151,11 +139,11 @@ "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 accounts for the contributions from the periodic images in the\n", - "constrained direction and its numerical cost grows linear with the number of\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 $3-D$\n", "Coulomb solver.\n", - "Furthermore, ELC can be extended straightforward to implement metallic boundary\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", @@ -165,7 +153,7 @@ "The total potential drop $\\Delta \\phi$ across the simulated system is readily\n", "obtained from the ion distribution and integrating twice over the one\n", "dimensional Poisson's equation:\n", - "$$-\\epsilon_{0}\\epsilon_{r}\\phi_\\mathrm{ind}(z)=\\iint_{0}^{z}\\rho(z^{'})(z-z^{'})dz^{'} $$\n", + "$$-\\epsilon_{0}\\epsilon_{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", @@ -181,13 +169,13 @@ "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", + "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 $P_z$, the constant potential correction requires no\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 ab iterative numerical scheme with\n", + "[2] which employs ab 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", @@ -196,10 +184,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## 1. System setup \n", "\n", @@ -209,44 +194,37 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "import espressomd\n", "import numpy as np\n", - "rng = np.random.default_rng(42)\n", + "import matplotlib.pyplot as plt\n", + "from scipy import constants as const\n", + "from tqdm import tqdm\n", + "from scipy.stats import sem\n", + "\n", + "import espressomd\n", + "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", "\n", "import espressomd.electrostatics\n", "import espressomd.electrostatic_extensions\n", - "from espressomd.interactions import *\n", "import espressomd.observables\n", "import espressomd.accumulators\n", "from espressomd import shapes\n", "\n", - "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", - "\n", - "import matplotlib.pyplot as plt\n", - "from scipy import constants as const\n", - "from tqdm import tqdm\n", - "from scipy.stats import sem" + "rng = np.random.default_rng(42)" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "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) and time $t=0.01 \\tau$.\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", @@ -256,43 +234,34 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# water at room temperature\n", - "EPSILON_R = 78.4 # Relative dielectric constant of the water\n", + "EPSILON_R = 78.4 # Relative dielectric constant of water\n", "TEMPERATURE = 300.0 # Temperature in Kelvin\n", "BJERRUM_LENGTH = const.elementary_charge**2 / (4*np.pi*const.epsilon_0*EPSILON_R*const.Boltzmann*TEMPERATURE) / const.nano\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 nanometers, not point-like\n", + "#Lennard-Jones parameters\n", + "LJ_SIGMA = 0.3 # Particle size nanometers\n", "LJ_EPSILON = 1.0\n", - "HS_ION_SIZE = 2**(1/6) * LJ_SIGMA\n", "\n", - "CONCENTRATION = 1e-2 # desired concentration 10 mmol/l\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 = np.array([1.0]) \n", + "q = 1.0 \n", "types = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", - "charges = {\"Cation\": q[0], \"Anion\": -q[0] }" + "charges = {\"Cation\": q, \"Anion\": -q }" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true, - "solution2": "hidden", - "solution2_first": true - }, + "metadata": {}, "source": [ "### 1.1 Setting up the box dimensions and create system\n", "\n", @@ -313,35 +282,42 @@ "\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 (`HS_ION_SIZE`) of about 0.33 nm is a\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.\n", + "ion size." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task\n", "\n", - "The first task now is to write a function \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*d^\\prime)$.\n" + "concentration is $\\rho = N/(A*d^\\prime)$." ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "source": [ + "```python\n", "def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):\n", " \"\"\" For a given number of particles, determine the lateral area of the box\n", " to match the desired concentration \"\"\"\n", @@ -352,29 +328,24 @@ " l_z = distance * debye_length\n", " \n", " box_volume = n_ionpairs / rho\n", - " area = box_volume / (l_z - 2*HS_ION_SIZE) # account for finite ion size in density calculation\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" + " return l_xy, l_z\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "box_l_xy, box_l_z = get_box_dimension(CONCENTRATION,DISTANCE,N_IONPAIRS)\n", @@ -386,10 +357,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "We now can create the **ESPResSo** system.\n", "\n", @@ -406,10 +374,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "ELC_GAP = 6*box_l_z\n", @@ -419,23 +384,26 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true, - "solution2": "hidden", - "solution2_first": true - }, + "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 \n", - "\n", - "First, we add two wall constraints at $z=0$ and $z=L_z$ to stop particles from\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", - "\n", "Refer to \n", "[espressomd.constraints.ShapeBasedConstraint](https://espressomd.github.io/doc/espressomd.html#espressomd.constraints.ShapeBasedConstraint)\n", "and its\n", @@ -445,15 +413,12 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "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", @@ -463,47 +428,42 @@ "ceil = espressomd.shapes.Wall(normal=[0, 0, -1],\n", " dist=-box_l_z) \n", "c2 = system.constraints.add(\n", - " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceil)" + " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceil)\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden", "solution2_first": true }, "source": [ "#### 1.2.2 Add particles for the ions\n", "\n", - "Now, place the ion pairs at random positions between the electrodes.\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": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "source": [ - "offset=HS_ION_SIZE # avoid unfavorable overlap at close to the walls\n", + "```python\n", + "offset=LJ_SIGMA # avoid unfavorable overlap at close to the walls\n", "Init_part_btw_z1=0+offset \n", "Init_part_btw_z2=box_l_z-offset\n", "ion_pos=np.empty((3),dtype=float)\n", @@ -518,69 +478,59 @@ " 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\"])" + " system.part.add(pos=ion_pos, type=types[\"Anion\"] , q=charges[\"Anion\"])\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden", "solution2_first": true }, "source": [ "#### 1.2.3 Add interactions:\n", - "For excluded volume interactions, we add a WCA potential. \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.\n", - "\n", - "Add the corresponding interaction parameters to the system." + "section." ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "source": [ - "for key, val in types.items():\n", - " for key1, val1 in types.items():\n", - " system.non_bonded_inter[val, val1].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)" + "```python\n", + "for val in types.values():\n", + " for val1 in types.values():\n", + " system.non_bonded_inter[val, val1].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden", "solution2_first": true }, @@ -598,20 +548,19 @@ "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", - "Write a function `setup_electrostatic_solver(potential_diff)` that\n", + "### Task\n", + "\n", + "* Write a function `setup_electrostatic_solver(potential_diff)` that\n", "returns the ELC instance." ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "source": [ + "```python\n", "def setup_electrostatic_solver(potential_diff):\n", "\n", " delta_mid_top = -1.0 #(Fully metallic case both -1) \n", @@ -632,25 +581,20 @@ " delta_mid_bot=delta_mid_bot,\n", " delta_mid_top=delta_mid_top)\n", " \n", - " return elc" + " return elc\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "Now add the solver to the system:" ] @@ -658,10 +602,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)" @@ -669,10 +610,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## 2. Equilibration\n", "\n", @@ -694,44 +632,27 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "# Relax the overlaps with steepest descent\n", - "\n", "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() # Switch bach to Velocity Verlet \n", - "\n", - "# Add thermostat \n", - "thermostat_seed = np.random.randint(np.random.randint(1000000))\n", - "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=thermostat_seed)" + "system.integrator.set_vv()\n", + "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ - "## Equilibrate the ion distribution\n", - "\n", - "Convergence after $t\\sim25$ time units, possible to run up to $t=100$ here...\n", - "This is a total of 10.000 time steps (~1 minute)." + "## Equilibrate the ion distribution" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# Equlibration parameters\n", @@ -754,10 +675,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# Plot the convergence of the total energy\n", @@ -770,11 +688,16 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convergence after $t\\sim50$ time units." + ] + }, { "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden", "solution2_first": true }, @@ -790,7 +713,9 @@ "The time average is obtained through a\n", "[espressomd.accumulators.MeanVarianceCalculator](espressomd.accumulators.MeanVarianceCalculator).\n", "\n", - "Write a function `setup_densityprofile_accumulators(bin_width)` that returns the\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", @@ -798,34 +723,20 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "deletable": true, - "editable": true, "solution2": "hidden" }, - "outputs": [], "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", - " Ion_id=[]\n", - " Cations = system.part.select(type=types[\"Cation\"])\n", - " Cations_id=[]\n", - " for i in Cations:\n", - " Cations_id.append(i.id)\n", - " Ion_id.append(i.id)\n", - " \n", - " Anions = system.part.select(type=types[\"Anion\"])\n", - " Anions_id=[]\n", - " for i in Anions:\n", - " Anions_id.append(i.id)\n", - " Ion_id.append(i.id)\n", " \n", " n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))\n", " \n", - " # Accumulator 1 : observable::Density_Profile\n", - " density_profile_cation = espressomd.observables.DensityProfile(ids=Cations_id,\n", + " density_profile_cation = espressomd.observables.DensityProfile(ids=cations.id,\n", " n_x_bins=1,\n", " n_y_bins=1,\n", " n_z_bins=n_z_bins,\n", @@ -836,10 +747,11 @@ " max_y=system.box_l[1],\n", " max_z=system.box_l[2] - ELC_GAP)\n", " \n", - " density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(obs=density_profile_cation, delta_N=20)\n", + " density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(\n", + " obs=density_profile_cation, delta_N=20)\n", " \n", " \n", - " density_profile_anion = espressomd.observables.DensityProfile(ids=Anions_id,\n", + " density_profile_anion = espressomd.observables.DensityProfile(ids=anions.id,\n", " n_x_bins=1,\n", " n_y_bins=1,\n", " n_z_bins=n_z_bins,\n", @@ -850,41 +762,35 @@ " max_y=system.box_l[1],\n", " max_z=system.box_l[2] - ELC_GAP)\n", " \n", - " density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(obs=density_profile_anion, delta_N=20)\n", + " density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(\n", + " obs=density_profile_anion, delta_N=20)\n", "\n", " zs = density_profile_anion.bin_centers()[0, 0, :, 2]\n", " \n", - " return zs, density_accumulator_cation, density_accumulator_anion" + " return zs, density_accumulator_cation, density_accumulator_anion\n", + "```" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(bin_width = DEBYE_LENGTH/10.)" + "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(\n", + " bin_width = DEBYE_LENGTH/10.)" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "### 3.2 Run the simulation\n", "\n", @@ -894,10 +800,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "N_SAMPLES_PROD = 10\n", @@ -921,24 +824,18 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "### Compare to analytical prediction\n", "\n", - "Since we assume pair-wise additivity, the total ion density follows from\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": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def gouy_chapman_potential(x, debye_length, phi_0):\n", @@ -954,10 +851,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10, 6))\n", @@ -965,14 +859,14 @@ "plt.plot(zs, anion_profile_mean, color='r', label='anion')\n", "plt.plot(zs, cation_profile_mean + anion_profile_mean, color='k', label='total')\n", "\n", - "x = np.linspace(HS_ION_SIZE, box_l_z-HS_ION_SIZE, 100)\n", + "x = np.linspace(LJ_SIGMA, box_l_z-LJ_SIGMA, 100)\n", "plt.plot(x, (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", - " + gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2., color='r', ls='--')\n", - "plt.plot(x, (gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2., color='r', ls='--')\n", + "plt.plot(x, (gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='b', ls='--')\n", "plt.plot(x, (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", - " + gouy_chapman_density(box_l_z-HS_ION_SIZE-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2 \\\n", - " + (gouy_chapman_density(box_l_z-HS_ION_SIZE-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", + " + (gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='k', ls='--', lw=2)\n", "\n", "plt.legend()\n", @@ -983,10 +877,14 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "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 Graham equation.\n", "To this end we calculate \n", @@ -996,15 +894,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "#Test sigma with Graham equation\n", - "sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:int(len(zs)/2.)]) * (zs[1]-zs[0])\n", - "sigma_right = np.sum((+cation_profile_mean-anion_profile_mean)[int(len(zs)/2.):]) * (zs[1]-zs[0])\n", + "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 graham_sigma(phi):\n", " return np.sinh(phi/4.) * np.sqrt(2*rho/(np.pi*BJERRUM_LENGTH))\n", @@ -1016,10 +910,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "The electric field is readily obtained from the integral \n", "$$E(z) = \\int_0^{z} \\frac{1}{\\epsilon_0 \\epsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$" @@ -1028,10 +919,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# plot the electric field\n", @@ -1055,21 +943,22 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ - "And the electric potential from $\\phi(z) = \\int_0^z -E(z^\\prime)\\,\\mathrm{d}z^\\prime$." + "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": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# plot the elecrostatic potential\n", @@ -1093,44 +982,35 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "measured_potential_difference = -(phi[-1]+phi[0])\n", + "measured_potential_difference = -(phi[-1]-phi[0])\n", "print('applied voltage', POTENTIAL_DIFF, 'measured voltage', measured_potential_difference,\n", " 'relative deviation:', 1-measured_potential_difference/POTENTIAL_DIFF)" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "## 4. Differential capacitance\n", "\n", - "With the above knowledge, we can now easily assess the \n", - "differential capacitance of the system, i.e. change the applied voltage\n", - "difference and determine the corresponding surface charge density." + "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": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sigma_vs_phi = []\n", "MIN_PHI = 0.5\n", - "MAX_PHI = 14\n", - "N_PHI = 10\n", + "MAX_PHI = 10\n", + "N_PHI = 7\n", "N_SAMPLES_EQUIL_CAP = 5\n", "N_SAMPLES_CAP = 5\n", "\n", @@ -1172,15 +1052,12 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "f, 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()*1.3)\n", + "x = np.linspace(0,sigma_vs_phi[:,0].max())\n", "phi_SI = sigma_vs_phi[:,0] / (const.elementary_charge / (const.Boltzmann * TEMPERATURE))\n", "plt.errorbar(-sigma_vs_phi[:,1]*const.elementary_charge/const.nano**2, phi_SI, xerr=sigma_vs_phi[:,2]*const.elementary_charge/const.nano**2, fmt='o',label='Sim')\n", "plt.plot(graham_sigma(x)*const.elementary_charge/const.nano**2,\n", @@ -1188,17 +1065,22 @@ "x = np.linspace(0,ax.get_ylim()[1])\n", "plt.plot(EPSILON_R*const.epsilon_0*x/2/(DEBYE_LENGTH*const.nano),x, label='linear PB', ls='--')\n", "plt.xlabel(r'$\\sigma\\,\\mathrm{[C/m^2]}$')\n", - "plt.ylabel(r'$\\Psi_0\\,\\mathrm{[V]}$')\n", + "plt.ylabel(r'$\\phi_\\mathrm{s}\\,\\mathrm{[V]}$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "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 deviations 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", @@ -1218,7 +1100,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -1232,7 +1114,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/testsuite/scripts/tutorials/test_electrodes_1.py b/testsuite/scripts/tutorials/test_electrodes_1.py index 6724cf8afe6..08419519e66 100644 --- a/testsuite/scripts/tutorials/test_electrodes_1.py +++ b/testsuite/scripts/tutorials/test_electrodes_1.py @@ -28,12 +28,13 @@ class Tutorial(ut.TestCase): def test_force(self): + # exclude largest distance (small force -> large relative error) msg = 'The force for particle 1 should agree with the analytical expression.' - np.testing.assert_allclose(tutorial.elc_forces_axial[:, 0], - tutorial.analytic_force_centered(tutorial.r / tutorial.BJERRUM_LENGTH, tutorial.box_l_z), rtol=1, err_msg=msg) + np.testing.assert_allclose(tutorial.elc_forces_axial[:-1, 0], + tutorial.analytic_force_centered(tutorial.r[:-1], tutorial.box_l_z), rtol=1e-1, err_msg=msg) msg = 'The force for particle 2 should agree with the analytical expression.' - np.testing.assert_allclose(-tutorial.elc_forces_axial[:, 1], - tutorial.analytic_force_centered(tutorial.r / tutorial.BJERRUM_LENGTH, tutorial.box_l_z), rtol=1, err_msg=msg) + np.testing.assert_allclose(-tutorial.elc_forces_axial[:-1, 1], + tutorial.analytic_force_centered(tutorial.r[:-1], tutorial.box_l_z), rtol=1e-1, err_msg=msg) if __name__ == "__main__": diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py index 355f69e4b09..e83393a40f8 100644 --- a/testsuite/scripts/tutorials/test_electrodes_2.py +++ b/testsuite/scripts/tutorials/test_electrodes_2.py @@ -34,10 +34,10 @@ class Tutorial(ut.TestCase): def test_potential_difference(self): # Test that the applied potential difference equals the one from - # integrating Poisson equation (within 30 %) + # 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.1, msg=msg) + 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!! @@ -45,7 +45,7 @@ def test_charge_profile(self): tutorial.cation_profile_mean + tutorial.anion_profile_mean) analytic = (tutorial.gouy_chapman_density(tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, -tutorial.POTENTIAL_DIFF / 2.) - + tutorial.gouy_chapman_density(tutorial.box_l_z - tutorial.HS_ION_SIZE - tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, tutorial.POTENTIAL_DIFF / 2.)) / 2. + + tutorial.gouy_chapman_density(tutorial.box_l_z - tutorial.LJ_SIGMA - tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, tutorial.POTENTIAL_DIFF / 2.)) / 2. msg = 'The density profile is not sufficiently equal to PB theory.' np.testing.assert_allclose( charge_profile, @@ -56,7 +56,7 @@ def test_charge_profile(self): def test_capacitance(self): # For low potentials the capacitance should be in line with Graham/DH - # equilibration limiting (2.5 minutes total) + # equilibration performance limiting graham = -tutorial.sigma_vs_phi[:, 0] / ( constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE)) msg = 'The capacitance at low potentials should be in line with Graham/DH.' From 58f088ef1693029381eefaf61fc42eef23c38bfe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 2 Oct 2023 22:10:47 +0200 Subject: [PATCH 3/3] Rework electrode tutorials Copy editing. Reduce code duplication. Improve matplotlib figures. --- doc/tutorials/electrodes/NotesForTutor.md | 7 +- .../electrodes/electrodes_part1.ipynb | 174 ++++----- .../electrodes/electrodes_part2.ipynb | 358 ++++++++---------- .../scripts/tutorials/test_electrodes_1.py | 21 +- .../scripts/tutorials/test_electrodes_2.py | 24 +- 5 files changed, 281 insertions(+), 303 deletions(-) diff --git a/doc/tutorials/electrodes/NotesForTutor.md b/doc/tutorials/electrodes/NotesForTutor.md index 3c87bb760b8..47497b59da4 100644 --- a/doc/tutorials/electrodes/NotesForTutor.md +++ b/doc/tutorials/electrodes/NotesForTutor.md @@ -8,9 +8,10 @@ * 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* +* Discrete image charges: ICC\* + +### Part 2 -## Part 2 * Nano-confined charged liquids as super-capacitors * Advanced Poisson-Boltzmann theory: Gouy-Chapman, Graham equation * Limits of PB: finite ion size, correlations, ... @@ -19,7 +20,7 @@ After the tutorial, students should be able to: -* Explain how ESPResSo implements 2d periodic electrostatics. +* 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. diff --git a/doc/tutorials/electrodes/electrodes_part1.ipynb b/doc/tutorials/electrodes/electrodes_part1.ipynb index 21e2f8861a3..acda4631df2 100644 --- a/doc/tutorials/electrodes/electrodes_part1.ipynb +++ b/doc/tutorials/electrodes/electrodes_part1.ipynb @@ -4,8 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Basic simulation of electrodes in ESPResSo part I:\n", - "# Ion-pair in a narrow metallic slit-like confinement using ICC $\\star$" + "# Basic simulation of electrodes in ESPResSo part I: ion-pair in a narrow metallic slit-like confinement using ICC$^\\star$" ] }, { @@ -18,10 +17,10 @@ "\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 respective tutorial\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", + " dipoles, surface and image charges\n", "- Reduced units, as described in the **ESPResSo** [user guide](https://espressomd.github.io/doc/introduction.html#on-units)" ] }, @@ -34,7 +33,7 @@ "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", + "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", @@ -50,63 +49,68 @@ "\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", - "$(\\epsilon_1\\vec{E}_1 - \\epsilon_2\\vec{E}_2).\\hat{n}=-\\sigma(\\vec{r})$.\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 $\\epsilon_1$ to another one,\n", - "$\\epsilon_2$.\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 ($\\epsilon_1 = \\epsilon_2 = 1$),\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 also the\n", - "polarization field $\\vec{P}$ will give a contribution. \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", + "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_{ind}$ on the surface. \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", + "(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", + "*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.\n", - "\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 ($\\epsilon_2 \\to\\infty$) the dielectric contrast is\n", - "$$ \\Delta = \\frac{\\epsilon_1 - \\epsilon_2} {\\epsilon_1 + \\epsilon_2} = -1 .$$\n", - "If the ions are placed in the center of the slab of width $w$ and a distance $r$\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 \\epsilon_0 \\epsilon_r w \\mathcal{G}(x) = \\sum_{n=-\\infty}^\\infty \\frac{-1^n}{\\sqrt{x^2+n^2}} ,$$\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 \\epsilon_0 \\epsilon_r w \\mathcal{G}(x) = \\frac{1}{x},$$\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 \\epsilon_0 \\epsilon_r w \\mathcal{G}(x) = \\sqrt{\\frac{8}{x}} e^{-\\pi x},$$\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." + "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", + "## 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", @@ -120,7 +124,7 @@ "[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 $3-D$\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", @@ -128,24 +132,23 @@ "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", + "We furthermore employ ICC$^\\star$ to solve the Poisson equation for an\n", "inhomogeneous dieletric:\n", - "$$ \\nabla (\\epsilon \\nabla \\phi)=-4\\pi \\rho$$\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_{ind} = \\frac{1}{4\\pi} \\oint\\, dA\\, \\cdot \\epsilon\\nabla \\phi = \\frac{A}{4\\pi}(\\epsilon_1\\vec{E}_1 \\cdot \\hat{n}-\\epsilon_2\\vec{E}_2 \\cdot\\hat{n})$$\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/\\epsilon_1\\sigma\\hat{n} $$\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{\\epsilon_1}{2\\pi} \\frac{\\epsilon_1-\\epsilon_2}{\\epsilon_1+\\epsilon_2}\\vec{E} \\cdot \\hat{n} =: \\frac{\\epsilon_1}{2\\pi} \\Delta \\, \\vec{E} \\cdot \\hat{n}$$\n", - "\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", + "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", @@ -177,10 +180,11 @@ "import numpy as np\n", "\n", "import espressomd\n", - "espressomd.assert_features(['ELECTROSTATICS'])\n", - "\n", "import espressomd.electrostatics\n", - "import espressomd.electrostatic_extensions" + "import espressomd.electrostatic_extensions\n", + "\n", + "espressomd.assert_features(['ELECTROSTATICS'])\n", + "plt.rcParams.update({'font.size': 18})" ] }, { @@ -190,16 +194,17 @@ "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", + "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\\epsilon_0\\epsilon_r k_\\mathrm{B}T}.$$ \n", - "In our case, if we choose the ion size ($\\sigma$) in simulations to represent a\n", - "typical value for mono-atomar salt, 0.3 nm in real units, then the\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." ] @@ -219,7 +224,7 @@ "# 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", + "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", @@ -246,8 +251,8 @@ "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])" + "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])" ] }, { @@ -267,9 +272,9 @@ "p3m = espressomd.electrostatics.P3M(\n", " prefactor=BJERRUM_LENGTH,\n", " accuracy=ACCURACY,\n", - " check_neutrality = False,\n", - " mesh = [100,100,150],\n", - " cao = 5\n", + " check_neutrality=False,\n", + " mesh=[100, 100, 150],\n", + " cao=5\n", " )" ] }, @@ -291,7 +296,7 @@ "solution2": "hidden" }, "source": [ - "```python \n", + "```python\n", "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)\n", "```" ] @@ -316,7 +321,7 @@ "metadata": {}, "outputs": [], "source": [ - "ICC_PARTCL_NUMBER_DENSITY = 1\n", + "ICC_PARTCL_NUMBER_DENSITY = 1.\n", "icc_partcls_bottom = []\n", "icc_partcls_top = []" ] @@ -331,7 +336,7 @@ "### 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", + "* 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``." ] }, @@ -341,7 +346,7 @@ "solution2": "hidden" }, "source": [ - "```python \n", + "```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", @@ -350,12 +355,12 @@ "# 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", + " 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", + " 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", "```" ] @@ -397,11 +402,10 @@ "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_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", - "\n", - "icc_normals = np.array([[0,0,1] for _ in range(n_icc_partcls//2)] + [[0,0,-1] for _ in range(n_icc_partcls//2)])\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", @@ -415,7 +419,7 @@ " areas=icc_areas,\n", " sigmas=icc_sigmas,\n", " epsilons=icc_epsilons\n", - " )\n", + ")\n", "system.electrostatics.extension = icc\n", "```" ] @@ -440,24 +444,25 @@ "metadata": {}, "outputs": [], "source": [ - "r = np.logspace(0,box_l_z/4.,10) \n", - "elc_forces_axial = np.empty((len(r), 2))\n", + "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.0, box_l_z/2.0]\n", + "p1.pos = [0., box_l_y / 2., box_l_z / 2.]\n", "\n", - "for i, x in enumerate(tqdm.tqdm(r)): \n", - " p2.pos = [x, box_l_y/2.0, box_l_z/2.0]\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)\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", + " part.q = 1. / n_icc_per_electrode\n", " for part in icc_partcls_bottom:\n", - " part.q = -1/n_icc_per_electrode" + " part.q = -1. / n_icc_per_electrode" ] }, { @@ -477,19 +482,18 @@ "outputs": [], "source": [ "def analytic_force_centered(r,w):\n", - " x = r/w\n", - " prefactor = BJERRUM_LENGTH / w**2\n", - "\n", " def summand(x, n):\n", - " return (-1)**n * x/(x**2+n**2)**(3./2)\n", + " return (-1)**n * x / (x**2 + n**2)**(3. / 2.)\n", " \n", " def do_sum(x):\n", - " limit = int(1e3)\n", - " sum_ = 0\n", - " for n in range(-limit+1,limit+1):\n", - " sum_ += summand(x,n)\n", - " return sum_\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", @@ -499,9 +503,9 @@ " return E\n", "\n", "def exponential_force(r,w):\n", - " x = 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", + " E = prefactor * np.sqrt(2.) * (1. / x)**(3. / 2.) * np.exp(-np.pi * x)\n", " return E" ] }, @@ -513,19 +517,19 @@ "source": [ "fig = plt.figure(figsize=(10, 6))\n", "\n", - "plt.plot(r/BJERRUM_LENGTH, -elc_forces_axial[:,1], color='red' , 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", - "x = np.logspace(-.25,1.45,100)\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(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='green', ls='--', label='Coulomb')\n", - "plt.plot(x/BJERRUM_LENGTH, exponential_force(x,box_l_z), color='red', ls='--', label='Exponential')\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()" + "plt.legend()\n", + "plt.show()" ] }, { diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb index 61a0e18b4c5..35f7723156a 100644 --- a/doc/tutorials/electrodes/electrodes_part2.ipynb +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -4,8 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Basic simulation of electrodes in ESPResSo part II:\n", - "# Electrolyte capacitor and Poisson-Boltzmann theory" + "# Basic simulation of electrodes in ESPResSo part II: Electrolyte capacitor and Poisson–Boltzmann theory" ] }, { @@ -18,10 +17,10 @@ "\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 respective tutorial\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", + " 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)." ] @@ -53,13 +52,13 @@ "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", + "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", + "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", @@ -79,7 +78,7 @@ "source": [ "## Theoretical Background \n", "\n", - "### Poisson-Boltzmann Theory\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", @@ -88,7 +87,7 @@ "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", - "withing the mean-field approximation of the Poisson Boltzmann (PB) equation. \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", @@ -101,7 +100,7 @@ "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{\\epsilon \\, k_{\\mathrm B} T}{\\sum_{j = 1}^N n_j^0 \\, q_j^2}\\right)^{1/2}, $$\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", @@ -121,8 +120,8 @@ "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 capactitance\n", - "$C=\\frac{\\mathrm{d}\\sigma}{\\mathrm{d}\\phi_\\mathrm{s}}$ within the mean-field\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." ] }, @@ -135,13 +134,13 @@ "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 two-dimensional a partially periodic system, we combine the efficient\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 $3-D$\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", @@ -151,9 +150,8 @@ "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\n", - "dimensional Poisson's equation:\n", - "$$-\\epsilon_{0}\\epsilon_{r}\\phi_\\mathrm{ind}(z)=\\iint_{0}^{z}\\rho(z^{\\prime})(z-z^{\\prime})dz^{\\prime} $$\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", @@ -163,9 +161,9 @@ "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 straightforward as the corresponding contribution from the induced\n", + "calculated straightforwardly, as the corresponding contribution from the induced\n", "charges is known:\n", - "$$ E_z^\\mathrm{(ind)} = \\frac{1}{\\epsilon_0 \\epsilon_r L^2 L_z} P_z$$\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", @@ -174,8 +172,8 @@ "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 ab iterative numerical scheme with\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", @@ -199,20 +197,20 @@ "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from scipy import constants as const\n", - "from tqdm import tqdm\n", - "from scipy.stats import sem\n", + "import scipy.constants as constants\n", + "import scipy.stats\n", + "import tqdm\n", "\n", "import espressomd\n", - "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", - "\n", "import espressomd.electrostatics\n", "import espressomd.electrostatic_extensions\n", "import espressomd.observables\n", "import espressomd.accumulators\n", - "from espressomd import shapes\n", + "import espressomd.shapes\n", "\n", - "rng = np.random.default_rng(42)" + "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", + "rng = np.random.default_rng(42)\n", + "plt.rcParams.update({'font.size': 18})" ] }, { @@ -224,7 +222,7 @@ "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", + "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", @@ -240,11 +238,12 @@ "# water at room temperature\n", "EPSILON_R = 78.4 # Relative dielectric constant of water\n", "TEMPERATURE = 300.0 # Temperature in Kelvin\n", - "BJERRUM_LENGTH = const.elementary_charge**2 / (4*np.pi*const.epsilon_0*EPSILON_R*const.Boltzmann*TEMPERATURE) / const.nano\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 nanometers\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", @@ -254,9 +253,9 @@ "POTENTIAL_DIFF = 5.0\n", "\n", "# Elementary charge \n", - "q = 1.0 \n", - "types = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", - "charges = {\"Cation\": q, \"Anion\": -q }" + "q = 1.\n", + "types = {\"Cation\": 0, \"Anion\": 1, \"Electrodes\": 2}\n", + "charges = {\"Cation\": q, \"Anion\": -q}" ] }, { @@ -286,7 +285,7 @@ "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", + "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." @@ -308,7 +307,7 @@ "\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*d^\\prime)$." + "concentration is $\\rho = N/(A \\cdot d^\\prime)$." ] }, { @@ -319,16 +318,18 @@ "source": [ "```python\n", "def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):\n", - " \"\"\" For a given number of particles, determine the lateral area of the box\n", - " to match the desired concentration \"\"\"\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 * (const.Avogadro / const.liter) * const.nano**3\n", - " debye_length = (4 * np.pi * BJERRUM_LENGTH * rho*2)**(-1./2) # desired Debye length in nm\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", + " 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", @@ -348,11 +349,11 @@ "metadata": {}, "outputs": [], "source": [ - "box_l_xy, box_l_z = get_box_dimension(CONCENTRATION,DISTANCE,N_IONPAIRS)\n", + "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" + "rho = N_IONPAIRS / (box_l_xy * box_l_xy * box_l_z) # in units of 1/nm^3" ] }, { @@ -365,7 +366,7 @@ "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 value of `ELC_GAP = 6*box_l_z` is a good choice.\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." @@ -377,8 +378,8 @@ "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", + "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" ] }, @@ -418,17 +419,16 @@ "solution2": "hidden" }, "source": [ - "```python \n", + "```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", - "ceil = espressomd.shapes.Wall(normal=[0, 0, -1],\n", - " dist=-box_l_z) \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=ceil)\n", + " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceiling)\n", "```" ] }, @@ -463,22 +463,22 @@ }, "source": [ "```python\n", - "offset=LJ_SIGMA # avoid unfavorable overlap at close to the walls\n", - "Init_part_btw_z1=0+offset \n", - "Init_part_btw_z2=box_l_z-offset\n", - "ion_pos=np.empty((3),dtype=float)\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", + " 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", + " 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", "```" ] }, @@ -515,9 +515,9 @@ }, "source": [ "```python\n", - "for val in types.values():\n", - " for val1 in types.values():\n", - " system.non_bonded_inter[val, val1].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)\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", "```" ] }, @@ -562,25 +562,21 @@ "source": [ "```python\n", "def setup_electrostatic_solver(potential_diff):\n", - "\n", - " delta_mid_top = -1.0 #(Fully metallic case both -1) \n", - " delta_mid_bot = -1.0\n", - "\n", + " delta_mid_top = -1. # (Fully metallic case both -1) \n", + " delta_mid_bot = -1.\n", " accuracy = 1e-3\n", - " check_accuracy = 1e-3\n", + " elc_accuracy = 1e-3\n", " p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,\n", - " accuracy=accuracy, \n", + " accuracy=accuracy,\n", " tune=True,\n", - " )\n", - " \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=check_accuracy,\n", + " maxPWerror=elc_accuracy,\n", " delta_mid_bot=delta_mid_bot,\n", " delta_mid_top=delta_mid_top)\n", - " \n", " return elc\n", "```" ] @@ -616,17 +612,11 @@ "\n", "### 2.1 Steepest descent\n", "\n", - "Before we can start the simulation, we need to remove the overlap between\n", - "particles to avoid large forces which would crash the simulation.\n", - "\n", - "For this, we use the steepest descent integrator with a relative convergence\n", - "criterion for forces and energies.\n", - "\n", - "After steepest descent, we switch to a Velocity Verlet integrator and set up a\n", - "Langevin thermostat.\n", - "Note, that we only analyze static properties, thus the damping and temperature\n", - "chosen here only determine the relaxation speed towards the equilibrium\n", - "distribution." + "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." ] }, { @@ -658,13 +648,13 @@ "# Equlibration parameters\n", "STEPS_PER_SAMPLE = 200\n", "N_SAMPLES_EQUIL = 50\n", - "N_PART = 2* N_IONPAIRS\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(range(N_SAMPLES_EQUIL)):\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", @@ -682,8 +672,8 @@ "plt.figure(figsize=(10, 6))\n", "plt.plot(times, e_total, label='total')\n", "plt.plot(times, e_kin, label='kinetic')\n", - "plt.xlabel('t')\n", - "plt.ylabel('E')\n", + "plt.xlabel('Simulation time')\n", + "plt.ylabel('Energy')\n", "plt.legend()\n", "plt.show()" ] @@ -728,45 +718,22 @@ "solution2": "hidden" }, "source": [ - "```python \n", + "```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", - " \n", " n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))\n", - " \n", - " density_profile_cation = espressomd.observables.DensityProfile(ids=cations.id,\n", - " n_x_bins=1,\n", - " n_y_bins=1,\n", - " n_z_bins=n_z_bins,\n", - " min_x=0,\n", - " min_y=0,\n", - " min_z=0,\n", - " max_x=system.box_l[0],\n", - " max_y=system.box_l[1],\n", - " max_z=system.box_l[2] - ELC_GAP)\n", - " \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", - " \n", - " \n", - " density_profile_anion = espressomd.observables.DensityProfile(ids=anions.id,\n", - " n_x_bins=1,\n", - " n_y_bins=1,\n", - " n_z_bins=n_z_bins,\n", - " min_x=0,\n", - " min_y=0,\n", - " min_z=0,\n", - " max_x=system.box_l[0],\n", - " max_y=system.box_l[1],\n", - " max_z=system.box_l[2] - ELC_GAP)\n", - " \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", - "\n", " zs = density_profile_anion.bin_centers()[0, 0, :, 2]\n", - " \n", " return zs, density_accumulator_cation, density_accumulator_anion\n", "```" ] @@ -785,7 +752,7 @@ "outputs": [], "source": [ "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(\n", - " bin_width = DEBYE_LENGTH/10.)" + " bin_width=DEBYE_LENGTH / 10.)" ] }, { @@ -810,13 +777,13 @@ "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(range(N_SAMPLES_PROD)):\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", + " times.append(system.time)\n", " energy = system.analysis.energy()\n", - " e_total.append( energy['total']) \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, :]" @@ -839,13 +806,13 @@ "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", + " 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)" + " return c0 / 2. * np.exp(-phi)" ] }, { @@ -854,24 +821,33 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(10, 6))\n", - "plt.plot(zs, cation_profile_mean, color='b', label='cation')\n", - "plt.plot(zs, anion_profile_mean, color='r', label='anion')\n", - "plt.plot(zs, cation_profile_mean + anion_profile_mean, color='k', label='total')\n", + "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", - "plt.plot(x, (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., color='r', ls='--')\n", - "plt.plot(x, (gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", - " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='b', ls='--')\n", - "plt.plot(x, (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", - " + (gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", - " + gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) )/2., color='k', ls='--', lw=2)\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", - "plt.legend()\n", - "plt.xlabel(r'$z\\,\\mathrm{[nm]}$')\n", - "plt.ylabel(r'$\\rho(z)\\,\\mathrm{[mol/l]}$')\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()" ] }, @@ -886,7 +862,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now check how well the surface charge agrees with Graham equation.\n", + "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 .$$" ] @@ -897,15 +873,16 @@ "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", + "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 graham_sigma(phi):\n", - " return np.sinh(phi/4.) * np.sqrt(2*rho/(np.pi*BJERRUM_LENGTH))\n", - "sigma_graham = graham_sigma(POTENTIAL_DIFF)\n", + "def grahame_sigma(phi):\n", + " return np.sinh(phi / 4.) * np.sqrt(2. * rho / (np.pi * BJERRUM_LENGTH))\n", "\n", - "# sigma in e/nm^2\n", - "print('simulation:', sigma_right, 'graham:', sigma_graham, 'relative:', sigma_right/sigma_graham)\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}%')" ] }, { @@ -913,7 +890,7 @@ "metadata": {}, "source": [ "The electric field is readily obtained from the integral \n", - "$$E(z) = \\int_0^{z} \\frac{1}{\\epsilon_0 \\epsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$" + "$$E(z) = \\int_0^{z} \\frac{1}{\\varepsilon_0 \\varepsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$" ] }, { @@ -923,18 +900,18 @@ "outputs": [], "source": [ "# plot the electric field\n", - "f, ax = plt.subplots(figsize=(10, 6))\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", - "dz_SI = (zs[1]-zs[0])*const.nano\n", - "chargedensity = (cation_profile_mean-anion_profile_mean)*const.elementary_charge/const.nano**3 \n", - "E_SI = 1/(EPSILON_R*const.epsilon_0)* np.cumsum(chargedensity*dz_SI)\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 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE) / const.nano)\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.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", @@ -962,20 +939,17 @@ "outputs": [], "source": [ "# plot the elecrostatic potential\n", - "f, ax = plt.subplots(figsize=(10, 6))\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 * (const.elementary_charge / (const.Boltzmann * TEMPERATURE))\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", - "ax2.plot(zs, phi)\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", - "ax.axhline(-5 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)))\n", "ax2.axhline(0, ls='--', color='k')\n", - "ax.axhline(0 / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)))\n", - "ax.set_xlim(0, 10*DEBYE_LENGTH)\n", + "ax.set_xlim(0, 10. * DEBYE_LENGTH)\n", "plt.show()" ] }, @@ -985,9 +959,10 @@ "metadata": {}, "outputs": [], "source": [ - "measured_potential_difference = -(phi[-1]-phi[0])\n", - "print('applied voltage', POTENTIAL_DIFF, 'measured voltage', measured_potential_difference,\n", - " 'relative deviation:', 1-measured_potential_difference/POTENTIAL_DIFF)" + "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}%')" ] }, { @@ -1015,38 +990,26 @@ "N_SAMPLES_CAP = 5\n", "\n", "# sample from high to low potential to improve sampling\n", - "for potential_diff in tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n", - "\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", - "\n", - " times=[]\n", - " e_total=[]\n", + " system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)\n", " sigmas = []\n", - " for tm in range(N_SAMPLES_EQUIL_CAP):\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", " for tm in range(N_SAMPLES_CAP):\n", - "\n", - " zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(bin_width = DEBYE_LENGTH/10.)\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", - "\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, :]\n", "\n", - " sigmas.append(np.sum((cation_profile_mean-anion_profile_mean)[:int(len(zs)/2.)]) * (zs[1]-zs[0]))\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), sem(sigmas)]) " + " sigma_vs_phi.append([potential_diff, np.mean(sigmas), scipy.stats.sem(sigmas)]) " ] }, { @@ -1055,15 +1018,18 @@ "metadata": {}, "outputs": [], "source": [ - "f, ax = plt.subplots(figsize=(10, 6))\n", + "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] / (const.elementary_charge / (const.Boltzmann * TEMPERATURE))\n", - "plt.errorbar(-sigma_vs_phi[:,1]*const.elementary_charge/const.nano**2, phi_SI, xerr=sigma_vs_phi[:,2]*const.elementary_charge/const.nano**2, fmt='o',label='Sim')\n", - "plt.plot(graham_sigma(x)*const.elementary_charge/const.nano**2,\n", - " x / (const.elementary_charge / (const.Boltzmann * TEMPERATURE)), label='Graham')\n", - "x = np.linspace(0,ax.get_ylim()[1])\n", - "plt.plot(EPSILON_R*const.epsilon_0*x/2/(DEBYE_LENGTH*const.nano),x, label='linear PB', ls='--')\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", @@ -1074,8 +1040,8 @@ "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 deviations which can be attributed to the fact that our ions are of finite size and thus the surface charge saturates." + "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." ] }, { @@ -1088,7 +1054,7 @@ "\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", + "[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", diff --git a/testsuite/scripts/tutorials/test_electrodes_1.py b/testsuite/scripts/tutorials/test_electrodes_1.py index 08419519e66..eaa490e7a9c 100644 --- a/testsuite/scripts/tutorials/test_electrodes_1.py +++ b/testsuite/scripts/tutorials/test_electrodes_1.py @@ -1,4 +1,5 @@ -# Copyright (C) 2019-2022 The ESPResSo project +# +# Copyright (C) 2023 The ESPResSo project # # This file is part of ESPResSo. # @@ -14,13 +15,14 @@ # # 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", + "@TUTORIALS_DIR@/electrodes/electrodes_part1.py", N_AXIAL_POINTS=6, script_suffix="@TEST_SUFFIX@") @@ -28,13 +30,14 @@ class Tutorial(ut.TestCase): def test_force(self): - # exclude largest distance (small force -> large relative error) - msg = 'The force for particle 1 should agree with the analytical expression.' - np.testing.assert_allclose(tutorial.elc_forces_axial[:-1, 0], - tutorial.analytic_force_centered(tutorial.r[:-1], tutorial.box_l_z), rtol=1e-1, err_msg=msg) - msg = 'The force for particle 2 should agree with the analytical expression.' - np.testing.assert_allclose(-tutorial.elc_forces_axial[:-1, 1], - tutorial.analytic_force_centered(tutorial.r[:-1], tutorial.box_l_z), rtol=1e-1, err_msg=msg) + 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__": diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py index e83393a40f8..5654db1ed88 100644 --- a/testsuite/scripts/tutorials/test_electrodes_2.py +++ b/testsuite/scripts/tutorials/test_electrodes_2.py @@ -1,4 +1,5 @@ -# Copyright (C) 2019-2022 The ESPResSo project +# +# Copyright (C) 2023 The ESPResSo project # # This file is part of ESPResSo. # @@ -14,6 +15,7 @@ # # 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 @@ -37,31 +39,33 @@ def test_potential_difference(self): # 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) + 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.gouy_chapman_density(tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, -tutorial.POTENTIAL_DIFF / 2.) - + tutorial.gouy_chapman_density(tutorial.box_l_z - tutorial.LJ_SIGMA - tutorial.zs, tutorial.CONCENTRATION, tutorial.DEBYE_LENGTH, tutorial.POTENTIAL_DIFF / 2.)) / 2. + 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, - analytic, + 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 Graham/DH + # For low potentials the capacitance should be in line with Grahame/DH # equilibration performance limiting - graham = -tutorial.sigma_vs_phi[:, 0] / ( + grahame = -tutorial.sigma_vs_phi[:, 0] / ( constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE)) - msg = 'The capacitance at low potentials should be in line with Graham/DH.' + msg = 'The capacitance at low potentials should be in line with Grahame/DH.' np.testing.assert_allclose( - graham, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg) + grahame, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg) if __name__ == "__main__":