diff --git a/Readme.md b/Readme.md index 874df4bf58c..95048e8c80b 100644 --- a/Readme.md +++ b/Readme.md @@ -1,9 +1,3 @@ -# Invitation to the ESPResSo Summer School 2023 - -[![CECAM Flagship School registration link](https://img.shields.io/badge/CECAM%20Flagship%20School-Register%20Now-blue?style=for-the-badge)](https://www.cecam.org/workshop-details/1229) - -The summer school "Simulating energy materials with ESPResSo and waLBerla" will take place on October 9-13, 2023, in Stuttgart. Registration is now open on [CECAM](https://www.cecam.org/workshop-details/1229). - # ESPResSo [![GitLab CI](https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/badges/python/pipeline.svg)](https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/-/commits/python) diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 200bf0cf707..8eeb22df13d 100644 --- a/doc/bibliography.bib +++ b/doc/bibliography.bib @@ -673,21 +673,14 @@ @Article{kolb99a doi = {10.1063/1.479208}, } -@PhdThesis{kruger11a, -author = {Kr\"{u}ger, Timm}, -title = {Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear}, -school = {Universit\"{a}t Bochum}, -year = {2011} -} - @Book{kruger12a, -author={Kr\"{u}ger, Timm}, -title={Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear}, -year={2012}, -publisher={Vieweg+Teubner Verlag}, -address={Wiesbaden}, -isbn={978-3-8348-2376-2}, -doi={10.1007/978-3-8348-2376-2}, + author = {Kr{\"u}ger, Timm}, + title = {Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear}, + year = {2012}, + publisher = {Vieweg+Teubner Verlag}, + isbn = {978-3-8348-2376-2}, + doi = {10.1007/978-3-8348-2376-2}, + address = {Wiesbaden}, } @Book{kruger17a, diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index a2ccc248511..39a0935b25c 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -272,7 +272,7 @@ Please contact the Biofluid Simulation and Modeling Group at the University of Bayreuth if you plan to use this feature. With the Immersed Boundary Method (IBM), soft particles are considered as an infinitely -thin shell filled with liquid (see e.g. :cite:`peskin02a,crowl10a,kruger11a`). When the +thin shell filled with liquid (see e.g. :cite:`peskin02a,crowl10a,kruger12a`). When the shell is deformed by an external flow, it responds with elastic restoring forces which are transmitted into the fluid. In the present case, the inner and outer liquid are of the same type and are simulated using @@ -281,7 +281,7 @@ lattice-Boltzmann. Numerically, the shell is discretized by a set of marker points connected by triangles. The marker points are advected with *exactly* the local fluid velocity, i.e., they do not possess a mass nor a -friction coefficient (this is different from the Object-in-Fluid method +friction coefficient (this is different from the :ref:`Object-in-Fluid` method below). We implement these marker points as virtual tracer particles which are not integrated using the usual velocity-Verlet scheme, but instead are propagated using a simple Euler algorithm with diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 57644a6d958..83b894d1b70 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -8,10 +8,10 @@ This chapter will describe how to get, compile and run the software. |es| releases are available as source code packages from the homepage [1]_. This is where new users should get the code. The code within release packages is tested and known to run on a number of platforms. -Alternatively, people that want to use the newest features of |es| or that -want to start contributing to the software can instead obtain the +Alternatively, people who want to use the newest features of |es| or +start contributing to the software can instead obtain the current development code via the version control system software [2]_ -from |es|'s project page at Github [3]_. This code might be not as well +from |es|'s project page at GitHub [3]_. This code might be not as well tested and documented as the release code; it is recommended to use this code only if you have already gained some experience in using |es|. @@ -93,7 +93,7 @@ To compile |es| on Ubuntu 22.04 LTS, install the following dependencies: .. code-block:: bash sudo apt install build-essential cmake cython3 python3-pip python3-numpy \ - libboost-all-dev openmpi-common fftw3-dev libhdf5-dev libhdf5-openmpi-dev \ + libboost-all-dev openmpi-common fftw3-dev libfftw3-mpi-dev libhdf5-dev libhdf5-openmpi-dev \ python3-scipy python3-opengl libgsl-dev freeglut3 Optionally the ccmake utility can be installed for easier configuration: @@ -221,7 +221,7 @@ To use Jupyter Notebook, install the following packages: .. code-block:: bash - pip3 install --user nbformat notebook 'jupyter_contrib_nbextensions==0.5.1' + pip3 install --user 'nbformat==5.1.3' 'nbconvert==6.4.5' 'notebook==6.4.8' 'jupyter_contrib_nbextensions==0.5.1' jupyter contrib nbextension install --user jupyter nbextension enable rubberband/main jupyter nbextension enable exercise2/main diff --git a/doc/sphinx/lb.rst b/doc/sphinx/lb.rst index 30828cd6055..7ceecf4e1a1 100644 --- a/doc/sphinx/lb.rst +++ b/doc/sphinx/lb.rst @@ -473,6 +473,13 @@ This allows the user to quickly set up a system with boundary conditions that simultaneously act on the fluid and particles. For a complete description of all available shapes, refer to :mod:`espressomd.shapes`. +When using shapes, keep in mind the lattice origin is offset by half a grid +size from the box origin. For illustration purposes, assuming ``agrid=1``, +setting a wall constraint with ``dist=1`` and a normal vector pointing along +the x-axis will set all LB nodes in the left side of the box as boundary +nodes with thickness 1. The same outcome is obtained with ``dist=1.49``, +but with ``dist=1.51`` the thickness will be 2. + .. _Prototyping new LB methods: Prototyping new LB methods 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/Readme.md b/doc/tutorials/Readme.md index e7499c39b4d..869ded3b331 100644 --- a/doc/tutorials/Readme.md +++ b/doc/tutorials/Readme.md @@ -59,6 +59,10 @@ physical systems. * **Electrokinetics** Modelling electrokinetics together with hydrodynamic interactions. [Guide](electrokinetics/electrokinetics.ipynb) +* **Electrodes** + Modelling electrodes and measuring differential capacitance with the ELC method. + [Part 1](electrodes/electrodes_part1.ipynb) | + Part 2 (work in progress) * **Constant pH method** Modelling an acid dissociation curve using the constant pH method. [Guide](constant_pH/constant_pH.ipynb) diff --git a/doc/tutorials/electrodes/CMakeLists.txt b/doc/tutorials/electrodes/CMakeLists.txt new file mode 100644 index 00000000000..8343dd1fbd4 --- /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) +# TODO: fix time step issues (#4798) before adding HTML_RUN back +nb_export(TARGET tutorial_electrodes SUFFIX "2" FILE "electrodes_part2.ipynb") diff --git a/doc/tutorials/electrodes/NotesForTutor.md b/doc/tutorials/electrodes/NotesForTutor.md new file mode 100644 index 00000000000..47497b59da4 --- /dev/null +++ b/doc/tutorials/electrodes/NotesForTutor.md @@ -0,0 +1,32 @@ +# Simulations of electrodes in ESPResSO + +## Physics learning goals + +### Part 1 + +* Give a short recap about image charges, dielectric media, ... +* Nano-confinement can exhibit a broad variety of interesting effects that can + be studied with computer simulations! +* Electrostatics in nano-confinement: concept of a Green's function +* Discrete image charges: ICC\* + +### Part 2 + +* Nano-confined charged liquids as super-capacitors +* Advanced Poisson-Boltzmann theory: Gouy-Chapman, Graham equation +* Limits of PB: finite ion size, correlations, ... +* Coarse-grained view: Surface charge density via ELC-IC +* How to apply a potential difference in the simulation. + +After the tutorial, students should be able to: + +* Explain how ESPResSo implements 2D periodic electrostatics. +* What are the limitations of the mean-field PB description. +* How to evaluate the differential capacitance from simulations. +* The basic idea of super-ionic states. + +## ESPResSo learning goals + +* Setting up ELC, ICC and ELC-IC simulations +* Why is choosing the ELC-gap a crucial parameter? +* Setting up observables and accumulators for the density profiles. diff --git a/doc/tutorials/electrodes/electrodes_part1.ipynb b/doc/tutorials/electrodes/electrodes_part1.ipynb new file mode 100644 index 00000000000..acda4631df2 --- /dev/null +++ b/doc/tutorials/electrodes/electrodes_part1.ipynb @@ -0,0 +1,572 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Basic simulation of electrodes in ESPResSo part I: ion-pair in a narrow metallic slit-like confinement using ICC$^\\star$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Prerequisites\n", + "\n", + "To work with this tutorial, you should be familiar with the following topics:\n", + "\n", + "- Setting up and running simulations in ESPResSo - creating particles,\n", + " incorporating interactions.\n", + " If you are unfamiliar with this, you can go through the tutorial\n", + " in the `lennard_jones` folder.\n", + "- Basic knowledge of classical electrostatics:\n", + " dipoles, surface and image charges\n", + "- Reduced units, as described in the **ESPResSo** [user guide](https://espressomd.github.io/doc/introduction.html#on-units)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "This tutorial introduces some basic concepts for simulating charges close to an\n", + "electrode interface using **ESPResSo**.\n", + "In this first part, we focus on the interaction of a single ion pair confined in\n", + "a narrow metallic slit pore using the ICC$^\\star$-algorithm\n", + "[1] for the computation of the surface polarization.\n", + "Here, we verify the strong deviation from a Coulomb-like interaction:\n", + "In metallic confinement, the ion pair interaction energy is screened\n", + "exponentially due to the presence of induced charges on the slit walls.\n", + "[2] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Theoretical Background \n", + "\n", + "The normal component of electric field across a surface dividing two dielectric\n", + "media in presence of a surface charge density $\\sigma$ is discontinuous and follows as\n", + "$(\\varepsilon_1\\vec{E}_1 - \\varepsilon_2\\vec{E}_2).\\hat{n}=-\\sigma(\\vec{r})$.\n", + "This expression describes the jump in the electric field across the material\n", + "interface going from a dielectric medium $\\varepsilon_1$ to another one,\n", + "$\\varepsilon_2$.\n", + "\n", + "While in the case of non-polarizable materials ($\\varepsilon_1 = \\varepsilon_2 = 1$),\n", + "this jump is only related to surface charges and the\n", + "potential is continuous across the interface, for polarizable materials, the\n", + "polarization field $\\vec{P}$ will also give a contribution. \n", + "In order to solve for the electric field in presence of a jump of the dielectric constant\n", + "across an interface, one must know the electric fields on both sides. \n", + "\n", + "Another approach is to replace this two-domain problem by an equivalent one\n", + "without the explicit presence of the dielectric jump.\n", + "This is achieved by introducing an additional fictitious charge, i.e. an induced\n", + "charge density $\\sigma_{\\mathrm{ind}}$ on the surface. \n", + "With this well known \"method of image charges\", it is sufficient to know the\n", + "electric field on one side of the interface. \n", + "**ESPResSo** provides the \"Induced Charge Calculation with fast Coulomb Solvers\"\n", + "(ICC$^\\star$) algorithm [1] which employs a numerical scheme for solving the boundary integrals and the induced charge. \n", + "\n", + "*Note*: Apart from ICC$^\\star$ that solves for image charges spatially resolved, **ESPResSo** offers the \"Electrostatic layer\n", + "correction with image charges\" (ELC-IC) method [3], for\n", + "planar 2D+h partially periodic systems with dielectric interfaces that accounts for laterally averaged surface charge.\n", + "The tutorial on *Basic simulation of electrodes in ESPResSo part II*\n", + "addresses this in detail." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Green's function for charges in a dielectric slab\n", + "\n", + "The Green's function for two point charges in a dielectric slab can be solved\n", + "analytically [2].\n", + "In the metallic limit ($\\varepsilon_2 \\to\\infty$) the dielectric contrast is\n", + "$$ \\Delta = \\frac{\\varepsilon_1 - \\varepsilon_2} {\\varepsilon_1 + \\varepsilon_2} = -1 .$$\n", + "If the ions are placed in the center of a slab of width $w$ and a distance $r$\n", + "away from each other, the Green's function accounting for all image charges\n", + "simplifies to [4] \n", + "$$ 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\sum_{n=-\\infty}^\\infty \\frac{-1^n}{\\sqrt{x^2+n^2}} ,$$\n", + "where we have introduced the scaled separation $x = r/w$.\n", + "\n", + "For $x\\to 0$ the term for $n = 0$ dominates and one recovers\n", + "$$ \\lim_{x\\to 0} 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\frac{1}{x},$$\n", + "which is the classical Coulomb law.\n", + "Contrary, for large distances $x \\to \\infty$ one finds\n", + "$$ \\lim_{x\\to \\infty} 4 \\pi \\varepsilon_0 \\varepsilon_r w \\mathcal{G}(x) = \\sqrt{\\frac{8}{x}} e^{-\\pi x},$$\n", + "i.e. the interaction decays exponentially.\n", + "Such screened electrostatic interactions might explain unusual features\n", + "concerning the nano-confined ionic liquids employed for supercapacitors referred\n", + "to 'super-ionic states'.\n", + "\n", + "We will explore this interaction numerically using ICC$^\\star$ in the following." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2D+h periodic systems, dielectric interfaces and Induced Charge Computation with ICC$^\\star$\n", + "\n", + "Partially periodic ionic systems with dielectric interfaces are very often\n", + "encountered in the study of energy materials or bio-macromolecular and membrane\n", + "studies. \n", + "These systems usually exhibit a confinement along one ($z$) direction, where the\n", + "confining boundary or interface imposes a dielectric discontinuity, while the\n", + "other $x$-$y$ directions are treated periodic. \n", + "To investigate such a partially periodic system, we combine the efficient scaling behavior\n", + "of three-dimensional mesh-based solvers (typically\n", + "$\\mathcal{O}(N \\log N)$ for P3M) with the Electrostatic Layer Correction (ELC)\n", + "[3].\n", + "The latter corrects for the contributions from the periodic images in the\n", + "constrained direction and its numerical cost grows linear with the number of\n", + "point charges $N$, hence the performance overall depends on the underlying 3D\n", + "Coulomb solver.\n", + "The method relies on an empty vacuum slab in the simulation box in the\n", + "$z$-direction perpendicular to slab.\n", + "While in theory this can become relatively small (typically 10% of the box\n", + "length), its choice in practice strongly affects the performance due to the\n", + "tuning of the P3M parameters to obtain the desired accuracy.\n", + "\n", + "We furthermore employ ICC$^\\star$ to solve the Poisson equation for an\n", + "inhomogeneous dieletric:\n", + "$$ \\nabla (\\varepsilon \\nabla \\phi)=-4\\pi \\rho$$\n", + "\n", + "The image charge formalism can be derived as follows:\n", + "- Integrate the latter expression at the boundary over an infinitesimally wide\n", + "pillbox, which will give the induced surface charge in this infinitesimal\n", + "segment as (Gauss law):\n", + "$$q_{\\mathrm{ind}} = \\frac{1}{4\\pi} \\oint\\, dA\\, \\cdot \\varepsilon\\nabla \\phi = \\frac{A}{4\\pi}(\\varepsilon_1\\vec{E}_1 \\cdot \\hat{n}-\\varepsilon_2\\vec{E}_2 \\cdot\\hat{n})$$\n", + "- The electric field in region 1 at the closest proximity of the interface, $\\vec{E}_{1}$,\n", + "can be written as a sum of electric field contributions from the surface charge\n", + "$\\sigma$ and the external electric field $\\vec{E}$:\n", + "$$ \\vec{E}_{1} =\\vec{E} + 2\\pi/\\varepsilon_1\\sigma\\hat{n} $$\n", + "- Combining this with the previous expression, the induced charge can be written in terms of the dielectric mismatch $\\Delta$ and the electric field as:\n", + "$$\\sigma = \\frac{\\varepsilon_1}{2\\pi} \\frac{\\varepsilon_1-\\varepsilon_2}{\\varepsilon_1+\\varepsilon_2}\\vec{E} \\cdot \\hat{n} =: \\frac{\\varepsilon_1}{2\\pi} \\Delta \\, \\vec{E} \\cdot \\hat{n}$$\n", + "\n", + "The basic idea of the ICC$^\\star$ formalism now is to employ a discretization\n", + "of the surface by means of spatially fixed ICC particles.\n", + "The charge of each ICC particle is not fixed but rather iterated using the\n", + "expressions for $\\vec{E}_{1}$ and $\\sigma$ above until a self-consistent\n", + "solution is found." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. System setup \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We first import all ESPResSo features and external modules." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import tqdm\n", + "import numpy as np\n", + "\n", + "import espressomd\n", + "import espressomd.electrostatics\n", + "import espressomd.electrostatic_extensions\n", + "\n", + "espressomd.assert_features(['ELECTROSTATICS'])\n", + "plt.rcParams.update({'font.size': 18})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We need to define the system dimensions and some physical parameters related to\n", + "length, time and energy scales of our system.\n", + "All physical parameters are defined in reduced units of length ($\\sigma=1$;\n", + "particle size), mass ($m=1$; particle mass), arbitrary time (we do not consider particle dynamics) and\n", + "elementary charge ($e=1$).\n", + "\n", + "Another important length scale is the Bjerrum Length, which is the length at\n", + "which the electrostatic energy between two elementary charges is comparable to\n", + "the thermal energy $k_\\mathrm{B}T$.\n", + "It is defined as\n", + "$$\\ell_\\mathrm{B}=\\frac{1}{4\\pi\\varepsilon_0\\varepsilon_r k_\\mathrm{B}T}.$$\n", + "\n", + "In our case, if we choose the ion size ($\\sigma$) in simulation units to represent a\n", + "typical value for mono-atomic salt, 0.3 nm in real units, then the\n", + "Bjerrum length of water at room temperature, $\\ell_\\mathrm{B}=0.71 \\,\\mathrm{nm}$ is\n", + "$\\ell_\\mathrm{B}\\sim 2$ in simulations units." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Box dimensions\n", + "# To construct a narrow slit Lz << (Lx , Ly)\n", + "box_l_x = 100.\n", + "box_l_y = 100.\n", + "box_l_z = 5.\n", + "\n", + "# Additional space for ELC\n", + "ELC_GAP = 6*box_l_z\n", + "\n", + "system = espressomd.System(box_l=[box_l_x, box_l_y, box_l_z + ELC_GAP])\n", + "\n", + "system.time_step = 0.01\n", + "system.cell_system.skin = 0.4\n", + "\n", + "# Elementary charge \n", + "q = 1.0 \n", + "\n", + "# Interaction parameters for P3M with ELC\n", + "BJERRUM_LENGTH = 2.0 # Electrostatic prefactor passed to P3M ; prefactor=lB KBT/e2 \n", + "ACCURACY = 1e-7 # P3M force accuracy \n", + "MAX_PW_ERROR = 1e-7 # maximum pairwise error in ELC\n", + "ICC_EPSILON_DOMAIN = 1. # epsilon inside the slit\n", + "ICC_EPSILON_WALLS = 1e5 # epsilon outside the slit. Very large to mimic metal\n", + "ICC_CONVERGENCE = 1e-3 # ICC numeric/performance parameters\n", + "ICC_RELAXATION = 0.95\n", + "ICC_MAX_ITERATIONS = 1000\n", + "\n", + "# Lennard-Jones parameters\n", + "LJ_SIGMA = 1.0\n", + "LJ_EPSILON = 1.0 \n", + "\n", + "# Particle parameters\n", + "TYPES = {\"Cation\": 0, \"Anion\": 1 ,\"Electrodes\": 2}\n", + "charges = {\"Cation\": q, \"Anion\": -q }\n", + "\n", + "# Test particles to measure forces\n", + "p1 = system.part.add(pos=[box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges[\"Cation\"], fix=3*[True])\n", + "p2 = system.part.add(pos=[3.0*box_l_x/4.0, box_l_y/2.0, box_l_z/2.0], q=charges[\"Anion\"], fix=3*[True])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup of electrostatic interactions\n", + "First, we define our 3D electrostatics solver (P3M)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p3m = espressomd.electrostatics.P3M(\n", + " prefactor=BJERRUM_LENGTH,\n", + " accuracy=ACCURACY,\n", + " check_neutrality=False,\n", + " mesh=[100, 100, 150],\n", + " cao=5\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task\n", + "\n", + "* Set up [ELC](https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatics.ELC) with ``p3m`` as its actor." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=ELC_GAP, maxPWerror=MAX_PW_ERROR)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we set up the ICC particles on both electrodes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ICC_PARTCL_NUMBER_DENSITY = 1.\n", + "icc_partcls_bottom = []\n", + "icc_partcls_top = []" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### TASK\n", + "\n", + "* Using the (area) density of ICC particles defined in the cell above, calculate the x/y positions of the particles for a uniform, quadratic grid. \n", + "* Add fixed particles on the electrodes. Make sure to use the correct ``type``. Give the top (bottom) plate a total charge of $+1$ ($-1$). \n", + "* Store the created particles in lists ``icc_partcls_bottom``, ``icc_partcls_top``." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "line_density = np.sqrt(ICC_PARTCL_NUMBER_DENSITY)\n", + "xs = np.linspace(0, system.box_l[0], num=int(round(system.box_l[0] * line_density)), endpoint=False)\n", + "ys = np.linspace(0, system.box_l[1], num=int(round(system.box_l[1] * line_density)), endpoint=False)\n", + "n_partcls_each_electrode = len(xs) * len(ys)\n", + "\n", + "# Bottom electrode\n", + "for x in xs:\n", + " for y in ys:\n", + " icc_partcls_bottom.append(system.part.add(pos=[x, y, 0.], q=-1. / n_partcls_each_electrode,\n", + " type=TYPES[\"Electrodes\"], fix=3*[True]))\n", + "# Top electrode\n", + "for x in xs:\n", + " for y in ys:\n", + " icc_partcls_top.append(system.part.add(pos=[x, y, box_l_z], q=1. / n_partcls_each_electrode,\n", + " type=TYPES[\"Electrodes\"], fix=3*[True]))\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task\n", + "\n", + "* Set ``elc`` as ``system.electrostatics.solver``\n", + "* Create an [ICC object]((https://espressomd.github.io/doc/espressomd.html#espressomd.electrostatic_extensions.ICC) and set it as ``system.electrostatics.extension``\n", + "\n", + "### Hints\n", + "\n", + "* ICC variables are defined in the second code cell from the top.\n", + "* Make sure to not mark our test particles ``p1`` and ``p2`` (with ids 0 and 1) as ICC particles." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "system.electrostatics.solver = elc\n", + "\n", + "n_icc_partcls = len(icc_partcls_top) + len(icc_partcls_bottom)\n", + "\n", + "# Common area, sigma and metallic epsilon\n", + "icc_areas = 1. / ICC_PARTCL_NUMBER_DENSITY * np.ones(n_icc_partcls)\n", + "icc_sigmas = np.zeros(n_icc_partcls)\n", + "icc_epsilons = ICC_EPSILON_WALLS * np.ones(n_icc_partcls)\n", + "icc_normals = np.repeat([[0, 0, 1], [0, 0, -1]], repeats=n_icc_partcls//2, axis=0)\n", + "\n", + "icc = espressomd.electrostatic_extensions.ICC(\n", + " first_id=min(system.part.select(type=TYPES[\"Electrodes\"]).id),\n", + " n_icc=n_icc_partcls,\n", + " convergence=ICC_CONVERGENCE,\n", + " relaxation=ICC_RELAXATION,\n", + " ext_field=[0, 0, 0],\n", + " max_iterations=ICC_MAX_ITERATIONS,\n", + " eps_out=ICC_EPSILON_DOMAIN,\n", + " normals=icc_normals,\n", + " areas=icc_areas,\n", + " sigmas=icc_sigmas,\n", + " epsilons=icc_epsilons\n", + ")\n", + "system.electrostatics.extension = icc\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Calculation of the forces" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N_AXIAL_POINTS = 10\n", + "r = np.logspace(0., box_l_z / 4., N_AXIAL_POINTS) \n", + "elc_forces_axial = np.empty((N_AXIAL_POINTS, 2))\n", + "n_icc_per_electrode = len(icc_partcls_top)\n", + "\n", + "p1.pos = [0., box_l_y / 2., box_l_z / 2.]\n", + "\n", + "for i in tqdm.trange(N_AXIAL_POINTS):\n", + " p2.pos = [r[i], box_l_y / 2., box_l_z / 2.]\n", + "\n", + " system.integrator.run(0, recalc_forces=True)\n", + " elc_forces_axial[i, 0] = p1.f[0]\n", + " elc_forces_axial[i, 1] = p2.f[0]\n", + " \n", + " # reset ICC charges to ensure charge neutrality \n", + " for part in icc_partcls_top:\n", + " part.q = 1. / n_icc_per_electrode\n", + " for part in icc_partcls_bottom:\n", + " part.q = -1. / n_icc_per_electrode" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Analysis and Interpretation of the data\n", + "\n", + "With this we can now compare the force between the two ions to the analytical prediction.\n", + "To evaluate the infinite series we truncate at $n=1000$, which already is well converged." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def analytic_force_centered(r,w):\n", + " def summand(x, n):\n", + " return (-1)**n * x / (x**2 + n**2)**(3. / 2.)\n", + " \n", + " def do_sum(x):\n", + " limit = 1000\n", + " accumulator = 0.\n", + " for n in range(-limit + 1, limit + 1):\n", + " accumulator += summand(x, n)\n", + " return accumulator\n", + "\n", + " x = r / w\n", + " prefactor = BJERRUM_LENGTH / w**2\n", + " F = do_sum(x) * prefactor\n", + " return F\n", + "\n", + "def coulomb_force(x):\n", + " prefactor = BJERRUM_LENGTH\n", + " E = prefactor / x**2\n", + " return E\n", + "\n", + "def exponential_force(r,w):\n", + " x = r / w\n", + " prefactor = BJERRUM_LENGTH\n", + " E = prefactor * np.sqrt(2.) * (1. / x)**(3. / 2.) * np.exp(-np.pi * x)\n", + " return E" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(10, 6))\n", + "\n", + "x = np.logspace(-0.25, 1.45, 100)\n", + "plt.plot(x / BJERRUM_LENGTH, analytic_force_centered(x, box_l_z), color='b', label=\"analytic\", marker='')\n", + "plt.plot(x / BJERRUM_LENGTH, coulomb_force(x), color='g', ls='--', label='Coulomb')\n", + "plt.plot(x / BJERRUM_LENGTH, exponential_force(x, box_l_z), color='r', ls='--', label='Exponential')\n", + "\n", + "plt.plot(r / BJERRUM_LENGTH, -elc_forces_axial[:,1], color='r', label=\"sim (p2)\", marker='o', ls='')\n", + "plt.plot(r / BJERRUM_LENGTH, +elc_forces_axial[:,0], color='k', label=\"sim (p1)\", marker='x', ls='')\n", + "\n", + "plt.xlabel(r'$r \\, [\\ell_\\mathrm{B}]$')\n", + "plt.ylabel(r'$F \\, [k_\\mathrm{B}T / \\sigma$]')\n", + "plt.loglog()\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] Tyagi, S.; Süzen, M.; Sega, M.; Barbosa, M.; Kantorovich, S. S.; Holm, C. An Iterative, Fast, Linear-Scaling Method for Computing Induced Charges on Arbitrary Dielectric Boundaries. J. Chem. Phys. 2010, 132 (15), 154112. https://doi.org/10.1063/1.3376011.\n", + " \n", + "[2] Kondrat, S.; Feng, G.; Bresme, F.; Urbakh, M.; Kornyshev, A. A. Theory and Simulations of Ionic Liquids in Nanoconfinement. Chem. Rev. 2023, 123 (10), 6668–6715. https://doi.org/10.1021/acs.chemrev.2c00728.\n", + "\n", + "[3] Tyagi, S.; Arnold, A.; Holm, C. Electrostatic Layer Correction with Image Charges: A Linear Scaling Method to Treat Slab 2D+h Systems with Dielectric Interfaces. J. Chem. Phys. 2008, 129 (20), 204102. https://doi.org/10.1063/1.3021064.\n", + "\n", + "[4] Loche, P.; Ayaz, C.; Wolde-Kidan, A.; Schlaich, A.; Netz, R. R. Universal and Nonuniversal Aspects of Electrostatics in Aqueous Nanoconfinement. J. Phys. Chem. B 2020, 124 (21), 4365–4371. https://doi.org/10.1021/acs.jpcb.0c01967." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb new file mode 100644 index 00000000000..35f7723156a --- /dev/null +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -0,0 +1,1088 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Basic simulation of electrodes in ESPResSo part II: Electrolyte capacitor and Poisson–Boltzmann theory" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Prerequisites\n", + "\n", + "To work with this tutorial, you should be familiar with the following topics:\n", + "\n", + "1. Setting up and running simulations in ESPResSo - creating particles,\n", + " incorporating interactions.\n", + " If you are unfamiliar with this, you can go through the tutorial\n", + " in the `lennard_jones` folder.\n", + "2. Basic knowledge of classical electrostatics:\n", + " dipoles, surface and image charges.\n", + "3. Reduced units and how to relate them to physical quantities, see the ESPResSo\n", + " [user guide](https://espressomd.github.io/doc/introduction.html#on-units)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "Ionic liquid (IL) based capacitors have long been established as promising\n", + "candidates in the area of efficient energy storage devices due to their\n", + "extraordinary capacitance, which is why they are also termed super-capacitors.\n", + "Typically, in these setups a fluid consisting of mobile charge carriers is\n", + "placed between two electrodes and thus gets polarized upon application of an\n", + "external potential formic an electric double layer at the interfaces\n", + "[1].\n", + "Electric double-layer capacitors (EDLCs) can be constructed from electrodes of\n", + "various geometries and materials where energy is stored by potential driven\n", + "adsorption of counterions on the surface of the electrodes and forming the\n", + "double layer. \n", + "Thus, conducting, high surface area electrode materials can further maximize the\n", + "energy per volume.\n", + "\n", + "In this tutorial we are going to investigate the ionic layer formation between\n", + "two conducting dielectric walls in presence of an applied voltage using **ESPResSo**'s\n", + "\"Electrostatic layer correction with image charges\" (ELC-IC) method\n", + "[2].\n", + "We employ a primitive model of an aqueous salt solution, where the solvent is\n", + "treated implicitly as a homogeneous dielectric background.\n", + "Thus, within the limits of a mean-field treatment and for not too small slit\n", + "widths, we can compare our findings with the analytical Gouy–Chapmann\n", + "solution for a single charged plane since additivity is assumed to hold if the\n", + "potential of both walls is screened sufficiently.\n", + "While this mean-field formalism properly describes the behavior of Coulomb\n", + "fluids composed of monovalent ions at low concentrations in the vicinity of\n", + "weakly charged interfaces, for strongly charged systems, where correlation and\n", + "finite size effects begin to dominate, Poisson–Boltzmann theory falls\n", + "inadequate.\n", + "Our goal in this tutorial is to demonstrate how coarse grained implicit solvent\n", + "simulations can corroborate on some of the approximations and hint on\n", + "extensions/deviations.\n", + "\n", + "The inclusion of dielectric inhomogeneities appearing at the conducting\n", + "interfaces demands to take into account image effects that involve the full\n", + "solution of the Poisson equation.\n", + "This is dealt with in a computationally cost-effective way using the ELC-IC method to\n", + "treat the image charge effects in the presence of 2D dielectric bounding\n", + "interfaces. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Theoretical Background \n", + "\n", + "### Poisson–Boltzmann Theory\n", + "\n", + "Charged surfaces in contact with a liquid containing free charges (ions) attract\n", + "oppositely charged ions that form a diffusive electric double layer. \n", + "The competition between electrostatic interactions and entropy of ions in\n", + "solution determines the exact distribution of mobile ions close to charged\n", + "membranes.\n", + "Gouy [3] and Chapman [4] derived in the\n", + "early 20th century the analytic solution for the case of a single planar wall\n", + "within the mean-field approximation of the Poisson–Boltzmann (PB) equation. \n", + "We will use it to describe our two-electrode system, which is justified if the electrodes are \n", + "so far apart that one surface does not influence the ion distribution in front of the other surface.\n", + "\n", + "In the case of a monovalent electrolyte, double integrating the PB equation and\n", + "employing the corresponding boundary conditions for the charged plane located at\n", + "$z=0$ yields the electrostatic potential:\n", + "$$\\phi(z) = -2\\ln\\left[\n", + " \\frac{1-\\tanh(\\phi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D} z}}\n", + " {1+\\tanh(\\phi_\\mathrm{s}/4)e^{-\\kappa_\\mathrm{D}z}} \\right].$$ \n", + "Here, $\\phi_\\mathrm{s}=\\phi(z=0)=$ const is the surface potential such\n", + "that $\\phi(z\\rightarrow \\infty)=0$.\n", + "$\\kappa_\\mathrm{D} = \\lambda_\\mathrm{D}^{-1}$ is the inverse Debye screening length given by\n", + "$$ \\lambda_\\mathrm{D} = \\left(\\frac{\\varepsilon \\, k_{\\mathrm B} T}{\\sum_{j = 1}^N n_j^0 \\, q_j^2}\\right)^{1/2}, $$\n", + "where $n_j^{(0)}$ and $q_j$ are the equilibrium number density and charge of the\n", + "$j$-th ion species.\n", + "For the monovalent salt this can conveniently be expressed in terms of the\n", + "Bjerrum length $\\ell_\\mathrm{B}$ and the equilibrium salt concentration\n", + "$\\rho^{(0)}=\\sum_j \\rho_j^{(0)}$,\n", + "$$ \\lambda_{\\mathrm D} = \\left(4 \\pi \\, \\ell_\\mathrm{B} \\, \\rho^{(0)}/e\\right)^{-1/2} . $$\n", + "\n", + "The cationic and anionic density profiles then follow from the Boltzmann\n", + "distribution as:\n", + "$$n_{\\pm}(z)=n_\\pm^{(0)}\\left(\\frac{1\\pm\\gamma e^{-\\kappa_\\mathrm{D}z}}\n", + " {1\\mp\\gamma e^{-\\kappa_\\mathrm{D}z}} \\right)^2$$\n", + "Here, $\\gamma$ is associated with the surface potential as\n", + "$\\phi_\\mathrm{s}=-4\\tanh^{-1}(\\gamma)$.\n", + "At large z, where the potential decays to zero, the ionic profiles tend to their\n", + "bulk (reservoir) densities, $n_\\pm(z\\to\\infty) = n_\\pm^{(0)}$\n", + "\n", + "The relation between the surface potential and the surface charge induced on the\n", + "electrode is given by the Grahame Equation [5] :\n", + "$$ \\sigma = \\sinh(\\phi_\\mathrm{s}/2) \\sqrt{\\frac{2 n_\\mathrm{b}}{\\pi \\ell_\\mathrm{B}}} $$\n", + "The latter expression thus yields the differential capacitance\n", + "$C=\\displaystyle\\frac{\\mathrm{d}\\sigma}{\\mathrm{d}\\phi_\\mathrm{s}}$ within the mean-field\n", + "solution for non-overlapping double layers." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ELC-IC for 2D+h periodic systems with dielectric interfaces\n", + "\n", + "In this tutorial we employ a parallel plate capacitor setup with constant\n", + "potential boundary conditions which needs to be treated appropriately by the\n", + "electrostatics solver.\n", + "To simulate a two-dimensional partially periodic system, we combine the efficient\n", + "scaling behavior of three-dimensional mesh-based solvers (typically\n", + "$\\mathcal{O}(N \\log N)$ for P3M) with the Electrostatic Layer Correction (ELC)\n", + "[1].\n", + "The latter removes the contributions from the periodic images in the\n", + "non-periodic direction and its numerical cost grows linear with the number of\n", + "point charges $N$, hence the performance overall depends on the underlying 3D\n", + "Coulomb solver.\n", + "Furthermore, ELC can be extended straightforwardly to metallic boundary\n", + "conditions (or any other dielectric contrast) by using the method of image charges,\n", + "which is referred to as the “Electrostatic Layer Correction with Image Charges”\n", + "(ELC-IC) approach used in this tutorial.\n", + "\n", + "A voltage difference can be applied between the electrodes by following\n", + "considerations:\n", + "The total potential drop $\\Delta \\phi$ across the simulated system is readily\n", + "obtained from the ion distribution and integrating twice over the one-dimensional Poisson equation:\n", + "$$-\\varepsilon_{0}\\varepsilon_{r}\\phi_\\mathrm{ind}(z)=\\iint_{0}^{z}\\rho(z^{\\prime})(z-z^{\\prime})dz^{\\prime}$$\n", + "Here, the subscript 'ind' indicates that this is the potential due to the\n", + "induced inhomogeneous charge distribution.\n", + "In order to set up a constant potential difference $\\Delta \\phi$, a homogeneous\n", + "electric field is superimposed such that\n", + "$$ \\Delta \\phi = \\Delta \\phi_\\mathrm{ind} + \\Delta \\phi_\\mathrm{bat},$$\n", + "where $\\Delta \\phi_\\mathrm{bat}$ corresponds to the potential of a (virtually)\n", + "applied battery.\n", + "In practice, the linear electric field in $E_z^\\mathrm{(bat)}=-\\phi_\\mathrm{bat}/L_z$\n", + "in the $z$-direction normal to the surface that one needs to apply can be\n", + "calculated straightforwardly, as the corresponding contribution from the induced\n", + "charges is known:\n", + "$$ E_z^\\mathrm{(ind)} = \\frac{1}{\\varepsilon_0 \\varepsilon_r L^2 L_z} P_z$$\n", + "Here, $L$ denotes the lateral system size, $L_z$ the distance between the plates\n", + "and $P_z = \\sum_i q_i z_i$ is the total dipole moment in $z$-direction due to\n", + "the charges $q_i$.\n", + "Then, to maintain $\\Delta \\phi$, a force $F_z^\\mathrm{bat} = q E_z^\\mathrm{(bat)}$\n", + "is applied on all charges.\n", + "Since ELC already calculates $P_z$, the constant potential correction requires no\n", + "additional computational effort.\n", + "\n", + "*Note*: Apart from ELC-IC, **ESPResSo** also provides the ICC$^\\star$ method\n", + "[2] which employs an iterative numerical scheme with\n", + "discretized surface particles to solve the boundary integrals at the dieletcric\n", + "interface.\n", + "The tutorial on *Basic simulation of electrodes in ESPResSo part I* addresses this\n", + "in detail." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. System setup \n", + "\n", + "First we import all ESPResSo features and external modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import scipy.constants as constants\n", + "import scipy.stats\n", + "import tqdm\n", + "\n", + "import espressomd\n", + "import espressomd.electrostatics\n", + "import espressomd.electrostatic_extensions\n", + "import espressomd.observables\n", + "import espressomd.accumulators\n", + "import espressomd.shapes\n", + "\n", + "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", + "rng = np.random.default_rng(42)\n", + "plt.rcParams.update({'font.size': 18})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We need to define system dimensions and some physical parameters related to\n", + "length, time and energy scales of our system.\n", + "As discussed in previous tutorials, all physical parameters are defined in terms\n", + "of a length $\\sigma$, mass $m$ and time $t$ and unit of charge $q$.\n", + "Since we are not explicitly interested in the dynamics of the system, we set the\n", + "mass to $m=1$ (particle mass).\n", + "For convenience, we choose the elementary charge as fundamental unit ($q=1e$)\n", + "and $\\sigma = 1 \\,\\mathrm{nm}$.\n", + "\n", + "With this we can now define the fundamental parameters of our system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# water at room temperature\n", + "EPSILON_R = 78.4 # Relative dielectric constant of water\n", + "TEMPERATURE = 300.0 # Temperature in Kelvin\n", + "BJERRUM_LENGTH = constants.elementary_charge**2 / constants.nano / \\\n", + " (4 * np.pi * constants.epsilon_0 * EPSILON_R * constants.Boltzmann * TEMPERATURE)\n", + "# BERRUM_LENGTH of water at room temperature is 0.71 nm; electrostatic prefactor passed to P3M KBT/e2 \n", + "\n", + "# Lennard-Jones parameters\n", + "LJ_SIGMA = 0.3 # Particle size in nanometers\n", + "LJ_EPSILON = 1.0\n", + "\n", + "CONCENTRATION = 1e-2 # desired concentration in mol/l\n", + "DISTANCE = 10 # 10 Debye lengths\n", + "N_IONPAIRS = 500\n", + "\n", + "POTENTIAL_DIFF = 5.0\n", + "\n", + "# Elementary charge \n", + "q = 1.\n", + "types = {\"Cation\": 0, \"Anion\": 1, \"Electrodes\": 2}\n", + "charges = {\"Cation\": q, \"Anion\": -q}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.1 Setting up the box dimensions and create system\n", + "\n", + "We want to make use of the optimal performance of **ESPResSo** in this tutorial,\n", + "which is roughly at 1000 particles/core.\n", + "Thus, we fixed above the number of ion pairs to `N_IONPAIRS = 500`.\n", + "\n", + "To be able to employ the analytical solution for the single plate also for the\n", + "double layer capacitor setup, the two electrodes need to be sufficiently far\n", + "away such that the additivity of the two surface potentials holds. In practice,\n", + "a separation of $d=10\\lambda_\\mathrm{D}$ is a good choice, represented by \n", + "`DISTANCE = 10`.\n", + "\n", + "Our choice of $c=10\\,\\mathrm{mmol}$ is a compromise between a sufficiently low\n", + "concentration for the PB theory to hold and not too large distances $d$ such\n", + "that the equilibration/diffusion of the ions is sufficiently fast\n", + "(`CONCENTRATION = 1e-2`).\n", + "\n", + "Note that in order to obtain results that we can interpret easily, we explicitly\n", + "set a unit system using nanometers as length-scale above.\n", + "The corresponding ion size of about 0.3 nm is a\n", + "typical value for a simple salt; this, however, is in sharp contrast to the\n", + "mean-field assumption of point-like ions.\n", + "The latter are not easily studied within Molecular Dynamics simulations due to\n", + "the required small time steps and are better suited for Monte Carlo type\n", + "simulations.\n", + "We instead focus here on analyzing deviations from PB theory due to the finite\n", + "ion size." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task\n", + "\n", + "* write a function \n", + "`get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS)`\n", + "that returns the lateral and normal box lengths `box_l_xy` and `box_l_z` (in\n", + "nanometers) for the given parameters.\n", + "\n", + "**Hint:** To account for the finite ion size and the wall interaction it is\n", + "useful to define the effective separation $d^\\prime = d-2\\sigma$, such that the\n", + "concentration is $\\rho = N/(A \\cdot d^\\prime)$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):\n", + " \"\"\"\n", + " For a given number of particles, determine the lateral area of the box\n", + " to match the desired concentration.\n", + " \"\"\"\n", + "\n", + " # concentration is in mol/l, convert to 1/sigma**3\n", + " rho = concentration * (constants.Avogadro / constants.liter) * constants.nano**3\n", + " debye_length = (4. * np.pi * BJERRUM_LENGTH * rho * 2.)**(-1. / 2.) # desired Debye length in nm\n", + " l_z = distance * debye_length\n", + " \n", + " box_volume = n_ionpairs / rho\n", + " area = box_volume / (l_z - 2. * LJ_SIGMA) # account for finite ion size in density calculation\n", + " l_xy = np.sqrt(area)\n", + "\n", + " return l_xy, l_z\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "box_l_xy, box_l_z = get_box_dimension(CONCENTRATION, DISTANCE, N_IONPAIRS)\n", + "\n", + "# useful quantities for the following calculations\n", + "DEBYE_LENGTH = box_l_z / DISTANCE # in units of nm\n", + "rho = N_IONPAIRS / (box_l_xy * box_l_xy * box_l_z) # in units of 1/nm^3" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now can create the **ESPResSo** system.\n", + "\n", + "Note that for ELC to work properly, we need to add a gap of `ELC_GAP` in the\n", + "non-periodic direction.\n", + "The precise value highly affects the performance due to the tuning of the P3M\n", + "electrostatic solver.\n", + "For $d=10\\lambda$ a gap value of $6 L_z$ is a good choice.\n", + "\n", + "We also set the time-step $dt = 0.01 \\tau$, which is limited by the choice of\n", + "$\\sigma$ and $\\tau$ in the repulsive WCA interaction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ELC_GAP = 6. * box_l_z\n", + "system = espressomd.System(box_l=[box_l_xy, box_l_xy, box_l_z + ELC_GAP])\n", + "system.time_step = 0.01" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Set up the double-layer capacitor\n", + "\n", + "We now set up an electrolyte solution made of monovalent cations and anions\n", + "between two metallic electrodes at constant potential. \n", + "\n", + "#### 1.2.1 Electrode walls " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "### Task\n", + "* add two wall constraints at $z=0$ and $z=L_z$ to stop particles from\n", + "crossing the boundaries and model the electrodes.\n", + "Refer to \n", + "[espressomd.constraints.ShapeBasedConstraint](https://espressomd.github.io/doc/espressomd.html#espressomd.constraints.ShapeBasedConstraint)\n", + "and its\n", + "[wall constraint](https://espressomd.github.io/doc/constraints.html?highlight=constraint#wall)\n", + "in the documentation to set up constraints and the `types` dictionary for the\n", + "particle type." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "# Bottom wall, normal pointing in the +z direction \n", + "floor = espressomd.shapes.Wall(normal=[0, 0, 1])\n", + "c1 = system.constraints.add(\n", + " particle_type=types[\"Electrodes\"], penetrable=False, shape=floor)\n", + "\n", + "# Top wall, normal pointing in the -z direction\n", + "ceiling = espressomd.shapes.Wall(normal=[0, 0, -1], dist=-box_l_z) \n", + "c2 = system.constraints.add(\n", + " particle_type=types[\"Electrodes\"], penetrable=False, shape=ceiling)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "#### 1.2.2 Add particles for the ions\n", + "\n", + "### Task\n", + "\n", + "* place ion pairs at random positions between the electrodes.\n", + "\n", + "Note, that unfavorable overlap can be avoided by placing the particles in the\n", + "interval $[\\sigma, d-\\sigma]$ in the $z$-direction only." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "offset = LJ_SIGMA # avoid unfavorable overlap at close distance to the walls\n", + "init_part_btw_z1 = offset \n", + "init_part_btw_z2 = box_l_z - offset\n", + "ion_pos = np.empty((3), dtype=float)\n", + "\n", + "for i in range (N_IONPAIRS):\n", + " ion_pos[0] = rng.random(1) * system.box_l[0]\n", + " ion_pos[1] = rng.random(1) * system.box_l[1]\n", + " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", + " system.part.add(pos=ion_pos, type=types[\"Cation\"], q=charges[\"Cation\"])\n", + " \n", + "for i in range (N_IONPAIRS):\n", + " ion_pos[0] = rng.random(1) * system.box_l[0]\n", + " ion_pos[1] = rng.random(1) * system.box_l[1]\n", + " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", + " system.part.add(pos=ion_pos, type=types[\"Anion\"], q=charges[\"Anion\"])\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "#### 1.2.3 Add interactions:\n", + "\n", + "### Task\n", + "\n", + "* For excluded volume interactions, add a WCA potential. \n", + "\n", + "Refer to the documentation to set up the\n", + "[WCA interaction](https://espressomd.github.io/doc/espressomd.html#espressomd.interactions.WCAInteraction) \n", + "under [Non-bonded](https://espressomd.github.io/doc/inter_non-bonded.html)\n", + "section." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "for t1 in types.values():\n", + " for t2 in types.values():\n", + " system.non_bonded_inter[t1, t2].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "For the (2D+h) electrostatic with dielectrics we choose the ELC-IC with P3M.\n", + "\n", + "Refer the documentation to set up\n", + "[ELCIC with P3M](https://espressomd.github.io/doc/electrostatics.html#electrostatic-layer-correction-elc)\n", + "under the [electrostatics](https://espressomd.github.io/doc/electrostatics.html)\n", + "section. \n", + "\n", + "As later we will study different potential drops between the electrodes, write a\n", + "function that sets up the electrostatic solver for a given value\n", + "`POTENTIAL_DIFF.`\n", + "This function will take care of tuning the P3M and ELC parameters.\n", + "For our purposes, an accuracy of $10^{-3}$ is sufficient.\n", + "\n", + "### Task\n", + "\n", + "* Write a function `setup_electrostatic_solver(potential_diff)` that\n", + "returns the ELC instance." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def setup_electrostatic_solver(potential_diff):\n", + " delta_mid_top = -1. # (Fully metallic case both -1) \n", + " delta_mid_bot = -1.\n", + " accuracy = 1e-3\n", + " elc_accuracy = 1e-3\n", + " p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,\n", + " accuracy=accuracy,\n", + " tune=True,\n", + " verbose=False)\n", + " elc = espressomd.electrostatics.ELC(actor=p3m,\n", + " gap_size=ELC_GAP,\n", + " const_pot=True,\n", + " pot_diff=potential_diff,\n", + " maxPWerror=elc_accuracy,\n", + " delta_mid_bot=delta_mid_bot,\n", + " delta_mid_top=delta_mid_top)\n", + " return elc\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now add the solver to the system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Equilibration\n", + "\n", + "### 2.1 Steepest descent\n", + "\n", + "Before we can start the simulation, we need to remove the overlap between particles.\n", + "For this, we use the steepest descent integrator.\n", + "Afterwards, we switch to a Velocity Verlet integrator and set up a Langevin thermostat.\n", + "Note, that we only analyze static properties, thus the damping and temperature chosen\n", + "here only determine the simulation time towards the equilibrium distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "system.integrator.set_steepest_descent(f_max=10, gamma=50.0,\n", + " max_displacement=0.02)\n", + "system.integrator.run(250)\n", + "system.integrator.set_vv()\n", + "system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Equilibrate the ion distribution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Equlibration parameters\n", + "STEPS_PER_SAMPLE = 200\n", + "N_SAMPLES_EQUIL = 50\n", + "N_PART = 2 * N_IONPAIRS\n", + "\n", + "times = np.zeros(N_SAMPLES_EQUIL)\n", + "e_total = np.zeros_like(times)\n", + "e_kin = np.zeros_like(times)\n", + "\n", + "for i in tqdm.trange(N_SAMPLES_EQUIL):\n", + " times[i] = system.time\n", + " energy = system.analysis.energy()\n", + " e_total[i] = energy['total']\n", + " e_kin[i] = energy['kinetic']\n", + " system.integrator.run(STEPS_PER_SAMPLE)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the convergence of the total energy\n", + "plt.figure(figsize=(10, 6))\n", + "plt.plot(times, e_total, label='total')\n", + "plt.plot(times, e_kin, label='kinetic')\n", + "plt.xlabel('Simulation time')\n", + "plt.ylabel('Energy')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convergence after $t\\sim50$ time units." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "## 3. Calculate and analyze ion profile\n", + "\n", + "### 3.1 Set up the density accumulators\n", + "\n", + "We now need to set up an \n", + "[espressomd.observables.DensityProfile](https://espressomd.github.io/doc/espressomd.html#espressomd.observables.DensityProfile)\n", + "observable to calculate the anion and cation density profiles.\n", + "\n", + "The time average is obtained through a\n", + "[espressomd.accumulators.MeanVarianceCalculator](espressomd.accumulators.MeanVarianceCalculator).\n", + "\n", + "### Task\n", + "\n", + "* Write a function `setup_densityprofile_accumulators(bin_width)` that returns the\n", + "`bin_centers` and the accumulators for both ion species in the $z$-range $[0,d]$.\n", + "Since we are not estimating errors in this tutorial, the choice of `delta_N` is\n", + "rather arbitrary and does not affect the results. In practice, a typical value is\n", + "`delta_N=20`." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def setup_densityprofile_accumulators(bin_width):\n", + " cations = system.part.select(type=types[\"Cation\"]) \n", + " anions = system.part.select(type=types[\"Anion\"])\n", + " n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))\n", + " density_profile_cation = espressomd.observables.DensityProfile(\n", + " ids=cations.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,\n", + " max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)\n", + " density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(\n", + " obs=density_profile_cation, delta_N=20)\n", + " density_profile_anion = espressomd.observables.DensityProfile(\n", + " ids=anions.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,\n", + " max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)\n", + " density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(\n", + " obs=density_profile_anion, delta_N=20)\n", + " zs = density_profile_anion.bin_centers()[0, 0, :, 2]\n", + " return zs, density_accumulator_cation, density_accumulator_anion\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(\n", + " bin_width=DEBYE_LENGTH / 10.)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.2 Run the simulation\n", + "\n", + "Now we take some measurement sampling the density profiles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N_SAMPLES_PROD = 10\n", + "\n", + "# Add the accumulators\n", + "system.auto_update_accumulators.clear()\n", + "system.auto_update_accumulators.add(density_accumulator_cation)\n", + "system.auto_update_accumulators.add(density_accumulator_anion)\n", + " \n", + "times = []\n", + "e_total = []\n", + "for tm in tqdm.trange(N_SAMPLES_PROD):\n", + " system.integrator.run(STEPS_PER_SAMPLE)\n", + " times.append(system.time)\n", + " energy = system.analysis.energy()\n", + " e_total.append(energy['total'])\n", + "\n", + "cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n", + "anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare to analytical prediction\n", + "\n", + "Since we assume additivity, the total ion density follows from\n", + "$$ \\rho (z) = \\rho_+(z) - \\rho_+ (d-z) + \\rho_-(z) - \\rho_-(d-z) .$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gouy_chapman_potential(x, debye_length, phi_0):\n", + " kappa = 1. / debye_length\n", + " return 2. * np.log((1. + np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))) \\\n", + " / (1. - np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))))\n", + "\n", + "def gouy_chapman_density(x, c0, debye_length, phi_0):\n", + " phi = gouy_chapman_potential(x, debye_length, phi_0)\n", + " return c0 / 2. * np.exp(-phi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, (ax1, ax2, ax3) = plt.subplots(figsize=(16, 4), nrows=1, ncols=3, sharey=True)\n", + "fig.subplots_adjust(wspace=0)\n", + "\n", + "x = np.linspace(LJ_SIGMA, box_l_z-LJ_SIGMA, 100)\n", + "anion_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2.\n", + "cation_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) \\\n", + " + gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.))/2.\n", + "\n", + "ax1.set_title('Cation')\n", + "ax2.set_title('Anion')\n", + "ax3.set_title('Total')\n", + "\n", + "ax1.plot(x, cation_profile_analytic, label='analytic')\n", + "ax2.plot(x, anion_profile_analytic, label='analytic')\n", + "ax3.plot(x, cation_profile_analytic + anion_profile_analytic, label='analytic')\n", + "\n", + "ax1.plot(zs[1:-1], cation_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n", + "ax2.plot(zs[1:-1], anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n", + "ax3.plot(zs[1:-1], cation_profile_mean[1:-1] + anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')\n", + "\n", + "ax1.legend(loc='upper center')\n", + "ax2.legend(loc='upper center')\n", + "ax3.legend(loc='upper center')\n", + "\n", + "ax2.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n", + "ax1.set_ylabel(r'$\\rho(z)\\,\\mathrm{[mol/l]}$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see good agreement between our simulation and the meanfield solution of Guy and Chapman. Low density and reasonably low potential make the assumptions of the analytical approach justified." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now check how well the surface charge agrees with Grahame's equation.\n", + "To this end we calculate \n", + "$$\\sigma = \\int_0^{d/2} \\rho(z) \\,\\mathrm{d}z .$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:len(zs)//2]) * (zs[1] - zs[0])\n", + "sigma_right = np.sum((cation_profile_mean-anion_profile_mean)[len(zs)//2:]) * (zs[1] - zs[0])\n", + "\n", + "def grahame_sigma(phi):\n", + " return np.sinh(phi / 4.) * np.sqrt(2. * rho / (np.pi * BJERRUM_LENGTH))\n", + "\n", + "sigma_grahame = grahame_sigma(POTENTIAL_DIFF)\n", + "print(f'simulation: {sigma_right:.3f} e/nm^2')\n", + "print(f'grahame: {sigma_grahame:.3f} e/nm^2')\n", + "print(f'relative deviation: {abs(1. - sigma_right/sigma_grahame) * 100.:.0f}%')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The electric field is readily obtained from the integral \n", + "$$E(z) = \\int_0^{z} \\frac{1}{\\varepsilon_0 \\varepsilon_r} \\rho(z^\\prime) \\,\\mathrm{d}z^\\prime .$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the electric field\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "\n", + "dz_SI = (zs[1] - zs[0]) * constants.nano\n", + "chargedensity = (cation_profile_mean - anion_profile_mean) * constants.elementary_charge / constants.nano**3 \n", + "E_SI = 1. / (EPSILON_R * constants.epsilon_0) * np.cumsum(chargedensity * dz_SI)\n", + "# integration constant: zero field in the center\n", + "E_SI -= E_SI.min()\n", + "E = E_SI / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE) / constants.nano)\n", + "ax2 = plt.twinx()\n", + "\n", + "ax.plot(zs, E_SI)\n", + "ax2.plot(zs, E)\n", + "ax.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n", + "ax.set_ylabel(r'$E_\\mathrm{ind}\\,\\mathrm{[V/m]}$')\n", + "ax2.set_ylabel(r'$E_\\mathrm{ind}\\,\\mathrm{[(k_\\mathrm{B}T/e)/nm]}$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the electric field reduces to 0 in the middle of the channel, justifying the assumption that the two electrodes are far enough apart to not influence each other." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The electric potential can be calculated from $\\phi(z) = \\int_0^z -E(z^\\prime)\\,\\mathrm{d}z^\\prime$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the elecrostatic potential\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "ax2 = ax.twinx()\n", + "phi_SI = -np.cumsum(E_SI * dz_SI)\n", + "phi = phi_SI * (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))\n", + "ax.plot(zs, phi_SI)\n", + "ax.set_xlabel(r'$z\\,\\mathrm{[nm]}$')\n", + "ax.set_ylabel(r'$\\phi\\,[V]$')\n", + "ax2.set_ylabel(r'$\\phi\\,[k_\\mathrm{B}T/e]$')\n", + "ax2.axhline(-5, ls='--', color='k')\n", + "ax2.axhline(0, ls='--', color='k')\n", + "ax.set_xlim(0, 10. * DEBYE_LENGTH)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "measured_potential_difference = -(phi[-1] - phi[0])\n", + "print(f'applied voltage: {POTENTIAL_DIFF:.2f} V')\n", + "print(f'measured voltage: {measured_potential_difference:.2f} V')\n", + "print(f'relative deviation: {abs(1. - measured_potential_difference / POTENTIAL_DIFF) * 100:.0f}%')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Differential capacitance\n", + "\n", + "With the above knowledge, we can now assess the \n", + "differential capacitance of the system, by changing the applied voltage\n", + "difference and determining the corresponding surface charge density." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sigma_vs_phi = []\n", + "MIN_PHI = 0.5\n", + "MAX_PHI = 10\n", + "N_PHI = 7\n", + "N_SAMPLES_EQUIL_CAP = 5\n", + "N_SAMPLES_CAP = 5\n", + "\n", + "# sample from high to low potential to improve sampling\n", + "for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n", + " system.electrostatics.solver = setup_electrostatic_solver(potential_diff)\n", + " system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)\n", + " sigmas = []\n", + "\n", + " for tm in range(N_SAMPLES_CAP):\n", + " zs, density_accumulator_cation, density_accumulator_anion = \\\n", + " setup_densityprofile_accumulators(bin_width=DEBYE_LENGTH / 10.)\n", + "\n", + " system.auto_update_accumulators.clear()\n", + " system.auto_update_accumulators.add(density_accumulator_cation)\n", + " system.auto_update_accumulators.add(density_accumulator_anion)\n", + " system.integrator.run(STEPS_PER_SAMPLE)\n", + "\n", + " cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n", + " anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]\n", + "\n", + " sigmas.append(np.sum((cation_profile_mean - anion_profile_mean)[:int(len(zs) / 2.)]) * (zs[1] - zs[0]))\n", + "\n", + " sigma_vs_phi.append([potential_diff, np.mean(sigmas), scipy.stats.sem(sigmas)]) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "sigma_vs_phi = np.array(sigma_vs_phi)\n", + "x = np.linspace(0, sigma_vs_phi[:,0].max())\n", + "phi_SI = sigma_vs_phi[:,0] / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))\n", + "plt.errorbar(-sigma_vs_phi[:,1] * constants.elementary_charge / constants.nano**2,\n", + " phi_SI, xerr=sigma_vs_phi[:,2] * scipy.constants.elementary_charge / scipy.constants.nano**2,\n", + " fmt='o',label='Simulation')\n", + "plt.plot(grahame_sigma(x) * constants.elementary_charge / constants.nano**2,\n", + " x / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE)), label='Grahame')\n", + "x = np.linspace(0, ax.get_ylim()[1])\n", + "plt.plot(EPSILON_R * constants.epsilon_0 * x / 2. / (DEBYE_LENGTH * constants.nano), x, label='linear PB',\n", + " ls='--')\n", + "plt.xlabel(r'$\\sigma\\,\\mathrm{[C/m^2]}$')\n", + "plt.ylabel(r'$\\phi_\\mathrm{s}\\,\\mathrm{[V]}$')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For small potential drops, one observes the expected Poisson–Boltzmann behavior. It also agrees with the linearized solution $\\sigma(\\phi_\\mathrm{s}) = \\varepsilon_r\\varepsilon_0 \\frac{\\phi_\\mathrm{s}}{2 \\lambda_\\mathrm{D}}$.\n", + "However, we observe already for potentials $\\sim 0.1\\,\\mathrm{V} = 4\\,k_\\mathrm{B}T / e$ a significant deviation which can be attributed to the fact that our ions are of finite size and thus the surface charge saturates." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] Conway, B. E. Electrochemical Supercapacitors; Springer US: Boston, MA, 1999. https://doi.org/10.1007/978-1-4757-3058-6.\n", + "\n", + "[2] Tyagi, S.; Arnold, A.; Holm, C. Electrostatic Layer Correction with Image Charges: A Linear Scaling Method to Treat Slab 2D+h Systems with Dielectric Interfaces. J. Chem. Phys. 2008, 129 (20), 204102. https://doi.org/10.1063/1.3021064.\n", + "\n", + "[3] Gouy, G. Constitution of the Electric Charge at the Surface of an Electrolyte. J. Phys. 1910, 9 (4), 457–467.\n", + "\n", + "[4] Chapman, D. L. A Contribution to the Theory of Electrocapillarity. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1913, 25 (148), 475. https://doi.org/10.1080/14786440408634187.\n", + "\n", + "[5] Grahame, D. C. The Electrical Double Layer and the Theory of Electrocapillarity. Chem. Rev. 1947, 41 (3), 441–501. https://doi.org/10.1021/cr60130a002.\n", + "\n", + "[6] Tyagi, S.; Süzen, M.; Sega, M.; Barbosa, M.; Kantorovich, S. S.; Holm, C. An Iterative, Fast, Linear-Scaling Method for Computing Induced Charges on Arbitrary Dielectric Boundaries. J. Chem. Phys. 2010, 132 (15), 154112. https://doi.org/10.1063/1.3376011." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/tutorials/electrokinetics/CMakeLists.txt b/doc/tutorials/electrokinetics/CMakeLists.txt index a3c419cb6fc..cf2099bada2 100644 --- a/doc/tutorials/electrokinetics/CMakeLists.txt +++ b/doc/tutorials/electrokinetics/CMakeLists.txt @@ -18,6 +18,6 @@ # configure_tutorial_target(TARGET tutorial_ek DEPENDS electrokinetics.ipynb - figures/schlitzpore_3d.png scripts/eof_analytical.py) + figures/schlitzpore_3d.png) nb_export(TARGET tutorial_ek SUFFIX "" FILE "electrokinetics.ipynb" HTML_RUN) diff --git a/doc/tutorials/electrokinetics/electrokinetics.ipynb b/doc/tutorials/electrokinetics/electrokinetics.ipynb index cde7cc5e494..28510f0d347 100644 --- a/doc/tutorials/electrokinetics/electrokinetics.ipynb +++ b/doc/tutorials/electrokinetics/electrokinetics.ipynb @@ -4,346 +4,200 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Electrokinetics" + "# Electrokinetics\n", + "### Table of contents\n", + "1. [Introduction](#1.-Introduction)\n", + "2. [Advection-Diffusion in 2D](#2.-Advection-Diffusion-equation-in-2D)\n", + "3. [Electroosmotic flow](#3.-Electroosmotic-flow)\n", + "4. [Reaction in turbulent flow](#4.-Reaction-in-turbulent-flow)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Table of Contents\n", - "1. [Introduction](#Introduction)\n", - "2. [Theoretical Background](#Theoretical-Background)\n", - " 1. [The Electrokinetic Equations](#The-Electrokinetic-Equations)\n", - " 2. [EOF in the Slit Pore Geometry](#EOF-in-the-Slit-Pore-Geometry)\n", - "3. [Simulation using ESPResSo](#Simulation-using-ESPResSo)\n", - " 1. [Mapping SI and Simulation Units](#Mapping-SI-and-Simulation-Units)\n", - " 2. [Setting up the slit pore system](#Setting-up-the-slit-pore-system)\n", - "4. [References](#References)\n", - " " + "## 1. Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Introduction\n", + "In this tutorial we're looking at the electrokinetics feature of ESPResSo, which allows us to describe the motion of potentially charged chemical species solvated in a fluid on a continuum level. The govering equations for the solvent are known as the Poisson-Nernst-Planck equations, which is the combination of the electrostatic Poisson equation and the dynamics of the chemical species described by the Nernst-Planck equation. For the advection we solve the incompressible Navier-Stokes equation. The total set of equations is given by\n", "\n", - "In recent years the lattice-Boltzmann method (LBM) has proven itself to be a viable way to introduce hydrodynamic interactions into coarse-grained MD simulations with moderate computational cost.\n", - "ESPResSo features such an algorithm, which can make use of the LBM and extend it to coarse-grain not only the solvent molecules but also ionic solutes. It is called EK and explicitly treats the ionic solutes in a continuum fashion and is valid for a wide range of salt concentrations [1-3].\n", - "\n", - "### Tutorial Outline\n", - "\n", - "To make our first steps using ELECTROKINETICS we will work on one of the few systems for which analytic solutions for the electrokinetic equations exist: the slip pore geometry with a counterion-only electrolyte.\n", - "The same slit pore system is also treated in the LBM tutorial, but there, the ionic species were modeled as explicit particles.\n", - "For this system, the two approaches lead to exactly the same results [4].\n", - "Differences became significant for multivalent ions, very high salt concentrations, and very high surface charge, since then the mean-field approach the EK employs, is basically solving the Poisson-Nernst-Planck formalism plus the Navier-Stokes equation on a lattice.\n", - "This leads to significantly different results from explicit ion approaches [5-7].\n", - "This tutorial is mainly divided into two sections.\n", - "* **Theoretical Background** introduces the electrokinetic equations and the analytical solution for the slit pore system.\n", - "* **Simulation using ESPResSo** deals exclusively with the simulation. \n", + "\\begin{aligned}\n", + "\\partial_{t} n_{i} &= - \\vec{\\nabla} \\cdot \\vec{j}_{i} \\\\\n", + "\\vec{j}_{i} &= - D_{i} \\vec{\\nabla} n_{i} - \\frac{z_{i} e}{k_{B} T} n_{i} \\vec{\\nabla} \\phi + n_{i} \\vec{u} \\\\\n", + "\\Delta \\phi &= \\frac{1}{4 \\pi \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i} \\\\\n", + "\\rho (\\partial_{t} \\vec{u} + (\\vec{u} \\cdot \\vec{\\nabla}) \\vec{u}) &= - \\vec{\\nabla} p + \\eta \\Delta \\vec{u} + \\sum_{i} \\frac{k_{B} T}{D_{i}} \\vec{j}_{i} + \\vec{f}_{\\mathrm{ext}} \\\\\n", + "\\vec{\\nabla} \\cdot \\vec{u} &= 0,\n", + "\\end{aligned}\n", "\n", - "If you already know about simple diffusion-migration-advection equations, continuum electrostatics, and Navier-Stokes, then you can skip the first section." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Theoretical Background\n" + "where $n_{i}$ denotes the ion density of species $i$, $\\vec{j}_{i}$ the density flux, $D_{i}$ the diffusion coefficient, $z_{i}$ the valency, $e$ the elementary charge, $k_{B}$ the Boltzmann constant, $T$ the temperature, $\\phi$ the electrostatic potential, $\\varepsilon_{0}$ the vacuum permittivity, $\\varepsilon_{\\mathrm{r}}$ the relative permittivity, $\\rho$ the fluid density, $\\vec{u}$ the fluid velocity, $p$ the hydrostatic pressure, $\\eta$ the dynamic viscosity, and $\\vec{f}_{\\mathrm{ext}}$ an external force density." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### The Electrokinetic Equations\n", - "\n", - "In the following, we will derive the equations modeling the time evolution of the concentrations of dissolved species as well as the solvent in the standard electrokinetic model.\n", - "We do so, neglecting the salt ions' contribution to the overall mass density, which allows us to treat the dynamics of the ions and the fluid separately [7].\n", - "The solvent fluid will be modeled using the Navier-Stokes equations while we use a set of diffusion-migration-advection equations for the ionic species.\n" + "# 2. Advection-Diffusion equation in 2D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Ionic Species\n", - "\n", - "The description starts with the ionic species' concentrations $c_{k}(\\vec{r}, t)$ (number density) and the associated flux densities $\\vec{j}_{k}(\\vec{r}, t)$, for which mass conservation holds\n", - "\n", - "\\begin{equation}\n", - "\\partial_{t} c_{k} = -\\nabla \\cdot\\vec{j}_{k} . \n", - "\\end{equation}\n", - "\n", - "Here $\\vec{r}$ denotes the spatial coordinate and $t$ the time, while $k$ enumerates the ionic species.\n", - "The fluxes are caused by diffusion (due to density variations and external forces) and advection.\n", - "\n", - "The advective contribution to the flux is given by\n", + "The first system that is simulated in this tutorial is the simple advection-diffusion of a drop of uncharged chemical species in a constant velocity field. To keep the computation time small, we restrict ourselves to a 2D problem, but the algorithm is also capable of solving the 3D advection-diffusion equation. Furthermore, we can also skip solving the electrostatic Poisson equation, since there are is no charged species present. The equations we solve thus reduce to\n", "\n", "\\begin{equation}\n", - "\\vec{j}_{k}^{\\mathrm{adv.}} = c_{k} \\vec{u} ,\n", + "\\partial_{t} n = D \\Delta n - \\vec{\\nabla} \\cdot (\\vec{v} n).\n", "\\end{equation}\n", "\n", - "where $\\vec{u}(\\vec{r}, t)$ denotes the fluid velocity (advective velocity).\n", - "This equation models advection as a simple co-movement of the dissolved ions with the surrounding fluid.\n", - "All inertial effects of the ions are neglected.\n", - "\n", - "The diffusive behavior of the ions is best described in a reference frame co-moving with the local advective velocity $\\vec{u}$.\n", - "We assume that the species' relative fluxes instantaneously relax to a local equilibrium.\n", - "This assumption allows us to derive the diffusive fluxes from a local free-energy density, which we define as\n", + "The fundamental solution of this partial diffential equation can be found analytically in the case of a constant velocity field $\\vec{v}$ and a constant diffusion coefficient $D$. For a $d$-dimensional system, the solution of an initally infinitessimaly small droplet at the origin can be written as\n", "\n", "\\begin{equation}\n", - "f \\big( c_{k}(\\vec{r}) \\big) = \\sum_{k} \\underbrace{k_{\\mathrm{B}}T c_{k}(\\vec{r}) \\left[ \\log \\left\\lbrace \\Lambda_{k}^{3} c_{k}(\\vec{r}) \\right\\rbrace - 1 \\right] }_{\\mathrm{ideal~gas~contribution}} + \\underbrace{z_{k} e c_{k}(\\vec{r}) \\Phi(\\vec{r})}_{\\mathrm{electrostatic~contribution}} ,\n", + "n(\\vec{x},t) = \\frac{1}{(4 \\pi D t)^{d/2}} \\exp \\left( - \\frac{(\\vec{x} - \\vec{v} t)^2}{4 D t} \\right).\n", "\\end{equation}\n", "\n", - "with the $\\Lambda_{k}$ the species' thermal de Broglie wavelengths, $z_{k}$ their valencies, and $\\Phi(\\vec{r})$ the electrostatic potential.\n", - "This free-energy density consists of only an ideal-gas and an electrostatic contribution.\n", - "The same assumptions form the basis of Poisson-Boltzmann (PB) theory.\n", - "Hence, the limitations of this model are the same as those of PB.\n", - "That means this model applies to monovalent ions at low to intermediate densities and surface charges [5,6,10,11].\n", - "\n", - "The species' chemical potentials $\\mu_{k}$ implied by the free-energy density read\n", - "\n", - "\\begin{equation}\n", - "\\mu_{k}(\\vec{r}) = \\delta_{c_k} f(c_{k}\\big( \\vec r ) \\big) = k_{\\mathrm{B}}T \\log(\\Lambda_{k}^{3} c_{k}(\\vec{r})) + z_{k} e \\Phi(\\vec{r}) .\n", - "\\end{equation}\n", + "After importing the necessary packages, we start by defining the necessary parameters for the simulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import espressomd\n", + "import espressomd.lb\n", + "import espressomd.electrokinetics\n", + "import espressomd.shapes\n", "\n", - "This in turn allows us to formulate the first-order approximation to the thermodynamic driving force as the gradient of the chemical potential, which we use to define an expression for the diffusive flux\n", + "espressomd.assert_features([\"WALBERLA\", \"WALBERLA_FFT\"])\n", "\n", - "\\begin{aligned}\n", - "\\vec{j}_{k}^{\\mathrm{diff}} &= \\xi_{k} \\left( -c_{k} \\nabla \\mu_{k} \\right) = -k_{\\mathrm{B}}T \\xi_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla\\Phi \\\\\n", - "&= -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi .\n", - "\\end{aligned}\n", + "import tqdm\n", + "import numpy as np\n", "\n", - "Here, $\\xi_{k}$ and $D_{k}$ denote the mobility and the diffusion coefficient of species $k$, which are related by the Einstein-Smoluchowski relation $D_{k} / \\xi_{k} = k_{\\mathrm{B}}T$ [11,12].\n", + "import scipy.optimize\n", "\n", - "Finally, the total number density flux combining effects of diffusion and advection reads\n", + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.animation as animation\n", + "import matplotlib.quiver\n", "\n", - "\\begin{equation}\n", - "\\vec{j}_{k} = \\vec{j}_{k}^{\\mathrm{diff}} + \\vec{j}_{k}^{\\mathrm{adv.}} = -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi + c_{k} \\vec{u} , \n", - "\\end{equation}\n", + "import tempfile\n", + "import base64\n", "\n", - "where the first term represents Fick's law of diffusion in the absence of an external potential, the second term gives the additional flux due to an external (in this case electrostatic) potential, and the last term introduces the influence of the motion of the underlying solvent." + "plt.rcParams.update({'font.size': 16})" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "#### Electrostatics\n", - "\n", - "The dynamics of the charged species in a typical micro- or nanofluidic system are slow compared to the relaxation of the electromagnetic fields.\n", - "This allows us to use stationary equations to model electromagnetic effects.\n", - "We further assume that the modeled species do not carry permanent magnetic dipoles and that electric currents in the system are small.\n", - "Under these conditions, the full set of Maxwell's equations reduces to the Poisson equation\n", - "\n", - "\\begin{equation}\n", - "\\nabla^2 \\Phi = - \\frac{1}{\\varepsilon} \\sum_{k} z_{k} e c_{k} = -4 \\pi l_\\mathrm{B} k_{\\mathrm{B}}T \\sum_{k} z_{k} c_{k} . \n", - "\\end{equation}\n", - "\n", - "Here $\\varepsilon = \\varepsilon_{0} \\varepsilon_r$ denotes the product of the vacuum permittivity $\\varepsilon_{0}$ with the relative permittivity of the solvent $\\varepsilon_r$.\n", - "We have also used the Bjerrum-length\n", - "\n", - "\\begin{equation}\n", - "l_\\mathrm{B} = \\frac{e^{2}}{4 \\pi \\varepsilon k_{\\mathrm{B}}T}.\n", - "\\end{equation}\n", - "\n", - "Finally, we have assumed that the permittivity is spatially homogeneous, since this will allow us to use efficient spectral methods to solve this equation.\n" + "BOX_L = [80, 80, 1] \n", + "AGRID = 1.0\n", + "DIFFUSION_COEFFICIENT = 0.06\n", + "TAU = 1.0\n", + "EXT_FORCE_DENSITY = [0, 0, 0]\n", + "FLUID_DENSITY = 1.0\n", + "FLUID_VISCOSITY = 1.\n", + "FLUID_VELOCITY = [0.04, 0.04, 0.0]\n", + "KT = 1.0\n", + "\n", + "RUN_TIME = 400" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "#### Hydrodynamics\n", - "\n", - "As said before, since the ionic species' contribute at most a few percent to the overall mass, we can safely approximate the overall fluid's mass by the mass of the solvent (typically water) and model the solvents velocity field $\\vec{u}(\\vec{r}, t)$ using the Navier-Stokes equations for an isotropic, incompressible Newtonian fluid\n", - "\n", - "\\begin{aligned}\n", - "\\rho \\big( \\partial_t \\vec{u} + \\left(\\vec{u} \\cdot \\nabla \\right) \\vec{u} \\big) &= -\\nabla p_H + \\eta \\nabla^{2} \\vec{u} + \\vec{f} ,\\\\\n", - "\\nabla \\cdot \\vec u &= 0 .\n", - "\\end{aligned}\n", - "\n", - "where $p_H$ denotes the hydrostatic pressure, $\\eta$ the shear viscosity, $\\rho$ the density of the fluid, and $\\vec{f}$ an external body force density.\n", - "For the assumption of incompressibility to hold, the Mach number needs to be small – a condition that is fulfilled for micro- and nanofluidic systems with flow velocities on the order of μm/s.\n", - "\n", - "Earlier we assumed that the ions' velocity relative to the fluid instantaneously relaxes to a stationary state and that this stationary velocity is given by the product of their mobility and the force exerted on them.\n", - "For this state to be stationary, all the momentum transferred into the ions by the external force needs to be dissipated into the fluid immediately.\n", - "From this we can conclude that the force density acting on the fluid must read\n", - "\n", - "\\begin{equation}\n", - "\\vec{f} = \\sum_{k} \\vec{j}^\\mathrm{diff}_k / \\xi_{k} = -\\sum_{k} (k_\\mathrm{B}T \\nabla c_{k} + z_{k} e c_{k} \\nabla \\Phi) .\n", - "\\end{equation}\n", - "\n", - "Summarizing, the set of electrokinetic equations we solve is given by\n", - "\n", - "\\begin{aligned}\n", - "\\vec{j}_{k} &= -D_{k} \\nabla c_{k} - \\xi_{k} z_{k} e c_{k} \\nabla \\Phi + c_{k} \\vec{u} ,\\\\\n", - "\\partial_{t} c_{k} &= -\\nabla \\cdot\\vec{j}_{k} ,\\\\\n", - "\\nabla^2 \\Phi &= -4 \\pi l_\\mathrm{B} k_\\mathrm{B}T \\textstyle\\sum_{k} z_{k} c_{k} ,\\\\\n", - "\\rho \\big( \\partial_t \\vec{u} + (\\vec{u} \\cdot \\nabla ) \\vec{u} \\big) &= -\\nabla p_H + \\eta \\nabla^{2} \\vec{u} - \\textstyle\\sum_{k} (k_\\mathrm{B}T \\nabla c_{k} + z_{k} e c_{k} \\nabla \\Phi) ,\\\\\n", - "\\nabla \\cdot \\vec{u} &= 0 .\n", - "\\end{aligned}" + "system = espressomd.System(box_l=BOX_L)\n", + "system.time_step = TAU\n", + "system.cell_system.skin = 0.4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### EOF in the Slit Pore Geometry\n", - "\n", - "The slit pore system depicted in Fig. 1 consists of two like charged parallel plates of infinite extent, confining a solution of water and the plates' counterions.\n", - "\n", - "
\n", - " missing\n", - "
\n", - "
Figure 1: Slit pore system and coordinate system used for the analytical calculations.
\n", - "
\n", - "
\n", - "\n", - "Due to the net neutrality of the system and the translational symmetry in directions parallel to the plates, the potential outside the two plates must be constant.\n", - "This means that using periodic or non-periodic boundary conditions makes no difference.\n", - "As the system is in equilibrium in the normal direction, the electrokinetic equations for this dimension reduce to the Poisson-Boltzmann equation for the electrostatic potential, which reads\n", - "\\begin{equation}\n", - "\\partial_x^2 \\Phi(x) = -4 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\, ze \\, c_0 \\cdot \\exp{\\left(-\\frac{ze\\Phi(x)}{k_\\mathrm{B}T}\\right)} \\; ,\n", - "\\end{equation}\n", - "where $x$ denotes the direction normal to the plates.\n", - "The constant $c_0$ has to be chosen such that charge neutrality is fulfilled.\n", - "Multiplying by $2 \\partial_x \\Phi(x)$ and applying the inverse chain rule further reduces the equation to first order.\n", - "Subsequent separation of variables yields the solution\n", - "\\begin{equation}\n", - "\\Phi(x) = -\\frac{k_\\mathrm{B}T}{ze} \\cdot \\log \\left[ \\frac{C^2}{8 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B}} \\cdot \\cos^{-2}\\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x\\right) \\right], \\quad \\left| \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right| < \\frac \\pi 2\\; .\n", - "\\end{equation}\n", - "Refer to [4] for details on this calculation.\n", - "Knowing that the counterion density $c$ resembles a Boltzmann distribution in the potential $ze \\Phi$ leads to the expression\n", - "\\begin{equation}\n", - "c(x) = \\frac{C^2}{8 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B}} \\cdot \\cos^{-2} \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right) \\; .\n", - "\\end{equation}\n", - "The constant $C$ is determined by fixing the number of counterions or requiring the E-field to vanish outside the volume contained by the plates.\n", - "Both yields\n", - "\\begin{equation}\n", - "C \\cdot \\tan \\left( \\frac{zed}{4 k_\\mathrm{B}T} \\cdot C \\right) = -4 \\pi \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\sigma \\; ,\n", - "\\end{equation}\n", - "where $d$ denotes the distance between the plates and $\\sigma$ their (constant) surface charge density.\n", - "\n", - "Applying an electric field along one of the directions parallel to the plates does not influence the charge distribution in the normal direction, which allows us to write down the hydrodynamic equations for the parallel direction.\n", - "After eliminating all terms from the Navier-Stokes Equations, which vanish due to symmetry, we are left with\n", - "\\begin{equation}\n", - "\\frac{\\partial_x^2 v_y(x)}{\\partial x^2} = -\\frac{q E C^2}{8 \\, k_\\mathrm{B}T \\, l_\\mathrm{B} \\, \\eta} \\cdot \\cos^{-2} \\left( \\frac{qC}{2 k_\\mathrm{B}T} \\cdot x \\right) \\; ,\n", - "\\end{equation}\n", - "which yields, by means of simple integration and the application of no-slip boundary conditions\n", - "\\begin{equation}\n", - "v_y(x) = \\frac{E}{2 \\pi \\, l_\\mathrm{B} \\, \\eta \\, ze} \\cdot \\log \\left[ \\frac{\\cos \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot x \\right)}{\\cos \\left( \\frac{zeC}{2 k_\\mathrm{B}T} \\cdot \\frac d 2 \\right)} \\right] \\; .\n", - "\\end{equation}\n", - "\n", - "With this tutorial comes a Python script eof_analytical.py, which evaluates all these expressions on the same grid as is used in the upcoming simulation." + "We use a lattice Boltzmann flow field with constant velocity for advection.\n", + "Note that we have to set ``kT=0.0`` here to avoid random fluctuations in the flow velocity." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "## Simulation using ESPResSo" + "lattice = espressomd.lb.LatticeWalberla(n_ghost_layers=1, agrid=AGRID)\n", + "lbf = espressomd.lb.LBFluidWalberla(\n", + " lattice=lattice, density=FLUID_DENSITY, kinematic_viscosity=FLUID_VISCOSITY,\n", + " tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=0.0, seed=42)\n", + "lbf[:, :, :].velocity = FLUID_VELOCITY\n", + "system.lb = lbf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Mapping SI and Simulation Units\n", - "\n", - "ESPResSo does not predefine any unit system.\n", - "This makes it more flexible but also requires us to spend some thought on the conversion from SI units to simulation units and back.\n", - "Since most first time users have trouble with this, we will go through that process in detail here.\n", - "\t\n", - "Important to note is that ESPResSo's unit system is nothing more than a rescaled variant of the SI system.\n", - "All physical formulas you are used to in the SI system remain valid and you can use them to find relations between your units.\n", - "Let's start by choosing a unit of length.\n", - "Since we are going to deal with Debye layers with extensions of nanometers, a sensible choice is\n", - "\n", - "\\begin{equation}\n", - "[x]=1\\mathrm{nm}.\n", - "\\end{equation}\n", - "\n", - "The involved energies are of the magnitude of $k_{\\mathrm{B}}T$.\n", - "We will simulate our system at room temperature ($300\\mathrm{K}$), hence we use as unit of energy\n", - "\\begin{equation}\n", - "[E]=k_\\mathrm{B}\\cdot 300\\mathrm{K}\\approx 4.14 \\cdot 10^{-21}\\mathrm{J}.\n", - "\\end{equation}\n", - "\n", - "By default ESPResSo has no concept for particle masses (but the feature can be activated).\n", - "That means all particle masses are assumed to be $1\\,[\\mathrm{m}]$, which forces us to use the particle mass as mass unit.\n", - "For this simulation we use the mass of sodium ions, which is\n", - "\\begin{equation}\n", - "[m]=23\\mathrm{u}\\approx 3.82\\cdot 10^{-26}\\mathrm{kg}.\n", - "\\end{equation}\n", - "\n", - "For the relation\n", - "\\begin{equation}\n", - "E=\\frac 1 2 mv^2\n", - "\\end{equation}\n", - "\n", - "to hold, the unit of time has to be defined so that\n", - "\\begin{equation}\n", - "[E]=[m]\\cdot\\frac{[x]^2}{[t]^2}.\n", - "\\end{equation}\n", - "\n", - "From that we get the missing unit of time\n", - "\\begin{equation}\n", - "[t]=[x]\\cdot\\sqrt{\\frac{[m]}{[E]}}=1\\mathrm{nm}\\cdot\\sqrt{\\frac{23\\mathrm{u}}{k_B\\cdot 300\\mathrm{K}}}\\approx 3.03760648\\cdot 10^{-12}\\mathrm{s}\\approx 3.04\\mathrm{ps}.\n", - "\\end{equation}\n", - "\n", - "The last unit we need is the one of charge.\n", - "We choose it to be the elementary charge\n", - "\\begin{equation}\n", - "[q]=e\\approx 1.60\\cdot 10^{-19}\\mathrm{C}.\n", - "\\end{equation}\n", - "\n", - "We now have all the units necessary to convert our simulation parameters." + "To use the electrokinetics-algorithm in ESPResSo, one needs to create an instance of the `EKContainer`-object and pass it a time step `tau` and Poisson solver `solver`.\n", + "Since our species is uncharged, we don't need to solve the electrostatic Poisson equation, so we can use the placeholder-class, which is called `EKNone`." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "|parameter |value (SI units) | value (simulation units)|\n", - "|:---------|----------------:|------------------------:|\n", - "|channel width $d$ | $50\\mathrm{nm}$ | $50\\mathrm{[x]}$|\n", - "|thermal energy $k_B T$ | $k_B\\cdot 300\\mathrm{K}$ | $1\\mathrm{[E]}$|\n", - "|Bjerrum length $l_B$ | $0.7095\\mathrm{nm}$ | $0.7095\\mathrm{[x]}$|\n", - "|counterion charge $q$ | $1e$ | $1\\mathrm{[q]}$|\n", - "|counterion diffusion coefficient $D$ | $2.0\\cdot 10^{-9}\\mathrm{m^2/s}$ | $0.006075\\mathrm{[x]^2/[t]}$|\n", - "|solvent density $\\rho$ | $1.0\\cdot 10^{3}\\mathrm{kg/m^3}$ | $26.18\\mathrm{[m]/[x]^3}$|\n", - "|solvent dynamic viscosity $\\eta$ | $1.0\\cdot 10^{-3}\\mathrm{Pa}\\mathrm{s}$ | $79.53\\mathrm{[m]/([x][t])}$|\n", - "|external electric field $E$ | $2.585\\cdot 10^{6}\\mathrm{V/m}$ | $0.1\\mathrm{[E]/([q][x])}$|\n" + "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n", + "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "ESPResSo determines the strength of the electrostatic interactions via the Bjerrum-length $l_\\mathrm{B}$.\n", - "That is the length for which the electrostatic interaction energy of two elementary charges equals the thermal energy\n", - "\n", - "\\begin{equation}\n", - "k_\\mathrm{B} T=\\frac{e^2}{4\\pi\\varepsilon_0\\varepsilon_r}\\cdot\\frac 1 {l_\\mathrm{B}}.\n", - "\\end{equation}\n", - "\n", - "This yields for water at $300K$ with $\\varepsilon_r = 78.54$, a Bjerrum length of $l_\\mathrm{B}\\approx 0.7095\\mathrm{nm}$." + "Now, we can add diffusive species to the container to integrate their dynamics." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, "source": [ - "### Setting up the slit pore system\n", - "\n", - "The script for this simulation comes with this tutorial and is called eof_electrokinetics.py.\n", - "All used commands are documented in the User's Guide in the section called **Electrokinetics**.\n", + "# Exercise:\n", + "- Create an instance of the [`espressomd.electrokinetics.EKSpecies`]() and add it to the system with [`system.ekcontainer.add()`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer.add). \n", "\n", - "We first set up a periodic simulation box of the desired dimensions.\n", - "Note that the dimensions are, of course, given in simulation units." + "# Hint:\n", + "- Use the variables `DIFFUSION_COEFFICIENT`, `KT` and `TAU` defined above.\n", + "- Enable both `advection` and `friction_coupling`.\n", + "- Make sure to initialize the `density` with 0.0, and disable electrostatics by setting `valency` to 0.0 as well." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "species = espressomd.electrokinetics.EKSpecies(\n", + " lattice=lattice, density=0.0, kT=KT,\n", + " diffusion=DIFFUSION_COEFFICIENT, valency=0.0,\n", + " advection=True, friction_coupling=True,\n", + " ext_efield=[0., 0., 0.], tau=TAU)\n", + "system.ekcontainer.add(species)\n", + "```" ] }, { @@ -351,38 +205,13 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "# Initializing espresso modules and the numpy package\n", - "import espressomd\n", - "import espressomd.lb\n", - "import espressomd.electrokinetics\n", - "import espressomd.shapes\n", - "\n", - "espressomd.assert_features([\"WALBERLA\", \"WALBERLA_FFT\"])\n", - "\n", - "import tqdm\n", - "import numpy as np\n", - "%matplotlib inline\n", - "import matplotlib.pyplot as plt\n", - "plt.rcParams.update({'font.size': 16})\n", - "\n", - "box_y = 6\n", - "box_z = 6\n", - "width = 50\n", - "\n", - "padding = 1\n", - "box_x = width + 2 * padding\n", - "\n", - "system = espressomd.System(box_l=[box_x, box_y, box_z])" - ] + "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We then store all the parameters we calculated earlier.\n", - "At this point, these parameters only reside in Python variables.\n", - "They will only be used by ESPResSo once they are being passed to the respective initialization functions." + "To compare our simulation to the fundamental solution of the advection-diffusion equation, we need to approximate a delta-droplet, which can be achieved by having a non-zero density only at the center of the domain." ] }, { @@ -391,28 +220,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Set the electrokinetic parameters\n", - "\n", - "agrid = 1.0\n", - "dt = 0.5\n", - "kT = 4.0\n", - "bjerrum_length = 0.7095\n", - "permittivity = 1. / (4 * np.pi * bjerrum_length)\n", - "D = 0.006075\n", - "valency = 1.0\n", - "viscosity_dynamic = 79.53\n", - "density_water = 26.15\n", - "sigma = -0.05\n", - "ext_force_density = [0.0, 0.1, 0.0]\n", - "\n", - "single_precision = False" + "system.ekcontainer[0][BOX_L[0] // 2, BOX_L[1] // 2, 0].density = 1.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Before we initialize the actual electrokinetics algorithm, we need to set the time step and some other parameters that are not actually used, but would otherwise lead to error messages." + "Now everything is set and we can finally run the simulation by running the integrator." ] }, { @@ -421,19 +236,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Set the simulation parameters\n", - "\n", - "system.time_step = dt\n", - "system.cell_system.skin = 0.2\n", - "system.thermostat.turn_off()\n", - "integration_length = 600" + "system.integrator.run(RUN_TIME)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can now set up the electrokinetics algorithm." + "For comparison, we prepare the analytical solution and show the 2D-density as well as a slice through the center of the droplet." ] }, { @@ -442,7 +252,20 @@ "metadata": {}, "outputs": [], "source": [ - "lattice = espressomd.lb.LatticeWalberla(agrid=agrid, n_ghost_layers=1)" + "def calc_gaussian(pos: np.ndarray, time: int, D: float):\n", + " dim = pos.shape[-1]\n", + " return (4 * np.pi * D * time)**(-dim / 2) * np.exp(-np.sum(np.square(pos), axis=-1) / (4 * D * time))\n", + "\n", + "pos = np.zeros((*BOX_L[:2], 2))\n", + "xpos = np.arange(-BOX_L[0] // 2, BOX_L[0] // 2)\n", + "ypos = np.arange(-BOX_L[1] // 2, BOX_L[1] // 2)\n", + "\n", + "pos[..., 1], pos[..., 0] = np.meshgrid(xpos, ypos)\n", + "\n", + "# add the advection shift\n", + "pos -= np.asarray(FLUID_VELOCITY[:2]) * RUN_TIME * TAU\n", + "\n", + "analytic_density = calc_gaussian(pos=pos, time=RUN_TIME * TAU, D=DIFFUSION_COEFFICIENT)" ] }, { @@ -451,11 +274,15 @@ "metadata": {}, "outputs": [], "source": [ - "viscosity_kinematic = viscosity_dynamic / density_water\n", - "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=density_water,\n", - " kinematic_viscosity=viscosity_kinematic,\n", - " tau=dt, single_precision=single_precision)\n", - "system.lb = lbf" + "fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(15, 7))\n", + "\n", + "ax1.imshow(system.ekcontainer[0][:, :, 0].density, origin=\"lower\", vmin=0, vmax=6e-3)\n", + "ax1.set_title(\"simulation\")\n", + "\n", + "imshow = ax2.imshow(analytic_density, origin=\"lower\", vmin=0, vmax=6e-3)\n", + "ax2.set_title(\"analytic\")\n", + "fig.colorbar(imshow, ax=[ax1, ax2], shrink=0.8)\n", + "plt.show()" ] }, { @@ -464,70 +291,120 @@ "metadata": {}, "outputs": [], "source": [ - "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice,\n", - " permittivity=permittivity,\n", - " single_precision=single_precision)\n", + "values_diagonal = np.diagonal(system.ekcontainer[0][:, :, 0].density)\n", + "analytic_diagonal = np.diagonal(analytic_density)\n", + "positions_diagonal = np.arange(len(values_diagonal)) * np.sqrt(2) * AGRID\n", + "\n", + "def gaussian_kernel(x, magnitude, mu, sigma):\n", + " return magnitude * np.exp(-(x - mu)**2 / (2 * sigma**2))\n", + "\n", + "popt, pcov = scipy.optimize.curve_fit(gaussian_kernel, positions_diagonal, analytic_diagonal, p0=[0.05,70,5.])\n", + "positions_analytic = np.concatenate([[positions_diagonal[0]],\n", + " np.linspace(popt[1] - 5 * popt[2], popt[1] + 5 * popt[2], 120),\n", + " [positions_diagonal[-1]]])\n", + "values_analytic = gaussian_kernel(positions_analytic, *popt)\n", + "\n", + "fig = plt.figure(figsize=(8, 5))\n", + "ax = fig.gca()\n", + "ax.plot(positions_diagonal, values_diagonal, \"o\", mfc=\"none\", label=\"simulation\")\n", + "ax.plot(positions_analytic, values_analytic, label=\"analytic\")\n", + "\n", + "ax.set_xlabel(\"position\")\n", + "ax.set_ylabel(\"density\")\n", "\n", - "system.ekcontainer = espressomd.electrokinetics.EKContainer(solver=eksolver, tau=dt)" + "plt.legend()\n", + "plt.show()" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "density_counterions = -2.0 * sigma / width\n", - "ekspecies = espressomd.electrokinetics.EKSpecies(\n", - " lattice=lattice, density=0.0, kT=kT, diffusion=D, valency=valency,\n", - " advection=True, friction_coupling=True, ext_efield=ext_force_density,\n", - " single_precision=single_precision, tau=dt)\n", - "system.ekcontainer.add(ekspecies)" + "From the plot one can see that the position of the density-peak matches well. However, one also sees that the droplet in the simulation has spread more than it should. The reason is that the discretization used for the advection term introduces an artifical, additional diffusion to the system. This is a fundamental limitation of the algorithm, which is why it cannot be applied to pure advection problems." ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "ekwallcharge = espressomd.electrokinetics.EKSpecies(\n", - " lattice=lattice, density=0.0, kT=kT, diffusion=0., valency=-valency,\n", - " advection=False, friction_coupling=False, ext_efield=[0, 0, 0],\n", - " single_precision=single_precision, tau=dt)\n", - "system.ekcontainer.add(ekwallcharge)" + "# 3. Electroosmotic flow" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "wall_left = espressomd.shapes.Wall(normal=[1, 0, 0], dist=padding)\n", - "wall_right = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(padding + width))" + "The next system in this tutorial is a simple slit pore, as shown in Figure 1. It consists of an infinite plate capacitor with an electrolytic solution trapped in between the plates. The plates of the capactior carry a constant surface charge and the counterions are solvated in the liquid. \n", + "\n", + "Charge attraction will cause the ions to accumulate near the surfaces, forming a characteristic ion density profile, which can be calculated analytically using the Poisson-Boltzmann equation. Since the system has translational symmetry in the directions parallel to the plates, the equations for parallel and orthogonal direction decouple. This means that applying an external electric field in a direction parallel to the plates will not change the distribution of the ions along the orthogonal direction. It will however cause motion of the ions and consequently the fluid: The characteristic flow profile of electroosmotic flow.\n", + "\n", + "
\n", + "\n", + "
\n", + "
Figure 1: Slit pore system and coordinate system used for the analytical calculations.
\n", + "
\n", + "
" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "ekspecies[padding:-padding, :, :].density = density_counterions" + "### Analytical solution\n", + "\n", + "Due to the symmetries of the system, it effectively reduces to a 1D problem along the orthogonal axis. The system can be described by the Poisson-Boltzmann equation:\n", + "\n", + "$$\n", + "\\partial_{x}^2 \\phi(x) = \\frac{1}{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sum_{i} z_{i} e n_{i}(x) \\exp \\left( -\\frac{z_{i} e \\phi(x)}{k_{\\mathrm{B}} T} \\right)\n", + "$$\n", + "\n", + "where $x$ is the normal-direction of the plates. Since we will only simulate a single ion species, the counterions, the sum only has a single summand. The solution for the potential is then given by:\n", + "\n", + "$$\n", + "\\phi(x) = -\\frac{k_{B}T}{z e} \\log \\left[ \\frac{C^2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}}{2 k_{B}T } \\cos^{-2} \\left( \\frac{z e C}{2 k_{B} T} x \\right) \\right], \\qquad \\text{with } \\left\\| \\frac{z e C}{2 k_{B} T} \\right\\| < \\frac{\\pi}{2},\n", + "$$\n", + "\n", + "where $C$ is an integration constant that is to be determined by the boundary conditions. The ion density follows then from the potential as\n", + "\n", + "$$\n", + "n(x) = \\frac{C^2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}}}{2 k_{B}T} \\cos^{-2} \\left( \\frac{z e C}{2 k_{B} T} x \\right).\n", + "$$\n", + "\n", + "To find the integration constant we use fact that the total system has to be charge neutral, i.e., the total charge on the plates is counterbalanced by the counterions. This leads to the following equation\n", + "\n", + "$$\n", + "C \\tan \\left( \\frac{z e d}{4 k_{B} T} C \\right) = - \\frac{e^2}{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}}} \\sigma,\n", + "$$\n", + "\n", + "where $\\sigma$ is the surface charge density of the plates. This is a transcendental equation, which must be solved numerically to find $C$. \n", + "\n", + "The electric field is applied in the $y$-direction, parallel to the plates.\n", + "Fluid flow is described by the incompressible Navier-Stokes equation, which due to the symmetries of the system reduces to the one-dimensional problem\n", + "\n", + "$$\n", + "\\frac{\\partial^2 v_{y}(x)}{\\partial x^2} = - \\frac{\\varepsilon_{0} \\varepsilon_{\\mathrm{r}} z e E C^2}{2 k_{B}T \\eta} \\cos^{-2}\\left( \\frac{q C}{2 k_{B} T} x \\right).\n", + "$$\n", + "\n", + "This equation can be solved analytically and the solution is given by\n", + "\n", + "$$\n", + "v_{y}(x) = \\frac{2 \\varepsilon_{0} \\varepsilon_{\\mathrm{r}} k_{B} T E}{\\eta z e} \\log \\left( \\frac{\\cos \\left( \\displaystyle\\frac{z e C}{2 k_{B} T} x \\right)}{\\cos \\left( \\displaystyle\\frac{z e C}{2 k_{B} T} \\frac{d}{2} \\right)} \\right),\n", + "$$\n", + "\n", + "where $d$ denotes the distance between the two plates. Finally, the shear-stress of this problem is given by\n", + "\n", + "$$\n", + "\\sigma(x) = \\mu \\frac{\\partial v_{y}(x)}{\\partial x}\n", + "$$" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "ekspecies[:padding, :, :].density = 0.0\n", - "ekspecies[-padding:, :, :].density = 0.0\n", + "### Numerical solution\n", "\n", - "ekwallcharge[:padding, :, :].density = -sigma / valency / padding\n", - "ekwallcharge[-padding:, :, :].density = -sigma / valency / padding" + "We start by resetting the system and defining the necessary parameters." ] }, { @@ -536,12 +413,8 @@ "metadata": {}, "outputs": [], "source": [ - "for shape_obj in (wall_left, wall_right):\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=[0., 0., 0.],\n", - " boundary_type=espressomd.electrokinetics.FluxBoundary)\n", - " ekspecies.add_boundary_from_shape(shape=shape_obj, value=0.0,\n", - " boundary_type=espressomd.electrokinetics.DensityBoundary)\n", - " lbf.add_boundary_from_shape(shape=shape_obj, velocity=[0., 0., 0.])" + "system.ekcontainer = None\n", + "system.lb = None" ] }, { @@ -550,101 +423,44 @@ "metadata": {}, "outputs": [], "source": [ - "# Integrate the system\n", - "for i in tqdm.trange(100):\n", - " system.integrator.run(integration_length)\n", - "\n", - "# Output\n", - "position_list = []\n", - "density_list = []\n", - "velocity_list = []\n", - "pressure_xy_list = []\n", - "\n", - "for i in range(int(box_x / agrid)):\n", - " if (i * agrid >= padding) and (i * agrid < box_x - padding):\n", - " position = i * agrid - padding - width / 2.0 + agrid / 2.0\n", - " position_list.append(position)\n", - " \n", - " node_idxs = (i, int(box_y / (2 * agrid)), int(box_z / (2 * agrid)))\n", + "AGRID = 1.0\n", + "TAU = 1.0\n", + "KT = 2.0\n", + "PERMITTIVITY = 0.28\n", + "DIFFUSION_COEFFICIENT = 0.25\n", + "VALENCY = 1.0\n", + "VISCOSITY_DYNAMIC = 0.5\n", + "DENSITY_FLUID = 1.0\n", + "SURFACE_CHARGE_DENSITY = -0.05\n", + "EXT_FORCE_DENSITY = [0.0, 0.01, 0.0]\n", + "\n", + "SINGLE_PRECISION = False\n", "\n", - " # density\n", - " density_list.append(ekspecies[node_idxs].density)\n", - "\n", - " # velocity\n", - " velocity_list.append(lbf[node_idxs].velocity[1])\n", + "padding = 1\n", + "WIDTH = 126\n", + "BOX_L = [(WIDTH + 2 * padding) * AGRID, 1, 1]\n", "\n", - " # xz component pressure tensor\n", - " pressure_xy_list.append(lbf[node_idxs].pressure_tensor[0, 1])\n", + "system.cell_system.skin = 0.4\n", + "system.box_l = BOX_L\n", + "system.time_step = TAU\n", "\n", - "np.savetxt(\"eof_simulation.dat\",\n", - " np.column_stack((position_list,\n", - " density_list,\n", - " velocity_list,\n", - " pressure_xy_list)),\n", - " header=\"#position calculated_density calculated_velocity calculated_pressure_xy\")" + "RUN_TIME = 200" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], + "cell_type": "markdown", + "metadata": {}, "source": [ - "from scripts import eof_analytical # executes automatically upon import\n", - "\n", - "# read analytical solution and simulation data\n", - "data_an = np.loadtxt(\"eof_analytical.dat\")\n", - "data_ek = np.loadtxt(\"eof_simulation.dat\")\n", - "\n", - "fig1 = plt.figure(figsize=(16, 4.5))\n", - "ax = fig1.add_subplot(131)\n", - "ax.plot(data_an[:, 0], data_an[:, 1], label=\"analytical\")\n", - "ax.plot(data_ek[:, 0], data_ek[:, 1], \"o\", mfc=\"none\", label=\"simulation\")\n", - "ax.set_xlabel(\"x-position\")\n", - "ax.set_ylabel(\"Counter-ion density\")\n", - "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", - "ax.legend(loc=\"best\")\n", - "\n", - "ax = fig1.add_subplot(132)\n", - "ax.plot(data_an[:, 0], data_an[:, 2], label=\"analytical\")\n", - "ax.plot(data_ek[:, 0], data_ek[:, 2], \"o\", mfc=\"none\", label=\"simulation\")\n", - "ax.set_xlabel(\"x-position\")\n", - "ax.set_ylabel(\"Fluid velocity\")\n", - "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", - "ax.legend(loc=\"best\")\n", - "\n", - "ax = fig1.add_subplot(133)\n", - "ax.plot(data_an[:, 0], data_an[:, 3], label=\"analytical\")\n", - "ax.plot(data_ek[:, 0], data_ek[:, 3], \"o\", mfc=\"none\", label=\"simulation\")\n", - "ax.set_xlabel(\"x-position\")\n", - "ax.set_ylabel(\"Fluid shear stress xz\")\n", - "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", - "ax.legend(loc=\"best\")\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" + "We can now set up the electrokinetics algorithm as in the first part of the tutorial, starting with the LB-method." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "## References\n", - "\n", - "[1] F. Capuani, I. Pagonabarraga and D. Frenkel *Discrete solution of the electrokinetic equations* The Journal of Chemical Physics, 2004 \n", - "[2] G. Rempfer *A Lattice based Model for Electrokinetics* Master's thesis, University of Stuttgart, 2013 \n", - "[3] G. Rempfer, G. B. Davies, C. Holm and J. de Graaf *Reducing spurious flow in simulations of electrokinetic phenomena* The Journal of Chemical Physics, 2016 \n", - "[4] G. Rempfer *Lattice-Boltzmann simulations in complex geometries* Bachelor's thesis, University of Stuttgart, Institute for Computational Physics, 2010 \n", - "[5] M. Deserno and C. Holm and S. May, *Fraction of Condensed Counterions around a Charged Rod: Comparison of Poisson-Boltzmann Theory and Computer Simulations* Macromolecules, 2000 \n", - "[6] C. Holm, P. Kékicheff and R. Podgornik *Electrostatic Effects in Soft Matter and Biophysics* Kluwer Academic Publishers, 2001 \n", - "[7] M. Deserno and C. Holm *Cell-model and Poisson-Boltzmann-theory: A brief introduction* Electrostatic Effects in Soft Matter and Biophysics, Kluwer Academic Publishers, 2001 \n", - "[8] J de Graaf., G. Rempfer and C. Holm *Diffusiophoretic Self-Propulsion for Partially Catalytic Spherical Colloids* IEEE T. Nanobiosci., 2014 \n", - "[9] M. Deserno *Counterion condensation for rigid linear polyelectrolytes* Universität Mainz, 2000 \n", - "[10] J. de Graaf, N Boon, M Dijkstra and R. van Roij *Electrostatic interactions between Janus particles* The Journal of Chemical Physics, 2012 \n", - "[11] A. Einstein *Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen* Annalen der Physik, 1905 \n", - "[12] M. von Smoluchowski *Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen* Annalen der Physik, 1906 \n" + "lattice = espressomd.lb.LatticeWalberla(agrid=AGRID, n_ghost_layers=1)" ] }, { @@ -652,12 +468,816 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "viscosity_kinematic = VISCOSITY_DYNAMIC / DENSITY_FLUID\n", + "lbf = espressomd.lb.LBFluidWalberla(lattice=lattice, density=DENSITY_FLUID,\n", + " kinematic_viscosity=viscosity_kinematic,\n", + " tau=TAU, single_precision=SINGLE_PRECISION)\n", + "system.lb = lbf" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since our species are going to carry a charge now, we need to solve the full electrostatic problem. For that, we have to specify an actual solver." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "# Exercise: \n", + "- Set up a Poisson solver for the electrostatic interaction and use it to create an instance of the [EKContainer](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer) \n", + "- Attach the container to the [`system.ekcontainer`](https://espressomd.github.io/doc/espressomd.html#espressomd.system.System.ekcontainer).\n", + "\n", + "# Hint:\n", + "- Use an [EKFFT](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKFFT)-object as the Poisson-solver." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "eksolver = espressomd.electrokinetics.EKFFT(lattice=lattice, permittivity=PERMITTIVITY,\n", + " single_precision=SINGLE_PRECISION)\n", + "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To simulate the system, we will use two different ion species: The counterions are propagated in the fluid. The second species will be used to describe the surface charge on the plates and therefore has to be stationary (i.e. no advection, no diffusion)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ekspecies = espressomd.electrokinetics.EKSpecies(lattice=lattice, density=0.0, kT=KT, diffusion=DIFFUSION_COEFFICIENT, valency=VALENCY, advection=True, friction_coupling=True, ext_efield=EXT_FORCE_DENSITY, single_precision=SINGLE_PRECISION, tau=TAU)\n", + "system.ekcontainer.add(ekspecies)\n", + "\n", + "ekwallcharge = espressomd.electrokinetics.EKSpecies(lattice=lattice, density=0.0, kT=KT, diffusion=0., valency=-VALENCY, advection=False, friction_coupling=False, ext_efield=[0, 0, 0], single_precision=SINGLE_PRECISION, tau=TAU)\n", + "system.ekcontainer.add(ekwallcharge)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we set the initial conditions for the ion densities. The counterions will be initialized with a homogeneous distribution, excluding the cells used as boundaries. The surface charge density is homogeneously distributed in the boundary cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "density_counterions = -2.0 * SURFACE_CHARGE_DENSITY / VALENCY / WIDTH\n", + "ekspecies[padding:-padding, :, :].density = density_counterions\n", + "\n", + "ekspecies[:padding, :, :].density = 0.0\n", + "ekspecies[-padding:, :, :].density = 0.0\n", + "\n", + "ekwallcharge[:padding, :, :].density = -SURFACE_CHARGE_DENSITY / VALENCY / padding\n", + "ekwallcharge[-padding:, :, :].density = -SURFACE_CHARGE_DENSITY / VALENCY / padding" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now have to specify the boundary conditions. For this, we use ESPResSo's`shapes`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "wall_left = espressomd.shapes.Wall(normal=[1, 0, 0], dist=padding)\n", + "wall_right = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(padding + WIDTH))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At both of them we specify no-flux and zero-density boundary conditions for the counterions. Furthermore, we set a no-slip boundary condition for the fluid." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "# Exercise\n", + "At both walls, set\n", + "\n", + "- No-flux boundary conditions for the counterions (``ekspecies``)\n", + "- Zero-density boundary conditions for the counterions\n", + "- No-slip boundary conditions for the fluid (``lbf``)\n", + "\n", + "# Hints\n", + "\n", + "- Use the shapes defined above and ``add_boundary_from_shape`` for [EK species](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKSpecies.add_boundary_from_shape) and [LB fluids](https://espressomd.github.io/doc/espressomd.html#espressomd.lb.LBFluidWalberla.add_boundary_from_shape)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "for wall in (wall_left, wall_right):\n", + " ekspecies.add_boundary_from_shape(shape=wall, value=[0., 0., 0.], boundary_type=espressomd.electrokinetics.FluxBoundary)\n", + " ekspecies.add_boundary_from_shape(shape=wall, value=0.0, boundary_type=espressomd.electrokinetics.DensityBoundary)\n", + " lbf.add_boundary_from_shape(shape=wall, velocity=[0., 0., 0.])\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can finally integrate the system and extract the ion density profile, the fluid velocity profile as well as the pressure-tensor profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "for i in tqdm.trange(80):\n", + " system.integrator.run(RUN_TIME)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mid_y = int(system.box_l[1] / (2 * AGRID))\n", + "mid_z = int(system.box_l[2] / (2 * AGRID))\n", + "density_eof = ekspecies[padding:-padding, mid_y, mid_z].density\n", + "velocity_eof = lbf[padding:-padding, mid_y, mid_z].velocity[:, 1]\n", + "pressure_tensor_eof = lbf[padding:-padding, mid_y, mid_z].pressure_tensor[:, 0, 1]\n", + "\n", + "positions = (np.arange(len(density_eof)) - WIDTH / 2 + 0.5) * AGRID" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For comparison, we calculate the analytic solution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def transcendental_equation(c, distance, kT, sigma, valency, permittivity) -> float:\n", + " elementary_charge = 1.0\n", + " return c * np.tan(valency * elementary_charge * distance / (4 * kT) * c) + sigma / permittivity\n", + "\n", + "solution = scipy.optimize.fsolve(func=transcendental_equation, x0=0.001, args=(WIDTH, KT, SURFACE_CHARGE_DENSITY, VALENCY, PERMITTIVITY))\n", + "\n", + "\n", + "def eof_density(x, c, permittivity, elementary_charge, valency, kT):\n", + " return c**2 * permittivity / (2 * kT) / (np.cos(valency * elementary_charge * c / (2 * kT) * x))**2\n", + "\n", + "def eof_velocity(x, c, permittivity, elementary_charge, valency, kT, ext_field, distance, viscosity):\n", + " return 2 * kT * ext_field * permittivity / (viscosity * elementary_charge * valency) * np.log(np.cos(valency * elementary_charge * c / (2 * kT) * x) / np.cos(valency * elementary_charge * c / (2 * kT) * distance / 2))\n", + "\n", + "def eof_pressure_tensor(x, c, elementary_charge, valency, kT, ext_field, permittivity):\n", + " return permittivity * ext_field * c * np.tan(valency * elementary_charge * c / (2 * kT) * x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "analytic_density_eof = eof_density(x=positions, c=solution, permittivity=PERMITTIVITY, elementary_charge=1.0, valency=VALENCY, kT=KT)\n", + "analytic_velocity_eof = eof_velocity(x=positions, c=solution, permittivity=PERMITTIVITY, elementary_charge=1.0, valency=VALENCY, kT=KT, ext_field=EXT_FORCE_DENSITY[1], distance=WIDTH, viscosity=VISCOSITY_DYNAMIC)\n", + "analytic_pressure_tensor_eof = eof_pressure_tensor(x=positions, c=solution, elementary_charge=1.0, valency=VALENCY, kT=KT, ext_field=EXT_FORCE_DENSITY[1], permittivity=PERMITTIVITY)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1 = plt.figure(figsize=(16, 4.5))\n", + "fig1.suptitle(\"electroosmotic flow\")\n", + "\n", + "ax = fig1.add_subplot(131)\n", + "ax.plot(positions, density_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_density_eof, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"Counter-ion density\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(132)\n", + "ax.plot(positions, velocity_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_velocity_eof, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"Fluid velocity\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(133)\n", + "ax.plot(positions, pressure_tensor_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_pressure_tensor_eof, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"Fluid shear stress xz\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the plots one can see that the analytic solution for the electroosmotic flow matches the simulation very well. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Comparison to pressure-driven flow\n", + "We can compare electroosmotic flow to pressure-driven flow. For this, we turn off the external electric field and enable a constant external force density on the fluid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "EXT_FORCE_DENSITY = [0.0, 0.000004, 0.0]\n", + "\n", + "ekspecies.ext_efield = [0.0, 0.0, 0.0]\n", + "lbf.ext_force_density = EXT_FORCE_DENSITY" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "for i in tqdm.trange(70):\n", + " system.integrator.run(RUN_TIME)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "density_pressure = ekspecies[padding:-padding, mid_y, mid_z].density\n", + "velocity_pressure = lbf[padding:-padding, mid_y, mid_z].velocity[:, 1]\n", + "pressure_tensor_pressure = lbf[padding:-padding, mid_y, mid_z].pressure_tensor[:, 0, 1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The analytic solution for pressure-driven flow between two infinite parallel plates is known as the Poiseuille flow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def pressure_velocity(x, distance, ext_field, viscosity):\n", + " return ext_field / (2 * viscosity) * (distance**2 / 4 - x**2)\n", + "\n", + "def pressure_pressure_tensor(x, ext_field):\n", + " return ext_field * x" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "analytic_velocity_pressure = pressure_velocity(x=positions, distance=WIDTH, ext_field=EXT_FORCE_DENSITY[1], viscosity=VISCOSITY_DYNAMIC)\n", + "analytic_pressure_tensor_pressure = pressure_pressure_tensor(x=positions, ext_field=EXT_FORCE_DENSITY[1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1 = plt.figure(figsize=(16, 4.5))\n", + "fig1.suptitle(\"pressure-driven flow\")\n", + "\n", + "ax = fig1.add_subplot(131)\n", + "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_density_eof, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"counter-ion density\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(132)\n", + "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_velocity_pressure, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"fluid velocity\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(133)\n", + "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"simulation\")\n", + "ax.plot(positions, analytic_pressure_tensor_pressure, label=\"analytic\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"fluid shear stress xz\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As one can again see, the body force on the fluid did non alter the ion-density profile.\n", + "However, because the force now applies homogeneously on the whole fluid, the flow profile looks parabolic." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To see the difference between the two types of flows, we plot the simulation data together in one plot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "fig1 = plt.figure(figsize=(16, 4.5))\n", + "fig1.suptitle(\"electroosmotic vs. pressure-driven flow comparison\")\n", + "\n", + "ax = fig1.add_subplot(131)\n", + "ax.plot(positions, density_eof, \"o\", mfc=\"none\", ms=4, markevery=0.015, label=\"eof\")\n", + "ax.plot(positions, density_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"counter-ion density\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(132)\n", + "ax.plot(positions, velocity_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"eof\")\n", + "ax.plot(positions, velocity_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"fluid velocity\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "ax = fig1.add_subplot(133)\n", + "ax.plot(positions, pressure_tensor_eof, \"o\", mfc=\"none\", markevery=0.015, label=\"eof\")\n", + "ax.plot(positions, pressure_tensor_pressure, \"o\", mfc=\"none\", markevery=0.015, label=\"pressure\")\n", + "ax.set_xlabel(\"x-position\")\n", + "ax.set_ylabel(\"fluid shear stress xz\")\n", + "ax.ticklabel_format(axis=\"y\", style=\"scientific\", scilimits=(0, 0))\n", + "ax.legend(loc=\"best\")\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looking at the fluid velocity plot, one can see that the electroosmotic flow profile flattens significantly faster towards the center of the channel when compared to the pressure driven flow. The reason for this is the accumulation of the counterion-density towards the oppositely charged plates. Here, the driving electric field causes the highest force on the fluid, which decays towards the center of the channel. In contrast, the Poiseuille-flow is driven by a constant, uniform driving force." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4. Reaction in turbulent flow" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To showcase the reaction feature of our electrokinetics algorithm, we simulate a simple reaction in complex flow.\n", + "For this, we choose a geometry of rigid cylinders.\n", + "At large flow velocities, a [Kármán vortex street](https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street), i.e., a repeating pattern of swirling vorticies behind the obstacle, develops.\n", + " \n", + "To this flow, we will add several species undergoing advection-diffusion, which is dominated by the downstream fluid flow in the channel.\n", + "The reaction will be included as a bulk-reaction, which means that the reaction can happen anywhere, the only requirement is that both species are present in the same lattice cell. When the reaction occurs, parts of the reactant species will turn into the product. How much of the species will transform within each timestep is determined by the respective reaction rate and the overall structure of the reaction.\n", + "\n", + "We start again by resetting the system and defining parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "system.ekcontainer = None\n", + "system.lb = None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "BOX_L = [80, 32, 1] \n", + "AGRID = 1.0\n", + "DIFFUSION_COEFFICIENT = 0.01\n", + "TAU = 0.03\n", + "EXT_FORCE_DENSITY = [0.6, 0, 0]\n", + "OBSTACLE_RADIUS = 6\n", + "\n", + "DENSITY_FLUID = 0.5\n", + "VISCOSITY_KINEMATIC = 2.0\n", + "KT = 1.0\n", + "\n", + "TOTAL_FRAMES = 50\n", + "\n", + "system.time_step = TAU\n", + "system.cell_system.skin = 0.4\n", + "system.box_l = BOX_L" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lattice = espressomd.lb.LatticeWalberla(n_ghost_layers=1, agrid=AGRID)\n", + "lbf = espressomd.lb.LBFluidWalberla(\n", + " lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=VISCOSITY_KINEMATIC,\n", + " tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=KT, seed=42)\n", + "system.lb = lbf\n", + "system.thermostat.set_lb(LB_fluid=lbf, seed=42)", + "\n", + "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n", + "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can focus on the reactions. In this tutorial we choose the simple case of $A + B \\rightarrow C$, which means that equal parts of the educt species $A$ and $B$ can turn into the product species $C$.\n", + "ESPResSo distinguishes between educts and products by the sign of their respective stoechiometric coefficients, where educts have negative coefficients and products positive coefficients. Intuitively this can be understood that when a reaction happens, the density of the educts will decrease, hence the stoechiometric coefficient is negative. \n", + "\n", + "The reaction rate constant $r$ is the rate at which the reaction happens. The order $O_i$ for a species $i$ specifies to which order the reaction depends on the density of that species. Positive orders mean that the reaction is faster the more density of this species is present, for negative orders the reaction slows down with higher density.\n", + "In general, this process can be written as $\\Gamma = r \\left[ A \\right]^{O_A} \\left[ B \\right]^{O_B} \\left[ C \\right]^{O_C}$, where $\\Gamma$ is known as the reaction rate. This is sometimes also called the [rate equation](https://en.wikipedia.org/wiki/Rate_equation).\n", + "\n", + "For our specific simulation this means that all stoechiometric coefficients are $-1$ for the educts and $+1$ for the product. We choose the order of the educts as $1$ and the order of the product as $0$. This means that the more amount of both educts is present, the more will react and the amount of product present won't have an influence." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "REACTION_RATE_CONSTANT = 2.5\n", + "EDUCT_COEFFS = [-1, -1]\n", + "EDUCT_ORDERS = [1,1]\n", + "PRODUCT_COEFFS = [1]\n", + "PRODUCT_ORDERS = [0]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create each involved species and directly specify their boundary-conditions for the domain-boundaries. We set the initial density of the species to 0 and also add Dirichlet boundary conditions of zero density at both the inlet and the outlet of the system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "educt_species = []\n", + "product_species = []\n", + "reactants = []\n", + "for coeff, order in zip(EDUCT_COEFFS, EDUCT_ORDERS):\n", + " species = espressomd.electrokinetics.EKSpecies(\n", + " lattice=lattice, density=0.0, kT=KT,\n", + " diffusion=DIFFUSION_COEFFICIENT, valency=0.0,\n", + " advection=True, friction_coupling=True,\n", + " ext_efield=[0., 0., 0.], tau=TAU)\n", + " system.ekcontainer.add(species)\n", + " reactants.append(\n", + " espressomd.electrokinetics.EKReactant(\n", + " ekspecies=species,\n", + " stoech_coeff=coeff,\n", + " order=order))\n", + " educt_species.append(species)\n", + " species[0,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n", + " species[-1,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n", + "\n", + "for coeff, order in zip(PRODUCT_COEFFS, PRODUCT_ORDERS):\n", + " species = espressomd.electrokinetics.EKSpecies(\n", + " lattice=lattice, density=0.0, diffusion=DIFFUSION_COEFFICIENT,\n", + " kT=KT, valency=0.0, advection=True, friction_coupling=True,\n", + " ext_efield=[0., 0., 0.], tau=TAU)\n", + " system.ekcontainer.add(species)\n", + " reactants.append(\n", + " espressomd.electrokinetics.EKReactant(\n", + " ekspecies=species,\n", + " stoech_coeff=coeff,\n", + " order=order))\n", + " product_species.append(species)\n", + " species[0,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)\n", + " species[-1,:,:].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "# Exercise:\n", + "- Create an instance of [`EKBulkReaction`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKBulkReaction) using the previously created `reactants` and activate the reaction by adding it to [`system.ekcontainer.reactions`](https://espressomd.github.io/doc/espressomd.html#espressomd.electrokinetics.EKContainer.reactions).\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "reaction = espressomd.electrokinetics.EKBulkReaction(\n", + " reactants=reactants, coefficient=REACTION_RATE_CONSTANT, lattice=lattice, tau=TAU)\n", + "\n", + "system.ekcontainer.reactions.add(reaction)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The next thing to add to the system is the cylindrical obstacles, which act as the boundaries for the Kármán vortices to form. These are placed close to the inlet of the system and also act as impenetrable boundaries for the species.\n", + "Since ESPResSo uses periodic boundary conditions, we need to add a total of three cylinders to the system, which will form two complete cylinders in the periodic system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cylinder_centers = [\n", + " [BOX_L[0] // 10, 0, 1],\n", + " [BOX_L[0] // 10, BOX_L[1] // 2, 1],\n", + " [BOX_L[0] // 10, BOX_L[1], 1], \n", + "]\n", + "\n", + "shape_cylinder = []\n", + "for cylinder_center in cylinder_centers:\n", + " shape_cylinder.append(espressomd.shapes.Cylinder(\n", + " center=cylinder_center,\n", + " axis=[0, 0, 1],\n", + " length=BOX_L[2],\n", + " radius=OBSTACLE_RADIUS,\n", + " ))\n", + "\n", + "for shape in shape_cylinder:\n", + " lbf.add_boundary_from_shape(shape)\n", + " for spec in (*educt_species, *product_species):\n", + " spec.add_boundary_from_shape(shape, value=[0,0,0], boundary_type=espressomd.electrokinetics.FluxBoundary)\n", + " spec.add_boundary_from_shape(shape, value=0., boundary_type=espressomd.electrokinetics.DensityBoundary)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Up to this point there is no species present anywhere in the system and also no way for it to enter the system. Since the reaction is irreversible in our setup, we need to introduce some density of both the educt species to the system.\n", + "For that we set two additional Dirichlet boundary conditions (sources) in the domain, where we fix the species' density to a constant, non-zero value. The sources are placed some distance apart such that the reaction happens further downstream when diffusion mixes the two species." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source_boundary = espressomd.electrokinetics.DensityBoundary(10.0)\n", + "source_x_pos = BOX_L[1] // 20\n", + "\n", + "educt_species[0][source_x_pos,BOX_L[1]//4-1:BOX_L[1]//4+1,:].density_boundary = source_boundary\n", + "educt_species[1][source_x_pos,3*(BOX_L[1]//4)-1:3*(BOX_L[1]//4)+1,:].density_boundary = source_boundary" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With this, the system is now finally complete and we can start the integration. To see the system evolve, we will render a movie from the timeseries of the system. For that we have to setup some helper functions for the plotting, which are beyond the scope of this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "VIDEO_TAG = \"\"\"\"\"\"\n", + "\n", + "# set ignore 'divide' and 'invalid' errors\n", + "# these occur when plotting the flowfield containing a zero velocity\n", + "np.seterr(divide='ignore', invalid='ignore')\n", + "\n", + "def anim_to_html(anim):\n", + " if not hasattr(anim, '_encoded_video'):\n", + " with tempfile.NamedTemporaryFile(suffix='.mp4') as f:\n", + " anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])\n", + " with open(f.name, \"rb\") as g:\n", + " video = g.read()\n", + " anim._encoded_video = base64.b64encode(video).decode('ascii')\n", + " plt.close(anim._fig)\n", + " return VIDEO_TAG.format(anim._encoded_video)\n", + "\n", + "animation.Animation._repr_html_ = anim_to_html" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "get_colormap = mpl.colormaps.get_cmap if hasattr(mpl.colormaps, \"get_cmap\") else mpl.cm.get_cmap\n", + "box_width = lattice.shape[1]\n", + "box_height = lattice.shape[0]\n", + "\n", + "boundary_mask = lbf[:, :, 0].boundary != None\n", + "\n", + "cmap = get_colormap(\"viridis\").copy()\n", + "cmap.set_bad(color=\"gray\")\n", + "cmap_quiver = get_colormap(\"binary\").copy()\n", + "cmap_quiver.set_bad(color=\"gray\")\n", + "\n", + "# setup figure and prepare axes\n", + "fig = plt.figure(figsize=(9.8, 5.5))\n", + "imshow_kwargs = {\"origin\": \"upper\", \"extent\": (0, BOX_L[1], BOX_L[0], 0)}\n", + "gs = fig.add_gridspec(1, 4, wspace=0.1)\n", + "ax1 = plt.subplot(gs[0])\n", + "ax2 = plt.subplot(gs[1], sharey=ax1)\n", + "ax3 = plt.subplot(gs[2], sharey=ax1)\n", + "ax4 = plt.subplot(gs[3], sharey=ax1)\n", + "ax1.set_yticks(np.arange(0, 80 + 1, 16))\n", + "ax1.set_xticks(np.arange(0, 32 + 1, 16))\n", + "ax2.set_xticks(np.arange(0, 32 + 1, 16))\n", + "ax3.set_xticks(np.arange(0, 32 + 1, 16))\n", + "ax4.set_xticks(np.arange(0, 32 + 1, 16))\n", + "\n", + "# set the background color for the quiver plot\n", + "bg_colors = np.copy(boundary_mask).astype(float)\n", + "bg_colors[boundary_mask] = np.NaN\n", + "ax4.imshow(bg_colors, cmap=cmap_quiver, **imshow_kwargs)\n", + "\n", + "for ax, title in zip(\n", + " [ax1, ax2, ax3, ax4],\n", + " [\"educt 1\", \"educt 2\", \"product\", \"fluid velocity\"]\n", + "):\n", + " ax.set_title(title)\n", + " ax.set_xlim((0, box_width))\n", + " ax.set_ylim((0, box_height))\n", + "\n", + "# create meshgrid for quiver plot subsampled by a factor 2\n", + "xs = np.arange(0, box_width, 2)\n", + "ys = np.arange(0, box_height, 2)\n", + "X, Y = np.meshgrid(xs, ys)\n", + "\n", + "flowfield = lbf[:, :, 0].velocity[::2, ::2, :]\n", + "quiver = ax4.quiver(X + 1., Y + 1., flowfield[..., 1], flowfield[..., 0], scale=100)\n", + "fig.subplots_adjust(left=0.025, bottom=0.075, right=0.975, top=0.925, wspace=0.0125)\n", + "\n", + "progress_bar = tqdm.tqdm(total=TOTAL_FRAMES)\n", + "\n", + "def draw_frame(frame):\n", + " system.integrator.run(50)\n", + " \n", + " flowfield = np.copy(lbf[:, :, 0].velocity)\n", + " e1 = np.copy(educt_species[0][:, :, 0].density)\n", + " e2 = np.copy(educt_species[1][:, :, 0].density)\n", + " p = np.copy(product_species[0][:, :, 0].density)\n", + " \n", + " # apply the mask for the boundary\n", + " e1[boundary_mask] = np.NaN\n", + " e2[boundary_mask] = np.NaN\n", + " p[boundary_mask] = np.NaN\n", + " flowfield[boundary_mask] = np.NaN\n", + "\n", + " ax1.imshow(e1, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n", + " ax2.imshow(e2, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n", + " ax3.imshow(p, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)\n", + " quiver.set_UVC((flowfield[::2, ::2, 1] + flowfield[1::2, 1::2, 1]) / 2.,\n", + " (flowfield[::2, ::2, 0] + flowfield[1::2, 1::2, 0]) / 2.)\n", + " \n", + " progress_bar.update()\n", + "\n", + "\n", + "animation.FuncAnimation(fig, draw_frame, frames=range(TOTAL_FRAMES), interval=300)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looking at the movie of the species densities one can see that the fluid flow advects the educt species from their source locations past the cylinders into the system. Here, they start to mix and react, such that the product forms.\n", + "The vortex street behind the obstacles enhances mixing through fluid turbulence. The density of the product then increases towards the outflow-location of the channel, where it is deleted because of our zero-density boundary condition." + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -671,7 +1291,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/doc/tutorials/electrokinetics/scripts/eof_analytical.py b/doc/tutorials/electrokinetics/scripts/eof_analytical.py deleted file mode 100644 index 83122d225a8..00000000000 --- a/doc/tutorials/electrokinetics/scripts/eof_analytical.py +++ /dev/null @@ -1,146 +0,0 @@ -# Copyright (C) 2010-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 . -# Set the slit pore geometry the width is the non-periodic part of the geometry -# the padding is used to ensure that there is no field outside the slit -import numpy as np - -box_x = 6 -box_y = 6 -width = 50 - -padding = 6 -box_z = width + 2 * padding - -# Set the electrokinetic parameters - -agrid = 1.0 -temperature = 1.0 -bjerrum_length = 0.7095 -valency = 1.0 -viscosity_dynamic = 79.53 -density_water = 26.15 -sigma = -0.05 -force = 0.1 - -viscosity_kinematic = viscosity_dynamic / density_water -density_counterions = -2.0 * sigma / width - -# Calculate the inverse length xi, which is a combination of various -# constants (xi = zeC/2kBT), with C a constant that needs to be -# solved for, or equivalently, xi needs to be solved for - -# root finding function - - -def solve(xi=None, d=None, bjerrum_length=None, sigma=None, valency=None): - el_char = 1.0 - return xi * np.tan(xi * d / 2.0) + 2.0 * np.pi * \ - bjerrum_length * sigma / (valency * el_char) - - -size = np.pi / (2.0 * width) - -pnt0 = 0.0 -pntm = pnt0 + size -pnt1 = pnt0 + 1.9 * size - -# the bisection scheme -tol = 1e-8 - -while size > tol: - val0 = solve(xi=pnt0, d=width, - bjerrum_length=bjerrum_length, sigma=sigma, valency=valency) - val1 = solve(xi=pnt1, d=width, - bjerrum_length=bjerrum_length, sigma=sigma, valency=valency) - valm = solve(xi=pntm, d=width, - bjerrum_length=bjerrum_length, sigma=sigma, valency=valency) - - if (val0 < 0.0) and (val1 > 0.0): - if valm < 0.0: - pnt0 = pntm - size = size / 2.0 - pntm = pnt0 + size - else: - pnt1 = pntm - size = size / 2.0 - pntm = pnt1 - size - elif (val0 > 0.0) and (val1 < 0.0): - if valm < 0.0: - pnt1 = pntm - size = size / 2.0 - pntm = pnt0 - size - else: - pnt0 = pntm - size = size / 2.0 - pntm = pnt0 - size - else: - raise Exception( - "Bisection method fails:\nTuning of domain boundaries may be required.") - - -# obtain the desired xi value -xi = pntm - -# function to calculate the density - - -def density(x=None, xi=None, bjerrum_length=None): - return (xi**2) / (2.0 * np.pi * bjerrum_length * np.cos(xi * x)**2) - -# function to calculate the velocity - - -def velocity(x=None, xi=None, d=None, bjerrum_length=None, force=None, - viscosity_kinematic=None, density_water=None): - return force * np.log(np.cos(xi * x) / np.cos(xi * d / 2.0)) / \ - (2.0 * np.pi * bjerrum_length * viscosity_kinematic * density_water) - -# function to calculate the nonzero component of the pressure tensor - - -def pressure_tensor_offdiagonal(x=None, xi=None, bjerrum_length=None, - force=None): - return force * xi * np.tan(xi * x) / (2.0 * np.pi * bjerrum_length) - - -position_list = [] -density_list = [] -velocity_list = [] -pressure_xy_list = [] - -for i in range(int(box_z / agrid)): - if (i * agrid >= padding) and (i * agrid < box_z - padding): - position = i * agrid - padding - width / 2.0 + agrid / 2.0 - position_list.append(position) - - # density - density_list.append(density(x=position, xi=xi, - bjerrum_length=bjerrum_length)) - - # velocity - velocity_list.append(velocity(x=position, xi=xi, - d=width, bjerrum_length=bjerrum_length, - force=force, viscosity_kinematic=viscosity_kinematic, - density_water=density_water)) - # xz component pressure tensor - pressure_xy_list.append(pressure_tensor_offdiagonal(x=position, xi=xi, - bjerrum_length=bjerrum_length, force=force)) - -np.savetxt( - "eof_analytical.dat", np.column_stack( - (position_list, density_list, velocity_list, pressure_xy_list)), - header="#position calculated_density calculated_velocity calculated_pressure_xy") diff --git a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb index bbc442c04d9..7e1c8587920 100644 --- a/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb +++ b/doc/tutorials/lattice_boltzmann/lattice_boltzmann_sedimentation.ipynb @@ -79,9 +79,9 @@ "espressomd.assert_features([\"LENNARD_JONES\", \"WALBERLA\"])\n", "\n", "# imports for data handling, plotting, and progress bar\n", + "import tqdm\n", "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import tqdm" + "import matplotlib.pyplot as plt" ] }, { diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 252296760c4..f7fd6f5abbf 100644 --- a/src/core/Observable_stat.cpp +++ b/src/core/Observable_stat.cpp @@ -58,8 +58,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size) constexpr std::size_t n_ext_fields = 1; // reduction over all fields constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions - auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb + - n_dipolar + n_vs + n_ext_fields; + auto const n_elements = n_kinetic + n_bonded + 2ul * n_non_bonded + + n_coulomb + n_dipolar + n_vs + n_ext_fields; m_data = std::vector(m_chunk_size * n_elements); // spans for the different contributions @@ -78,8 +78,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size) } Utils::Span -Observable_stat::non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const { +Observable_stat::get_non_bonded_contribution(Utils::Span base_pointer, + int type1, int type2) const { auto const offset = static_cast(Utils::upper_triangular( std::min(type1, type2), std::max(type1, type2), max_seen_particle_type)); return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size}; diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index b34157004bc..b84bcfdda00 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -38,8 +38,9 @@ class Observable_stat { std::size_t m_chunk_size; /** Get contribution from a non-bonded interaction */ - Utils::Span non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const; + Utils::Span + get_non_bonded_contribution(Utils::Span base_pointer, int type1, + int type2) const; public: explicit Observable_stat(std::size_t chunk_size); @@ -91,29 +92,30 @@ class Observable_stat { return {bonded.data() + offset, m_chunk_size}; } - void add_non_bonded_contribution(int type1, int type2, + void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, Utils::Span data) { assert(data.size() == m_chunk_size); - auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter; - auto const dest = non_bonded_contribution(source, type1, type2); + auto const span = (molid1 == molid2) ? non_bonded_intra : non_bonded_inter; + auto const dest = get_non_bonded_contribution(span, type1, type2); boost::transform(dest, data, dest.begin(), std::plus<>{}); } - void add_non_bonded_contribution(int type1, int type2, double data) { - add_non_bonded_contribution(type1, type2, {&data, 1}); + void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, + double data) { + add_non_bonded_contribution(type1, type2, molid1, molid2, {&data, 1}); } /** Get contribution from a non-bonded intramolecular interaction */ Utils::Span non_bonded_intra_contribution(int type1, int type2) const { - return non_bonded_contribution(non_bonded_intra, type1, type2); + return get_non_bonded_contribution(non_bonded_intra, type1, type2); } /** Get contribution from a non-bonded intermolecular interaction */ Utils::Span non_bonded_inter_contribution(int type1, int type2) const { - return non_bonded_contribution(non_bonded_inter, type1, type2); + return get_non_bonded_contribution(non_bonded_inter, type1, type2); } /** MPI reduction. */ diff --git a/src/core/accumulators.cpp b/src/core/accumulators.cpp index 6f707afd0a3..bc2b8306033 100644 --- a/src/core/accumulators.cpp +++ b/src/core/accumulators.cpp @@ -60,24 +60,33 @@ int auto_update_next_update() { }); } +namespace detail { +struct MatchPredicate { + AccumulatorBase const *m_acc; + template bool operator()(T const &a) const { + return a.acc == m_acc; + } +}; +} // namespace detail + void auto_update_add(AccumulatorBase *acc) { - assert(acc); - assert(std::find_if(auto_update_accumulators.begin(), - auto_update_accumulators.end(), [acc](auto const &a) { - return a.acc == acc; - }) == auto_update_accumulators.end()); + assert(not auto_update_contains(acc)); auto_update_accumulators.emplace_back(acc); } + void auto_update_remove(AccumulatorBase *acc) { - assert(std::find_if(auto_update_accumulators.begin(), - auto_update_accumulators.end(), [acc](auto const &a) { - return a.acc == acc; - }) != auto_update_accumulators.end()); + assert(auto_update_contains(acc)); + auto const beg = auto_update_accumulators.begin(); + auto const end = auto_update_accumulators.end(); auto_update_accumulators.erase( - boost::remove_if( - auto_update_accumulators, - [acc](AutoUpdateAccumulator const &au) { return au.acc == acc; }), - auto_update_accumulators.end()); + std::remove_if(beg, end, detail::MatchPredicate{acc}), end); +} + +bool auto_update_contains(AccumulatorBase const *acc) noexcept { + assert(acc); + auto const beg = auto_update_accumulators.begin(); + auto const end = auto_update_accumulators.end(); + return std::find_if(beg, end, detail::MatchPredicate{acc}) != end; } } // namespace Accumulators diff --git a/src/core/accumulators.hpp b/src/core/accumulators.hpp index 68a9da09a05..bc8199b05d6 100644 --- a/src/core/accumulators.hpp +++ b/src/core/accumulators.hpp @@ -31,6 +31,7 @@ namespace Accumulators { */ void auto_update(boost::mpi::communicator const &comm, int steps); int auto_update_next_update(); +bool auto_update_contains(AccumulatorBase const *) noexcept; void auto_update_add(AccumulatorBase *); void auto_update_remove(AccumulatorBase *); diff --git a/src/core/constraints/Constraints.hpp b/src/core/constraints/Constraints.hpp index a8c9b1ecd22..edc46cc20cc 100644 --- a/src/core/constraints/Constraints.hpp +++ b/src/core/constraints/Constraints.hpp @@ -50,23 +50,21 @@ template class Constraints { container_type m_constraints; public: + bool contains(std::shared_ptr const &constraint) const noexcept { + return std::find(begin(), end(), constraint) != end(); + } void add(std::shared_ptr const &constraint) { auto const &box_geo = *System::get_system().box_geo; if (not constraint->fits_in_box(box_geo.length())) { throw std::runtime_error("Constraint not compatible with box size."); } - assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) == - m_constraints.end()); - + assert(not contains(constraint)); m_constraints.emplace_back(constraint); on_constraint_change(); } void remove(std::shared_ptr const &constraint) { - assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) != - m_constraints.end()); - m_constraints.erase( - std::remove(m_constraints.begin(), m_constraints.end(), constraint), - m_constraints.end()); + assert(contains(constraint)); + m_constraints.erase(std::remove(begin(), end(), constraint), end()); on_constraint_change(); } diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 7043bbb31f2..450544ef39e 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -164,7 +164,9 @@ void ShapeBasedConstraint::add_energy(const Particle &p, } } // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) - if (part_rep.type() >= 0) - obs_energy.add_non_bonded_contribution(p.type(), part_rep.type(), energy); + if (part_rep.type() >= 0) { + obs_energy.add_non_bonded_contribution( + p.type(), part_rep.type(), p.mol_id(), part_rep.mol_id(), energy); + } } } // namespace Constraints diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index c693051f83d..40cbe84dbc1 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -184,7 +184,7 @@ inline void add_non_bonded_pair_energy( if (do_nonbonded(p1, p2)) #endif obs_energy.add_non_bonded_contribution( - p1.type(), p2.type(), + p1.type(), p2.type(), p1.mol_id(), p2.mol_id(), calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist, coulomb_kernel)); diff --git a/src/core/event.cpp b/src/core/event.cpp index 37a1e3bad28..7fa4acaf646 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -157,6 +157,12 @@ void on_particle_charge_change() { partCfg().invalidate(); } +void on_particle_local_change() { + cells_update_ghosts(global_ghost_flags()); + + recalc_forces = true; +} + void on_particle_change() { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; diff --git a/src/core/event.hpp b/src/core/event.hpp index 3fda162f9ad..bec60711a08 100644 --- a/src/core/event.hpp +++ b/src/core/event.hpp @@ -66,6 +66,9 @@ void on_particle_change(); /** called every time the charge of a particle has changed. */ void on_particle_charge_change(); +/** called every time that local properties of a particle have changed. */ +void on_particle_local_change(); + /** called every time the Coulomb parameters are changed. all Coulomb methods have a short range part, aka near field diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index a761cdd7e2f..319872fc633 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -71,10 +71,8 @@ inline void add_non_bonded_pair_virials( .f + calc_non_central_force(p1, p2, ia_params, d, dist).f; auto const stress = Utils::tensor_product(d, force); - - auto const type1 = p1.mol_id(); - auto const type2 = p2.mol_id(); - obs_pressure.add_non_bonded_contribution(type1, type2, flatten(stress)); + obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(), + p2.mol_id(), flatten(stress)); } #ifdef ELECTROSTATICS diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 7c5a8194d26..f1f2e6c7c01 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -157,6 +157,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, auto const random_index = i_random(number_of_particles_with_type(type)); return get_random_p_id(type, random_index); }; + auto only_local_changes = true; for (int i = 0; i < std::min(n_product_types, n_reactant_types); i++) { auto const n_product_coef = reaction.product_coefficients[i]; auto const n_reactant_coef = reaction.reactant_coefficients[i]; @@ -164,6 +165,11 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, // particles of reactant_types(i) to product_types(i) auto const old_type = reaction.reactant_types[i]; auto const new_type = reaction.product_types[i]; +#ifdef ELECTROSTATICS + if (charges_of_types.at(new_type) != charges_of_types.at(old_type)) { + only_local_changes = false; + } +#endif for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) { auto const p_id = get_random_p_id_of_type(old_type); on_particle_type_change(p_id, old_type, new_type); @@ -175,7 +181,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, } bookkeeping.changed.emplace_back(p_id, old_type); } - on_particle_change(); // create product_coefficients(i)-reactant_coefficients(i) many product // particles iff product_coefficients(i)-reactant_coefficients(i)>0, // iff product_coefficients(i)-reactant_coefficients(i)<0, hide this number @@ -188,7 +193,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); bookkeeping.created.emplace_back(p_id); } - on_particle_change(); + only_local_changes = false; } else if (delta_n < 0) { auto const type = reaction.reactant_types[i]; for (int j = 0; j < -delta_n; j++) { @@ -197,7 +202,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); hide_particle(p_id, type); } - on_particle_change(); + only_local_changes = false; } } // create or hide particles of types with noncorresponding replacement types @@ -212,7 +217,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); hide_particle(p_id, type); } - on_particle_change(); } else { // create additional product_types particles auto const type = reaction.product_types[i]; @@ -221,9 +225,14 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); bookkeeping.created.emplace_back(p_id); } - on_particle_change(); } } + // determine which fine-grained event to trigger + if (n_product_types == n_reactant_types and only_local_changes) { + on_particle_local_change(); + } else { + on_particle_change(); + } } std::unordered_map diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 4b5b711dae7..2739b0cceb8 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -114,6 +114,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { BOOST_REQUIRE_GE(get_particle_node_parallel(pid2), rank ? -1 : 1); BOOST_REQUIRE_GE(get_particle_node_parallel(pid3), rank ? -1 : 1); } + set_particle_property(pid1, &Particle::mol_id, type_a); + set_particle_property(pid2, &Particle::mol_id, type_b); + set_particle_property(pid3, &Particle::mol_id, type_b); auto const reset_particle_positions = [&start_positions]() { for (auto const &kv : start_positions) { @@ -226,10 +229,10 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { on_non_bonded_ia_change(); // matrix indices and reference energy value - auto const max_type = type_b + 1; - auto const n_pairs = Utils::upper_triangular(type_b, type_b, max_type) + 1; - auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, max_type); - auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, max_type); + auto const size = std::max(type_a, type_b) + 1; + auto const n_pairs = Utils::upper_triangular(type_b, type_b, size) + 1; + auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, size); + auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, size); auto const frac6 = Utils::int_pow<6>(sig / r_off); auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift); @@ -237,6 +240,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { auto const obs_energy = calculate_energy(); if (rank == 0) { for (int i = 0; i < n_pairs; ++i) { + // particles were set up with type == mol_id auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10); diff --git a/src/python/espressomd/electrokinetics.py b/src/python/espressomd/electrokinetics.py index 1fb2a59ae9a..336cea4780d 100644 --- a/src/python/espressomd/electrokinetics.py +++ b/src/python/espressomd/electrokinetics.py @@ -31,7 +31,16 @@ class EKFFT(ScriptInterfaceHelper): """ A FFT-based Poisson solver. + Intrinsically assumes periodic boundary conditions. + Parameters + ---------- + lattice : :obj:`espressomd.lb.LatticeWalberla ` + Lattice object. + permittivity : :obj:`float` + permittivity of the fluid :math:`\\epsilon_0 \\epsilon_{\\mathrm{r}}`. + single_precision : :obj:`bool`, optional + Use single-precision floating-point arithmetic. """ _so_name = "walberla::EKFFT" _so_features = ("WALBERLA_FFT",) @@ -44,6 +53,12 @@ class EKNone(ScriptInterfaceHelper): The default Poisson solver. Imposes a null electrostatic potential everywhere. + Parameters + ---------- + lattice : :obj:`espressomd.lb.LatticeWalberla ` + Lattice object. + single_precision : :obj:`bool`, optional + Use single-precision floating-point arithmetic. """ _so_name = "walberla::EKNone" _so_features = ("WALBERLA",) @@ -576,16 +591,61 @@ def __repr__(self): @script_interface_register class EKReactant(ScriptInterfaceHelper): + """ + Reactant-object which specifies the contribution of a species to a reaction. + + Parameters + ---------- + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` + EK species to react + stoech_coeff: :obj:`float` + Stoechiometric coefficient of this species in the reaction. + Products have positive coefficients whereas educts have negative ones. + order: :obj:`float` + Partial-order of this species in the reaction. + + """ _so_name = "walberla::EKReactant" _so_creation_policy = "GLOBAL" class EKBulkReaction(ScriptInterfaceHelper): + """ + Reaction type that is applied everywhere in the domain. + + Parameters + ---------- + lattice : :obj:`espressomd.electrokinetics.LatticeWalberla ` + Lattice object. + tau : :obj:`float` + EK time step, must be an integer multiple of the MD time step. + coefficient : :obj:`float` + Reaction rate constant of the reaction. + reactants: array_like of :obj:`~espressomd.electrokinetics.EKReactant` + Reactants that participate this reaction. + + """ _so_name = "walberla::EKBulkReaction" _so_creation_policy = "GLOBAL" class EKIndexedReaction(ScriptInterfaceHelper): + """ + Reaction type that is applied only on specific cells in the domain. + Can be used to model surface-reactions. + + Parameters + ---------- + lattice : :obj:`espressomd.electrokinetics.LatticeWalberla ` + Lattice object. + tau : :obj:`float` + EK time step, must be an integer multiple of the MD time step. + coefficient : :obj:`float` + Reaction rate constant of the reaction. + reactants: array_like of :obj:`~espressomd.electrokinetics.EKReactant` + Reactants that participate this reaction. + + """ _so_name = "walberla::EKIndexedReaction" _so_creation_policy = "GLOBAL" @@ -648,21 +708,64 @@ def __setitem__(self, key, values): @script_interface_register class EKReactions(ScriptObjectList): + """ + Container object holding all EK-reactions that are considered. + + Methods + ------- + clear() + Remove all reactions. + + """ _so_name = "walberla::EKReactions" _so_creation_policy = "GLOBAL" _so_bind_methods = ("clear",) def add(self, reaction): + """ + Add a reaction to the container. + + Parameters + ---------- + reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction` + Reaction to be added. + + """ if not isinstance(reaction, (EKBulkReaction, EKIndexedReaction)): raise TypeError("reaction object is not of correct type.") self.call_method("add", object=reaction) def remove(self, reaction): + """ + Remove a reaction from the container. + + Parameters + ---------- + reaction : :obj:`~espressomd.electrokinetics.EKBulkReaction` or :obj:`~espressomd.electrokinetics.EKIndexedReaction` + Reaction to be removed. + + """ self.call_method("remove", object=reaction) @script_interface_register class EKContainer(ScriptObjectList): + """ + Container object holding the :obj:`~espressomd.electrokinetics.EKSpecies`. + + Parameters + ---------- + tau : :obj:`float` + EK time step, must be an integer multiple of the MD time step. + solver : :obj:`~espressomd.electrokinetics.EKNone` or :obj:`~espressomd.electrokinetics.EKFFT` + Solver defining the treatment of the electrostatic Poisson-equation. + + Methods + ------- + clear() + Removes all species. + + """ _so_name = "walberla::EKContainer" _so_creation_policy = "GLOBAL" _so_features = ("WALBERLA",) @@ -675,8 +778,36 @@ def __init__(self, *args, **kwargs): else: super().__init__(**kwargs) + @property + def reactions(self): + """ + Returns + ------- + Reactions-container of the current reactions (:obj:`~espressomd.electrokinetics.EKReactions`). + + """ + return self._getter("reactions") + def add(self, ekspecies): + """ + Add an :obj:`~espressomd.electrokinetics.EKSpecies` to the container. + + Parameters + ---------- + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` + Species to be added. + + """ self.call_method("add", object=ekspecies) def remove(self, ekspecies): + """ + Remove an :obj:`~espressomd.electrokinetics.EKSpecies` from the container. + + Parameters + ---------- + ekspecies : :obj:`~espressomd.electrokinetics.EKSpecies` + Species to be removed. + + """ self.call_method("remove", object=ekspecies) diff --git a/src/script_interface/ObjectList.hpp b/src/script_interface/ObjectList.hpp index 72b7bb8f9b7..f9ad6f61093 100644 --- a/src/script_interface/ObjectList.hpp +++ b/src/script_interface/ObjectList.hpp @@ -52,6 +52,8 @@ class ObjectList : public ObjectContainer { virtual void add_in_core(const std::shared_ptr &obj_ptr) = 0; virtual void remove_in_core(const std::shared_ptr &obj_ptr) = 0; + virtual bool + has_in_core(const std::shared_ptr &obj_ptr) const = 0; public: ObjectList() { @@ -74,6 +76,12 @@ class ObjectList : public ObjectContainer { * @param element The element to add. */ void add(std::shared_ptr const &element) { + if (has_in_core(element)) { + if (Base::context()->is_head_node()) { + throw std::runtime_error("This object is already present in the list"); + } + throw Exception(""); + } add_in_core(element); m_elements.push_back(element); } @@ -84,6 +92,12 @@ class ObjectList : public ObjectContainer { * @param element The element to remove. */ void remove(std::shared_ptr const &element) { + if (not has_in_core(element)) { + if (Base::context()->is_head_node()) { + throw std::runtime_error("This object is absent from the list"); + } + throw Exception(""); + } remove_in_core(element); m_elements.erase(std::remove(m_elements.begin(), m_elements.end(), element), m_elements.end()); diff --git a/src/script_interface/accumulators/AutoUpdateAccumulators.hpp b/src/script_interface/accumulators/AutoUpdateAccumulators.hpp index 842f2e899bb..a5ada2cbd98 100644 --- a/src/script_interface/accumulators/AutoUpdateAccumulators.hpp +++ b/src/script_interface/accumulators/AutoUpdateAccumulators.hpp @@ -29,6 +29,11 @@ namespace ScriptInterface { namespace Accumulators { class AutoUpdateAccumulators : public ObjectList { + bool + has_in_core(std::shared_ptr const &obj_ptr) const override { + return ::Accumulators::auto_update_contains(obj_ptr->accumulator().get()); + } + void add_in_core(std::shared_ptr const &obj_ptr) override { ::Accumulators::auto_update_add(obj_ptr->accumulator().get()); } diff --git a/src/script_interface/constraints/Constraints.hpp b/src/script_interface/constraints/Constraints.hpp index f2b9d537c37..a8b304c872a 100644 --- a/src/script_interface/constraints/Constraints.hpp +++ b/src/script_interface/constraints/Constraints.hpp @@ -32,6 +32,9 @@ namespace ScriptInterface { namespace Constraints { class Constraints : public ObjectList { + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return ::Constraints::constraints.contains(obj_ptr->constraint()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { ::Constraints::constraints.add(obj_ptr->constraint()); } diff --git a/src/script_interface/shapes/Union.hpp b/src/script_interface/shapes/Union.hpp index 1227ec0b504..475e899ef24 100644 --- a/src/script_interface/shapes/Union.hpp +++ b/src/script_interface/shapes/Union.hpp @@ -39,6 +39,9 @@ class Union : public ObjectList { Union() : m_core_shape(std::make_shared<::Shapes::Union>()) {} private: + bool has_in_core(const std::shared_ptr &obj_ptr) const override { + return m_core_shape->contains(obj_ptr->shape()); + } void add_in_core(const std::shared_ptr &obj_ptr) override { m_core_shape->add(obj_ptr->shape()); } diff --git a/src/script_interface/tests/ObjectList_test.cpp b/src/script_interface/tests/ObjectList_test.cpp index 2645ffdd87e..4298d37fa1b 100644 --- a/src/script_interface/tests/ObjectList_test.cpp +++ b/src/script_interface/tests/ObjectList_test.cpp @@ -43,6 +43,10 @@ struct ObjectListImpl : ObjectList { std::vector mock_core; private: + bool has_in_core(const ObjectRef &obj_ptr) const override { + return std::find(mock_core.begin(), mock_core.end(), obj_ptr) != + mock_core.end(); + } void add_in_core(const ObjectRef &obj_ptr) override { mock_core.push_back(obj_ptr); } @@ -82,8 +86,8 @@ BOOST_AUTO_TEST_CASE(removing_elements) { BOOST_AUTO_TEST_CASE(clearing_elements) { // A cleared list is empty. ObjectListImpl list; - list.add(ObjectRef{}); - list.add(ObjectRef{}); + list.add(std::make_shared()); + list.add(std::make_shared()); list.clear(); BOOST_CHECK(list.elements().empty()); BOOST_CHECK(list.mock_core.empty()); diff --git a/src/script_interface/walberla/EKContainer.hpp b/src/script_interface/walberla/EKContainer.hpp index 5f33f50c28e..26a6c617ae4 100644 --- a/src/script_interface/walberla/EKContainer.hpp +++ b/src/script_interface/walberla/EKContainer.hpp @@ -61,6 +61,9 @@ class EKContainer : public ObjectList { std::shared_ptr<::EK::EKWalberla::ek_container_type> m_ek_container; bool m_is_active; + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return m_ek_container->contains(obj_ptr->get_ekinstance()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { context()->parallel_try_catch( [this, &obj_ptr]() { m_ek_container->add(obj_ptr->get_ekinstance()); }); diff --git a/src/script_interface/walberla/EKReactions.hpp b/src/script_interface/walberla/EKReactions.hpp index 84104049454..12f0fe81695 100644 --- a/src/script_interface/walberla/EKReactions.hpp +++ b/src/script_interface/walberla/EKReactions.hpp @@ -40,6 +40,9 @@ namespace ScriptInterface::walberla { class EKReactions : public ObjectList { std::shared_ptr<::EK::EKWalberla::ek_reactions_type> m_ek_reactions; + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return m_ek_reactions->contains(obj_ptr->get_instance()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { m_ek_reactions->add(obj_ptr->get_instance()); } diff --git a/src/shapes/include/shapes/Union.hpp b/src/shapes/include/shapes/Union.hpp index 492a7c42cc4..c719ac37826 100644 --- a/src/shapes/include/shapes/Union.hpp +++ b/src/shapes/include/shapes/Union.hpp @@ -33,6 +33,10 @@ namespace Shapes { class Union : public Shape { public: + bool contains(std::shared_ptr const &shape) const noexcept { + return std::find(m_shapes.begin(), m_shapes.end(), shape) != m_shapes.end(); + } + void add(std::shared_ptr const &shape) { m_shapes.emplace_back(shape); } diff --git a/src/shapes/unit_tests/Union_test.cpp b/src/shapes/unit_tests/Union_test.cpp index eb2357f5386..dfacb5ecfa8 100644 --- a/src/shapes/unit_tests/Union_test.cpp +++ b/src/shapes/unit_tests/Union_test.cpp @@ -42,8 +42,12 @@ BOOST_AUTO_TEST_CASE(dist_function) { wall2->d() = -10.0; Shapes::Union uni; + BOOST_CHECK(not uni.contains(wall1)); + BOOST_CHECK(not uni.contains(wall2)); uni.add(wall1); uni.add(wall2); + BOOST_CHECK(uni.contains(wall1)); + BOOST_CHECK(uni.contains(wall2)); auto check_union = [&wall1, &wall2, &uni](Utils::Vector3d const &pos) { double wall1_dist; diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index ce792497bc6..a3b20673a03 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -252,6 +252,7 @@ python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1) python_test(FILE script_interface.py MAX_NUM_PROC 4) python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 2) python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4) +python_test(FILE reaction_trivial.py MAX_NUM_PROC 2) python_test(FILE reaction_complex.py MAX_NUM_PROC 1) python_test(FILE reaction_bookkeeping.py MAX_NUM_PROC 1) python_test(FILE widom_insertion.py MAX_NUM_PROC 1) diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py index cbe2e829e16..9b1eeff8063 100644 --- a/testsuite/python/analyze_energy.py +++ b/testsuite/python/analyze_energy.py @@ -49,8 +49,8 @@ def setUpClass(cls): cls.system.bonded_inter[5] = cls.harmonic def setUp(self): - self.system.part.add(pos=[1, 2, 2], type=0) - self.system.part.add(pos=[5, 2, 2], type=0) + self.system.part.add(pos=[1, 2, 2], type=0, mol_id=6) + self.system.part.add(pos=[5, 2, 2], type=0, mol_id=6) def tearDown(self): self.system.part.clear() @@ -86,12 +86,14 @@ def test_non_bonded(self): self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7) self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7) self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7) + self.assertAlmostEqual(energy["non_bonded_intra"], 1., delta=1e-7) + self.assertAlmostEqual(energy["non_bonded_inter"], 0., delta=1e-7) # Test the single particle energy function self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum( [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7) # add another pair of particles - self.system.part.add(pos=[3, 2, 2], type=1) - self.system.part.add(pos=[4, 2, 2], type=1) + self.system.part.add(pos=[3, 2, 2], type=1, mol_id=7) + self.system.part.add(pos=[4, 2, 2], type=1, mol_id=7) energy = self.system.analysis.energy() self.assertAlmostEqual(energy["total"], 3., delta=1e-7) self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7) @@ -101,6 +103,18 @@ def test_non_bonded(self): self.assertAlmostEqual(energy["non_bonded", 0, 0] + energy["non_bonded", 0, 1] + energy["non_bonded", 1, 1], energy["total"], delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 0, 0], 1., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 1, 1], 1., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 0, 1], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 0, 0], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 1, 1], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 0, 1], 1., delta=1e-7) # Test the single particle energy function self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum( [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7) diff --git a/testsuite/python/reaction_trivial.py b/testsuite/python/reaction_trivial.py new file mode 100644 index 00000000000..7229928b11f --- /dev/null +++ b/testsuite/python/reaction_trivial.py @@ -0,0 +1,78 @@ +# +# Copyright (C) 2023 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import espressomd +import espressomd.reaction_methods +import unittest as ut +import numpy as np + + +class Test(ut.TestCase): + + """Test the trivial reaction 1A <-> 1B""" + + system = espressomd.System(box_l=[1., 1., 1.]) + system.time_step = 0.1 + system.cell_system.skin = 0.1 + np.random.seed(42) + + def test_equilibration(self): + system = self.system + N0 = 40 + c0 = 0.001 + types = {"A": 0, "B": 1} + system.box_l = np.ones(3) * np.cbrt(N0 / c0) + RE = espressomd.reaction_methods.ReactionEnsemble( + seed=42, kT=1., exclusion_range=1., search_algorithm="parallel") + RE.set_non_interacting_type(type=max(types.values()) + 1) + system.part.add( + pos=np.random.random((N0, 3)) * system.box_l, + type=N0 * [types["A"]]) + RE.add_reaction( + gamma=1., + reactant_types=[types["A"]], + reactant_coefficients=[1], + product_types=[types["B"]], + product_coefficients=[1], + default_charges={types["A"]: 0, types["B"]: 0}) + + # warmup + system.thermostat.set_langevin(kT=1., gamma=3., seed=42) + system.integrator.run(200) + RE.reaction(steps=5 * N0) + + # sampling + average_NA = 0.0 + num_samples = 100 + for _ in range(num_samples): + RE.reaction(steps=10) + system.integrator.run(20) + average_NA += system.number_of_particles(type=types["A"]) + average_NA /= num_samples + + alpha = average_NA / float(N0) + rate0 = RE.get_acceptance_rate_reaction(reaction_id=0) + rate1 = RE.get_acceptance_rate_reaction(reaction_id=1) + self.assertAlmostEqual(alpha, 0.50, delta=0.01) + self.assertAlmostEqual(rate0, 0.85, delta=0.05) + self.assertAlmostEqual(rate1, 0.85, delta=0.05) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/script_interface.py b/testsuite/python/script_interface.py index 7665a682514..8848e45e18c 100644 --- a/testsuite/python/script_interface.py +++ b/testsuite/python/script_interface.py @@ -102,6 +102,22 @@ def test_autoparameter_exceptions(self): with self.assertRaisesRegex(AttributeError, "Object 'HarmonicBond' has no attribute 'unknown'"): bond.unknown + def test_objectlist_exceptions(self): + """Check ObjectList framework""" + wall = espressomd.shapes.Wall(normal=[-1, 0, 0]) + constraint = espressomd.constraints.ShapeBasedConstraint(shape=wall) + constraints = espressomd.constraints.Constraints() + self.assertEqual(len(constraints), 0) + constraints.add(constraint) + self.assertEqual(len(constraints), 1) + with self.assertRaisesRegex(RuntimeError, "This object is already present in the list"): + constraints.add(constraint) + self.assertEqual(len(constraints), 1) + constraints.remove(constraint) + self.assertEqual(len(constraints), 0) + with self.assertRaisesRegex(RuntimeError, "This object is absent from the list"): + constraints.remove(constraint) + def test_feature_exceptions(self): """Check feature verification""" all_features = set(espressomd.code_info.all_features()) diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index e93755f95ab..7991e706561 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) # TODO: unstable, see issue #4798 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..eaa490e7a9c --- /dev/null +++ b/testsuite/scripts/tutorials/test_electrodes_1.py @@ -0,0 +1,44 @@ +# +# Copyright (C) 2023 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import unittest as ut +import importlib_wrapper +import numpy as np + +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/electrodes/electrodes_part1.py", N_AXIAL_POINTS=6, + script_suffix="@TEST_SUFFIX@") + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + + def test_force(self): + ref_force = tutorial.analytic_force_centered( + tutorial.r, tutorial.box_l_z) + msg = "The force for particle 1 doesn't agree with analytical expression" + np.testing.assert_allclose(np.log(tutorial.elc_forces_axial[:, 0]), + np.log(ref_force), rtol=0.05, err_msg=msg) + msg = "The force for particle 2 doesn't agree with analytical expression" + np.testing.assert_allclose(np.log(-tutorial.elc_forces_axial[:, 1]), + np.log(ref_force), rtol=0.05, err_msg=msg) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/scripts/tutorials/test_electrodes_2.py b/testsuite/scripts/tutorials/test_electrodes_2.py new file mode 100644 index 00000000000..5654db1ed88 --- /dev/null +++ b/testsuite/scripts/tutorials/test_electrodes_2.py @@ -0,0 +1,72 @@ +# +# Copyright (C) 2023 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import unittest as ut +import importlib_wrapper +import numpy as np +from scipy import constants + +params = {'N_SAMPLES_EQUIL': 25, 'N_SAMPLES_PROD': 5, + 'N_SAMPLES_EQUIL_CAP': 5, 'N_SAMPLES_CAP': 1, + 'MIN_PHI': 1, 'MAX_PHI': 2.5, 'N_PHI': 4} + +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/electrodes/electrodes_part2.py", + script_suffix="@TEST_SUFFIX@", **params) + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + + def test_potential_difference(self): + # Test that the applied potential difference equals the one from + # integrating Poisson equation + msg = 'The potential difference is not equal to the one from integrating Poisson equation.' + self.assertAlmostEqual( + tutorial.measured_potential_difference / tutorial.POTENTIAL_DIFF, + 1., delta=0.25, msg=msg) + + def test_charge_profile(self): + # Roughly test the profile, deviations are expected!! + charge_profile = ( + tutorial.cation_profile_mean + + tutorial.anion_profile_mean) + analytic = ( + tutorial.cation_profile_analytic + + tutorial.anion_profile_analytic) + msg = 'The density profile is not sufficiently equal to PB theory.' + np.testing.assert_allclose( + charge_profile[1:-1], + analytic[1:-1], + rtol=5e-2, + atol=5e-2, + err_msg=msg) + + def test_capacitance(self): + # For low potentials the capacitance should be in line with Grahame/DH + # equilibration performance limiting + grahame = -tutorial.sigma_vs_phi[:, 0] / ( + constants.elementary_charge / (constants.Boltzmann * tutorial.TEMPERATURE)) + msg = 'The capacitance at low potentials should be in line with Grahame/DH.' + np.testing.assert_allclose( + grahame, tutorial.sigma_vs_phi[:, 1], atol=.015, err_msg=msg) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/scripts/tutorials/test_electrokinetics.py b/testsuite/scripts/tutorials/test_electrokinetics.py index 611f8dafceb..e138fe209ad 100644 --- a/testsuite/scripts/tutorials/test_electrokinetics.py +++ b/testsuite/scripts/tutorials/test_electrokinetics.py @@ -1,4 +1,5 @@ -# Copyright (C) 2019-2022 The ESPResSo project +# +# Copyright (C) 2019-2023 The ESPResSo project # # This file is part of ESPResSo. # @@ -14,35 +15,109 @@ # # 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 as iw import numpy as np +import scipy.optimize tutorial, skipIfMissingFeatures = iw.configure_and_import( - "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", integration_length=400) + "@TUTORIALS_DIR@/electrokinetics/electrokinetics.py", TOTAL_FRAMES=0) @skipIfMissingFeatures class Tutorial(ut.TestCase): system = tutorial.system - def normalize_two_datasets(self, a, b): - offset = min(np.min(a), np.min(b)) - a -= offset - b -= offset - scale = max(np.max(a), np.max(b)) - a /= scale - b /= scale - - def test_simulation(self): - for varname, tol in zip(["density", "velocity"], [2, 5]): - sim = np.array(tutorial.__dict__[varname + "_list"]) - ana = np.array(tutorial.eof_analytical.__dict__[varname + "_list"]) - self.normalize_two_datasets(sim, ana) - accuracy = np.max(np.abs(sim - ana)) - # expecting at most a few percents deviation - self.assertLess(accuracy, tol / 100.) + def test_analytic_solutions(self): + # check advection-diffusion + mu_simulated = tutorial.positions_diagonal[np.argmax( + tutorial.values_diagonal)] + mu_analytic = tutorial.positions_analytic[np.argmax( + tutorial.values_analytic)] + tol = 2. * tutorial.AGRID + self.assertAlmostEqual(mu_simulated, mu_analytic, delta=tol) + # check electroosmotic flow + np.testing.assert_allclose( + np.copy(tutorial.density_eof), + tutorial.analytic_density_eof, + rtol=0.005) + np.testing.assert_allclose( + np.copy(tutorial.velocity_eof), + tutorial.analytic_velocity_eof, + rtol=0.05) + np.testing.assert_allclose( + np.copy(tutorial.pressure_tensor_eof), + tutorial.analytic_pressure_tensor_eof, + rtol=0.05) + # check pressure-driven flow + np.testing.assert_allclose( + np.copy(tutorial.density_pressure), + tutorial.analytic_density_eof, + rtol=0.005) + np.testing.assert_allclose( + np.copy(tutorial.velocity_pressure), + tutorial.analytic_velocity_pressure, + rtol=0.05) + np.testing.assert_allclose( + np.copy(tutorial.pressure_tensor_pressure), + tutorial.analytic_pressure_tensor_pressure, + rtol=0.05) + + def test_karman_vortex_street(self): + """ + Check for the formation of a Kármán vortex street. Due to turbulence, + a wavy pattern emerges. + """ + def get_frequency_species(species): + """ + Compute the principal wavevector of a chemical species. + """ + fdata = np.zeros(16) + for i in range(32, 80): + rdata = species[i, :, 0].density + fdata += np.abs(np.fft.fft(rdata - np.mean(rdata)))[:16] + return np.argmax(fdata) + + def get_phase_karman(species): + """ + Compute the time-dependent phase of a turbulent flow profile. + """ + phase = [] + k = 2 # wavevector for product species + for i in range(36, 68): + rdata = species[i, :, 0].density + phase.append(np.angle(np.fft.fft(rdata - np.mean(rdata))[k])) + return np.array(phase) + + def cosine_kernel(x, magnitude, freq, phase): + return magnitude * np.cos(x * freq + phase) + + tutorial.system.integrator.run(2000) + # check for finite values + for species in (*tutorial.educt_species, *tutorial.product_species): + assert np.all(np.isfinite(species[:, :, :].density)) + assert np.all(species[:, :, :].density >= 0) + assert np.all(np.isfinite(tutorial.lbf[:, :, :].velocity)) + # there is only one inlet per educt, thus wavelength == box width + self.assertEqual(get_frequency_species(tutorial.educt_species[0]), 1) + self.assertEqual(get_frequency_species(tutorial.educt_species[1]), 1) + # reaction happens in the mixing region, thus the frequency doubles + self.assertEqual(get_frequency_species(tutorial.product_species[0]), 2) + # check for turbulence onset + ref_params = np.array([2., 0.12, 1. / 2. * np.pi]) + sim_phase = get_phase_karman(tutorial.product_species[0]) + xdata = np.arange(sim_phase.shape[0]) + popt, _ = scipy.optimize.curve_fit( + cosine_kernel, xdata, sim_phase, p0=ref_params, + bounds=([-4., 0.08, 0.], [4., 0.24, 2. * np.pi])) + fit_phase = cosine_kernel(xdata, *popt) + rmsd = np.sqrt(np.mean(np.square(sim_phase - fit_phase))) + self.assertAlmostEqual(popt[2], ref_params[2], delta=0.20) + self.assertAlmostEqual(popt[1], ref_params[1], delta=0.02) + self.assertAlmostEqual(popt[0], ref_params[0], delta=0.80) + self.assertLess(rmsd / abs(popt[0]), 0.2) if __name__ == "__main__":