From 18a58174f654d630c2bdea943a3adaba0b063fa5 Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Mon, 6 Nov 2023 19:24:16 +0100 Subject: [PATCH 1/2] Add tutorial on GCMC --- doc/tutorials/CMakeLists.txt | 1 + doc/tutorials/Readme.md | 4 +- doc/tutorials/constant_pH/constant_pH.ipynb | 4 +- .../CMakeLists.txt | 25 + .../NotesForTutor.md | 14 + .../figures/schematic.svg | 257 ++++++ .../grand_canonical_monte_carlo.ipynb | 783 ++++++++++++++++++ testsuite/scripts/tutorials/CMakeLists.txt | 1 + .../test_grand_canonical_monte_carlo.py | 43 + 9 files changed, 1129 insertions(+), 3 deletions(-) create mode 100644 doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt create mode 100644 doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md create mode 100644 doc/tutorials/grand_canonical_monte_carlo/figures/schematic.svg create mode 100644 doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb create mode 100644 testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt index 479f3491e53..fc0ee979ba6 100644 --- a/doc/tutorials/CMakeLists.txt +++ b/doc/tutorials/CMakeLists.txt @@ -115,6 +115,7 @@ add_subdirectory(ferrofluid) add_subdirectory(constant_pH) add_subdirectory(widom_insertion) add_subdirectory(electrodes) +add_subdirectory(grand_canonical_monte_carlo) 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 869ded3b331..d70dd75b54a 100644 --- a/doc/tutorials/Readme.md +++ b/doc/tutorials/Readme.md @@ -69,7 +69,9 @@ physical systems. * **Widom particle insertion method** Measuring the excess chemical potential of a salt solution using the Widom particle insertion method. [Guide](widom_insertion/widom_insertion.ipynb) - +* **Grand-Canonical Monte Carlo** + Simulating a polyelectrolyte solution coupled to a reservoir of salt. + [Guide](grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb) [comment]: # (End of tutorials landing page) diff --git a/doc/tutorials/constant_pH/constant_pH.ipynb b/doc/tutorials/constant_pH/constant_pH.ipynb index 95a88f258b9..965941fee82 100644 --- a/doc/tutorials/constant_pH/constant_pH.ipynb +++ b/doc/tutorials/constant_pH/constant_pH.ipynb @@ -1119,7 +1119,7 @@ "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" }, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -1133,7 +1133,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.10.6" } }, "nbformat": 4, diff --git a/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt b/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt new file mode 100644 index 00000000000..895254c5e43 --- /dev/null +++ b/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt @@ -0,0 +1,25 @@ +# +# 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 . +# + +configure_tutorial_target(TARGET tutorial_grand_canonical_monte_carlo DEPENDS + grand_canonical_monte_carlo.ipynb + figures/schematic.svg) + +nb_export(TARGET tutorial_grand_canonical_monte_carlo SUFFIX "" FILE + "grand_canonical_monte_carlo.ipynb" HTML_RUN) diff --git a/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md b/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md new file mode 100644 index 00000000000..098aeae3a31 --- /dev/null +++ b/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md @@ -0,0 +1,14 @@ +# Notes for Tutors: Grand-Canonical Monte Carlo + +## Physics learning goals + +After the tutorial, students should be able to explain: + +* what the grand-canonical ensemble is and where it is applicable +* how to simulate in grand-canonical ensemble using GCMC + +## ESPResSo learning goals + +In the course of this tutorial, students should learn to: + +* how to set up a GCMC simulation in ESPResSo diff --git a/doc/tutorials/grand_canonical_monte_carlo/figures/schematic.svg b/doc/tutorials/grand_canonical_monte_carlo/figures/schematic.svg new file mode 100644 index 00000000000..1f605c7dab9 --- /dev/null +++ b/doc/tutorials/grand_canonical_monte_carlo/figures/schematic.svg @@ -0,0 +1,257 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb b/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb new file mode 100644 index 00000000000..56376895e15 --- /dev/null +++ b/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb @@ -0,0 +1,783 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Simulations in the Grand-Canonical Ensemble\n", + "\n", + "## Table of Contents\n", + "* [Introduction](#Introduction)\n", + "* [Grand-Canonical Ensemble](#Grand-Canonical-Ensemble)\n", + "* [How to simulate the grand-canonical ensemble?](#How-to-simulate-a-grand-canonical-ensemble?)\n", + "* [References](#References)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", + "Often, in soft-matter physics and chemistry we are interested in systems where the particle numbers are not fixed. This is for instance the case in systems where chemical reactions occur or when particles can be exchanged with a reservoir. The case of changing particle numbers is best described by the grand canonical ensemble.\n", + "\n", + "Canonical Monte-Carlo or molecular dynamics simulations are not suitable for this task, since they conserve the particle numbers. However, the Metropolis Monte Carlo method can be generalized to the Grand-Canonical ensemble in a straightforward manner, yielding the so-called Grand-Canonical Monte Carlo method (GCMC), which is introduced in this tutorial.\n", + "\n", + "This tutorial builds on the tutorial [Widom insertion](https://espressomd.github.io/tutorials/widom_insertion/widom_insertion.html), in which you measured the chemical potential of a monovalent salt solution by using a combination of MD and MC techniques. In this part you will make use of the results (i. e. the chemical potentials) of the previous tutorial to simulate the partitioning of salt ions between a reservoir and a polymer solution. Therefore we introduce and use the grand-canonical Monte Carlo method." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Grand-Canonical Ensemble\n", + "The case of changing particle numbers is best described by the grand canonical ensemble [1].\n", + "In the grand canonical ensemble, the volume $V$, the temperature $T$ and the chemical potentials $\\mu_i$ of the exchangeable species $i$ rather than the particle numbers $N_i$ are fixed. The chemical potentials $\\mu_i$ correspond to the free energy (Helmholtz free energy $F$, will be defined later) cost of adding a particle of type $i$ while keeping everything else fixed:\n", + "\\begin{align}\n", + "\\mu_i = \\left(\\frac{\\partial F}{\\partial N_i}\\right)_{T, V, \\left\\{N_j\\right\\}_{j\\neq i}}.\n", + "\\end{align}\n", + "The grand canonical ensemble is thus also referred to as the $(\\mu, V, T)$-ensemble.\n", + "In the grand canonical ensemble, the probability for the system to contain $N_i$ particles of type $i$, i.e. a total of $N=\\sum_{i=1}^{N_\\mathrm{s}}N_i$ ($N_\\mathrm{s}$ is the number of different types) and to be in a specific microstate $\\xi\\in\\Gamma_{N}$ is given by the probability distribution\n", + "\\begin{align}\n", + "p^\\mathrm{G}\\left(\\left\\{N_i\\right\\},\\xi\\right) = \\frac{\\exp\\left(-\\beta \\left(H(\\xi)-\\sum_{i=1}^{N_\\mathrm{s}}\\mu_iN_i\\right)\\right)}{Z^\\mathrm{G}h^{3N}\\prod_{i=1}^{N_\\mathrm{s}} N_i!}. \\quad \\tag{1}\n", + "\\end{align}\n", + "Here, $H$ denotes the Hamiltonian of the system with $N_i$ particles of type $i$, $h$ is the planck constant, $\\beta=1/k_\\mathrm{B}T$ the inverse temperature and $Z^\\mathrm{G}$ is the grand-canonical partition function which is given by\n", + "\\begin{align}\n", + "Z^\\mathrm{G}\\left(\\left\\{\\mu_i\\right\\},V,T\\right) = \\sum_{N_1=0}^{\\infty}...\\sum_{N_\\mathrm{s}=0}^{\\infty}Z\\left(\\left\\{N_i\\right\\},V,T\\right)\\exp\\left(\\beta \\sum_{i=1}^{N_\\mathrm{s}}\\mu_iN_i\\right),\n", + "\\end{align}\n", + "with the canonical partition function\n", + "\\begin{align}\n", + "\tZ\\left(\\left\\{N_i\\right\\},V,T\\right) = \\frac{1}{h^{3N} \\prod_{i=1}^{N_s}N_i!}\\int_{\\Gamma_N} \\mathrm{d}\\xi \\ \\exp(-\\beta H(\\xi)).\n", + "\\end{align}\n", + "The partition function is connected to a thermodynamic potential, called the grand potential (Landau free energy):\n", + "\\begin{align}\n", + " \\Omega\\left(\\left\\{\\mu_i\\right\\},V,T\\right) = -k_\\mathrm{B}T\\ln\\left(Z^\\mathrm{G}\\left(\\left\\{\\mu_i\\right\\},V,T\\right)\\right).\n", + "\\end{align}\n", + "In a similiar way, the Helmholtz free energy is defined as\n", + "\\begin{align}\n", + "F\\left(\\left\\{N_i\\right\\},V,T\\right) = -k_\\mathrm{B}T\\ln\\left(Z\\left(\\left\\{N_i\\right\\},V,T\\right)\\right).\n", + "\\end{align}\n", + "\n", + "When considering a system that can exchange particles with a reservoir, the chemical potential $\\mu_i$ for each exchangeable species $i$ has to be equal in the reservoir and in the system:\n", + "\\begin{equation}\n", + "\\mu_i^{\\mathrm{sys}} = \\mu_i^{\\mathrm{res}}.\n", + "\\end{equation}\n", + "This condition is referred to as a chemical equilibrium between the two phases and is equivalent to a minimization of the free energy of the total system (system + reservoir).\n", + "In the case where one requires electroneutrality of the phases, due to the additional constraint, the electrochemical potential $\\overline{\\mu}_i=\\mu_i+z_ie\\psi$ rather than the chemical potential has to be equal. Here, $z_i$ is the valency of type $i$, $e$ is the elementary charge and $\\psi$ is the local electrostatic potential." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## How to simulate a grand-canonical ensemble?\n", + "### Metropolis-Hastings Algorithm\n", + "The grand-canonical distribution (see eq. (1)) can be sampled using the Metropolis Monte Carlo method, similiar to the canonical ensemble. Before we derive the acceptance criterion for the GCMC method, let us shortly recall the generic Metropolis-Hastings algorithm.\n", + "\n", + "We want to calculate averages of observables $A$ which are distributed according to some distribution $p^\\mathrm{eq}_s$ (e.g. the canonical or grand-canonical distribution)\n", + "\\begin{equation}\n", + "\\langle A\\rangle = \\frac{\\sum_{s\\in\\Gamma} A_s p^\\mathrm{eq}_s}{\\sum_{s\\in\\Gamma} p^\\mathrm{eq}_s}. \\quad \\tag{2}\n", + "\\label{eq:can_average_mc}\n", + "\\end{equation}\n", + "Here, $s$ labels the different states of the systems, $A_s$ is the value of the observable $A$ in the state $s$ and the sum runs over the set of all states $\\Gamma$.\n", + "The sum can also be an integral for systems with continuous degrees of freedom, for this derivation it does not matter.\n", + "For almost all systems, the ensemble average (2) can not be evaluated exactly and only numerical approximations can be obtained.\n", + "The basic idea of all MC methods is to randomly sample only a subset of states.\n", + "In the most naive approach (simple sampling), states are simply selected with a uniform probability from the space $\\Gamma$ of states $s$, resulting in a sample $\\gamma\\subset\\Gamma$.\n", + "The canonical average can then be approximated by a sum over $\\gamma$:\n", + "\\begin{equation}\n", + "\\langle A\\rangle \\approx \\frac{\\sum_{s\\in\\gamma} A_s p^\\mathrm{eq}_s}{\\sum_{s\\in\\gamma} p^\\mathrm{eq}_s} \\quad \\tag{3}\n", + "\\end{equation}\n", + "This method is highly inefficient, as most states will typically not contribute significantly to the average.\n", + "To get a more efficient sampling method, let us select the states $s$ according to a non-uniform probability distribution $p_s$.\n", + "Then, the average can be written as \n", + "\\begin{equation}\n", + "\\langle A\\rangle^{\\mathrm{can}} \\approx \\frac{\\sum_{s\\in\\gamma} A_s p^\\mathrm{eq}_s/p_s}{\\sum_{s\\in\\gamma} p^\\mathrm{eq}_s/p_s}, \\quad \\tag{4}\n", + "\\label{eq:can_average_mc_prob}\n", + "\\end{equation}\n", + "where we have to divide in the numerator and in the denominator by $p_s$ to compensate for the fact that the samples are distributed non-uniformly.\n", + "It is obvious that we recover (3) from (4) as the special case of a uniform distribution. \n", + "To sample the averages more efficiently, we simply pick\n", + "\\begin{equation}\n", + "p_s \\propto p^\\mathrm{eq}_s \\quad \\tag{5}\n", + "\\label{eq:mc_prob}\n", + "\\end{equation}\n", + "such that the selected samples are representative of the desired equilibrium distribution.\n", + "In this case, the estimators of the averages are thus simply given by\n", + "\\begin{equation}\n", + "\\langle A\\rangle \\approx \\frac{1}{N}\\sum_{s\\in\\gamma}A_s.\n", + "\\end{equation}\n", + "The remaining problem is now to find a method which allows us to sample states according to (5).\n", + "This can be done in terms of so-called Markov chains. \n", + "A Markov chain is a stochastic process without memory.\n", + "In our case, the Markov chain is simply a sequence $\\{s_i\\}_{i=1,...N}$ of $N$ states $s_i$ of the system and the fact that the Markov chain has no memory means that the transition probability from a state $s_i$ to a state $s_j$ can be given in terms of a transition matrix $W_{i\\rightarrow j}$ which is constant.\n", + "The equilibrium probability distribution has to be invariant under $W$. Mathematically this condition can be expressed in the form\n", + "\\begin{align}\n", + "\\sum_{i}W_{i\\rightarrow j}p_i^\\mathrm{eq} = p_j^\\mathrm{eq}\n", + "\\end{align}\n", + "A sufficient constraint on $W_{i\\rightarrow j}$ which ensures the validity of this condition\n", + "is the so-called criterion of detailed balance:\n", + "\\begin{align}\n", + "\\frac{W_{i\\rightarrow j}}{W_{j\\rightarrow i}} = \\frac{p_j^\\mathrm{eq}} {p_i^\\mathrm{eq}}.\n", + "\\end{align}\n", + "Typically, the transition matrix elements are factorized into a proposal probability $g_{i\\rightarrow j}$ and acceptance probability $P_{i\\rightarrow j}^\\mathrm{acc}$:\n", + "\\begin{align}\n", + "W_{i\\rightarrow j} = g_{i\\rightarrow j} P_{i\\rightarrow j}^\\mathrm{acc}. \n", + "\\end{align}\n", + "For the symmetric choice $g_{i\\rightarrow j}=g_{j\\rightarrow i}$, the condition of detailed balance reduces to \n", + "\\begin{align}\n", + "\\frac{P_{i\\rightarrow j}^\\mathrm{acc}}{P_{j\\rightarrow i}^\\mathrm{acc}} = \\frac{p_j^\\mathrm{eq}}{p_i^\\mathrm{eq}}. \n", + "\\quad \\tag{6}\n", + "\\end{align}\n", + "By checking the cases $p_j^\\mathrm{eq}>p_i^\\mathrm{eq}$ and $p_j^\\mathrm{eq} \\leq p_i^\\mathrm{eq}$ one can see that the acceptance probability \n", + "\\begin{align}\n", + "P_{i\\rightarrow j}^\\mathrm{acc,\\text{ }Metropolis} = \\min\\left(1,\\frac{p_j^\\mathrm{eq}}{p_i^\\mathrm{eq}}\\right) \\quad \\tag{7}\n", + "\\end{align}\n", + "fulfills equation (6). With this, we can now formulate the Metropolis-Hastings algorithm [2].\n", + "We start with an arbitrary initial state $s_1$ and generate a sequence $\\{s_i\\}_{i=1,...N}$ of $N$ states.\n", + "Each new state (n) is generated from the previous state (o) in the following way. \n", + "First, we propose a new state according to the symmetric proposal probability $g_{\\mathrm{o}\\rightarrow\\mathrm{n}}$.\n", + "Then, we accept the new configuration with a probability given by equation (7).\n", + "If the proposed new state is rejected, the old state is retained.\n", + "Applying this algorithm many times results in the desired Markov chain of representative configurations.\n", + "\n", + "### GCMC method\n", + "To simulate a system in a grand-canonical ensemble the described algorithm can be applied the grand-canonical distribution [1]. When considering the possible states of the system, it becomes clear that every state can be reached by a combination of two kinds of moves:\n", + "1. Displacement moves that leave the particle numbers fixed and only change positions. This is equivalent to sampling in the canonical ensemble with the current particle number. Due to the ergodic hypothesis this operation can be performed either using canonical Monte Carlo methods or canonical molecular dynamic simulations (e. g. Langevin dynamics). In this tutorial we choose the latter.\n", + "\n", + "2. Insertion and Deletion moves that add or delete a single particle while keeping all other particle positions fixed.\n", + "\n", + "In the GCMC method the insertion and deletion moves are implemented in the following fashion. First of all it is determined randomly with equal probability if an insertion or deletion move is performed. In the case of an insertion the additional particle is placed at a uniformly chosen random position in the simulation box. Then the proposed new state is accepted according to the criterion\n", + "\\begin{align}\n", + "\\begin{split}\n", + "P_{N_i\\rightarrow N_i+1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{1}{N_i+1}\\left(\\frac{V}{\\Lambda_i^3}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i+1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)+\\beta \\mu_i\\right)\\right).\n", + "\\end{split}\n", + "\\end{align}\n", + "In the case of a deletion a randomly selected particle is removed from the system. \n", + "We have the analogous result\n", + "\\begin{align}\n", + "\\begin{split}\n", + "P_{N_i\\rightarrow N_i-1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,N_i\\left(\\frac{\\Lambda_i^3}{V}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i-1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)-\\beta \\mu_i\\right)\\right).\n", + "\\end{split}\n", + "\\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise**\n", + "\n", + "* Derive the acceptance criteria for the GCMC method by combining the grand-canonical distribution with the Metropolis-Hastings acceptance criteria.\n", + "\n", + "**Hint**\n", + "\n", + "* Use the grand-canonical distribution with all the momenta integrated out:\n", + " \\begin{align}\n", + "p^\\mathrm{eq}_{N} \\propto \\exp\\left(-\\beta U(\\mathbf{q}) \\right ) \\prod_i \\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp \\left (\\beta \\mu_i N_i\\right)\n", + "\\end{align}\n", + "This alternative expression is valid if we are only interested in observables which do not depend on the momenta." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "For the case ($N_i\\rightarrow N_i + 1$) one gets\n", + "\\begin{align}\n", + "P_{N_i\\rightarrow N+1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{p^\\mathrm{eq}_{N_i+1}}{p^\\mathrm{eq}_{N_i}}\\right) = \\min\\left(1,\\frac{\\frac{1}{(N_i+1)!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i+1}\\times \\exp\\left(-\\beta U^{N_i+1}(\\mathbf{q})+\\beta \\mu (N_i+1)\\right)}{\\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp\\left(-\\beta U^{N_i}(\\mathbf{q})+\\beta \\mu N_i\\right)}\\right)\\\\\n", + "&= \\min\\left(1,\\frac{1}{N_i+1}\\left(\\frac{V}{\\Lambda_i^3}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i+1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)+\\beta \\mu_i\\right)\\right).\n", + "\\end{align}\n", + "For the other case ($N_i\\rightarrow N_i-1$) we have the analogous result\n", + "\\begin{align}\n", + "P_{N_i\\rightarrow N_i-1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{p^\\mathrm{eq}_{N_i-1}}{p^\\mathrm{eq}_{N_i}}\\right) = \\min\\left(1,\\frac{\\frac{1}{(N_i-1)!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i-1}\\times \\exp\\left(-\\beta U^{N_i-1}(\\mathbf{q})+\\beta \\mu_i (N_i-1)\\right)}{\\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp\\left(-\\beta U^{N_i}(\\mathbf{q})+\\beta \\mu_i N_i\\right)}\\right)\\\\\n", + "&= \\min\\left(1,N_i\\left(\\frac{\\Lambda_i^3}{V}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i-1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)-\\beta \\mu_i\\right)\\right).\n", + "\\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The Simulated System: Polyelectrolyte solution coupled to a reservoir\n", + "\n", + "
\n", + "
\n", + "
Fig. 1: Schematic representation of the considered two phase setup.
\n", + "
\n", + "
\n", + "Here we want to investigate a polyelectrolyte solution that is coupled to a reservoir containing a monovalent salt, as shown schematically in figure 1. Similar systems are useful to investigate for example the swelling behavior of hydrogels. The surrounding of the gel can be considered as the reservoir and the inner of the gel is the system. The small salt ions (X$^+$ and X$^-$) can be exchanged between the two phases while the polyelectrolyte is confined to the system. Experimentally such a setup could be realized by seperating the two phases with a semi-permeable membrane. Similiar setups also occur in the study of polyelectrolyte hydrogels.\n", + "\n", + "Due to the macroscopic electroneutrality of both phases the salt ions obey an electrochemical equilibrium, i. e. there is a difference in the electrostatic potential between the two phases $$\\mu_i^{\\mathrm{sys}}+z_ie\\psi^{\\mathrm{Donnan}} = \\mu_i^{\\mathrm{res}}.$$\n", + "This effect is called the Donnan effect [3] and leads to an unequal partioning of salt between the reservoir and the system. The partitioning can be quantified using the partition coeffiecents $\\xi_i \\equiv c^{\\mathrm{sys}}_i/c^{\\mathrm{res}}_i$, where $c^{\\mathrm{sys}}_i$ and $c^{\\mathrm{res}}_i$ are the particle concentrations of type $i$ in the system respectively in the system. Defining the \"universal\" partition coefficient,\n", + "\\begin{align*}\n", + " \\xi\\equiv\\left\\{\n", + " \\begin{array}{ll}\n", + " \\xi_+ - \\frac{c_{\\text{M}^-}^{\\text{sys}}}{c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}& \\text{for positive ions}\\\\\n", + " \\xi_- & \\text{for negative ions,}\\\\\n", + " \\end{array}\\right.\n", + "\\end{align*}\n", + "one can obtain the following solution for an ideal system without any interaction between the particles:\n", + "\\begin{align}\n", + " \\xi^{\\text{Donnan}} = -\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}} + \\sqrt{\\left(\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}\\right)^2+1}.\n", + " \\end{align}\n", + " This solution shows that salt is rejected by the polyelctrolyte solution. In the following we want to investigate this partitioning for a system with electrostatic interactions.\n", + " \n", + "Because of the electroneutrality condition the two definitions of the \"universal\" partition coefficients for positive and negative ions are equal. The additional term in the definition for positive ions is a correction due to the counterions in the system, which neutralize the charged polyelectrolytes." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the folllowing, we plot the universal parition coefficient over the ratio of the concentrations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "plt.rcParams.update({\"font.size\": 18})\n", + "\n", + "ratios = np.logspace(-1, 2, 1000)\n", + "plt.figure(figsize=(10, 6))\n", + "plt.loglog(ratios, -ratios + np.sqrt((ratios)**2 + 1), label=r\"$\\xi^{\\mathrm{Donnan}}$\")\n", + "plt.xlabel(r\"$\\dfrac{c_{\\mathrm{M}^-}^{\\mathrm{sys}}}{2c_{\\mathrm{X}^+,\\mathrm{X}^-}^{\\mathrm{res}}}$\")\n", + "plt.ylabel(r\"$\\xi$\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation Setup\n", + "### Load data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us first add required modules." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import tqdm # module for progress bar\n", + "import scipy.interpolate\n", + "import pint # module for unit conversions\n", + "\n", + "import espressomd\n", + "import espressomd.interactions\n", + "import espressomd.electrostatics\n", + "import espressomd.reaction_methods\n", + "import espressomd.polymer\n", + "\n", + "espressomd.assert_features([\"ELECTROSTATICS\", \"P3M\", \"WCA\"])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the next step we define the [Lennard-Jones](https://espressomd.github.io/tutorials4.2.0/lennard_jones/lennard_jones.html) units. We use the module *pint* in order to make unit conversions between reduced simulation units and SI units easier." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "##### Units\n", + "ureg = pint.UnitRegistry()\n", + "sigma = 0.355 * ureg.nm # Sigma in SI units.\n", + "# This is corresponds to the half of the Bjerrum length in water (eps_r approx 80) and room temperature (293 K)\n", + "T = 298.15 * ureg.K # 25 degree celcius corresponds to approx room temperature\n", + "CREF_IN_MOL_PER_L = 1.0 * ureg.molar #in mol/l\n", + "\n", + "# Define the reduced units\n", + "ureg.define(f\"reduced_length = {sigma}\")\n", + "ureg.define(f\"reduced_charge = 1* e\")\n", + "ureg.define(f\"reduced_energy = {T} * boltzmann_constant\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We load the resulting excess chemical potentials for different salt concentrations of the reservoir from the [Widom insertion tutorial](https://espressomd.github.io/tutorials/widom_insertion/widom_insertion.html). The loaded excess chemical potentials are in units of $k_\\mathrm{B}T$ (i.e. in reduced units)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "salt_concentration_magnitudes_si = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3]\n", + "salt_concentration_si = np.array(salt_concentration_magnitudes_si) * ureg.molar\n", + "salt_concentration_sim = salt_concentration_si.to(\"molecules/reduced_length^3\")\n", + "excess_chemical_potential_data = [\n", + " -0.0672291585811938, -0.123732052028611, -0.218687259792629,\n", + " -0.326207586236404, -0.510091568808112, -0.657928029835777\n", + "]\n", + "excess_chemical_potential_data_error = [\n", + " 0.000994798843587, 0.00152160176511698, 0.00220162667953136,\n", + " 0.0029299206854553, 0.00422309102221265, 0.00603137482524256\n", + "]\n", + "excess_chemical_potential_monovalent_pairs_in_bulk = scipy.interpolate.interp1d(\n", + " salt_concentration_sim.magnitude, excess_chemical_potential_data)\n", + "fig = plt.figure(figsize=(10, 6))\n", + "plt.errorbar(salt_concentration_si.magnitude,\n", + " excess_chemical_potential_data,\n", + " yerr=excess_chemical_potential_data_error,\n", + " linestyle=\"-\",\n", + " marker=\"o\",\n", + " markersize=5.,\n", + " capsize=5.)\n", + "plt.xlabel(r\"salt concentration $c_{\\mathrm{salt}}$ (mol/l)\")\n", + "plt.ylabel(r\"excess chemical potential $\\mu_{\\mathrm{ex}}$ ($k_{\\mathrm{B}}T$)\")\n", + "plt.xscale(\"log\")\n", + "plt.xlim((5e-4, 0.5))\n", + "plt.ylim((-0.8, 0));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Initializing the system\n", + "Now we are ready to set up the system. In the first step we define the system parameters such as charges and system size as well as simulation paramters such as the time step. For small salt concentrations we need to produce more samples to get good statistics for $\\xi$ due to the finite system size. For a real production run it is recommended to use larger sample sizes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "##### Particle types and charges\n", + "types = {\n", + " \"Xplus\": 1, # Positive ions\n", + " \"Xminus\": 2, # Negative ions\n", + " \"Monomer\": 3, # Negative monomers\n", + "}\n", + "\n", + "charges = {\n", + " \"Xplus\": 1.0,\n", + " \"Xminus\": -1.0,\n", + " \"Monomer\": -1.0,\n", + "}\n", + "# Initialize concentrations\n", + "c_salt_res_sim = salt_concentration_sim[0]\n", + "c_monomer = 0.1 * ureg.molar\n", + "C_MONOMERS_SIM = c_monomer.to(\"molecules/reduced_length^3\")\n", + "\n", + "##### Number of monomers\n", + "N_CHAINS = 10\n", + "N_MONOMERS_PER_CHAIN = 5\n", + "\n", + "##### System size\n", + "BOX_LENGTH = np.power(\n", + " N_CHAINS * N_MONOMERS_PER_CHAIN / C_MONOMERS_SIM.magnitude, 1.0 / 3.0)\n", + "\n", + "##### Initial number of salt ion pairs\n", + "n_salt = int(c_salt_res_sim.magnitude * BOX_LENGTH**3)\n", + "\n", + "##### Weeks-Chandler-Andersen interaction\n", + "LJ_EPSILON = 1.0 # in reduced units, i.e. 1 sigma\n", + "LJ_SIGMA = 1.0 # in reduced units, i.e. 1 k_B T\n", + "\n", + "##### Electrostatic interaction\n", + "L_BJERRUM = 2.0 # in reduced units, i.e. 2 sigma\n", + "\n", + "##### Langevin-Thermostat\n", + "KT = 1.0 # in reduced units\n", + "GAMMA = 1.0 # in reduced units\n", + "\n", + "##### Integrator parameters\n", + "DT = 0.01 # Integrator time step\n", + "SKIN = 0.4\n", + "\n", + "warmup_loops = 50\n", + "number_of_loops_for_each_concentration = [\n", + " 750, 750, 250, 100, 100, 100\n", + "] # set it higher for better accuracy\n", + "steps_per_loop = 1000\n", + "samples_per_loop = 100\n", + "\n", + "##### Seeding\n", + "langevin_seed = 42\n", + "reaction_seed = 42" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we define our box and add our particles. We create the polymer chains as in the corresponding [tutorial](https://espressomd.github.io/tutorials/polymers/polymers.html). In addition, we enable the electrostatic interactions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "##### Create an instance of the system class\n", + "system = espressomd.System(box_l=[BOX_LENGTH] * 3)\n", + "system.time_step = DT\n", + "system.cell_system.skin = SKIN\n", + "print(\"Created system class\")\n", + "\n", + "#### Add the FENE interaction\n", + "fene = espressomd.interactions.FeneBond(k=30, d_r_max=1.5)\n", + "system.bonded_inter.add(fene)\n", + "\n", + "#### Add the polymer chains\n", + "polymer_positions = espressomd.polymer.linear_polymer_positions(\n", + " n_polymers=N_CHAINS,\n", + " beads_per_chain=N_MONOMERS_PER_CHAIN,\n", + " bond_length=0.9,\n", + " seed=42)\n", + "\n", + "for positions in polymer_positions:\n", + " monomers = system.part.add(pos=positions,\n", + " type=[types[\"Monomer\"]] * N_MONOMERS_PER_CHAIN,\n", + " q=[charges[\"Monomer\"]] * N_MONOMERS_PER_CHAIN)\n", + " previous_part = None\n", + " for part in monomers:\n", + " if not previous_part is None:\n", + " part.add_bond((fene, previous_part))\n", + " previous_part = part\n", + "\n", + "##### Add the particles\n", + "N_plus = n_salt + N_CHAINS * N_MONOMERS_PER_CHAIN\n", + "system.part.add(type=[types[\"Xplus\"]] * N_plus,\n", + " pos=np.random.rand(N_plus, 3) * BOX_LENGTH,\n", + " q=[charges[\"Xplus\"]] * N_plus)\n", + "system.part.add(type=[types[\"Xminus\"]] * n_salt,\n", + " pos=np.random.rand(n_salt, 3) * BOX_LENGTH,\n", + " q=[charges[\"Xminus\"]] * n_salt)\n", + "print(\"Added particles\")\n", + "\n", + "##### Add non-bonded interactions\n", + "for i in types:\n", + " for j in types:\n", + " system.non_bonded_inter[types[i],\n", + " types[j]].wca.set_params(epsilon=LJ_EPSILON,\n", + " sigma=LJ_SIGMA)\n", + "\n", + "print(\"Added WCA interaction\")\n", + "\n", + "p3m = electrostatics.P3M(prefactor=L_BJERRUM * KT, accuracy=1e-3)\n", + "system.electrostatics.solver = p3m\n", + "\n", + "print(\"Added electrostatics\")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Now we want to add a coupling to the reservoir. We can formally represent the insertion and deletion of ion pairs as a chemical reaction:\n", + "\n", + "\\begin{align}\n", + "\\emptyset &\\Leftrightarrow \\mathrm{X}^- + \\mathrm{X}^+\n", + "\\end{align}\n", + "It is important to note that the reaction can run in both directions, i. e. insertion and deletion of salt ions. These reactions are described by a concentration based equilibrium constant:\n", + "\\begin{align}\n", + " K = c_{\\mathrm{salt,res}}^2 \\exp \\left ( \\mu_{\\mathrm{salt,res}}^{\\mathrm{ex}} \\right )\n", + "\\end{align}\n", + "Here $\\mu_{\\mathrm{salt,res}}^{\\mathrm{ex}}$ denotes the excess chemical potential of a salt ion pair as we measured in the [Widom insertion tutorial](https://espressomd.github.io/tutorials/widom_insertion/widom_insertion.html).\n", + "\n", + "**Exercise**\n", + "\n", + "* Set up the reaction defined above with the corresponding equilibrium constant.\n", + "\n", + "**Hint**\n", + "\n", + "* You can find further explanation for setting up a reaction [here](https://espressomd.github.io/doc/reaction_methods.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "K_XX = c_salt_res_sim.magnitude**2 * np.exp(\n", + " excess_chemical_potential_monovalent_pairs_in_bulk(\n", + " c_salt_res_sim.magnitude))\n", + "RE = espressomd.reaction_methods.ReactionEnsemble(kT=KT,\n", + " exclusion_range=1.0,\n", + " seed=reaction_seed)\n", + "RE.add_reaction(gamma=K_XX,\n", + " reactant_types=[],\n", + " reactant_coefficients=[],\n", + " product_types=[types[\"Xplus\"], types[\"Xminus\"]],\n", + " product_coefficients=[1, 1],\n", + " default_charges={\n", + " types[\"Xplus\"]: charges[\"Xplus\"],\n", + " types[\"Xminus\"]: charges[\"Xminus\"]\n", + " })\n", + "# Set the non interacting type at the smallest integer.\n", + "# This may speed up the simulation (see Espresso docummentation section 19. Reaction methods)\n", + "RE.set_non_interacting_type(type=len(types) + 1)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we perform the initial warm up." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "##### Warmup\n", + "system.integrator.set_steepest_descent(f_max=0,\n", + " gamma=10.0,\n", + " max_displacement=0.001)\n", + "print(\"Removing overlaps with steepest descent...\")\n", + "for i in tqdm(range(warmup_loops),\n", + " bar_format='{l_bar}{bar:60}{r_bar}{bar:-60b}'):\n", + " system.integrator.run(steps_per_loop)\n", + "print(\"Done.\")\n", + "\n", + "system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=langevin_seed)\n", + "system.integrator.set_vv()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Now we are ready to perform the actual sampling.\n", + "\n", + "**Exercise**\n", + "* In order to sample for different reservoir salt concentrations, set up a function ``perform_sampling(salt_concentration)`` that takes as input the salt concentration of the reservoir.\n", + "* In this function perform a warm up loop in order to equilibrate the system to the new salt concentration before the production run.\n", + "* Measure the particle numbers of positive and negative salt ions in each iteration in order to return the partition coefficients $\\xi_+$ and $\\xi_-$ at the end.\n", + "\n", + "**Hint**\n", + "* Make sure to reset the equilibrium constant according to the salt concentration." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def perform_sampling(salt_concentration, number_of_loops):\n", + " K_XX = salt_concentration**2 * np.exp(\n", + " excess_chemical_potential_monovalent_pairs_in_bulk(salt_concentration))\n", + " RE.change_reaction_constant(gamma=K_XX, reaction_id=0)\n", + " for i in tqdm.trange(warmup_loops, bar_format=\"{l_bar}{bar:60}{r_bar}{bar:-60b}\"):\n", + " system.integrator.run(steps_per_loop)\n", + " RE.reaction(samples_per_loop)\n", + "\n", + " particle_numbers_positive = []\n", + " particle_numbers_negative = []\n", + "\n", + " for i in tqdm.trange(number_of_loops, bar_format=\"{l_bar}{bar:60}{r_bar}{bar:-60b}\"):\n", + " system.integrator.run(steps_per_loop)\n", + " RE.reaction(samples_per_loop)\n", + "\n", + " particle_numbers_positive.append(\n", + " system.number_of_particles(type=types[\"Xplus\"]))\n", + " particle_numbers_negative.append(\n", + " system.number_of_particles(type=types[\"Xminus\"]))\n", + " partition_coefficient_positive = np.mean(np.asarray(particle_numbers_positive))\\\n", + " / (BOX_LENGTH**3 * salt_concentration)\n", + " partition_coefficient_negative = np.mean(np.asarray(particle_numbers_negative))\\\n", + " / (BOX_LENGTH**3 * salt_concentration)\n", + " return partition_coefficient_positive, partition_coefficient_negative\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can perform the actual simulations for the different salt concentrations and measure the partition coefficients." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "partition_coefficients_positives_array = np.zeros_like(salt_concentration_sim)\n", + "partition_coefficients_negatives_array = np.zeros_like(salt_concentration_sim)\n", + "\n", + "for i, salt_concentration in enumerate(salt_concentration_sim.magnitude):\n", + " print(f\"Salt concentration {i+1}/{len(salt_concentration_sim)}\")\n", + " partition_coefficients_positives_array[\n", + " i], partition_coefficients_negatives_array[i] = perform_sampling(\n", + " salt_concentration, number_of_loops_for_each_concentration[i])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To compare the results of our simulations we define a function for the analytical solution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def analytical_solution(ratio):\n", + " return -ratio + np.sqrt(ratio**2 + 1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we can measure the partition coefficients derived from the simulations and compare them to the analytical results for an ideal system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(10, 6))\n", + "plt.loglog(c_monomer.magnitude/(2*salt_concentration_si.magnitude), partition_coefficients_negatives_array,\\\n", + " label = r'$\\xi$ from $\\xi_-$',linestyle = 'none',marker = 's', markersize=10)\n", + "universal_partion_coefficient_positive = \\\n", + " partition_coefficients_positives_array - c_monomer.magnitude/salt_concentration_si.magnitude\n", + "plt.loglog(c_monomer.magnitude/(2*salt_concentration_si.magnitude),universal_partion_coefficient_positive ,\\\n", + " label = r'$\\xi$ from $\\xi_+$',marker = 'd', markersize=7, linestyle = 'none', color=\"lime\")\n", + "ratios = np.logspace(-1, 2, num=200)\n", + "plt.loglog(ratios, analytical_solution(ratios), color=\"black\")\n", + "plt.xlabel(r'$\\frac{c_{\\mathrm{M}^-}}{2c_{\\mathrm{salt}}^{\\mathrm{res}}}$')\n", + "plt.ylabel(r'$\\xi$')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "**Exercise**\n", + "* Interpret the deviation of the simulation results from the ideal prediction." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "First, we note that the data points for negative and positive ions\n", + "always collapse onto the same value for ξ, i.e. ξ is indeed a universal value that\n", + "is always the same for negative and positive ions. When comparing these results\n", + "to the analytical prediction for an ideal gas, it is easy to see that the partition\n", + "coefficient is higher for the case with electrostatic interactions. This is essentially\n", + "due to the fact that the overall charge density inside the system is higher than in\n", + "the reservoir, and thus the activity coefficient of an ion pair inside the system is\n", + "smaller than in the reservoir, i.e. the free energy cost of putting additional ion\n", + "pairs into the system is overall lowered as compared to an ideal gas. Still, we see\n", + "that ξ < 1 even for the interacting case, i.e. the salt concentration is still smaller\n", + "inside the system than in the reservoir." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[1] D. Frenkel, Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. 2nd edition, chapter 5.6: Grand-Canonical Ensemble. Academic Press, 2002, ISBN: 9780122673511.\n", + "\n", + "[2] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21, 1087 (1953); https://doi.org/10.1063/1.1699114. \n", + "\n", + "[3] F. Donnan. The theory of membrane equilibria. Chemical reviews 1.1, 73-90 (1924). " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index 7991e706561..1779def85d6 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -65,6 +65,7 @@ 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) +tutorial_test(FILE test_grand_canonical_monte_carlo.py) add_custom_target( check_tutorials diff --git a/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py b/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py new file mode 100644 index 00000000000..8abf7f3815d --- /dev/null +++ b/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py @@ -0,0 +1,43 @@ +# +# 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 pint +import numpy as np + +ureg = pint.UnitRegistry() +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/grand_canonical_monte_carlo/grand_canonical_monte_carlo.py", + salt_concentration_magnitudes_si=[0.003, 0.01], number_of_loops_for_each_concentration=[250, 100], + excess_chemical_potential_data=[-0.123732052028611, -0.218687259792629], + excess_chemical_potential_data_error=[0.00152160176511698, 0.00220162667953136]) + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + + def test(self): + sim_xi_minus = tutorial.partition_coefficients_negatives_array + self.assertLess(np.abs(sim_xi_minus[0] - 0.05), 0.02) + self.assertLess(np.abs(sim_xi_minus[1] - 0.15), 0.1) + + +if __name__ == "__main__": + ut.main() From 9ac5c463e5e9f128d22f341840cbe760a6cd58a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 6 Nov 2023 19:38:14 +0100 Subject: [PATCH 2/2] Modify GCMC tutorial Reduce sampling times. Improve testing. Redesign plots. --- .../CMakeLists.txt | 9 +- .../NotesForTutor.md | 4 +- .../grand_canonical_monte_carlo.ipynb | 247 ++++++++++-------- .../test_grand_canonical_monte_carlo.py | 15 +- 4 files changed, 158 insertions(+), 117 deletions(-) diff --git a/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt b/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt index 895254c5e43..a1981aabe09 100644 --- a/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt +++ b/doc/tutorials/grand_canonical_monte_carlo/CMakeLists.txt @@ -17,9 +17,10 @@ # along with this program. If not, see . # -configure_tutorial_target(TARGET tutorial_grand_canonical_monte_carlo DEPENDS - grand_canonical_monte_carlo.ipynb - figures/schematic.svg) +configure_tutorial_target( + TARGET tutorial_grand_canonical_monte_carlo DEPENDS + grand_canonical_monte_carlo.ipynb figures/schematic.svg) nb_export(TARGET tutorial_grand_canonical_monte_carlo SUFFIX "" FILE - "grand_canonical_monte_carlo.ipynb" HTML_RUN) + "grand_canonical_monte_carlo.ipynb" HTML_RUN VAR_SUBST + "\"p3m_params={'mesh':10,'cao':6,'r_cut':8.22}\"") diff --git a/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md b/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md index 098aeae3a31..618949fff21 100644 --- a/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md +++ b/doc/tutorials/grand_canonical_monte_carlo/NotesForTutor.md @@ -5,10 +5,10 @@ After the tutorial, students should be able to explain: * what the grand-canonical ensemble is and where it is applicable -* how to simulate in grand-canonical ensemble using GCMC +* how to simulate in the grand-canonical ensemble using GCMC ## ESPResSo learning goals In the course of this tutorial, students should learn to: -* how to set up a GCMC simulation in ESPResSo +* to set up a GCMC simulation in ESPResSo diff --git a/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb b/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb index 56376895e15..41552acb3a0 100644 --- a/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb +++ b/doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb @@ -32,30 +32,42 @@ "## Grand-Canonical Ensemble\n", "The case of changing particle numbers is best described by the grand canonical ensemble [1].\n", "In the grand canonical ensemble, the volume $V$, the temperature $T$ and the chemical potentials $\\mu_i$ of the exchangeable species $i$ rather than the particle numbers $N_i$ are fixed. The chemical potentials $\\mu_i$ correspond to the free energy (Helmholtz free energy $F$, will be defined later) cost of adding a particle of type $i$ while keeping everything else fixed:\n", - "\\begin{align}\n", + "\n", + "$$\n", "\\mu_i = \\left(\\frac{\\partial F}{\\partial N_i}\\right)_{T, V, \\left\\{N_j\\right\\}_{j\\neq i}}.\n", - "\\end{align}\n", + "$$\n", + "\n", "The grand canonical ensemble is thus also referred to as the $(\\mu, V, T)$-ensemble.\n", "In the grand canonical ensemble, the probability for the system to contain $N_i$ particles of type $i$, i.e. a total of $N=\\sum_{i=1}^{N_\\mathrm{s}}N_i$ ($N_\\mathrm{s}$ is the number of different types) and to be in a specific microstate $\\xi\\in\\Gamma_{N}$ is given by the probability distribution\n", - "\\begin{align}\n", + "\n", + "$$\n", "p^\\mathrm{G}\\left(\\left\\{N_i\\right\\},\\xi\\right) = \\frac{\\exp\\left(-\\beta \\left(H(\\xi)-\\sum_{i=1}^{N_\\mathrm{s}}\\mu_iN_i\\right)\\right)}{Z^\\mathrm{G}h^{3N}\\prod_{i=1}^{N_\\mathrm{s}} N_i!}. \\quad \\tag{1}\n", - "\\end{align}\n", + "$$\n", + "\n", "Here, $H$ denotes the Hamiltonian of the system with $N_i$ particles of type $i$, $h$ is the planck constant, $\\beta=1/k_\\mathrm{B}T$ the inverse temperature and $Z^\\mathrm{G}$ is the grand-canonical partition function which is given by\n", - "\\begin{align}\n", + "\n", + "$$\n", "Z^\\mathrm{G}\\left(\\left\\{\\mu_i\\right\\},V,T\\right) = \\sum_{N_1=0}^{\\infty}...\\sum_{N_\\mathrm{s}=0}^{\\infty}Z\\left(\\left\\{N_i\\right\\},V,T\\right)\\exp\\left(\\beta \\sum_{i=1}^{N_\\mathrm{s}}\\mu_iN_i\\right),\n", - "\\end{align}\n", + "$$\n", + "\n", "with the canonical partition function\n", - "\\begin{align}\n", + "\n", + "$$\n", "\tZ\\left(\\left\\{N_i\\right\\},V,T\\right) = \\frac{1}{h^{3N} \\prod_{i=1}^{N_s}N_i!}\\int_{\\Gamma_N} \\mathrm{d}\\xi \\ \\exp(-\\beta H(\\xi)).\n", - "\\end{align}\n", + "$$\n", + "\n", "The partition function is connected to a thermodynamic potential, called the grand potential (Landau free energy):\n", - "\\begin{align}\n", + "\n", + "$$\n", " \\Omega\\left(\\left\\{\\mu_i\\right\\},V,T\\right) = -k_\\mathrm{B}T\\ln\\left(Z^\\mathrm{G}\\left(\\left\\{\\mu_i\\right\\},V,T\\right)\\right).\n", - "\\end{align}\n", + "$$\n", + "\n", "In a similiar way, the Helmholtz free energy is defined as\n", - "\\begin{align}\n", + "\n", + "$$\n", "F\\left(\\left\\{N_i\\right\\},V,T\\right) = -k_\\mathrm{B}T\\ln\\left(Z\\left(\\left\\{N_i\\right\\},V,T\\right)\\right).\n", - "\\end{align}\n", + "$$\n", + "\n", "\n", "When considering a system that can exchange particles with a reservoir, the chemical potential $\\mu_i$ for each exchangeable species $i$ has to be equal in the reservoir and in the system:\n", "\\begin{equation}\n", @@ -74,64 +86,83 @@ "The grand-canonical distribution (see eq. (1)) can be sampled using the Metropolis Monte Carlo method, similiar to the canonical ensemble. Before we derive the acceptance criterion for the GCMC method, let us shortly recall the generic Metropolis-Hastings algorithm.\n", "\n", "We want to calculate averages of observables $A$ which are distributed according to some distribution $p^\\mathrm{eq}_s$ (e.g. the canonical or grand-canonical distribution)\n", - "\\begin{equation}\n", + "\n", + "$$\n", "\\langle A\\rangle = \\frac{\\sum_{s\\in\\Gamma} A_s p^\\mathrm{eq}_s}{\\sum_{s\\in\\Gamma} p^\\mathrm{eq}_s}. \\quad \\tag{2}\n", "\\label{eq:can_average_mc}\n", - "\\end{equation}\n", + "$$\n", + "\n", "Here, $s$ labels the different states of the systems, $A_s$ is the value of the observable $A$ in the state $s$ and the sum runs over the set of all states $\\Gamma$.\n", "The sum can also be an integral for systems with continuous degrees of freedom, for this derivation it does not matter.\n", "For almost all systems, the ensemble average (2) can not be evaluated exactly and only numerical approximations can be obtained.\n", "The basic idea of all MC methods is to randomly sample only a subset of states.\n", "In the most naive approach (simple sampling), states are simply selected with a uniform probability from the space $\\Gamma$ of states $s$, resulting in a sample $\\gamma\\subset\\Gamma$.\n", "The canonical average can then be approximated by a sum over $\\gamma$:\n", - "\\begin{equation}\n", + "\n", + "$$\n", "\\langle A\\rangle \\approx \\frac{\\sum_{s\\in\\gamma} A_s p^\\mathrm{eq}_s}{\\sum_{s\\in\\gamma} p^\\mathrm{eq}_s} \\quad \\tag{3}\n", - "\\end{equation}\n", + "$$\n", + "\n", "This method is highly inefficient, as most states will typically not contribute significantly to the average.\n", "To get a more efficient sampling method, let us select the states $s$ according to a non-uniform probability distribution $p_s$.\n", "Then, the average can be written as \n", - "\\begin{equation}\n", + "\n", + "$$\n", "\\langle A\\rangle^{\\mathrm{can}} \\approx \\frac{\\sum_{s\\in\\gamma} A_s p^\\mathrm{eq}_s/p_s}{\\sum_{s\\in\\gamma} p^\\mathrm{eq}_s/p_s}, \\quad \\tag{4}\n", "\\label{eq:can_average_mc_prob}\n", - "\\end{equation}\n", + "$$\n", + "\n", "where we have to divide in the numerator and in the denominator by $p_s$ to compensate for the fact that the samples are distributed non-uniformly.\n", "It is obvious that we recover (3) from (4) as the special case of a uniform distribution. \n", "To sample the averages more efficiently, we simply pick\n", - "\\begin{equation}\n", + "\n", + "$$\n", "p_s \\propto p^\\mathrm{eq}_s \\quad \\tag{5}\n", "\\label{eq:mc_prob}\n", - "\\end{equation}\n", + "$$\n", + "\n", "such that the selected samples are representative of the desired equilibrium distribution.\n", "In this case, the estimators of the averages are thus simply given by\n", - "\\begin{equation}\n", + "\n", + "$$\n", "\\langle A\\rangle \\approx \\frac{1}{N}\\sum_{s\\in\\gamma}A_s.\n", - "\\end{equation}\n", + "$$\n", "The remaining problem is now to find a method which allows us to sample states according to (5).\n", "This can be done in terms of so-called Markov chains. \n", "A Markov chain is a stochastic process without memory.\n", "In our case, the Markov chain is simply a sequence $\\{s_i\\}_{i=1,...N}$ of $N$ states $s_i$ of the system and the fact that the Markov chain has no memory means that the transition probability from a state $s_i$ to a state $s_j$ can be given in terms of a transition matrix $W_{i\\rightarrow j}$ which is constant.\n", "The equilibrium probability distribution has to be invariant under $W$. Mathematically this condition can be expressed in the form\n", - "\\begin{align}\n", + "\n", + "$$\n", "\\sum_{i}W_{i\\rightarrow j}p_i^\\mathrm{eq} = p_j^\\mathrm{eq}\n", - "\\end{align}\n", + "$$\n", + "\n", "A sufficient constraint on $W_{i\\rightarrow j}$ which ensures the validity of this condition\n", "is the so-called criterion of detailed balance:\n", - "\\begin{align}\n", + "\n", + "$$\n", "\\frac{W_{i\\rightarrow j}}{W_{j\\rightarrow i}} = \\frac{p_j^\\mathrm{eq}} {p_i^\\mathrm{eq}}.\n", - "\\end{align}\n", + "$$\n", + "\n", "Typically, the transition matrix elements are factorized into a proposal probability $g_{i\\rightarrow j}$ and acceptance probability $P_{i\\rightarrow j}^\\mathrm{acc}$:\n", - "\\begin{align}\n", + "\n", + "$$\n", "W_{i\\rightarrow j} = g_{i\\rightarrow j} P_{i\\rightarrow j}^\\mathrm{acc}. \n", - "\\end{align}\n", + "$$\n", + "\n", "For the symmetric choice $g_{i\\rightarrow j}=g_{j\\rightarrow i}$, the condition of detailed balance reduces to \n", - "\\begin{align}\n", + "\n", + "$$\n", "\\frac{P_{i\\rightarrow j}^\\mathrm{acc}}{P_{j\\rightarrow i}^\\mathrm{acc}} = \\frac{p_j^\\mathrm{eq}}{p_i^\\mathrm{eq}}. \n", "\\quad \\tag{6}\n", - "\\end{align}\n", + "$$\n", + "\n", "By checking the cases $p_j^\\mathrm{eq}>p_i^\\mathrm{eq}$ and $p_j^\\mathrm{eq} \\leq p_i^\\mathrm{eq}$ one can see that the acceptance probability \n", - "\\begin{align}\n", + "\n", + "$$\n", "P_{i\\rightarrow j}^\\mathrm{acc,\\text{ }Metropolis} = \\min\\left(1,\\frac{p_j^\\mathrm{eq}}{p_i^\\mathrm{eq}}\\right) \\quad \\tag{7}\n", - "\\end{align}\n", + "$$\n", + "\n", "fulfills equation (6). With this, we can now formulate the Metropolis-Hastings algorithm [2].\n", "We start with an arbitrary initial state $s_1$ and generate a sequence $\\{s_i\\}_{i=1,...N}$ of $N$ states.\n", "Each new state (n) is generated from the previous state (o) in the following way. \n", @@ -147,18 +178,21 @@ "2. Insertion and Deletion moves that add or delete a single particle while keeping all other particle positions fixed.\n", "\n", "In the GCMC method the insertion and deletion moves are implemented in the following fashion. First of all it is determined randomly with equal probability if an insertion or deletion move is performed. In the case of an insertion the additional particle is placed at a uniformly chosen random position in the simulation box. Then the proposed new state is accepted according to the criterion\n", - "\\begin{align}\n", + "\n", + "$$\n", "\\begin{split}\n", "P_{N_i\\rightarrow N_i+1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{1}{N_i+1}\\left(\\frac{V}{\\Lambda_i^3}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i+1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)+\\beta \\mu_i\\right)\\right).\n", "\\end{split}\n", - "\\end{align}\n", + "$$\n", + "\n", "In the case of a deletion a randomly selected particle is removed from the system. \n", "We have the analogous result\n", - "\\begin{align}\n", + "\n", + "$$\n", "\\begin{split}\n", "P_{N_i\\rightarrow N_i-1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,N_i\\left(\\frac{\\Lambda_i^3}{V}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i-1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)-\\beta \\mu_i\\right)\\right).\n", "\\end{split}\n", - "\\end{align}" + "$$" ] }, { @@ -175,9 +209,9 @@ "**Hint**\n", "\n", "* Use the grand-canonical distribution with all the momenta integrated out:\n", - " \\begin{align}\n", + "$$\n", "p^\\mathrm{eq}_{N} \\propto \\exp\\left(-\\beta U(\\mathbf{q}) \\right ) \\prod_i \\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp \\left (\\beta \\mu_i N_i\\right)\n", - "\\end{align}\n", + "$$\n", "This alternative expression is valid if we are only interested in observables which do not depend on the momenta." ] }, @@ -188,11 +222,14 @@ }, "source": [ "For the case ($N_i\\rightarrow N_i + 1$) one gets\n", + "\n", "\\begin{align}\n", "P_{N_i\\rightarrow N+1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{p^\\mathrm{eq}_{N_i+1}}{p^\\mathrm{eq}_{N_i}}\\right) = \\min\\left(1,\\frac{\\frac{1}{(N_i+1)!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i+1}\\times \\exp\\left(-\\beta U^{N_i+1}(\\mathbf{q})+\\beta \\mu (N_i+1)\\right)}{\\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp\\left(-\\beta U^{N_i}(\\mathbf{q})+\\beta \\mu N_i\\right)}\\right)\\\\\n", "&= \\min\\left(1,\\frac{1}{N_i+1}\\left(\\frac{V}{\\Lambda_i^3}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i+1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)+\\beta \\mu_i\\right)\\right).\n", "\\end{align}\n", + "\n", "For the other case ($N_i\\rightarrow N_i-1$) we have the analogous result\n", + "\n", "\\begin{align}\n", "P_{N_i\\rightarrow N_i-1}^\\mathrm{acc,\\text{ }GCMC} &= \\min\\left(1,\\frac{p^\\mathrm{eq}_{N_i-1}}{p^\\mathrm{eq}_{N_i}}\\right) = \\min\\left(1,\\frac{\\frac{1}{(N_i-1)!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i-1}\\times \\exp\\left(-\\beta U^{N_i-1}(\\mathbf{q})+\\beta \\mu_i (N_i-1)\\right)}{\\frac{1}{N_i!}\\left(\\frac{V}{\\Lambda_i^3}\\right)^{N_i}\\times \\exp\\left(-\\beta U^{N_i}(\\mathbf{q})+\\beta \\mu_i N_i\\right)}\\right)\\\\\n", "&= \\min\\left(1,N_i\\left(\\frac{\\Lambda_i^3}{V}\\right)\\times \\exp\\left(-\\beta \\left(U^{N_i-1}(\\mathbf{q})-U^{N_i}(\\mathbf{q})\\right)-\\beta \\mu_i\\right)\\right).\n", @@ -210,31 +247,41 @@ "
Fig. 1: Schematic representation of the considered two phase setup.
\n", "\n", "\n", + "\n", "Here we want to investigate a polyelectrolyte solution that is coupled to a reservoir containing a monovalent salt, as shown schematically in figure 1. Similar systems are useful to investigate for example the swelling behavior of hydrogels. The surrounding of the gel can be considered as the reservoir and the inner of the gel is the system. The small salt ions (X$^+$ and X$^-$) can be exchanged between the two phases while the polyelectrolyte is confined to the system. Experimentally such a setup could be realized by seperating the two phases with a semi-permeable membrane. Similiar setups also occur in the study of polyelectrolyte hydrogels.\n", "\n", - "Due to the macroscopic electroneutrality of both phases the salt ions obey an electrochemical equilibrium, i. e. there is a difference in the electrostatic potential between the two phases $$\\mu_i^{\\mathrm{sys}}+z_ie\\psi^{\\mathrm{Donnan}} = \\mu_i^{\\mathrm{res}}.$$\n", - "This effect is called the Donnan effect [3] and leads to an unequal partioning of salt between the reservoir and the system. The partitioning can be quantified using the partition coeffiecents $\\xi_i \\equiv c^{\\mathrm{sys}}_i/c^{\\mathrm{res}}_i$, where $c^{\\mathrm{sys}}_i$ and $c^{\\mathrm{res}}_i$ are the particle concentrations of type $i$ in the system respectively in the system. Defining the \"universal\" partition coefficient,\n", + "Due to the macroscopic electroneutrality of both phases the salt ions obey an electrochemical equilibrium, i. e. there is a difference in the electrostatic potential between the two phases\n", + "\n", + "$$\n", + "\\mu_i^{\\mathrm{sys}}+z_ie\\psi^{\\mathrm{Donnan}} = \\mu_i^{\\mathrm{res}}.\n", + "$$\n", + "\n", + "This effect is called the Donnan effect [3] and leads to an unequal partioning of salt between the reservoir and the system. The partitioning can be quantified using the partition coefficients $\\xi_i \\equiv c^{\\mathrm{sys}}_i/c^{\\mathrm{res}}_i$, where $c^{\\mathrm{sys}}_i$ and $c^{\\mathrm{res}}_i$ are the particle concentrations of type $i$ in the system and reservoir, respectively. Defining the \"universal\" partition coefficient,\n", + "\n", "\\begin{align*}\n", " \\xi\\equiv\\left\\{\n", " \\begin{array}{ll}\n", - " \\xi_+ - \\frac{c_{\\text{M}^-}^{\\text{sys}}}{c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}& \\text{for positive ions}\\\\\n", - " \\xi_- & \\text{for negative ions,}\\\\\n", + " \\displaystyle\\xi_+ - \\frac{c_{\\text{M}^-}^{\\text{sys}}}{c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}& \\text{for positive ions}\\\\\n", + " \\displaystyle\\xi_- & \\text{for negative ions,}\\\\\n", " \\end{array}\\right.\n", "\\end{align*}\n", + "\n", "one can obtain the following solution for an ideal system without any interaction between the particles:\n", - "\\begin{align}\n", - " \\xi^{\\text{Donnan}} = -\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}} + \\sqrt{\\left(\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}\\right)^2+1}.\n", - " \\end{align}\n", - " This solution shows that salt is rejected by the polyelctrolyte solution. In the following we want to investigate this partitioning for a system with electrostatic interactions.\n", + "\n", + "$$\n", + "\\xi^{\\text{Donnan}} = -\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}} + \\sqrt{\\left(\\frac{c_{\\text{M}^-}^{\\text{sys}}}{2c_{\\text{X}^+,\\text{X}^-}^{\\text{res}}}\\right)^2+1}.\n", + "$$\n", + "\n", + "This solution shows that salt is rejected by the polyelctrolyte solution. In the following we want to investigate this partitioning for a system with electrostatic interactions.\n", " \n", - "Because of the electroneutrality condition the two definitions of the \"universal\" partition coefficients for positive and negative ions are equal. The additional term in the definition for positive ions is a correction due to the counterions in the system, which neutralize the charged polyelectrolytes." + "Because of the electroneutrality condition, the two definitions of the \"universal\" partition coefficients for positive and negative ions are equal. The additional term in the definition for positive ions is a correction due to the counterions in the system, which neutralize the charged polyelectrolytes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the folllowing, we plot the universal parition coefficient over the ratio of the concentrations." + "In the folllowing, we plot the universal partition coefficient over the ratio of the concentrations." ] }, { @@ -245,6 +292,7 @@ "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "np.random.seed(42)\n", "plt.rcParams.update({\"font.size\": 18})\n", "\n", "ratios = np.logspace(-1, 2, 1000)\n", @@ -329,16 +377,16 @@ "metadata": {}, "outputs": [], "source": [ - "salt_concentration_magnitudes_si = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3]\n", + "salt_concentration_magnitudes_si = [0.001, 0.003, 0.01, 0.03, 0.1]\n", "salt_concentration_si = np.array(salt_concentration_magnitudes_si) * ureg.molar\n", "salt_concentration_sim = salt_concentration_si.to(\"molecules/reduced_length^3\")\n", "excess_chemical_potential_data = [\n", " -0.0672291585811938, -0.123732052028611, -0.218687259792629,\n", - " -0.326207586236404, -0.510091568808112, -0.657928029835777\n", + " -0.326207586236404, -0.510091568808112\n", "]\n", "excess_chemical_potential_data_error = [\n", " 0.000994798843587, 0.00152160176511698, 0.00220162667953136,\n", - " 0.0029299206854553, 0.00422309102221265, 0.00603137482524256\n", + " 0.0029299206854553, 0.00422309102221265\n", "]\n", "excess_chemical_potential_monovalent_pairs_in_bulk = scipy.interpolate.interp1d(\n", " salt_concentration_sim.magnitude, excess_chemical_potential_data)\n", @@ -418,10 +466,10 @@ "\n", "warmup_loops = 50\n", "number_of_loops_for_each_concentration = [\n", - " 750, 750, 250, 100, 100, 100\n", - "] # set it higher for better accuracy\n", - "steps_per_loop = 1000\n", - "samples_per_loop = 100\n", + " 200, 200, 150, 150, 100\n", + "] # set these higher for better accuracy\n", + "steps_per_loop = 200\n", + "samples_per_loop = 25\n", "\n", "##### Seeding\n", "langevin_seed = 42\n", @@ -487,8 +535,22 @@ "\n", "print(\"Added WCA interaction\")\n", "\n", - "p3m = electrostatics.P3M(prefactor=L_BJERRUM * KT, accuracy=1e-3)\n", + "##### Minimize forces\n", + "system.integrator.set_steepest_descent(f_max=N_plus, gamma=10.0,\n", + " max_displacement=LJ_SIGMA / 100.)\n", + "n_steps = system.integrator.run(steps_per_loop)\n", + "assert n_steps < steps_per_loop\n", + "\n", + "print(\"Removed overlaps with steepest descent...\")\n", + "\n", + "##### Add thermostat and long-range solvers\n", + "system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=langevin_seed)\n", + "system.integrator.set_vv()\n", + "\n", + "p3m_params = {\"mesh\": 10}\n", + "p3m = espressomd.electrostatics.P3M(prefactor=L_BJERRUM * KT, accuracy=1e-3, **p3m_params)\n", "system.electrostatics.solver = p3m\n", + "system.integrator.run(2 * steps_per_loop)\n", "\n", "print(\"Added electrostatics\")" ] @@ -502,13 +564,16 @@ "source": [ "Now we want to add a coupling to the reservoir. We can formally represent the insertion and deletion of ion pairs as a chemical reaction:\n", "\n", - "\\begin{align}\n", - "\\emptyset &\\Leftrightarrow \\mathrm{X}^- + \\mathrm{X}^+\n", - "\\end{align}\n", + "$$\n", + "\\emptyset \\rightleftharpoons \\mathrm{X}^- + \\mathrm{X}^+\n", + "$$\n", + "\n", "It is important to note that the reaction can run in both directions, i. e. insertion and deletion of salt ions. These reactions are described by a concentration based equilibrium constant:\n", - "\\begin{align}\n", - " K = c_{\\mathrm{salt,res}}^2 \\exp \\left ( \\mu_{\\mathrm{salt,res}}^{\\mathrm{ex}} \\right )\n", - "\\end{align}\n", + "\n", + "$$\n", + "K = c_{\\mathrm{salt,res}}^2 \\exp \\left ( \\mu_{\\mathrm{salt,res}}^{\\mathrm{ex}} \\right )\n", + "$$\n", + "\n", "Here $\\mu_{\\mathrm{salt,res}}^{\\mathrm{ex}}$ denotes the excess chemical potential of a salt ion pair as we measured in the [Widom insertion tutorial](https://espressomd.github.io/tutorials/widom_insertion/widom_insertion.html).\n", "\n", "**Exercise**\n", @@ -555,33 +620,6 @@ "outputs": [], "source": [] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we perform the initial warm up." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "##### Warmup\n", - "system.integrator.set_steepest_descent(f_max=0,\n", - " gamma=10.0,\n", - " max_displacement=0.001)\n", - "print(\"Removing overlaps with steepest descent...\")\n", - "for i in tqdm(range(warmup_loops),\n", - " bar_format='{l_bar}{bar:60}{r_bar}{bar:-60b}'):\n", - " system.integrator.run(steps_per_loop)\n", - "print(\"Done.\")\n", - "\n", - "system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=langevin_seed)\n", - "system.integrator.set_vv()" - ] - }, { "cell_type": "markdown", "metadata": { @@ -696,18 +734,21 @@ "metadata": {}, "outputs": [], "source": [ - "fig = plt.figure(figsize=(10, 6))\n", - "plt.loglog(c_monomer.magnitude/(2*salt_concentration_si.magnitude), partition_coefficients_negatives_array,\\\n", - " label = r'$\\xi$ from $\\xi_-$',linestyle = 'none',marker = 's', markersize=10)\n", "universal_partion_coefficient_positive = \\\n", - " partition_coefficients_positives_array - c_monomer.magnitude/salt_concentration_si.magnitude\n", - "plt.loglog(c_monomer.magnitude/(2*salt_concentration_si.magnitude),universal_partion_coefficient_positive ,\\\n", - " label = r'$\\xi$ from $\\xi_+$',marker = 'd', markersize=7, linestyle = 'none', color=\"lime\")\n", + " partition_coefficients_positives_array - c_monomer.magnitude / salt_concentration_si.magnitude\n", + "fig = plt.figure(figsize=(10, 6))\n", "ratios = np.logspace(-1, 2, num=200)\n", - "plt.loglog(ratios, analytical_solution(ratios), color=\"black\")\n", - "plt.xlabel(r'$\\frac{c_{\\mathrm{M}^-}}{2c_{\\mathrm{salt}}^{\\mathrm{res}}}$')\n", - "plt.ylabel(r'$\\xi$')\n", - "plt.legend();" + "plt.loglog(ratios, analytical_solution(ratios), color=\"black\", label=r\"$\\xi^{\\mathrm{Donnan}}$\")\n", + "plt.loglog(c_monomer.magnitude / (2 * salt_concentration_si.magnitude), partition_coefficients_negatives_array,\n", + " marker=\"o\", markersize=10, markeredgewidth=1.5, markeredgecolor=\"k\", markerfacecolor=\"none\",\n", + " linestyle=\"none\", label=r\"$\\xi_-$\")\n", + "plt.loglog(c_monomer.magnitude / (2 * salt_concentration_si.magnitude), universal_partion_coefficient_positive,\n", + " marker=\"o\", markersize=4.5, markeredgewidth=0., markeredgecolor=\"k\", markerfacecolor=\"k\",\n", + " linestyle=\"none\", label=r\"$\\xi_+ - c_{\\mathrm{M}^-}^{\\mathrm{sys}}/c_{\\mathrm{salt}}^{\\mathrm{res}}$\")\n", + "plt.xlabel(r\"$\\dfrac{c_{\\mathrm{M}^-}^{\\mathrm{sys}}}{2c_{\\mathrm{salt}}^{\\mathrm{res}}}$\")\n", + "plt.ylabel(r\"$\\xi$\")\n", + "plt.legend()\n", + "plt.show()" ] }, { @@ -751,11 +792,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "[1] D. Frenkel, Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. 2nd edition, chapter 5.6: Grand-Canonical Ensemble. Academic Press, 2002, ISBN: 9780122673511.\n", - "\n", - "[2] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21, 1087 (1953); https://doi.org/10.1063/1.1699114. \n", - "\n", - "[3] F. Donnan. The theory of membrane equilibria. Chemical reviews 1.1, 73-90 (1924). " + "[1] Daan Frenkel, Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. 2nd edition, chapter 5: Monte Carlo Simulations in Various Ensembles, section 5.6: Grand-Canonical Ensemble, pp. 126–135. Academic Press, 2002, ISBN: 978-0-12-267351-1, doi:[10.1016/B978-012267351-1/50007-9](https://doi.org/10.1016/B978-012267351-1/50007-9). \n", + "[2] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller. Equation of state calculations by fast computing machines. *J. Chem. Phys.* 21, pp. 1087–1092 (1953), doi:[10.1063/1.1699114](https://doi.org/10.1063/1.1699114). \n", + "[3] F. Donnan. The theory of membrane equilibria. *Chemical Reviews* 1(1), pp. 73–90 (1924), doi:[10.1021/cr60001a003](https://doi.org/10.1021/cr60001a003). " ] } ], diff --git a/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py b/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py index 8abf7f3815d..7b4408d4643 100644 --- a/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py +++ b/testsuite/scripts/tutorials/test_grand_canonical_monte_carlo.py @@ -19,24 +19,25 @@ import unittest as ut import importlib_wrapper -import pint import numpy as np -ureg = pint.UnitRegistry() tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/grand_canonical_monte_carlo/grand_canonical_monte_carlo.py", - salt_concentration_magnitudes_si=[0.003, 0.01], number_of_loops_for_each_concentration=[250, 100], - excess_chemical_potential_data=[-0.123732052028611, -0.218687259792629], - excess_chemical_potential_data_error=[0.00152160176511698, 0.00220162667953136]) + p3m_params={"mesh": 10, "cao": 6, "r_cut": 8.22}) @skipIfMissingFeatures class Tutorial(ut.TestCase): def test(self): + ratios = tutorial.c_monomer.magnitude / \ + (2. * tutorial.salt_concentration_si.magnitude) + ref_xi = tutorial.analytical_solution(ratios) sim_xi_minus = tutorial.partition_coefficients_negatives_array - self.assertLess(np.abs(sim_xi_minus[0] - 0.05), 0.02) - self.assertLess(np.abs(sim_xi_minus[1] - 0.15), 0.1) + sim_xi_plus = tutorial.universal_partion_coefficient_positive + np.testing.assert_allclose( + sim_xi_minus, sim_xi_plus, rtol=1e-5, atol=1e-5) + np.testing.assert_allclose(sim_xi_minus / ref_xi, 2., rtol=0., atol=2.) if __name__ == "__main__":