From 54611500971accee63bf132134c4e122b9d5f53a Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 15 Oct 2019 13:57:29 +0200 Subject: [PATCH 01/11] core: Move ParticleList declaration to separate file --- src/core/ParticleList.hpp | 24 ++++++++++++++++++++++++ src/core/particle_data.hpp | 21 +-------------------- 2 files changed, 25 insertions(+), 20 deletions(-) create mode 100644 src/core/ParticleList.hpp diff --git a/src/core/ParticleList.hpp b/src/core/ParticleList.hpp new file mode 100644 index 00000000000..5973689728b --- /dev/null +++ b/src/core/ParticleList.hpp @@ -0,0 +1,24 @@ +#ifndef ESPRESSO_CORE_PARTICLE_LIST_HPP +#define ESPRESSO_CORE_PARTICLE_LIST_HPP + +#include "Particle.hpp" + +#include + +/** List of particles. The particle array is resized using a sophisticated + * (we hope) algorithm to avoid unnecessary resizes. + * Access using \ref realloc_particlelist, ... + */ +struct ParticleList { + ParticleList() : part{nullptr}, n{0}, max{0} {} + /** The particles payload */ + Particle *part; + /** Number of particles contained */ + int n; + /** Number of particles that fit in until a resize is needed */ + int max; + + Utils::Span particles() { return {part, static_cast(n)}; } +}; + +#endif diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index f5cac4eb0ad..0f11856bdb9 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -37,6 +37,7 @@ #include "config.hpp" #include "Particle.hpp" +#include "ParticleList.hpp" #include #include @@ -77,26 +78,6 @@ enum { #endif -/************************************************ - * data types - ************************************************/ - -/** List of particles. The particle array is resized using a sophisticated - * (we hope) algorithm to avoid unnecessary resizes. - * Access using \ref realloc_particlelist, ... - */ -struct ParticleList { - ParticleList() : part{nullptr}, n{0}, max{0} {} - /** The particles payload */ - Particle *part; - /** Number of particles contained */ - int n; - /** Number of particles that fit in until a resize is needed */ - int max; - - Utils::Span particles() { return {part, static_cast(n)}; } -}; - /************************************************ * exported variables ************************************************/ From 5ddda5686eb66a770c61b3f11b8c2f2d18b38967 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 15 Oct 2019 15:28:46 +0200 Subject: [PATCH 02/11] core: Narrow includes --- src/core/Cell.hpp | 3 +- src/core/PartCfg.hpp | 2 ++ src/core/ParticleList.hpp | 36 +++++++++++++++++++ src/core/cells.cpp | 1 + src/core/collision.cpp | 5 +-- src/core/cuda_common_cuda.cu | 1 + src/core/debug.cpp | 1 + src/core/domain_decomposition.cpp | 1 + .../magnetic_non_p3m_methods.cpp | 1 + .../mdlc_correction.cpp | 1 + src/core/forces_inline.hpp | 1 + src/core/ghosts.cpp | 1 + src/core/global.cpp | 2 ++ .../lb_particle_coupling.cpp | 1 + .../immersed_boundary/ImmersedBoundaries.cpp | 1 + src/core/integrators/steepest_descent.cpp | 1 + src/core/integrators/velocity_verlet_npt.cpp | 5 +-- src/core/io/mpiio/mpiio.cpp | 1 + .../object-in-fluid/oif_global_forces.cpp | 1 + src/core/particle_data.cpp | 35 ++++-------------- src/core/particle_data.hpp | 8 ----- src/core/rattle.cpp | 1 + 22 files changed, 68 insertions(+), 42 deletions(-) diff --git a/src/core/Cell.hpp b/src/core/Cell.hpp index d53905b80ab..fe1c68bf15d 100644 --- a/src/core/Cell.hpp +++ b/src/core/Cell.hpp @@ -19,7 +19,8 @@ #ifndef CORE_CELL_HPP #define CORE_CELL_HPP -#include "particle_data.hpp" +#include "Particle.hpp" +#include "ParticleList.hpp" #include diff --git a/src/core/PartCfg.hpp b/src/core/PartCfg.hpp index cc05ecf73fa..89784a14418 100644 --- a/src/core/PartCfg.hpp +++ b/src/core/PartCfg.hpp @@ -23,6 +23,8 @@ #include "ParticleCache.hpp" #include "cells.hpp" #include "grid.hpp" +#include "particle_data.hpp" + #include "serialization/Particle.hpp" #include diff --git a/src/core/ParticleList.hpp b/src/core/ParticleList.hpp index 5973689728b..8693249ce57 100644 --- a/src/core/ParticleList.hpp +++ b/src/core/ParticleList.hpp @@ -4,6 +4,7 @@ #include "Particle.hpp" #include +#include /** List of particles. The particle array is resized using a sophisticated * (we hope) algorithm to avoid unnecessary resizes. @@ -19,6 +20,41 @@ struct ParticleList { int max; Utils::Span particles() { return {part, static_cast(n)}; } + + int resize(int size) { +/** granularity of the particle buffers in particles */ + constexpr int PART_INCREMENT = 8; + + assert(size >= 0); + int old_max = max; + Particle *old_start = part; + + if (size < max) { + if (size == 0) + /* to be able to free an array again */ + max = 0; + else + /* shrink not as fast, just lose half, rounded up */ + max = + PART_INCREMENT * + (((max + size + 1) / 2 + PART_INCREMENT - 1) / PART_INCREMENT); + } else + /* round up */ + max = PART_INCREMENT * ((size + PART_INCREMENT - 1) / PART_INCREMENT); + if (max != old_max) + part = Utils::realloc(part, sizeof(Particle) * max); + return part != old_start; + } }; +/** Allocate storage for local particles and ghosts. This version + does \em not care for the bond information to be freed if necessary. + \param plist the list on which to operate + \param size the size to provide at least. It is rounded + up to multiples of \ref PART_INCREMENT. + \return true iff particle addresses have changed */ +inline int realloc_particlelist(ParticleList *l, int size) { + return assert(l), l->resize(size); +} + #endif diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 1638cbd2435..4fb7e100003 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -38,6 +38,7 @@ #include "layered.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "nsquare.hpp" +#include "particle_data.hpp" #include #include diff --git a/src/core/collision.cpp b/src/core/collision.cpp index bf7045d6cfa..90666eaf540 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -18,6 +18,8 @@ */ #include "collision.hpp" + +#ifdef COLLISION_DETECTION #include "Particle.hpp" #include "cells.hpp" #include "communication.hpp" @@ -27,6 +29,7 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "rotation.hpp" #include "virtual_sites/VirtualSitesRelative.hpp" +#include "particle_data.hpp" #include #include @@ -37,8 +40,6 @@ #include -#ifdef COLLISION_DETECTION - /// Data type holding the info about a single collision typedef struct { int pp1; // 1st particle id diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index b34db0295fb..bc29fde74f6 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -27,6 +27,7 @@ #include "cuda_utils.hpp" #include "errorhandling.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "particle_data.hpp" #include diff --git a/src/core/debug.cpp b/src/core/debug.cpp index b9e6b186511..2132ec1fb1b 100644 --- a/src/core/debug.cpp +++ b/src/core/debug.cpp @@ -29,6 +29,7 @@ #include "cells.hpp" #include "errorhandling.hpp" #include "grid.hpp" +#include "particle_data.hpp" void check_particle_consistency() { int n, c; diff --git a/src/core/domain_decomposition.cpp b/src/core/domain_decomposition.cpp index dfd08cb0e8c..f8ce2a35a5b 100644 --- a/src/core/domain_decomposition.cpp +++ b/src/core/domain_decomposition.cpp @@ -28,6 +28,7 @@ #include "communication.hpp" #include "errorhandling.hpp" #include "grid.hpp" +#include "particle_data.hpp" #include "serialization/ParticleList.hpp" #include diff --git a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp index 7e239a6f35d..7cfe5c4597f 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp @@ -29,6 +29,7 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "particle_data.hpp" #include diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp index c59496e8531..b1d9c919c23 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp @@ -44,6 +44,7 @@ #include "errorhandling.hpp" #include "global.hpp" #include "grid.hpp" +#include "particle_data.hpp" DLC_struct dlc_params = {1e100, 0, 0, 0, 0}; diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index c6fc53365d4..06ebfc8e297 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -61,6 +61,7 @@ #include "object-in-fluid/oif_local_forces.hpp" #include "object-in-fluid/out_direction.hpp" #include "thermostat.hpp" +#include "particle_data.hpp" #ifdef DIPOLES #include "electrostatics_magnetostatics/dipole_inline.hpp" diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 057699a9e57..93c44266cbb 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -28,6 +28,7 @@ #include "Particle.hpp" #include "communication.hpp" #include "errorhandling.hpp" +#include "particle_data.hpp" #include #include diff --git a/src/core/global.cpp b/src/core/global.cpp index e749e70a344..194dd7bb992 100644 --- a/src/core/global.cpp +++ b/src/core/global.cpp @@ -37,6 +37,8 @@ #include "rattle.hpp" #include "thermostat.hpp" #include "tuning.hpp" +#include "particle_data.hpp" + #include #include diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index c4633855772..d399187e4b0 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -28,6 +28,7 @@ #include "lb_interpolation.hpp" #include "lbgpu.hpp" #include "random.hpp" +#include "particle_data.hpp" #include #include diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 7a80d56254f..f9d58a8b89f 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -25,6 +25,7 @@ #include "communication.hpp" #include "errorhandling.hpp" #include "grid.hpp" +#include "particle_data.hpp" #include diff --git a/src/core/integrators/steepest_descent.cpp b/src/core/integrators/steepest_descent.cpp index 581fd0ced08..c6d1f4d2b97 100644 --- a/src/core/integrators/steepest_descent.cpp +++ b/src/core/integrators/steepest_descent.cpp @@ -25,6 +25,7 @@ #include "event.hpp" #include "integrate.hpp" #include "rotation.hpp" +#include "particle_data.hpp" #include diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 3e7a0615c3b..1302c8be4a8 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -20,7 +20,6 @@ #include "config.hpp" #ifdef NPT -#include "Particle.hpp" #include "ParticleRange.hpp" #include "cells.hpp" #include "communication.hpp" @@ -29,7 +28,9 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" #include "thermostat.hpp" -#include "utils/math/sqr.hpp" +#include "particle_data.hpp" + +#include void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0; diff --git a/src/core/io/mpiio/mpiio.cpp b/src/core/io/mpiio/mpiio.cpp index 687a7b17a5c..a22ff999b03 100644 --- a/src/core/io/mpiio/mpiio.cpp +++ b/src/core/io/mpiio/mpiio.cpp @@ -57,6 +57,7 @@ #include "event.hpp" #include "integrate.hpp" #include "mpiio.hpp" +#include "particle_data.hpp" #include diff --git a/src/core/object-in-fluid/oif_global_forces.cpp b/src/core/object-in-fluid/oif_global_forces.cpp index 4571f8a6850..b3956b29d4a 100644 --- a/src/core/object-in-fluid/oif_global_forces.cpp +++ b/src/core/object-in-fluid/oif_global_forces.cpp @@ -25,6 +25,7 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "grid_based_algorithms/lb_interface.hpp" +#include "particle_data.hpp" #include using Utils::angle_btw_triangles; diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 98f4801d4cc..dba226c32b1 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -21,11 +21,11 @@ /** \file * Particles and particle lists. * - * The corresponding header file is Particle.hpp. + * The corresponding header file is particle_data.hpp. */ +#include "particle_data.hpp" #include "Particle.hpp" - #include "PartCfg.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "cells.hpp" @@ -58,9 +58,6 @@ * defines ************************************************/ -/** granularity of the particle buffers in particles */ -#define PART_INCREMENT 8 - /** my magic MPI code for send/recv_particles */ #define REQ_SNDRCV_PART 0xaa @@ -542,10 +539,12 @@ void clear_particle_node() { particle_node.clear(); } \param part the highest existing particle */ void realloc_local_particles(int part) { + constexpr auto INCREMENT = 8; + if (part >= max_local_particles) { - /* round up part + 1 in granularity PART_INCREMENT */ + /* round up part + 1 in granularity INCREMENT */ max_local_particles = - PART_INCREMENT * ((part + PART_INCREMENT) / PART_INCREMENT); + INCREMENT * ((part + INCREMENT) / INCREMENT); local_particles = Utils::realloc(local_particles, sizeof(Particle *) * max_local_particles); @@ -555,28 +554,6 @@ void realloc_local_particles(int part) { } } -int realloc_particlelist(ParticleList *l, int size) { - assert(size >= 0); - int old_max = l->max; - Particle *old_start = l->part; - - if (size < l->max) { - if (size == 0) - /* to be able to free an array again */ - l->max = 0; - else - /* shrink not as fast, just lose half, rounded up */ - l->max = - PART_INCREMENT * - (((l->max + size + 1) / 2 + PART_INCREMENT - 1) / PART_INCREMENT); - } else - /* round up */ - l->max = PART_INCREMENT * ((size + PART_INCREMENT - 1) / PART_INCREMENT); - if (l->max != old_max) - l->part = Utils::realloc(l->part, sizeof(Particle) * l->max); - return l->part != old_start; -} - void update_local_particles(ParticleList *pl) { Particle *p = pl->part; int n = pl->n, i; diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index 0f11856bdb9..6122b68f726 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -112,14 +112,6 @@ void free_particle(Particle *part); /* Functions acting on Particle Lists */ /************************************************/ -/** Allocate storage for local particles and ghosts. This version - does \em not care for the bond information to be freed if necessary. - \param plist the list on which to operate - \param size the size to provide at least. It is rounded - up to multiples of \ref PART_INCREMENT. - \return true iff particle addresses have changed */ -int realloc_particlelist(ParticleList *plist, int size); - /** Append a particle at the end of a particle List. reallocates particles if necessary! This procedure does not care for \ref local_particles. diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index 4ea17e8be7c..736d922d3de 100644 --- a/src/core/rattle.cpp +++ b/src/core/rattle.cpp @@ -34,6 +34,7 @@ int n_rigidbonds = 0; #include "grid.hpp" #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "particle_data.hpp" #include From e2347b75dbb1de299b44fde4ccd65139973fb16f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 16:20:32 +0200 Subject: [PATCH 03/11] Remove tutorial 10 reaction ensemble Has been superseded by tutorial 12 constant pH. --- doc/sphinx/introduction.rst | 1 - .../10-reaction_ensemble.tex | 382 ------------------ .../10-reaction_ensemble/CMakeLists.txt | 18 - .../figures/alpha_vs_C.pdf | Bin 12577 -> 0 bytes .../figures/alpha_vs_C.png | Bin 7307 -> 0 bytes .../10-reaction_ensemble/figures/qdistrib.pdf | Bin 16799 -> 0 bytes .../10-reaction_ensemble/figures/qdistrib.png | Bin 10338 -> 0 bytes .../figures/titration.pdf | Bin 14164 -> 0 bytes .../figures/titration.png | Bin 8737 -> 0 bytes .../10-reaction_ensemble/references.bib | 10 - .../10-reaction_ensemble/scripts/RE_vs_cpH.py | 96 ----- .../scripts/RE_vs_cpH_poly.py | 229 ----------- doc/tutorials/CMakeLists.txt | 2 - doc/tutorials/Readme.rst | 1 - 14 files changed, 739 deletions(-) delete mode 100644 doc/tutorials/10-reaction_ensemble/10-reaction_ensemble.tex delete mode 100644 doc/tutorials/10-reaction_ensemble/CMakeLists.txt delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.pdf delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.png delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/qdistrib.pdf delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/qdistrib.png delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/titration.pdf delete mode 100644 doc/tutorials/10-reaction_ensemble/figures/titration.png delete mode 100644 doc/tutorials/10-reaction_ensemble/references.bib delete mode 100644 doc/tutorials/10-reaction_ensemble/scripts/RE_vs_cpH.py delete mode 100644 doc/tutorials/10-reaction_ensemble/scripts/RE_vs_cpH_poly.py diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index b7818b0773b..cbb296a69f7 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -297,7 +297,6 @@ Currently, the following tutorials are available: * :file:`06-active_matter`: Modelling of self-propelling particles. * :file:`07-electrokinetics`: Modelling electrokinetics together with hydrodynamic interactions. * :file:`08-visualization`: Using the online visualizers of |es|. -* :file:`10-reaction_ensemble`: Modelling chemical reactions by means of the reaction ensemble. * :file:`11-ferrofluid`: Modelling a colloidal suspension of magnetic particles. * :file:`12-constant_pH`: Modelling an acid dissociation curve using the constant pH method diff --git a/doc/tutorials/10-reaction_ensemble/10-reaction_ensemble.tex b/doc/tutorials/10-reaction_ensemble/10-reaction_ensemble.tex deleted file mode 100644 index a525ec1a7a6..00000000000 --- a/doc/tutorials/10-reaction_ensemble/10-reaction_ensemble.tex +++ /dev/null @@ -1,382 +0,0 @@ -% Copyright (C) 2010-2019 The ESPResSo project -% Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 -% Max-Planck-Institute for Polymer Research, Theory Group -% -% 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 . -% - - -\documentclass[ -a4paper, % paper size -11pt, % font size -twoside, % two sided -footsepline, % add a line to separate the footer -headsepline, % add a line to separate the header -headexclude, % header does not belong to the text -footexclude, % footer does not belong to the text -pagesize, % set the pagesize in a DVI document -]{scrartcl} - -\usepackage[version=4]{mhchem} -\usepackage{hyperref} -\usepackage{graphicx} -\usepackage{fancyvrb} - -\include{common} - -\begin{document} -\title{Tutorial 10: Reaction Ensemble and Constant pH ensemble} -\author{Andrea Tagliabue, Jonas Landsgesell} - -\maketitle -\tableofcontents - -\section{Introduction} - -This tutorial introduces the basic features for simulating titratable systems via two different methods, the Constant pH method and the Reaction Ensemble method, pointing out differences and similarities between these two approaches and explaining the reasons to prefer the latter in some particular cases. -The first part (example one) of this tutorial will consider a homogeneous aqueous solution of a titratable acidic species HA that can dissociate as follows -\begin{equation}\label{disso} -\ce{HA <--> A- + H+}\text{,} -\end{equation} -while the second part (example two) will consider a polyelectrolyte chain formed by a series of titratable monomers bonded together. - -\noindent Both constant pH and Reaction Ensemble methods are implemented via a Monte Carlo algorithm with the classical Metropolis-Hastings acceptation rule. Here we compare both methods based on the paper published by Landsgesell et al.\cite{landsgesell2017simulation}. - -\section{Theoretical Background} - -In this tutorial we deal with chemical reactions of the following kind: -\begin{equation}\label{disso} -\ce{HA <--> A- + H+}\text{.} -\end{equation} -If $N_0 = N_{\text{HA}} + N_{\text{A}^-}$ is the number of titratable groups in solution, a degree of dissociation $\alpha$ can be defined: -\begin{equation}\label{alpha} -\alpha = \frac{N_{\text{A}^-}}{N_0}. -\end{equation} -Another observable is the degree of association which is given by: -\begin{equation*} -\overline{n}=\frac{N_\text{HA}}{N_0}=1-\alpha. -\end{equation*} -The equilibrium concentration of each species is defined by the equilibrium reaction constant -\begin{equation} -K =\exp(\beta \sum_i \nu_i(\mu_i-\mu_i^0) )= \frac{a(\text{A}^-)a(\text{H}^+)}{a(\text{HA})}; -\end{equation} -for each species $i$ $a(i)=e^{\beta (\mu_i-\mu_i^0)} = (c_i/c_{i}^{0})\gamma_i$ denotes the relative activity, $\mu_i$ is the chemical potential, $\mu_i^0$ is the chemical potential under some standard reference conditions, $\nu_i$ is the stoichiometric coefficient, $\gamma_i$ is the activity coefficient, $c_{i}$ and $c_{i}^{0}$ are the concentration of the species $i$ and its standard concentration, respectively. In the case of non-interacting particles (ideal gas) or very dilute solutions ($\gamma_i \approx 1$) activities are equal to $c_i/c_{i}^{0}$. In chemical equilibrium, which is defined by $\Delta G=\sum_i \nu_i \mu_i=0$, we obtain that $K=\exp(\beta \Delta (G-G^0))=\exp(-\beta \sum_i \nu_i \mu_i^0)$. -Therefore we have (in the case of no interactions): -\begin{equation} -\label{(eq5)} -\exp( \beta [ (\mu_{\text{H} ^+}- \mu_{\text{H}^+}^0) + (\mu_{\text{A}^-} - \mu_{\text{A}^-}^0) - (\mu_\text{HA} - \mu_\text{HA}^0) ])= -\prod_i \left( \frac{c_i}{c_i^0} \gamma_i \right)^{\nu_i} = \text{const} \text{.} -\end{equation} -Equivalently, for reaction (\ref{disso}) and in case of a dilute solutions ($\gamma_i \approx 1$), equation (\ref{(eq5)}) can be approximated as follows: -\begin{equation} -K_{\text{c}} = \frac{c_{\text{H}}c_{\text{A}}}{c_{\text{HA}}}=\text{const}, -\end{equation} -where $K_{\text{c}}$ carries the dimension 1/volume. - - -\subsection{Constant pH method} - -In the constant pH method, the acceptance probability for a reaction is - -\begin{equation} -P_{\text{acc}} = \text{min}\left\lbrace 1, e^{\beta \Delta E_\text{pot} \pm \ln_{10} (\text{pH - pK}_\text{a})}\right\rbrace \text{,} -\end{equation} -and the acceptance probability of a reaction is $P_\text{acc}=\frac{N_\text{HA}}{N0}$ for a dissociation and $P_\text{acc}=\frac{N_\text{A}}{N0}$ for an association reaction \cite{landsgesell2017simulation}. Here $\Delta E_\text{pot}$ is the potential energy change due to the reaction, while $\text{pH - p}K_\text{a}$ is an input parameter composed by two terms, pH and -p$K_\text{a}$, that represent, respectively, the concentration of \ce{H+} ions in the solution and the logarithm in base 10 of the thermodynamic dissociation constant for HA when the system is approximated as a dilute solution ($\gamma_i \approx 1$): -\begin{equation} -\text{pH} = -\log_{10} \frac{c_{\text{H}^+}}{c_{\text{H}^+}^0}\text{;} -\end{equation} -\begin{equation} -K_\text{a} = \frac{(c_{\text{H}^+}/c_{\text{H}^+}^0)(c_{\text{A}^+-}/c_{\text{A}^+-}^0)}{(c_{\text{HA}}/c_{\text{HA}}^0)} \text{.} -\end{equation} -The chemical prefactor $\pm 1$ defines the direction of the reaction (+1 dissociation, -1 association). -When a dissociation move is attempted, the titratable molecule HA is charged and a counterion \ce{H+} is randomly placed into the simulation box, while when an association move is attempted, a \ce{A-} is neutralized and a random counterion \ce{H+} is removed from the cell. - - -\subsection{Reaction Ensemble method} -A chemical reaction involving ${n}$ species of type ${s_i}$ and with stoichiometric coefficients $\nu_i$ can be written as -\begin{equation} -\sum_{i = 1}^n \nu_{i} s_i = 0\text{.} -\end{equation} -For such a general reaction, the acceptance probability in the Reaction Ensemble method is defined as -\begin{equation}\label{RE} -P_{\text{acc}} = \text{min} \left\lbrace 1, K_{c}^\xi V^{\bar{\nu}\xi} \prod_{i = 1}^{n} \left[ \frac{N_{i}^0!}{(N_{i}^{0} + \xi \nu_i)} \right] e^{-\beta \Delta E_\text{pot}} \right\rbrace\text{.} -\end{equation} -Here $K_c$ is the ideal reacting gas quantity introduced above, $V$ the volume, $\bar{\nu}$ the total change in the number of particles, $N_0^i$ the number of particles prior to the reaction and $\xi$ is the extent of the reaction, which could assume value +1 (forward reaction) or -1 (backward reaction). Each reaction is proposed with uniform probability. - -\noindent For reaction (\ref{disso}) eq. (\ref{RE}) can be simplified: -\begin{equation} -P_{\text{acc}} = \text{min} \left\lbrace 1, ( K_c \frac{N_{\text{HA}}}{N_\text{A} - N_{\text{H}+}/V)} e^{-\beta \Delta E_\text{pot}} \right\rbrace\text{.} -\end{equation} -Notice that in this case you can also define $K_\text{a} = K_c/c_0$ . - -The main difference in the two methods consists in the fact that in the Reaction Ensemble the system pH is determined via the actual proton concentration in the simulation box, while in Constant pH method it represents an input parameter which remains constant during all the simulation. - - -\section{Example one: homogeneous aqueous solution of acidic species} - -Compile Espresso with the following features in your -\texttt{myconfig.hpp} to be set throughout the whole tutorial: - -\begin{verbatim} -# define LENNARD_JONES -# define ELECTROSTATICS -\end{verbatim} - -\subsection{System Setup} - -We start by defining some important input parameters and setting the particles randomly inside the box: -\begin{verbatim} -import espressomd -from espressomd import reaction_ensemble - -box_l = 15 #side length of our cubic simul box -system = espressomd.System(box_l=[box_l]*3) - -system.set_random_state_PRNG() -np.random.seed(seed=system.seed) - - -# Particle setup -############################################################# -# type 0 = HA -# type 1 = A- -# type 2 = H+ - -N0 = 50 #number of titratable particles -K_diss = 8.3**(-4) #dissociation costant - -# titratable particles (HA) -for i in range(N0): - system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=1) -# counterions (H+) -for i in range(N0, 2 * N0): - system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=2) -\end{verbatim} - -For each negatively charged \ce{A-} particle (type = 1) we put in the simulation box the respective counterion \ce{H+} (type = 2) in order to maintain the electroneutrality of the system. -Notice that the implementation in Espresso requires that the dimension of the equilibrium constant is consistent with its internal unit of volume; rules for a correct conversion of $K_a$ (experimental) to $K_c$ (in internal units) are explained in the user guide: \url{http://espressomd.org/html/doc/advanced_methods.html}. The next step is to define the reaction system and seed the pseudo random number generator which is used for the Monte Carlo steps: - -\begin{verbatim} -mode="reaction_ensemble" -#mode="constant_pH_ensemble" - -RE=None -mode=="reaction_ensemble" -if(mode=="reaction_ensemble"): - RE = reaction_ensemble.ReactionEnsemble(temperature=1, exclusion_radius=1, seed=42) -elif(mode=="constant_pH_ensemble"): - RE = reaction_ensemble.ConstantpHEnsemble(temperature=1, exclusion_radius=1, seed=42) - RE.constant_pH=7 - -RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], -product_types=[1, 2], product_coefficients=[1, 1], default_charges= -{0: 0, 1: -1, 2: +1}) - -system.setup_type_map([0, 1, 2]) -\end{verbatim} -You can switch from one method to the other simply by changing the \texttt{mode} parameter from \texttt{"reaction\_ensemble"} to \texttt{"constant\_pH\_ensemble"}. Now the system is ready to be simulated. - -\begin{verbatim} -Nreac = 10**4 #number of association/dissociation attempt to be performed -for i in range(Nreac): - RE.reaction() -\end{verbatim} -Notice that, at this level, we're not taking into account electrostatic interactions between charged species. Therefore, if you compare the input equilibrium constant $K_c$ and the effective one which can be calculated at the end of a simulation run $K_c^{\text{(eff)}} = \frac{\langle c_{\text{A}^-} \rangle \langle c_{\text{H}^+} \rangle }{\langle c_\text{HA} \rangle}$ you will obtain $K_c \simeq K_c^\text{(eff)}$. This is not true when e.g. electrostatic interactions are enabled due to the fact the latter introduces an excess chemical potential. - - -\subsection{Dissociation degree versus concentration} - -Performing several simulations with constant $N_0$ but varying the dimension of the simulation box it is possible to obtain the dissociation degree $\alpha$ (eq. \ref{alpha}) as a function of the concentration of titratable units. Results are shown in figure \ref{alpha_vs_C}. -\begin{figure}[tb] - \centering - \includegraphics[scale=0.6]{figures/alpha_vs_C.pdf} - \caption{Dissociation degree $\alpha$ as function of the concentration of titratable groups. Yellow squares: Constant pH method; light blue triangles: Reaction Ensemble method; blue solid line: ideal behavior.} - \label{alpha_vs_C} -\end{figure} -As you can see, only the curve obtained with the Reaction Ensemble method fits the ideal behavior described by the dilution law -\begin{equation} -K_{a} = \frac{\alpha^2}{1 - \alpha}c_\text{titr}\text{,} -\end{equation} -where $c_\text{titr}$ is the concentration of titratable units in the solution, and this is due to the fact that the acceptance rule in Constant pH method does not depend on the volume of the system. As $\alpha$ differs, so will the number of counterions in the cell. This could have a strong impact on screening effects when electrostatic interactions are taken into account. Moreover, in the Constant pH method the real chemical nature of counterions is unknown. In fact, when a HA molecule dissociates at high pH value, the generated \ce{H+} would react instantaneously with an \ce{OH-} ions, so the positively charged ions that remains in solution must represent a different species (e.g. a \ce{Na+} cation). - - - - -\section{Example two: linear weak polyelectrolytes} -Weak polyelectrolytes are a family of responsive materials whose application spans the range from flocculation-induced water purification to tissue targeted drug delivery of expensive or cytotoxic medicines. Their properties, however, depend on their specific chemical environment, with details such as pH, salt concentration and valency, and their topology markedly modifying the chemical behavior of their ionizing groups. -Here we're going to study the titration curves of a linear weak polyelectrolyte composed by 50 titratable units with both Reaction Ensemble and Constant pH methods. In this case, we'll enable electrostatic interactions between charged species in order to correctly describe their properties. - -\subsection{System setup} - -First of all, we start to define some important variables and create a dictionary containing the type and the charge for all the species in solution (this could be very useful when you have to handle with a lot of different species in your simulations): -\begin{verbatim} -# system parameters -box_l = 56.134 -l_bjerrum = 2.0 -temperature = 1.0 - -system = espressomd.System(box_l=[box_l]*3) - -# particle setup -N_P = 1 #number of chains -MPC = 50 #monomers per chain -N0 = N_P*MPC #total number of monomers -nNaOH = 0 #number of initial Na+OH- -nHCl = 0 #number of initial H+Cl- (additional H+'s) - -type_HA = 0 # type 0 = HA -type_A = 1 # type 1 = A- -type_H = 2 # type 2 = H+ -type_OH = 3 # type 3 = OH- -type_Na = 4 # type 4 = Na+ -type_Cl = 5 # type 5 = Cl- - -charges={} -charges[type_HA] = 0 -charges[type_A] = -1 -charges[type_H] = 1 -charges[type_OH] = -1 -charges[type_Na] = 1 -charges[type_Cl] = -1 -\end{verbatim} -Here, we've defined the type and the charge of polymer titratable units (\ce{HA <--> A-}), counterions \ce{H+}, and additional ionic species (\ce{Na+OH-} and \ce{H+Cl-}) that can be inserted inside the simulation box. Now we have to set up all these particles and to define all the interactions. -\begin{verbatim} -# bonding interactions -bond_l = 1.2 #bond length -kbond = 100 #force constant for harmonic bond -harmonic_bond = interactions.HarmonicBond(k=kbond, r_0=bond_l) -system.bonded_inter.add(harmonic_bond) - -# non-bonding interactions parameters -lj_eps = 1.0 -lj_sig = 1.0 -lj_cut = 1.12246 -lj_shift = 0.0 - -# setting up the polymer -positions = polymer.positions( - n_polymers=N_P, beads_per_chain=MPC, bond_length=bond_l, seed=13) -for polymer in positions: - for i, pos in enumerate(polymer): - id = len(system.part) - system.part.add(id=id, pos=pos, type=type_A, q=charges[type_A]) - if i > 0: - system.part[id].add_bond((harmonic_bond, id - 1)) - -# setting up counterions -for i in range(N0): - system.part.add(pos=np.random.random(3) * system.box_l, type=type_H, - q=charges[type_H]) - -# setting up other background ions -# - Na+ and OH- -for i in range(nNaOH): - system.part.add(pos=np.random.random(3) * system.box_l, type=type_OH, - q=charges[type_OH]) -for i in range(nNaOH): - system.part.add(pos=np.random.random(3) * system.box_l, type=type_Na, - q=charges[type_Na]) -# - (additional) H+ and Cl- -for i in range(nHCl): - system.part.add(pos=np.random.random(3) * system.box_l, type=type_H, - q=charges[type_H]) -for i in range(nHCl): - system.part.add(pos=np.random.random(3) * system.box_l, type=type_Cl, - q=charges[type_Cl]) -\end{verbatim} -As already explained, we need to enable electrostatic interactions to correctly simulate the physics of this type of system. -\begin{verbatim} -# setting up electrostatics -from espressomd import electrostatics -p3m = electrostatics.P3M(prefactor = l_bjerrum*temperature, accuracy=1e-3) -system.actors.add(p3m) -\end{verbatim} -To titrate a weak polyelectrolytic chain, which is composed by $N_0$ titratable units with a certain $\text{p}K_a$, in the Reaction Ensemble, we need to modulate the pH of the environment. To do this, we have to add additional \ce{H+} ions, increasing the pH, or \ce{OH-} ions, to decrease it. We have also to introduce the respective counterions (\ce{Cl-} or \ce{Na+}) to preserve the electroneutrality of the system. Finally, we have to take into account also the autoprotolysis of water (\ce{H2O <--> H+ + OH-}) beside the main dissociation reaction. - -\begin{verbatim} -K_diss = 1.0*0.002694 #eq constant HA <--> A- + H+ -# 0.00269 is the conversion constant from mol/L to internal units -# when sigma is 3.55 angstrom - -K_w = 10.0**(-14)*0.02694**2 #eq constant for autoprotolysis -# notice that here we converted the value from (mol/L)^2 to -# internal units - -#HA <--> A- + H+ -RE.add_reaction(gamma=K_diss, reactant_types=[type_HA], - reactant_coefficients=[1], product_types=[type_A, type_H], - product_coefficients=[1,1], default_charges={type_HA: charges[type_HA], - type_A: charges[type_A], type_H: charges[type_H]}) - -#H2O autoprotolysis -RE.add_reaction(gamma=(1/K_w), reactant_types=[type_H, type_OH], -reactant_coefficients=[1,1], product_types=[], product_coefficients=[], -default_charges={type_H: charges[type_H], type_OH: charges[type_OH]}) -\end{verbatim} -Please notice that the order in which the species are written in reactants and products lists is very important because, when a reaction is performed, the first species in the reactants list is replaced by the first species in the product lists, the second reactant species is replaced by the second product species, and so on. Moreover, if the reactant list has more species than the products list, reactant molecules in excess are deleted from the system, while if the products list has more species than the reactants list, product molecules in excess are created and randomly placed inside the simulation box. For example, reversing the order of products in our reaction (i.e. from \verb|product_types=[type_H, type_A]| to \verb|product_types=[type_A, type_H]|), a neutral monomer would be positively charged and a negative monovalent counterion would be randomly placed inside the cell. - -Due to the fact that, at a certain pH, the ionization degree of a weak polyelectrolytes depends on its spatial conformation, in order to obtain correctly averaged values we have also to couple the reaction algorithm with MD simulations. -\begin{verbatim} -# Integration parameters -system.time_step = 0.01 -system.cell_system.skin = 10. #only for tutorial purposes -system.cell_system.max_num_cells = 2744 -system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42) -\end{verbatim} -Finally, we are ready to run the simulation! -\begin{verbatim} -for i in range(12000): - RE.reaction(50) - system.integrator.run(500) - print(i,") HA", system.number_of_particles(type=type_HA), "A-", - system.number_of_particles(type=type_A), "H+", - system.number_of_particles(type=type_H), 'OH-', - system.number_of_particles(type=type_OH), 'Cl-', - system.number_of_particles(type=type_Cl), 'NA+', - system.number_of_particles(type=type_Na)) -\end{verbatim} - - - -\subsection{Titration curves} - -\begin{figure}[h] - \centering - \includegraphics[scale=0.6]{figures/titration.pdf} - \caption{Titration curves for a linear weak polyelectrolyte chain composed by $N_0=50$ titratable monomers. Yellow squares: Constant pH method; light blue triangles: Reaction Ensemble method; blue solid line: ideal behavior.} - \label{titration} -\end{figure} -For a solution of weak acidic molecules with a certain $\text{p}K_\text{a}=-\log_{10}(K_c/c^0)$, the trend of $\alpha$ as function of pH is described by the Henderson-Hasselbalch equation: -\begin{equation} -\text{pH} = -\log_{10}(c_{\text{H}^+})/c^0). -\end{equation} -However, when titratable units are bonded together in a polyelectrolyte chain their effective $\text{p}K_\text{a}^\text{(eff)}$ differs from the ideal one ($\text{p}K_\text{a}$); this can be explained by the fact that the charge carried by a dissociated monomer tends to partially inhibit the dissociation of its neighbors, and this results into a lower total degree of dissociation with respect to the non-bonded acidic units case. Anyway, this effect can be partially compensated by the presence of counterions, which are able to screen repulsive interactions between dissociated monomers. As you can observe in figure \ref{titration}, Constant pH and Reaction Ensemble method results are very similar at high pH values, but they show very pronounced differences at low pH values. More in details, $\text{p}K_\text{a}^\text{(eff)}$ tends to the ideal one when the concentration of \ce{H+} is high. This depends on the fact that with Reaction Ensemble method we have to inject a strong acid (\ce{H+}\ce{Cl-}) in order to titrate the polymer, and this results in a more salty solution with a strong screening power. This behavior would be reversed in case of a weak poly-base, with superimposable curves at low pH values and the Reaction Ensemble one approaching the ideal one when the amount of \ce{OH-} ions becomes relevant. - - -\section{Charge distribution along the chain} - -\begin{figure}[h] - \centering - \includegraphics[scale=0.6]{figures/qdistrib.pdf} - \caption{Mean charges of monomers as function of their position along the chain (calculated via Reaction Ensemble method) at $\alpha \simeq 0.4$ Notice that the curve is noisy because of the short duration of the simulation. Make sure to perform longer simulations to achieve a better convergence.} - \label{qdistrib} -\end{figure} -Figure \ref{qdistrib} shows the mean charge assumed by each monomer during the simulation as a function of its position along the chain. As you can observe, monomers lying at the extremities tend to be more charged than those lying in the innermost regions. This could be easily explained thinking that the ends of the chain can better arrange in space in order to minimize repulsive interactions between charges. This results do not depend on the method, i.e. the shape of the curve at a certain dissociation degree would be te same also with the Constant pH; however, as previously discussed, at a certain value $\text{pH - p}K_\text{a}$, the total dissociation degree is method-dependent, so the Constant pH curve would be shifted to lower mean charge values for a weak poly-acid at high pH values. - -\bibliographystyle{unsrt} -\bibliography{references} - -\end{document} diff --git a/doc/tutorials/10-reaction_ensemble/CMakeLists.txt b/doc/tutorials/10-reaction_ensemble/CMakeLists.txt deleted file mode 100644 index b333ad94adf..00000000000 --- a/doc/tutorials/10-reaction_ensemble/CMakeLists.txt +++ /dev/null @@ -1,18 +0,0 @@ -file(COPY "scripts" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) - -get_filename_component(BASENAME ${CMAKE_CURRENT_SOURCE_DIR} NAME) - -# we assume the tex filename is the same as the directory name -add_custom_command( - OUTPUT ${BASENAME}.pdf - COMMAND sh ../../latexit.sh - ${CMAKE_CURRENT_SOURCE_DIR}:${CMAKE_CURRENT_SOURCE_DIR}/../common - ${BASENAME} - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${BASENAME}.tex - ${CMAKE_CURRENT_SOURCE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR}/../common -) - - -add_custom_target(tutorials_10 DEPENDS ${BASENAME}.pdf) - diff --git a/doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.pdf b/doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.pdf deleted file mode 100644 index bedb55fdc09843816fcce1bb751d60e734486f09..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12577 zcmdsd2UJsAx9$c+JfMOKhy}qeO^SpfRk5NJX`x9d(yI`V5{kzICu@QK}%H$so^4c~kjXZdEn%^^Hz zCU$?@^(b-IV*#;4W1izjcOQGcQ|ue^k!abl7dz{I$WWMjp3~Ot+7Ghe7c3Ag1^hpd zLtYMk4L4-|l$I&qSzq*deKe3wSXmkF0DqEJQlKF$sc}kJT3FJ?3Os?t;Klk=(pbmZ z$<_wV9OH`+M!>ql@4&wwJmmi5AuMU&=V2=>scU6#%jBnP*hW0>F(|E7tzMVZ)jNM09`Zw1JzFn#Z%-Rr zFJahjdVkRoRvWhGU%amHr&z!8I?QHWydMT+>ZI>&jb>WtoV%whv*+}`j3Fl_B`k}O z`l0^6WRd?V%bzehd`Eg~h%nCvYCv6(dRdL&u)mu6S=YM-dbi|MQJWvz9U;9)-=^{K ziCy#QH-}o!@~(tM#O{uI_F?19U7^BJH9=?2?FIXWch=Ku2w~)>tz@XFkXPrQh=Miz zfhYy}pWcS0S+@jJ{{Jna6onOG_5K}EaM$0gEe*rvKM%71#@rWOe_6}%aV0mU!!`;Q*Q{ThN24MK#@F95-VA0dV9S-Y1ktBKA5$&aG>#5=H>f!DN z<6Pn2f_=Jg&>LG(A@>RWK+ilPvW&g#&H}aKA@6lg4-rg7ODD#~WrxebLyOId?e8oV zH{APbt8!_(wE*FoxV*5 zw(6UN;Uh@DdUD>!h2i^;{^ulY&C>2>G$JJ?v&A(~b>qK!PT#|sw*QLwI@9j`sN&qL zSDOR*vi|ZMlSR^FdW(=#%U&JJ|UBI6&)=nwnZd(_zabf&6z$ z_i&Kq>U32_pBqch+mN&UH!j?^x`Ne$2}2Lh+_{6-#!xLSj; zmNZ^XsQ<)t|NC0A!jSLF!^>-hg$L{8GiCJ9(2`YJuZ6y|<=a-X-Q&xw34|Z9H}g&+34h`9L1t#Qn^)|DB*F zxxRky*%{Omn{}d8y>7D9nD`o;rI>xugN^aEBo1!Gw>0dNl( zmj%xSB&3j@H-^L^O!jny8r%HO~C+f!&cMtL1i1MjY~ahvVj7hwq`lsV@~;0 zaBCF1_oYvnxAcIBVzsl(7uo}+I5h=#X%a&_Xki0KWY}*lHs%VA*77b+<c)hnYlP0Aohe^fl+qYbO5T|JBABYhj>WLT9 zQa}VOdRZKvsx1_XX0rDXCZ7A)+OS$Abfh8QDiX6a2N;lv&qh8K>k!)(tD10VNz3~; zgi$Pe-xaRjSFIOXu2r6TZ^9^GFhwCdi1=SUp3o~&YW9|L2o*@1Wymhs*J^(8n&+#z znLCSU#9{h_m!{&>VlYCb@iB5fu38N$mtt`iqB4{3kr2oHIYwLUkV|Xu(!~7^q3PP& z#dze$HS8(8s0N-@;a(^iN$QYK4qlflw|5zp%k;UwUUt=CSnVaW8? zRn;==5>v}_`I9kg?^|XLwyz%ayB}eH)Kv0@v6iGGv1=*n9@Bf#8E+{_uF7cwz3uZG z^+6+vj?qBt@Plj*ftKhmu`%DPQpZNaC9W4qr#hC)JwI{@FgZE*zG0P&32JvInc_`J zhz9KEYcTvdI0NsN#^dG}^S}a7Zf@NiCwhvBv@gfz&9jxSh@KUX>NsLx1I9vUj<7wH z8Wn5ff>_u$Lql3+g;Znv`8*-_0W+{rs2jJi)}W`0yH2_k_SKJcd+rPH#;0M}KV{(R z&hzlAjH{{@6I&{F=2h55E@w@LM=ldnb?v=z@ki7j+` zuX0r~1X)V7ERIEZpD?Z^q$(b*u3M>0++ zpsKmBBfRc$4Kele36(xFW?wL7uf?vOlyR}(RB2*tx!O=QqC;xB_(IQkDP1i_g-5+(ZiV@Au7h{Qa8{J(Q;P@7H>O7MdyAVw#*cAf?+MTL7WMBFO7@>k zWO*^tlL~MeHQdcZUsZloeFt=GH3J*am^^0(A18``10QLLV3^LM4^KQbJz6 zt|#;7r|7S*Gm6am?bq&`5MHn9bJ=4gQ4;30bF)euSv^lSL7^b?1+IG0u-}6gy6!aq z=n&fpySbBgrfO2nk0(@VG;4D=&{I>sS7VAjw$u`k_{6;5`NARH+h7M|^&4|OPyGLb zCPwQso*%MAnleu7)zQJ)2%>Yt^CihMc-D~oB08)WPkcLphs+7RX{GE0>!rFe7&|Nm zZMrmxVoL8w;{t{g7S9Xb4<`SS@qckLiy32bi9(~LaHirjr|m0gYir5cSw&`)#1fQB zBy7sF&NO{J%>vP}jPd2mjaR{r$fq)6dJ*9&d8i;9mRihoQtgta;`zSWcu(Y3q)E3n zt_pnsRxhHRg4?VORG2Gq$YTWh3ANWO`KzR8*K@Fh$e(ehAxE);Hx7ND6yk2+VJL$D zhok|wfk}d9Pa0S$`>PzFbSqTMXZ(w+kAvBvhx~a-Y()&BWYqyV>Rx+IMnuAUsH{^vfCuWjjHe zPCrAvlSPdKH0sJy)~$W+hqb@0TYL8pYd`zps>O9zMg4G9%MVymLNALM0fGM4WSD8^ zzm~z!WidD?sCjH*Wo5;8OBQL|r?Rag)s^Y8ymJiy>db3#_)K+#MeT^^#XH!>B^Pbj zR89zrh&jj`A6{A|7DxZ=A}y%)sLH}wz$u3zTK&KR=c$qu3$j|S&nRldBhDiGXS-sB zJ5@Vtud?1N%^sT=BMv)0>Qanp2OBhXE)w-F*V@~PMSZX78WWsNl&4H<8YwOGgjcuk z{Fv1oi|5ZIWpOo;YFpjh&Ew#dg(z07-8vaLTFJYRoz-(f(R^evu4G_S5EEstfh|l- zI#R7XD%KM5^Am5j!7thj7;Ruu=8#U7Z^m@E{TilV6H|P9dYkHGWPien)t2GagfUdL z{G1qNp@M*$ zYvT$M7=}TYem}V|sf%`t$Zv5$QhQbiWAn8d7QM~s`@+|0s@f+-u(9kQB-U4rYnl&T zqbe(`4pKJXf`=GP5pob3A`Lmp*+Eu39{!9kdlYq+IwE$}82VHlF?pgl28zaJtC<(0@mcEb+eV= znQL-)c^<>mia064%mYrvOMdf)l*FAbI;e%aD{h&V>)}vBA9PvXTOn$$#nMF9e9IaV zU7U~qme|2PW^zhgniJ~7b_RPhj*Ts9uF#M`6f|yz@|3jyjj}=~jq;S)L7_Yqu^p*u z(J$N0&x~=r&X==XZ<}*U4I-^D+;&4#9O|dC?h7`IEwh;n>pfpnjj^xm$ zTqChGp_yd7ks@_#`>@zHE~vtwfH5&Rye#-~%)ffHFbmv|462ST`1u$*jwLid;F@0| zFksdyB5ThAcOStvvD7mp`pM(C= zxdEZz2)f^Nnug6}WIR3<3>|IzvJye7h8cC*s9kQSiypy>Cb~4;!s4Zfm%Z&5lfCOc`$8VX*37#*KqOL?2-+$0X#?ZYU zudoLjO)1n}1tgj`d_Yo-({Rh>`yfksrPA4<@qGgp9as9n*-}`1I$(YA4_wEKQ}HY} z@~l8mM^^XJ$q0CxnQJ;ttzDP_@O5i9aIx2KhC*Bz^)&R!5pV7p74y%u^mc<>u&wx} zr%ywwFlS`$Ug?=`|Gb1b4F`~@+{_T`IO+TArV+hK|B}N@*v$%Eo@b-aBZeEcHSf#C zn8QJAh`pYsYV21P##@k@*W;l8_qDxQTfVLp_OR@99#P8pDGTKP`7}7S<8z;xnee~~ zCoN(D_p-yeAHTMn{`hyppSUhIUGABdi_#Nu%MiOGM+ofI2co1pbgf5_M8JOBD3*bt zm3~`-pvwC=PCo>Sd=^rd_i<})AgN&RUOLi48NYu6D|9M4*533WEcCn~(^PHmD4R__ zEbR)rh~HiSrLz^!8%F32QJy{vnT(89g14A%MRRMf5jde{EIY7A`W^}8U=Q9oax;mx zWgDxykkbfOnRj(zAy(fA?4$L|DuiT}kF6)+eF?0F5ym3=8O|3_UcEWq((7xY+o)3U zXe{QiHA~Ht_sGzy1RNTydK@`1+@rep8#pQ-V-Qt_1;Ad5Prj}#Kg|B%o<_U9|1|VX zT_2ndZyOFear!EuGPVFFTR;9>BkSl`YMa0!YS)1cAzL2>8@>zR!*X#!Cv~qYAD0}Q zp}+3Z(OQVR%?_Ra{1KDYj(LRJ2n_D=(+q}Us6Mw+(ztA&Hu3?!U$Wan*K7n%3tO{{ z@QKYOyVIsH<>vLDywPdH@Y3Mn(AGf!SFVl$QzhMl{_}Z}etNXTT}zH>EBt1@V_!FK zW3BmaNJJ`5Px+|$@!jsdD%Tfugq^Ur9xRq&%p_`)t*Yv)nkJcWDM%XEsk>M;jwY#M zi&891I@MyXwl^Yg0URo%S1k09y0X68CEqq|o+=M7518m`7FivUDhbBRoY@35@VN$t zQRh{$Hy8^8p-fopFjJ1O_ovS*6EY$z`G zWZJCW8^E{F3f~&a4mBJ=@x3yvYeYgG3e*^Gx+7S|vK*Kwwk@9jKqTr|Op8jUrQ%0J zd-H<5gZIL26EB``tIas$^)1H(!R}%VqipY`Jj~KmRXvB$zD8Kk3jj*vx8e%dJ|%M^ z^=n4Sny~fG=3%_mdPKRmj)tq=ug~Qm>n#??-uz{R&59kJ1$ZF$7pX=65fRI}T zCgA;uEIgjO{U|Hce+0#+_tF}kI0PnwmoPFMVlajm#xb_Q+lhz4fFkCACNRL~9kL=H zN(VRP`!kqptsNjY6^olU0jN~hhey0HhO^BPz`Jk0t%wV|xQitO0n^Sl*7yVZBn)1t zFppOb5kLbsSOaxucw#g14v0jDKR|CPGr$sIF9b6k5_J#wu$;MtRRUvkhOq@)&P07S z=me02ip>2n@j@OZ1GqmY7RQ;b2Rooh*b-mym*ZAR&xof8aE=l;n z^>~4tJl{6$bE#*1avbfCux`ugbbOoR$pWbtVd8!#$%=r}dLzxP@Kf9YBTZ|`lo!@a zgH}8p!Scp>;tWb;6a=MFcZ?z3OGOxORKJ?Da46OU*AZ~Eo8*$WQG^RxK7^9!P~R#F zj_D|smeiTFM+ZToavR2a(vT|_C*Oo;9_0gxGI=5@$a~u;N(*OIwHf3IAJXZ~ouu5O z6Q1SDd#OoFMuI0{Hu;78)-|b~dWkQC5!21PorlUMQ6E249MfuG4H1;x&bIszoW>pn zmNd5&$?3iIFX`DGWiUlMmJbqD24XuKg`{1ln(S~C;7W1&ZAi1MVZi}@Xev8f;8Se7 zs1ewT^%e26RkT>vMF3tg8*S8ux^p|j#Zn)C-wO3Oo!E#@{}>Ti@*)?|+rw;jMO-y5 zMS+Q0vU}K~Uyek;%DjlTsd<91fN4J@`EFKtxh&?R_4dIl=*yPrhLBq)&%oq$-|V-p z*(R7_z4Ji3WiAiY@A*)`LxCg57-(POoi5-Z-oX4Vht3u75Dv#O`#sA9Jb2INcE)YH z)W{lQd;Tsh6Hqx1ANXv6aGa(Hqi(sJ7wPe8x6KS zw?;-LPx>xQAgbO9ZSR;tdrfNu zJgP0gZdU&_sRL{zOPfYyW|kKg7A8zsji&bSAF;h?4zhQVz`X~pj}kI_PK+%mj$nJH z<4Z^#r17GAZ}&wQ6{M|t6^LKw3bDDv&!*>*cq9X?Jmmwf@tJyWH7`zO>^O_4QhBj>7Nwt~u-l^{JLU>tfTlBJN!v9lrzT%GJ`vWR+WIF-stlYi}d05St|dxp26 zX;{Y_?Z$fw1(w4c_{=VqHyZ<`O4;wSC>jN-m_7>3jS5#eL^`yb6~brMqN>pk$cixy zO`BFO)#ca;5pTn0Hu#xwvHTU!Ebj0>jaA5R5*eb>spu1&b-sN8K_2pc5 zH2p>Avw-&Ows-FBz6aEsb~}nwo>^BG@D#n3fc7~(_x=_WF>OdS5!+D;i8|$+{dOl} z%Tro<9PoF5Q~v&Mw<5NDJ7(f6mAXmSI7x^us9SJf&>Fg~*~{8#R4BB1KFj)j{FpOh zPX&-#r-s;wLxTUTcaJ8XJ^ho&?$wc zT@<~=whiW0^KGafx1y%=T;|1TSWMO?Rqm^M>hpJHZ&p17lfRi)l>9dHxr_h|Fk?&m z^LMOH=W&@-HXgpongXJJUkM=da#9tqY;fAMdCEEc2{cmOIXGPP*Sn0Z|S1 zS=Ur-Ph{sJKy;>(b7mMJQ0H5~w0cYAkvP|kMfaS*Q4~vl($!lXC)0LnikjKP=BVEL z^Njbbbn>Wa9c%ur#=)ukamalVae27Ms>45T9;w~Fl|fGp#i40z(p#fMi@EFiUC1hd zdp2zyvy|z=y0Z=4s|)gU5Sm=-vP$(YR}3$rP}b_ZMEM~LL#0RwzhdFgoMQ@nTu|G4 zWXhES@G$|#3j2CoibH=KhON#bEY8FUWW_W)__V-uj=C|^h0q&jBK`kEMdi1aRi#B~i{`y(+C)mWijS5u$tpMe(#VFvn%OplshnfptV$))7{I`Rt zRH{#9RIpahc-)d&ze(M)L2_!N%%b5hpRe+$iGJBHPhvY$e_e&MUMK?&+FDEhH60>c z{L`M9{ai?>nZOEJ+A+Uv`;F62L6F_gm4Y9F+~2LAO%nKbvj3_YY3qc4bu{*Yc3M6u zJuaalAAFBiN+U`r*Uc}y0w3di;$!eJn*kyITeM(3iD8ZTY# zbX{99qp2+TnanV71_j8$g|8pzQAqLTG22}u+?^HfCB#UQHl9)?o`52yh{=83gC~V` z6O(CG;)RiF6prp(u9t2P(o9iB=Kx{niL2EvBN;Ed@+bjzwXfF5d4fRXyv7HrR~E*M z`{aAZ>Kw-i{b>mqp*OQyN*GQP%a-nRIVC`63KghLNJ5S6vB(>na2O~Z%A73~ zRh7KtkBTDPpBTE8w$na-7`tKC#eACgD90AiC|e*wBu~&GBVa+>Jp<$ zoJ{z`rYCR_6W4D#^to(pU~cH#zffmW*Ygp-4 z_uwbESlsew%Vpnlc>h&Ru7gwxE)zFgPH9rZCw&ai6pKc?JiysibMi54UbW#eLoW&@ zt)SWjXL|E@1l(P03qJoCQ=FA~t~%jyFrr6d>|O-80A9e^=}XSEQ4p?l?5+vgLf5qi z#N4#3ujgMI5fUrB*kl%P{62f=7w?)`e69RI=qBYpOAl?nXK!RNk=U0*u$L3Lv~=ob z3v6khcUuPROqgbW9F=Z3wLk`5W5j%(Htd{8X@4_TGiK{3igKmjY(KrR*~0ZZwlpq0 zexRPzCN-V$-E*#T)ce55+DZ#JU~0M%p1zk? zjDD~iCK|pBh&lyPjf5QhsNJm;q*!Fh9H->_xV5%^nZgzR_77Hu-?7RP#Hp z4AzdWSzQks(0P8YA4JedSB$}TtLvUND}D8?lr)A@kg4e*OzaKQC3)kL_qfdaIw2#G zCXAc$YVp{C?30C=4qC-=zY5xa_MvAdC$YGXaw2MbiSwKSB<=Qk~Q`gPT`*_F$|S^=5vQ4M=Ja-KX&dpX)I!o+l16=p{_e+_385NV0THdsvRo=k-GP_= zedoU0H(0$6#VqCOHoO*uPhueue|5?1tESUXKy(=xHASd&) zfEO;t{#jhRUZD(^YkM&ZU@tjd^8%&Wa7DM?pJmwIa0T~YYN4h7NbcTal}5hfp9PxouyHlR}e|9slN z>dmG9Rn-2o##+M)ZRO%_zfK$0=10l-zb5&5+S+kS2}_^gl=|~eSPp?Wfe^M6Uhk8Z zR#X5*^zfgs+mAjeDQN{zjr_MhDQQp-{w_rHyk zmji14txw?}nrTUC8JQDOoQDse)j7xcA5VVhM*si- diff --git a/doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.png b/doc/tutorials/10-reaction_ensemble/figures/alpha_vs_C.png deleted file mode 100644 index fc7530c482838fda0ada43edeb6e2dbeef480c96..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 7307 zcmb_>c|4TwzxOl-Lt%=n;U-&%knGI(k|kp)SwgaJSt3hW$C8AMqOwb(2$6l4wd`a! z_C3qk*D=pMzs`AH=X=igIp_TGyo}de_jP^t&-?m(uIru%U2QcQYF26p1VW>su5ud! zAp;-~s1ckDoXPv-t_*>kk=E7HS0xgOkW>#PC8f#9Npo}ay1F_}PR_KnwA$L*bRytK zM0VDHZ6G$Jmv-Wb_|CFaBC(!G3?kOo6Vu&@=?{qMNMaB&CNF@ ztCt7?|A^sTzafxL2qcJw#SH?9gg}TQk%run+z?_Q#2S2p5QiWI>+8hQhIL%<6p@I( zi_aTgz^-HQ;^Lr3^Yil^9Ub-c_20gIOGrov3=FijwMC=R5)u-utgK)JkW)RN-ya19 z1Bt&tZ;(MG9oHTN)dvwFK}cdffO42_YlGAyA&~X;dQcxYPx^wSr>BGLpb8ZQDk^Yv z=J9RqI}qf;ur2sN<)&`z0fEpqk^axb3(>McAgrqzDky!QM|DLcQ-wa`~sMx{OFhHb^7e7FsxPylaK1AT=x_Fih+2vWKI5zPGuvwZ9=ii!Z* zqBQHP<|m63(-Ed(6qvYc6gqRSP!pG|d>{1dP(oKFN6bw?b-lYV=p{!6OuBHtttL*C zIPOxMe45AlTKeJ8a8{D*&Zy(|B^R02HesSt;O*1PL^v#4o_gsR@8s}QRmuE!tX6Id zZO7veYXfmcQcEa>Zn2fei$pDYR-$Y2IZ$W}XO>B+47G}+pr`(t;$O?l&vfc>_N;kH zi!wL!SWGo@BV^Dk#n!OkI}6Ai3&4z1vdeb)A$^-T=6xB#G2@uh%*~w~GZm-AR4AJw zPh?&^+v#||9hAZRHam+~P_8;)Trnleg(_6JHr!m*^`3;hVV{(D(nO}HMRcUjtaMTs zzqeKL&tetNOKpRJ5z6LiC-MdGIVI>J_aLFt{IsHOOnxJclVNjq-vhgwRnWN_iqxgq z|K`WwR({dRtgW|O(UIL+PTiiJM-JL><<6O0^0a!K2bbpWu5Q#gZPfavnEm#sR3qYC zrl*Zpt|E9SB=tbiC1nNqb7_RG*E!_ z{N}eElR^WhWl@Y4Z+=d=^JisdvrMg z2q<6JHk2qUt|Ju#Z+u(`q?-=F5Ue>zXYDcZ!V#J@0f!B=akMQ|$LBwi^SW=C z1D~TDDbLeUvtwvc*n^IcIdtOXU7I;wzK#)e0HviA&O^>{fkBqGg?C~~N9!w3U|S&* zhq1&*a~3I>hc9LKJh1DYtgtPUDBo_6V`=ZfqgkDWx-_RR+n6issCCEnyug&=equ*V zA8zNfXAdX)@i?(0-|dOH$X3zOs;F�U9&5K(3xC4M(Gd9 z6Ysx~8uY-H8anoeOe6=?C&5sumQMV+K# ztN-WD@DD9xoB8erPb4i2+-?OI6>`x)i zCK%AevQ`>h?7#MQp?>ym*xzFK9?`>bAwW^zei1IxkFhu z;>5R4&N-V9>g8`Gfrg^4uFEh+B6X?$Q-A!MXJs&P1Y8jHw>tK}X|kY@sZyi>OtSn6 z!j|h3*?(=re^UcSc8k)oc)*s!zbK6llTJHAwhVSN`2;T|&`|aw%DzO4F_r9r05EW} z)^tsa={ws>tn}~xkQRzIdZCxVh{&39sP~6)CT+O1WJ8gOuEh z7fLN92UmK=QtPD9A{#GR$4fXvjo~?@#N@Mw!GZe>Dp! z#U8?@z7<7@(+i|-^r{@SiA2D)7`~$k|LgA^F2EhJ{%Y^Ngs!Y?e!?4AKv9xK8v@0@ zEj)*Y362Vwbk-4X1{Y}I>}-`8m<&hjNxCAMrrlqX1$0l;IoJ-8KN&82@hY^~vRjWP z?jC)Pc?O5gN1AKIjQ=&9I9*Ye>1k5v^?MiORH5uEdDMDEu#Qs&{a5azlv*b(H!$x)JVF>8O0)X$v_NmniSeDc3kb z)r*lYRZXy6>blSBL_1gadKU8b)`Xt)ssG-0J`Vug8P%B|o>11sn9fs$+6e0d)AkeE zGyw&Z%egTZ;5p~(5{9mij<&#FcKLgNR&3$QIKiXJ8Bk@J*RAoCD5TnGQ(}t13T=|7adT0-bJoP zbloM9c#0-kUS$t!uT3Z;)wHkza#te3mq8g2Fx^AmBVov?{w(G|84e4(Ra2C4=!C2j zSE}o_^XUvYXZ?;kmh3Rxm$3Wol-$8}=JRRojRG^e+}0egw9W3yjcVd4Cx!d|EwcE7 z*#V#XE4C9#1?#m`O>1W3l9J>zY;f|e9|at*lJ0Ly`5y9F*R96j6jn+;2gY<-tdBe- zlL<&fw#8#qoJ(-EN^J-+%8@T2@>x4us#nVlF+CJ7;{4e&J?Y*}{^${XQ*%FRljV6@ zLJtKsW8KB}Cg} zli?o-Vk(0-<(f^JUD9GAX;DLH`Dj-w=3KwEoT%oT>y%GARExA<8^bjj>u4>Xz29h! z$98}Si4&e)qoeTt*fcJ*Y?{WTT6OrK8H3a( z1P?I#4IqXxofj0I{uOek{}S1I1cRU~wv~4lEC4%BQm|lUn>wFeHWL1XcGu{NIKN1L zjS35HHm|OXsQUMX`K!al1I80)6G31@Z)_id?Lk%WGT0sj^!L0g{Ov%8ZrcQgfA-G$ zA(Ht+m49K6@9V5Q{<=|%a%_?7LP?TE72jiPjir8YA|jtgg?!R|QGn#7_8tL9LXO|N zg&%IhsRkIgC|hLYtESX)T7RS6de_o!qDGlb@;uI8N!_W1j&M}t_JOhRUt+;74g&F> zg7^)@n?l{x)ninP>}qk^wX3AZ4DfLxm9O1MK$j)F>X>N57`zonPfFr*Bm`v&c#V8~ z)R5#OLR26JT~gQcAQ^*3VZaP+sfbwSnM$ESDalJR>Wh;WKbP9~LckSh^M*#C>C>H8 zkvBTL5pe1|P0hJ)sz&AT?D4F9nHJ?>>nYP^AEseH6HcnCPXS7mz@J82G7cN#c!CCMyI}$H3B758^`4 zku(=7PR0U<*$CVq(NXK#ljfeelY!d{8pleC%lj#$FlNfDnL?c8Y1}1nZajY+Ndavj z{5giy_7=FRXnQn(|yVi7Hoql}$Gnr+^)d;?$DE;Zp=sZgdWi zC4^1cj*#`PM|<5sCxh2ww}Ms7@Q}q_r@pL;fXuO9x46r@=P^cOeA*rab?Sy1mxPd} zx@p<$>{$^#(SF|{6otsEuq3l@7jpxJ4JrcDXh0DkRxtwqYJzEFbm9sJo3*Jt!E$qx zuO}Sj*7zU)Rkj@+>QF2H4kQeYX-?!+$m7?8uj>YYmm;>i2zc0(q9jb2s~mk1b1nua z`6Xwz-;rckq%72!+CuDX+FcEh()ZT|kdko#R&*s)1<$oh;y|fl|(8c3X{Zjj8O4Q&=S&09wf+dYpQ;Ary_ z6I!gV$CG1p;raT3K(30>T{V|UMBCTu5$-VHTLIiUN&U2Uv#OQ{YsbwHoLtL?xeE|02{2?LG(>YmB7@%q8h1@ery z61_V1&e-z3BisXKSdS^<9~?9zpB>z!V5|t3(4m;&QBmeqTD2o%>0`Bs{jElD2ej|; z=x=!hC}6oNK%?``m4OMi@Nm5!Isi5So%Z@`+b#ia-F#25fe};vRP+L4MJ)pcRx<>T z_v3`h?6~jRYDvF znQyUgsT8lNsSFHJsa%+(iU;Qu=To#l&3VDLES3|#wt4bHE7TnM_1)>A6?axFnd1E} zVHf~sn^aB32;|S$KB`DlPGunM#mG%9-l5CbB3@f%;fnL>=$9TNhRWi`2zCGOx6Hfm zojwnf5}9>6}Gd!DsoT|}8N zaDZ}ffWmI)J~yp({vqrO?AaJW58a8_#H@|i~8@4%O zopC?wOUTrP&%KwbZS!>&z$=Zv5IP76NDpbRlZ2iGAS73rN-eXgR(TA!A1<1TlkhUC zeN461Yk3I+e2^N8& zJGA@|gX_$_5DQp|NTwm(k6(*-O4se&IT#1zW1GafDEbmNLmXabu`Ip~piW~8<$opJ zGJ-5-RT4mnOFZK{I;Jk4Zp^>y#54Pk zPISzoIvVVnUQVB;;OAkK4g@Z*-tanhKSQ=$dRhcE?FG-Tljb^^P7ghf)MjfFn`#D< zFrh<6r^sObX}-W2_r|3-Lw%p9Qp4H(9zS<1o;gJB2LEY06`!GF6qd%A6DWbBm zz&7xP!whp4dtGjx%)UycnAmHt>V8{bS@Mk%cJLh4Wa$Uctf;3@;Z828CtqN zUYZLWR=*x^nrr)Di^2Jwnkt7bnwcifa3+q+oA!7-5WpDc)wJpbOq0Zv(FINV%^clS z2Ewf}3I#8yqyA1M&h>i1w zAB26t11UsD_*@5Q(WjwK)M$YEp?m0|#SXL?1)O@j>}*7R;D1%bbF~9C^R}PXuET1R zs(XGAscA~7=e)BOzoVjE5F&pgszs4yIzlDz4lN_0Do*U>(5DNvc})2sK7s+~E%qxo zX)+}nefg(d6d(Ct=E8=E^$lt?hJcYRvka?!7Y9ZK+?Ov-*El&EiQk_SQ`k6o){a>@ zW>fNvra94?UycozDaoE3^NcU@kFE(&VM1H&(a-o_m2%GEBQ*})zlUVzg@m^1&Mtay z-nyY?j{X@VBY$&7^Qbq*RvkxI#A1J74{W4`R{3R^;nc^?grd(A*w4|57i*Xs`a~Tn zb)Cx?%``ij=6_XAc084_D7L_+Qe;4KC{#; z&sBhnRLS9~k3T{)FJoW8@}2Q516J5zq>I(gWA5Eq+1X0p@qiku3cJdt#45>@w-a4N zoACM%vE-QXX0zcQv-RGkE_18Knw(=>&sO&{Tk{;8!OpmNbW{~{+xWBNYjoRw*WCeZ zRR_~W2M#8g$4qkiUjE@MS!c1L`&yF1LlEQqv4@6@@|@9^;Bz9|p2;Tz&sz+4wr)>O z6ux|L*cxIsxzB974nG}wH^m@y*=NCQ} z{PAe_RzT&;oGMAh`w6ywQ(NZS*I=)hFyjk0-mTkVHIw3?oD&-W diff --git a/doc/tutorials/10-reaction_ensemble/figures/qdistrib.pdf b/doc/tutorials/10-reaction_ensemble/figures/qdistrib.pdf deleted file mode 100644 index 8b3e0fcaa25c923e63fc4e001a383b3c13a61e45..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16799 zcmbt*2|Sc*`|zk!PEKgkV#!icWH+`f6vzxSocFxvUH<>?`{vg#Gxu`c*LB_3c3<~>NgmffB&#Tg=aI}G zdp5_TET$mlahAxVr6neB?BMR?dtOWhG@0>;iHXS{a(4A~@CJWf?R*{d9qc{MIq>M{ z@c8(8JJ`AL1SZSScuwx7DIHf(%x%q6Q9PlL|MKLARo)%IqBS9+nc^{dkMaSLm*E5zS}S3UM4kjGxKl;C}8WSn!RF5E~2jUx0_2 zD!ZHgL+MX!Irtd_`W`m-1-OZ6X|Wr?N8U~W{1B7ZJ0PYgCU0*C{s$O?j^+R5PaQi; zbg=j3j0wbx;n}pZ`@p}SYpDLYhM4@x3!V;Q^2hC*960M8xAS&z_XQyTV33K!IcK|r z9)V(3ph-atuc0BQp@dgZ&=6Bm0bf;=)$nR!Hk`S79`3$ivX2;>*{}^_3w6_z05LtKN-lO$=vU(FNZ^iJiOgFEeHR>3{?dM zF%`VRPxAkPi`pN!{3%9~Hx&6UiLLu(?qX%g{-W0+v>@I625!ISnq5`Zf!8KF?rM4Z zV7=a*hmMU0-|lQS*fbLs5i5A}QT^)i-$nDj=!&@TY{=bC+wy+#b$J+;){Nom@-`j& zSBYX%_*0_P)&AIRHa5$g;NbuND^VI^8f^0ZS)$lYf8(|yTU`FLlKl_tKI-;|L_*Hz zfX1pTC~&Yp;^6Fb-d7B8p7T@B!_~vv+|$k;$nZbOhx}0>i{73dY_7ynDENQ)t@GadP-)F%%~pA&0loo4j)ERM1dz zQc{e=S6+SHL%(X@$NdkV4;d2^69|Ink{8!5@kjT|PS4F_$t}=m(!V6Q^vZ9zPdtr76!FvA+oL^sWUSvvp5?+N zD23E|%zM3hJ^xOJ-K|5(V>^#zLcMqF3_-0sYPth~?`L(3*UE65UG6$Qm zn!SzVo03Zxv;N&EoomOx4^3n&IjU_DsX9j=YxRU@WCJddGylzCkw%xZ9rw=D{EN+% zdowML1z%v8_WFjilD1Co5`2|1-fA3>MKy<0i_O(~2ec~qa3abdx7c-XJ+R+3E3tMT z&bMA4zGlaLDJt=6q{1lfgW|^E^Cgg?3zuaUvk}+Rs}l|FH?s1!G0gM zd|JD4q5d*yYGEO`q}Q}<1+)X^0s^yyIk9!3mJwYwNMB_sB$Lu*^{Zp3k&?yqjaDZ&y51@BG=C5mLpIx02CWA?5)cGeT zapYwP6C;fTH9$=pL|!rL2409{c(Y|@?gEs0d{zuLK#Hfhw7&d@=Gx#4$K}y{Stb`> zjal|5F$wFT@qKuXUN!j~ck8O-ZJ%x18;3ZyQ`NM?xUCBYPd2i}ShwjwzB7{$5yXOQ zR)()&a3ucz(1R`&OxmCfqddE;$K}V777iEzCLtBAXApze@_X)nV9HtbD1>%B<_{jC zIK4;frT@8WmF`o6O$&62+0bp83>KluAy!J@!j7=(fZBc+<7%O%5x?Zw<IikOo}x;8aIkdSA=@}_)C%xtTH_7LmMJqX02O)RsOjdC0^wgEC| zl-<-C2yZRznaSm}g>5qy)?!(H?LabToBMIw&9&fH>X!xn96uY$ohXlAQoke{HNZ7R zqeQM%FX7k$r)nt2CkXwPpI*@-YRLrGo$B1r;BXOC-q;mBSc&;1g<`WuDhd z@f;C9hSt6^qBTKIo4+3FB;k3XE!#Z(71`EzGIq5Xcm)d@IHBLd4$q^ z8L{HElVU7_BF3U(4w}z<-mUVP0t0N82hb>O>5Wi}rzGr`J(|8EPh?cW{wKfe0Kb?t z3X3lM)N7C>B@=^O5o8h;gPgCkEoI1E5>*qE4l)23#gYt`H`m@`F(iQGV(^6&XNTNE zB5om1lR1OyYY)^<<_L5$saP6cYY`p$Dv4Z{PL=at zrXmtwZhgkGj!;$JaoBJOW$a%h#ranK_>t? z0#~)aU;XW4RYkiWpF5K}TfC(ST0e4AqH=Sl0!+_L3A6S|)4hd+{epj-IaxI@8R%?L zbj1aFSN!Gq{X4V4>0v;xq*zUugMUtU5UQ#0WvX=Z2FX4%soQOJh*L-A8-rc5ljb3V z(RHgVEJV{<>wbHANc^nG`KH!sDsmL%;X!$CE5jm`vy!l+V6H1BJte)TR0#9MqzwnE z7$!FOIOQDfZlqi=L_o*jM{d2%jtdQTx?o~scrt5hR|C|&r5in)HZZFWNs+d7WPTY}yYO@@jq<=14N4Unl4DU5SEaSaJnf6Mm~Lk6 zdrmS~02>khM>wXwD{reReU>s$qwGgDpQn-HPBWA*;)f_Q!l50XUu+Ei)u2u&PglyG z$sB8)^nau`+vUyuBSCK8sFQ1QU9oOa4a99>ta27?_v)8-k%Cc!%EPHa9J9dfs^Muf z@X2B)Y~tdo2&e2}o$||-)-K4YL)R(MZe76`^Gep&I4W`Js?~7?xIxhKMKUXL%rbYF z6z7z)TXlTQYc4lN^7^5sPozk3{w zVxM0J4HIuzIj1sWzEFyK`w3_I!^M`Ze!yha;1N{wmBqUnH9b_brz;t`GJkWvD2-&m z%drqP!aau)S=8FrL~&SBy8IQ%FR)Te{#yObo);FpSk~a9bwV7^Qc?{)lE8=)2;ZIf zym=WKC3eJ-viqY&RQJrv3z-vt8MXuLrF5>Td^Ob6a<%g~AdY8&W61+1~-Z-WC_h~s0UM zRgHW5^-)6`Xd@==vxX^?;oKG)gD4MCn-SMtLVj3=);IA$Oum6=swNpTEKu zEi{rO0IZ?auInqiCH+vt!udpMNj-EzLLY=iLlM3eoE;kPDD9!&mL;Q_<@LFtGS-Xr z43Qzq?oj=HPvFWps}&l3-{Kxgqg<)Ge|q7RjntCac?no#_hB?`*+cp>sYlQVl$842 zdlx)rBh$vP;+fEaQ-&Rp!{W;uQy7aJ0#^UMvqSr+xfQBeP@C&_u>ry&yYih<-7TpT zO=h|rq1uNL?`XhrmP$>`DhrJQcPC8xa}C8Q>GtnbV2G7Dn>a?Je4^hGhed1`JEW`H z#*|=zJ+i%vo6Bw$?YXm=vwz(+(3?8b8Yn!0oO+jQ2FQ4&bk0MP0EPD@C&1}{_Jy@U zBZo#Q;_V|iwqVlTmfCd6*Ai>Vx$c@UN#aTqP80#qz7oJ|bO%Y3C+&@Gy%2Q-l_r2l zhfwfK66dt!;X8Wl2#Dkiq_&uri8P(N+q&U<|jaTyW%p~jI};&7#`@<#wsqd)9(Ly;Jkb^QgF^>B;jBFXPbv#*7!^>+=G zN=p+O<*`qID6r}$CU~Mtw$wm*s;wVlD7V@EUll|Rak~WW>+P$CxPh>KjzKmghEcls zTzX@Xh~24(<%Ne}%c@8OPq~ifON+o8GMids(r!7m2D_F;ubYtjK7W?R8Sz!S5u$GZ;(8arc;K8LtRdA`T0jIMTd9Dufnn{Zsq69QSe&C)Z?<5ber_$BkEPdhaq#pGO54sB0F(zXau;=j_Xp%b)i?{9WyrtYws zp-~dcFh&zB->=t0Fk0D~wpSACnfr80Rn{%G06hwnT}R zAggCS3)I?J{4HSAK>&T05Tpb`f5gV=_U;Yp3A( zAnuv+*&n3SnNkbUzqKmw<*Z4)Ckw);&kiRX0)9V*|55`DUpd} z@QK%li;tj&Gwqj21Yq1s<3Om-aOPJAXJfxFlq_3!5W;q{@dp_;8Pa||4+p(>t! ztErlV8zjG?_pT)~R3?E?HNRL58MMYD1$%FqXSat;lJNeIyN;elE#)s*2Htrp?2=m) z=Xky*QLv;2`dH{1CK6LYqe##wez}4^`6Rywn=#_W0VaVs;Xer2`TecENXps#a7jQM zkwuV{SHl?lg&6fvW{EtUE|XKu_3#X87`!&ysqV>nej|XZz!%ERy2tK%J@MIh^600e zfoCIHSTc-b><$XoB_K+`IPtZbJ)XaQIOCK_t%ZS1&2hjAaI|qNm@PWjEXZ10&5mE1 z&+>kn790c$msY--P62^Z_QFn^gp6<rqlH5Pz*kYOTRk8NKw%Y@ncj}8#x=X{XG*%~eyo8EtW{y{rMVULAAY}r5xW|a zSk%iFZt~BgNq*itOd~bRjlUb_;lFsM10R;4I!`T>x<{mH%XL zcCtNzs`x~NAmao?I{Flt zKQ%=8Djz*(tQ&u8Y^xAO0&FE!12{a(=>jP~!(%p|`fT+{DW;Jz zNZ=mUbJ2?$wCtd1siV!Ld`g7!*y%5n$5_Uv10JVGOGKNt6iq^2HIRlu6J%x4R4N}Q zJ79cy8=r$v-EAZ=^9b5iLi(Oi@!pUH&YhC3@CZ%i_RGv0QKn34m3QObb!G2eWi_m) zYT~P*^}lJ&A3TzP5Tn@m#tl*IDWx)xPrkgjxoEMzZv(US;u!%&@jr}S#b$u@>bl>=-(Xy%*oV@>$VXSWc(_P8O3k{9)fBmw$y zF*x@baI%3pMUwi8ZtpKAydve-)qJ#-_nbo*cXVM9ISbS9kfS#|`*BmBpfwHspr_OFnZB zYw)7K0KAew80#ZX3G$zt#~8guHd1hC44twKFn$bpr2*;bECo2{j<5`=7NUJL1WNDz zrJTFnpaweO-FSIE791sd_cN*cnAFQ%bV?^Xl6fb1Y_q^fbXzhTdL6W_5qqqBZ(O+^ zzaB`WVL(od@5&lzq8hsBqN`gD-9w1oyX!0rqJ8!(HX{rVY#Y$yFecS!iO;&5j~r5* zfJEVUfRl}GDlo;jxihcGra1x%iSr1pw{TZnW|iw;o+ zHt?vqyt*37BY+HuzMqJke{dy{e-+q7a4Z`75?u6lf#iB26=b)LuPbXVbMcEs8n>jj z)S|cDqrkb|U|Vdo!wUk#Z?x9E;6ZoG5jiesAUY|pD%8KMlfT3 zr!E3aRgqL-zi=B!Ez#q<$1V)}Vufrl9CpAJ%(BZJ!rBV=Q z75uK({Od_8paNa-QX^6TXT2%_n~Hb6eF~i7vJM>uC`^&^_ZAe0T-@b@*5(;6zMWD1 z1La<<>1MeJqZUl55)xtCc3??%zj~`HoaQ0}r|4exk0DJ(0k!ANWDz>;p%N!YQ^QS( zL?QN2J|I5KEW`~lwX;_%l#i8L;fA?rR{*DJ8-p;_2JKNzxKeR_w6xpNB!sa|229am zQd#BeWlji8T@T#3Hx-GB6VBZ>6iM_&tcD{=aZk5Vfz#mK-B5YGZLF%5X4%!u9 z$=hxn7lYk4uh)iiL|Jso9{$ zrbZ`Ml~h*Qf2tu7NMC6n&aMyHcNYy;-$|Kl&h}i@{_MIWjF%oi3e0H9IE&CIxI?Z7 z%rwQKFKgKBgm2zw3X!elp4+}1ou?XKKgHdpoF6=GVv){G}*_W9PG1*b`qJhrKLwF#Zl5WQ9 zr!o>|t?phY8CWw3C7%y%!H7G6=ec;)$(B)ISIxn!wyfsnhDlIamRL-d{>#O zhA{l3_g81Beo31uBoP6EUg2P*_du}w+qe&D;WpAWHNM}5fErqQuWvcvUc94e!~M`Z zNK)X82^s{hpwbV$IEg0N%IOr?CFRfVg7anE-A|+Q8|f5HsXI2BjnD;RkJ*V4`#IA^ zpD=FV7Jm^bKZ7Q<&5e6i(J1}~$MiuYf}Q=1u@`n2k0cS6>g;+HHGjCAPk7887<6Hg zI8j*t(g-Cc0b1GGkS5P0Y;+2j+xqN^FMEikE!pprT*_DslC+C`0yoH1tUW6mgAn62 zVS#N;FYb$$_PP)3TNi!~ua(+6lxNQ*w8V2WPJxWb!~JMhA}Q`ZiIB%xtyl@ma`xhH zZFyw)&ZNx+W!O*OMwqb(aTmNkl0<#^m74*SBzlMcG4p6rT)%5QWz4a`%}5$>E*@Dd z4G5(?r%|5ZlaRDcR<}q3>F=ZKd=*ON>7p!qX%AaXU{n`R$i^Y>DV-I(-Qi0W0u$wA zKeB8*8}!82^g77IcsIyKoFHnwGA^;apmDOV)HlX|-V7riIbIigvO>9n)kKC*Y)|Eo;J^8n?I)IH<#yJ5GFF2EfQVh&x@ z&brGsAY^?qvqc2liiVxS>wF~?km8Et>mv5)?$W&XvGbY)oPtcwD~&Nm>96Gs_(U8g zw8ta4Y7r8?MTL_5BvsWOCRDUIi5k-oFkInIOTT8(LF2l z5L{}csc$NKojQFjwRi_iZ*N(w-~t(xL5DFTgRa;_5|s;+E?=CLfa$G(6J>HSh@2^y zw16=>d#^Rau+k7^iV-B(m!Nr3<2m$2_1QvJJh8srll*a# zl3cQ6IdpctH%VK+im32h07$r^8SiK4gT$9uqRgt zkF^tPp#6?a0WJAkNC=snae{c=O3%o;e|w~!q{0sPjq>1>@|Fv_&s3QLw-QOKE|aK= z7*Dk^dSyS#xJz1fZWnyuLjcOt4r3uh`ohwew_!Y2v=Yt;lG zC7YjXgyiMVBV)7E%kNqqN>(flntW+b=Fwdc3cBwKJBSZbmQig+dXwr8=OJcrN<|z@;2?&gB=SgurBtJpUUv#vV$ZL{yuZL#l7A6C;<@3^; zH;E?uKK%v6NoS-#r@~DO>)u*j2Cs6&9e9_*`A6Jdi5@f)Ldt_QVp`tiH)1kkY;ZBQ zR|wx0iaZyZ12ZxW-PJU~Q<-&%ZMwQ!1Ig>2-Y2b628d{eE?>hObE+(MHAYteV|C^b zTD#|6mZycAzo7}bqDvBf^E-(O_N#h={)W8+4#U%Ab_S29+5<#C&it)li@~b~niTvd zhD43^tr$P4`z?yJ3Y_YXnGI3&Lsyyrs7L=CWR^A)Y)EHc-gS>zG%e7SW06bqG?~E6 z3=;ktrD`yDp66|3hx#aYK@s2>03<(k2a)qmWr%iIi7KHW;z}P{7k) z`hdiyr!omYyoUtQ53ksm{Tw#%+oE+pYKpFq2)b~@>%%1w%vFwoo2(JSej z3eP{4Wg!8|T4>9P{DZF^utw#czDr=u9ZU8Jt0$%s6fBRsgR;dv@e)O5w`|cxAk40=D0I#-Opn@>Usd$=f1RuVRMn_#U7v< zyFb6aIrT{sT*JyY!MsgHba9Q^(U)<|r2UNOy`%^=1jCf&z@GCTEa(Hgqs{``B zLs})g3w{#;4h$@h2jFUvyeJq${T8V8h+1|NcwGdS$`r_E?0{`LF;f0aj^3pqGXJKi zt-uWBj~k;wGgje-XMSTdL<|m{5S?ck-kt{bhYN&u-%Q^HT;T}``8FA7O>-YJIsr09 zCDHGlML)=6JQm=Zn=NK3u^5BC0${iqEOTBP4h{~snWjC_fTr5H zyElL%tW!U^!m1j2o_3WqH8n+_oen5)g#jn_qgMoShGg0c5qa=#>IQqn6v^*4>*NGT zI=NS_9hCqSo+>==8GLa)z`|7#;2h~S>-phQ0ywc%A|Cv(f(Sb6S(jkQUuRs!2q|xp z6B4TJv@Id6x+7uzrWo8fltgy3UP@Wph&K0rpWFvd7!ThJD(k6h`7FAX-RtHLZjls` zy(~+UIo0lqN{@t;A8R_z8y)Us5zn8|1w7SF0}U$w=;!NPBQ#2dBSa2uF`w8RqP1Q8 zcatgxZOZ-z$TUpj+6NwR2E>mGZ}~fW+W8#k=Ki0hiZV$r55w7eYi2w$3UYU`OXH z=GzdA=g497O7t|kUv4f~MVVvFtXtmUtpYEuOfFwx2sC`h5E7Lb9O8N1c;vq5+9&q- zlA_?C#J0^|AQMHn!{>NX$R*Kd@va@#)Z`IXBf>oS8;4C)L)HxkMm#zpqeR>p| zFaolx;IZ6+DCOs&B&smRgOiqm@rThBmuhREhT914fj5)j1Kv+QP5xVV$}0nK?!QJ_ zg#$o=G_1U$eT=g@?!zxw zmN-)Yr*Z%#UIze8ZXtqffE#I)FaYphy4C$7S92^q&JsNQIS7!o0C8iW#wgme4JW-5 z&IfM^=pJivz?2e!&BK916&T31Uxv^6HRm>qYrh_ivWpFR9hL>IO>rvvDDesy2_6e2 z0onOWH=}gfW!s$#TMDm|Pk7ZoL0)t@6?jpAJ?-3ttRT2+0FGaL5ARab_NE^9Lf}9u z^4YLAV~pQHB&F?y*V>+XS`MiCYj#?X6dY`FwdLCnx)3?95~LY$4`i9tZNMdK=nn)} zbdME;B(_ak^7UL4#Z@0dHSO(DU}gp*BX1ASX&7;A$fC{zjgs0@IIE zF+nY7F>KrC58swXC+ooz9^mRvzpX_IZUJnZnbhWmko&25g4_9<_HKzqaF-HjlyvZz z+F4U@dpRk_og_qhd#-B2bF zJNu}=^f9o)dvAUWdS8wXy(XML1C5NWv2yBC#XSgVAx#a%y=ag*I|x!2G8OqQJ1gFb zPMEI7wbo#i^v)q}(deqWCjQ1}XZW{R2Lp0iNcq9CxCa}UR5Nuu{LchL%| zU)5aT+it5@ zp(V?vOym}lrL=Z#c4jjd+pJ|NeW5h-a^a%DSGGX{ugdfybm&Js<`;)zB1sxA>U^Z6G^efhFY!B(uAB&?wD0I9y{4!eDQbD2)=t zUXmNj+76cdegMr12TKMAUid9w8D_CKOkm5Sdo>ithAG4(NU~vukqB!5%>UAT)ovpj zGV^=P^ws}7z~I6QCy#U;q#7#kbYYf z;E-|~X#;`Q+6&G<$kTzlu2F8GVO7yp_7~r=)k2v`Ee3DX*uEWJolyn1I9W9MHn7** zh7o%KSYTIw41Nm?bbQq@=6~^mtqX{aaGE>d94o%+I&q!}p|BAoY$@aR=#HQ5-FK0L z*@E&dT|&9wnfiBqGznuhvE?wq#iUE;)$_;9UaJ*rxzx(=3PXCQ}lO#C-80_TxoiRiHMa0!w+a0t}1hy(9af4u1FUv z8Gl~^On&7y|4YeXqrT<+>}2awm5gTCbCt4Bs=wrvg}w%FMjD^5q)i5Dh3S81njgud zIO(RA3x>v?(aX(~PosaPQ9iHTw|{$H^yOm>e3cbvOZXZm+;0`BguKX(A$QMJ%PlV- z?YT(w_o(MSO6U-VoGiI8l^y<3hHqCiWZY;0~>B~0mbnB~h($p8=ik6-oF9O+#^K${UNUH5$D z$s_3Jomkdh&?SyV#L_4)0G6RN8Pgb~R%sMB<&*6w&*{d2Y#Y1fgXn+#d4E#IoBbjN zlorI02)9510!Xr~+V=nXGmb=v0HAC^Mvo2ozyDk=Y50%z5H|j4ON)z(v$KN5y2~&0 z!FG@N{}u^s(j+!8d4EZ8bBFE4jTL?RI|^G;l7>z5gBW$@0gNa6M1K5fu#JC(#DjUG zfhF9U8#WZ~#P=WSLe?rSzX;VE^mWH0keEA{i@% zScL0@K?lngkhljy_dpdVvHfqMM!r9W-PD!a68Ng^`@)-bxwi455^a~KaX6c88d za+q>E%!T_dbuCSoth-Ri-1X@@>uO77*g25>DeU3e{_o0%#O{GuYc8zT2DY)e>Nx{) zE5Bk33%xxycKv=CA7#PIi@x5j!pa=hvwgDeDbBFtVE^os#DZ6qCdzoA<0|L9Mt*nf zb@m#`E_^kdqh}3j5~jx@+}7-4G1q;c4`^8&o{?zyGQ*SOn^7q%fVL>BT>L?vO)&?+ z;>J}9{allFUGs3e7hH41$I9vDsA;5?Q+k7mhe(=>*G#Hl6y@Q9{IMU#OQfya%jaz; zva;D=@N+C#VVjD<&4rDI0ZHAOc@-fe_Z-vrHMk{Zelu~)+o;pvmh$lh8{s%rwKBi! z?$d7VW0tF9PPab^@0kq-t30{)U^r(ZfTp2p61qM5@ZIc;!&m(a8{@}S?9v;?+(pcO zl~SlEv{YLNnN+&vdf|zsJV@hPF(#qFp6_jIp7c5gja)NYe6L!FIs@{3Tdt{wSG7!f z_>Pp0EYon)U1Y;H|Cq_*&~G7L8!!Gao58ZCWbrTkwO;E8^=md61a$Rk4y@=MF4=10 zk>n6y*`qR0S{1nKF}Sd6PGOUF&#FiB+{Wqn<{b4jkk)Hyix=cOoHw#AsZXeQ9-}3IsnuL*P-GWraqFA}*U0Wi1b8_@dZ;tc#!r;i1`1$QQHSOZnH{DYLdQTJ2 zj{L5PE)p{rO;|&5&)f5h{49!D-UFpMk!M_eBCW!G9ovUZ`XxQ;s+#yiw1$qyIRpsy znhk8!s~ZyP)to49>1pd+{MfilwSmog6(e8*UC)lx8<(_vZaeMl!)+CQ-I6pnm!5F> zoWD#=n9kGYswqguqva)gAu##$E4!1A?m5yIAvd?A^wVO_?j!ZSqy++d;@de9xS`Ku znw~f`S!l<<#A19I@Ut8S{A4%;ykBW5YUPyNa<5ynWp&g0@>eD-Qh&;3wTAjpBnPmo zb+4Ripf}a)bgWi(<;oY)V1-SJ)Tee5*0M0rDa^1T`5PQCHP(=?Xtvwc(&Z-s zrh2cjklJ{Cc3*BcB|6Ahf)S>=0L7G%F`9oE|DU@tc)|VBK9!w&+Ubnb5x`~q3X<6~fb#&D5 z@`P$bYj&TkWK*e7xIG}`LjA5M3A%hsO?__?QHN>wDSH32Ib+>Qj#&|Aex!%yM z%M+bw+>vu%x${Na_o{>pecNX*=Eq91AWM8Ovv)3RV!P=0=+}^P{#T!8kj|JPnWP2A zJVuE#QTF>ZPkrwN!o{bZ@lW#>lG2BxtJ;zle~dFy=+brlC?zPJFp<>}gYX4R_py2_o{mjvAeVU+|iJ^AbqYaMe2< zkeMFht~wo(POt*jf3CC-Ycw_Dqpa2D&5xEMH$t0ugmtS5nb?~Ek>rJ<7tT-Ko|V_t z{qc-;a}XTDkLTSyG1S`-)duMnv?-Zhi}ZY=cCN#Fkx{@}ga)O_eK1*IgG>YsMY>ko6^65r=|`n?Ef}S zO%;&&w>I^E8KrEs_Cz-#hoa@7vk0r+TWZrb72rS9QcWL%lQw^+gJZ{z@$m4V;%j{Q zM4@0_NJxZGdQoXKLjUkMdNH&il!6J>pfo`6wXb9{ra=RPSzK&D`Jw%P8qBk2&(Ps0 zi}n{es^}e+kBOlvMq_61DhklvAXs^0FpN!q{!~vT84qJH5Jq1cK7T7^A(gf)69O@! z?z7||FUDF5vxv)am#k|5P;a)>k?#B|`6WA5{N6mKfn49qjv;dcRGPbiM%tBI>`Q6--|UZs2Hs%Hanji;ENOpXkcIr@DEMwYG8H|RRT@i#{ zFYtjI&}T%Z7!dqH*y{jnvZ_r=fl0j(p&uf+q;dU0;sHMJ0Q!dz$3sZg0k{-J3WKb| zpj?<#@P9Btyeo)Tz{R51e?!O5MEt>w@E^GYaic9vs)0#O5Mdl5IDmK$5HA-7am7mG zViW8bfe)O3{#gjVA;p4=!N z(Ip?CsqD!2x)STyi)0)#cWD) zP+~H9V0FNn_x1fFZ3{#wk#+rSfI{jik68}eud%99R_?7aP@S0;$#Qh|Lcyt4|CZC+ z<^~jJHo$Yr9C54=)S}H`V7`3BvC=N=4s;DBwZo)0dxv0oo08`mq^wiScyACdo-UDk zgnXDE>;B>`-FrHUIwaE<`rFnB~NgbGvm9`xpAUI z5#49=7lJdcRyc%t3d3Hfz8%&D7amc;qm}W{kQ@0N^4&oj#ntjYx`-!)-zS-3{dHOJ z@7H-40sW635Pum6o9K3$}e-lPu`c~yZ}I7UkpAjn1=JZn@k%! z6n@Hrgvn*_d@AQ{z0c!*gy3iWoe4MDQ`zRW?y5n0k|IWUM)i7s3KuJKTVoL-atVJ5 z9SXLwZU1?pf(K<-PFH?PMx^Yno_;8KJRtNl=bBMVq^@p^!0%18j zkll9KxoXk4!ig!|L3HK8?U!HlkyHCvxe)5<^?ncAK+SMc0jdVR_3&l{eSFn|vUWx8#Yz6BX z(DybhO!l~A1LA%nQ-i55n41!6GMKXC3dK$&82MUc3qN03sf% zg8d7Dr(TI9Et?#A3z@tq+#UCXe{NcdVXj}O(Ok@H}kJ|Jih}!o#93x zHP2e0BnX}_f}+HUz9qtl%})@4>^#?*{9Z5QOz!rKy(s~@Z5r=EIEM&VCa0RgE*SLY z*>zHpAqL3-FnGr1A(T#;;*(K?&*vEVAh#yE6Jih3v0UF(qIKBS?KA?L9nQc$V|cfn z>lJsP?z#p;_9jxR1i3S|>Je5Ee)A3Yac%nsJ)EWA)zI%elauC-WuDjsn{19y^*MYz zQxJTj;VAXZ_vQBSnZ;*0IckmB+8F7(oT9+K9P=-E65AS#X`5aCTV+yZl`D}o_GM2 z$sfcTcj$mXNy|eK;^YN*mtc7SA&1!yg+Xba(>EG5zWXVbJY09apn3<8`}s}}>)Dyk z8e&hH;FJ3_(@k4S32=H>0LZm-?8RtHO!SVF9t%CQ0&Xf=x8-xI^S6-VQ0>i@s5g$J~ZS%x0OQ*Q_@p#V52#0>3k$-A*MB`G-&sLq1eQrGU+< z+}Iq&%Gve#a%+K(Dc)7kVhq)=)yymCM`l`Gq+zQBKpPv>aye*5l4xV|wtnd3PWDIT z-l4ssnQL+|pzp$p_+$!v#XgE*HU^B>hF=Wzm&NqzJ`f=cZw#6N$~DnXf^27jnbz%u zV+X;eQqvuJhg+9LDG><%H~?2=A_=Ssil$FQ`l2>o3ocI1h zsYWVBvjo{D?0+%9$2ZxUriBn*JH4c$Z*Ku+KV#cgE4n)S!1^Xy%`=TEc1DQ+EbM0zSN)w*w;Ns8#ee&#&`s*(XIH2TT?;j*Q8t}e3Z)^Woj zQO1RLQ99F!tzQdlHfUv#yCpolF2JvU9l*V4)&Ak%ehgGuz zJ$=9PES9}apk_qvRp)BECQG+cmdV$}x?%P4nkFh79=6F4;Lxu0Ze?$J@iE(k2R?tt z$0OV2vhCLk11Zjt_S=fTNT7joXl;|)pWb0M_MAhdaeA1XyOOq@YX$F=&Elzf_JDveW-|PvA z*CZ`88FuIFKYMlSJ%2z*tu;`xm{g{#4Z%&*3ajl=y6jD!cdk{pA-uts*(&mN32J1n&mU!gMuY)|wQvX?yb12Q2F zs@{=ncKs4WfB!SlS;+Ug4r!gWzB2cOzCCp*vfqV1KQme(C;8wBYYsNW<_Y#J6)W)m znZS+u+2c0nK;IDsB9zC22b*d*>WXE=i8gL$4lGx%>LOzG{&;w!f6$&)`@-O*D+}xl zc(S`zN1V0gaE0sWX>e0(HTY|)R%dYg7i5r4z29bBDsTKVuXi}BWBG%!;cUsU_9#Bo zPtTLa*ym-4LZ988L*NGNIvbwgXVvaLtpS=UKZJ^wss#o)Y%Iac^n`mer|)GVgQgTS*oO}iEF93A zFn4fjDFm#3(l28ljhw@k(t4c()4$T1y2Q^c$#QZ|Z82VP!x!X8xrqToKsUXE@`>s& zsBJTC>|HV6J;*2%n)KpZE-j%kQ!8e6NViJJ_lQ+^7$7n;rz<76UI>mfoy z)t=)q#)TG+rX}Gy{AxAN$(u04{olgMGO~ncPpeovwshF+L091Rjn=-IxMLph6e-$X zeSi}$m4g#<3<=Ox(gR)3iNZEj9`1dEg&zc1@i$mMrnxUo=kdq2@$Qf4ZL;8d^F|&f zTOPl_`n`3(kC z={A3==FmF~+nk1}*IucrdvPGVmSSJZIJBzVp5u@w9lJ)R(X3&jIzilU`|YX=;L=YE z4()-u5hK1L1! z?10S9<P2IdI&yVAe zA3|{p38_$$NsE5*=p5GAUGmT+?6r}8AAtL1S*E412?2;;rGc-D`E`3NKuYD(#5iY(cCo#J7j4WBCas(-T- z9UBHYNgdH+OFbh=kV`oz-tuK)72=!=OR%ObrM}O@S#2DszRV0a+D$#O1J97(AT_Ga zt*du~eytoce$uu&U-eT}nO3N9EuVfe-kT5IuzHHaQ%|Fz38ni$gr|Ls73fz*CM&KNI%AOa|RV?jFNvftUn3bEHu^(-C zPtYFTs0%#4aKQs^H!ACI{lt4qm~Jw-*9an>NWMU-%h{K)dC5;MI^qOxj6bSNet$eG zjMhPR-Ns;h+NRl)$*&69x$XiFH`iaB5+`~Jy8_Fd{2HLygBpN#@8=;`va=p(p;;(+ zN{smGQL!LVUKbBX^U))_m*~0eNF_`w-g#Y2z)*zO6*&4t91AX^efDVNy7S16@m!wa>|n<^l5iCE1E10X{bNT%9yI8Y%-6+q712ffjur!844RX~_Qc`F zpNB>1R%~iEjDEzbFr7E^WuKb5-BbL!thxa3lw|kHsUtyo>sg#xygXLxZN-L$L$aa- z5qK24+{Gd^&qb(JcLzeL5!Rlv%d_J(w-n|@gCB} zF(KzWk*%zQ+`t^I2tm&7#dk01{!hgeI{}JIaH$|cbLQJ>p~dUMr(R!XnR?&kGxU8^ zJd9+aL^Wx9+<#lVV#yk;P6S{R1X6ttY^NMF9yeBc;8zHU+p7#o(Oljf_v(g8{NHjc z9u5a%uko98vm4ewuKN)d*i$GF*~pAPg6p|Fcei5CFkoRzW$wLL|BnG}XNPG3VYtiE zl=>$#>NCU0@Q34p9@@WSY9T^2MAp!R3n0R2wfxP+0ZDSY0KC{HMVJgglg<}80q)8I z=%)Hc>j79-Mv5@{^X5TV7X`_iHg$cbA`T|)0A}5Wr;8fTr%MY_?xLn3z zEpq}VÎUjH~6$_ay5sep*aXsBzQfrvxS+USxPVM72G(ZrcN-Cvr}c6SgYCqiUx z44hOoCwt6aN0q76<#Nfi!(Suk1Rf2x)NsttPP&%aTn8N+Pj$}wWZc-{WI)riUAZ~C zld|MloQaO=<8f(4E|j%8r|c+a2;8DZ{Ja!4Xyr)RZ$5Rx79H7QPbvy{6J5{b;B^n9;g>0xj8FBL8WC=m$ zqN%QT4Qb%3GQFr_K5LK?C4GBZPoR2I!^vMpjZ*TsQvq~Z-aX#7J1I%j60OZ0sj@QM z)c^RL=5-|p7uXc zT)?XKT?vKFvJK4!JAN+=6foX#X=#l+N$&+IRJ?Qf?Ue>CNZeT8IyF!p@W%IYdo(0| zQ6WShPmojXDyMMIrSCF7L~baO#UV21`#hW^`O|%ST9djs{j!^HLU#9YvDARoD_izF zSbRy1EZVkKeLmTT7kx#DzEcTdV92!0KiIymAFdYg*8fC}6{zdqg@HLz(8-H_sX5*; zt-;_F1L*kf8nVXuJ9?Hj9$}lF+1oj(a#(^Ew=#VSUi(u}nQHr8r9rX|gXlq^?A0(Q zgu8y}tt{CJWL&>?5xub{RYY7A=-;IV?noDGVsX^gmP_=0BilGHSIea{|4=x>02?j- z%)mwJzw1hnHC&Ov=>=VeCugWWg+dr=95wXnbFUs4`j0iaI}@IE0M!y&wBe64`#IvRw7{J+H!>jKC2eC_kunH3)e?ccX=CtRcb>mm4sfn&s%5cEElZ5^ z#a=B2f;H4qgb5>Mhv5Tx3u65kL_@9p_O}Y}wElzJc-#0dIOkNKZCYHrUA$Q%SpY5uAD2C4j| zGa)T_QiSTprTu#kA&l!ks`}qM%?FYoTtj#b@YTVN{9qJ_NTz9E{+}wImT5^^V5l@| zjei-1^cb?q#|=G-hdfO&N@0In>OT#}rK|~t3QYc$SPHa@cQdNJRRX}!nG8CUzP?`` zUyOZ{IZxp2XfBy5Zo4_}3vxkmdB!XTOC@m{1QaV^gs4Mvx~=xi$&s?}Q}5J}V+ z+Z~2B|0VRlUYOW7JH7LmX#{?_*!>9=`XqbY*LAX%NKQ#?LbcbQg#=m=KgftEyrsJjDzVP1Z4}-Mi4HyEFAUzKc&sGF}fbblDt-<5*(!0h5d3Ny{AUXc~zw%!uy&cu}`U=jVWZG^>sQ&n)^Lpw0U^Z8k z-YAG9*NrVf)h>sl{wkgn{L;%=Hc)XXQ@s6*BR!(woL4H7|69KAcY#H9YLFZVG_K^S zzdx$w@%n<>_>6Hq%lQunfjh5yO|@JGUo~G{<;aJB=LqKnhAQB^zh2_QU=CFMISW92 z$-nP?WpSSl+vgMG{4l6Zh3qM)r9SklAR=l+S=`Ddmu^Y88evrxIU4E90BmO-Um7Qv z6(?D8W_b(>;+mc%c7@sa9~QLdyhOv`i|1S3H0o=h-w|5MF>b*^w zPZ(>uVm{G~IGR412T*#|3^o=qeoD92WPh|sSEsyVMYJP4Uqj=_n-D#X zm9E4eyQN7IlXK^2GA5O|l_P5g1wWSk>xZKqs2mN?!3WHPNS?=CI~yaIwoLz=TL{uE z`54TWuYQ9eXl?r(^)FE+p+BqP+&$;rh?zp2%iqH=HpyB}J^5ApWd%FN1;H6&@v{R| zzkx@^&D0eq{%gAYw>gP+lAc4=n~PQFo}^1Q^2C_Nn}Wr|HQC4V)D#ktt==hx>CHmI zH>*d_84EoOvLQ8vZ5hdaycX=|Vc_#Dw^jO8ox+wgm_GcZ*rDRcP`-#a7?Xr5A}xbo zrs3Di54=q|BCDld2EU&{SC%WYZ~40Z>W0+b&aE=VtM2wXGEN#lmeAIjfzl3tH$G(l zn_nR@3}DtS(ho;|j#tm??jAbN>rN+yY<;q5cr6kI4HT2}x%&(34dP%%d0EK;ydJ!< zIOdzVzkOSYc++hs=fbUiM}|qO#~GCDwL&u3lrOY@h}mt{>ajS{$YAEmbdsYOA#{F* zhHlb2dFHGAV>zO|!Eo7HF`26=M@WYfv>-|T`<3Ha9WUNqDmqMsT`uUjWuNF|U^j7K zWjX1t3bJ+Hv7)}yf~5Qx$zpNhL~S+IE@bysJltUJ4<`$Qs?Z=XJu+~SCd9uF71H#v z|4`dHOg&m%&D9*h$1iqm-}RX4@%DTXQ&e5Y->?Qj zp?5pbgmsuIYKW`uFg_#f?<9e_9rIo00 zRa~s>$J>PpDl3~qT586NIAkhn<&i*=ZFO0k7AXjv0@rT{tnM07?hee&eOY%3P<4_i zOHv8HHY-YFxWpdgCpOBcm>F+deMnK=#IGVJAed{^osYlyi4(5R;;LD#yUak)gl)6! zm6{_sFJf@9f8Tz3>usj|aM8;Tj_tOz2Xy8PZi-L3S+!rvTmHba%Cc{&mDK!Mc_jUY z>HGQawJuLF1DzdZgW8`T+dT^To$P978EJEFT6kJ;d(a}h89_^S&09`^?tQ=B40=|L zbZ)mVZti~-pIq3(JuKU1+mn&2eZ94rd@=i_fljz`hs5Y!|M^EnyK26!e77Talb6w| zo8WMBx}S5xo2bNk?{JaNYe4-H<7TRYv~iq)+}^>#8lZU+u#46ng{jlS>07^8+FCWPKc?F*A_ z9uCfI_dpddcl=uWXDeMrf;=cXcY@zol14WRxPE~W$ zZ~+~xh6;spF|~*iNcv9M>zt#k(daeLm!Pp5w^35OWwgVgOq1QOJ2#VyocVn(iS?{y zny!lVv$<>?v&cFu5~=((!Izs84U)lp)sT~gzf!kVoJONYF5SK94-BOp%&nV+$9LV` z=oI4$+9vywBBi>*0<{c0g6=0Jd*nE^S3QVAclWWx>hZnD$zQNY7Cr-tnQNtMphwT6 pXu~^_x~J3AIYa;P#MvE+(13C5kBbYV=rK2pzK)@Gk>*9@KLE`bO>qDK diff --git a/doc/tutorials/10-reaction_ensemble/figures/titration.pdf b/doc/tutorials/10-reaction_ensemble/figures/titration.pdf deleted file mode 100644 index d542afae3b0d78ef9591725f73b11601300bf2f0..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 14164 zcmdUWc{J4D|L{<@K`ADYvJ_(MyOgaIDP$XjX|nHQ-^NmwY7z<+35klqkTCX)Eqke< zvM-T6CNvnHJL>aYKHu+ie!p{`zaFO(_rC6HU+(t4!ba!LOUgg(s@fbxI_ zrNbworY9%mjAUd>K9wrSUXdw#WBQ)KhYM>ABf?e)&`LHKTO`N1S>EqFkJETk55rGP z`^<>)p3LEfe@8qPCrtV>&C?0h@FR_q0Sy7^vswbO0@8LU@IR0k0Dk{3ZDDxb)xpkhcTA9?fFey- z8VvmV!$avGJOrdogV7EG(ncsJhh2U~C?5w;KS1&?16_8ocR^`;1qoaQT`~fSs;W|| za*8rCssaiM;7pc*;jyO@!#gX;N#_w{uj}IbKC9xOB)j)gAdBn z7fs`A7yMVZo`9c^|E^3lmCu4v=NtlD>>P|QXw!!L(bd?&*UR6>&cRoJW;f%1(UGP$ z&7S|_b(MdJ_20ZsWAj_QKL*^@$;AJ<->!wud--_m_Gtgp7)mlS0t$*Uf7Jh{EXse$ z@{ceQz9YLYQh;OcW@vM`CgJVT$q)@q9gnyY<0vI{BwMPZ$j4V%2hJuxbL`e8p6J!# zT93sh@!osh$+#HDUpl36)D3#DFMbq$4M5!qMHlA9=J_BCszL%6Vv{g0;{$Q@idi8|MoNFDl;9YSWWZ%j)X zkCIs+|HIEh$mwH9B$A(>-)V^%qRC!YU$3E|@ex6RLauQGKtk*Bp+r0cvL=kAgG6ln zs^1Z%R$%rY1@K=c1|~LHyj#+HR$!3(VH`>5vhQ{#b*DBjGe8W4ttB@;r{M!eI0MEC zSvszA_Emay<#A3rr2k?hy;Lh&?~CAHA@P?}hc`U-1MtHNkpktBG)|M&^jdM>sx(`D z)y>IV%Fc_5)v;frVBFfq)4#^GRh;L;W_=mw<4IMTC&};JP^k}Ve-bg#wcXu1z02%I zq@ScBN$HyH*R_M&I}0Nf=j*=OT&^dlEcZo!D^z{_$@IqPk~YJAm(o$ME`Kq1i{Z55 z8pZ09_tR<=pQR6+5aUhVCNGcJ**4Q{k+`yp-_ivn7&j?9kIwtIzn`&r!BsQ~-_9w1 zTbCnpszX-9z&_*-rkNMYTBDbCU2Ci6jd8r$#7_4RZk;k9aSnBZ{hU}I^Rd*2T@2zB ztRmMTi?dndX8lmpxsRqtOY$d{DmY>o82Oea!X_@k#IM(WeDVNqaTY~)1ACMjti8Eh z{buRD`Yz(v6Bn;&N1wfMFI~BeJ7Ov3q=bS_R(86{gA|2Tk%(>3HZ^>a{U-3M`H}c8 zJrgxKE+D)Syaa1k6MX4Mk9=EhUUnP{820Jhi}7nq{G5CRp#qnWcK$NDP1I&+5|6LY zAGnxy>gkPaj+hwzq&$`8%}h?NS2tyx0Uh*}>#&z6lMB{{^!nsOXdRNg=IP(5JuLC^ zDHwk#J6%C2>#}MsmN~%CoVv7F@ua3j*BB7CxDL~Dce|>}9K;fvtj%?Z{z2YtyR(D` z1yW5vXz-KQVb5w(3$Tkvs%5$XU5;pbos-gKg$e~l9nL1ViNu)`@N7%|q7yir{1lC^ zR2$a0?*zOgpZ9zQpw}&Ad7vw;JjLFOXTlK3(CVw17`e@cxNerhrUjJc(|axSBE~)! zWv0C63JsoEaNFgXhp=VYTeKWTs6jxy!~a_#pEnBC+oT`Uyey{Ca$9>DJ3 zYfP~xN$fEN*qaXun2RzCn#9>m0PT6**Fvwu2QI!AlDB|m7Z4fF-i$2CTTB71Z2J3PBN+MrQP&-(3+fLPz1>yej>#uS- zqu#jS7k(E<_WlgLbvGzNvo{z^A$Hc8&8YX{`%R9 z(kj=S_ap}KIb_3TG~7Uz2Pe{R!(H?l5!GlDZ^m!JTDp$AFJzvv&$D6xH#Yml9)*fW zH(Q!d!P~GyJj>847ViMsOONJdwh-Vyl_yV`H#(x%mX<`}V_?8K7{DuPeV`$TXAM(t zuj3M1r_J^7p`Q8P4;m6G&K4{PE(Zw6d6$Q19gYsxTTswuqqUQ0Z8y+Pc05dLCxUiF zEd0RT-neMM=d!nw+Ow@Ml?O{sc3rc$P3+w_hBY#`H8ZNqJ~@FBgz6tZ0oN|kR?|8R zW4&TwC18Ev9w5Nj0+`RTAkLKz)|_&9$h5X1OcwAu=Rn&DcFa@CqN_r1ZbdjG+Wxh&<$q9UIp%=jCeL5@v{zqr$aRSy460! z3d~8?>QHv^r*i}$kuEq0U|>*oakep zI2xA@AP8WthsUFuoPGekD;37D??8f>kg+mxnE=`k#7UikJ-Ws$37{td^w2I;H!hoo z3*`Hfh89L1h`#ZLBjy`68ZN|HMBuP<5rlSpbGAK@dOXVdZV3%!1wg*`i5;@(KsxK& zXdN``*lb>~j_&zSKsVL{MYk`Jv2eLpG-%8e>^=*^I$5g>ARdf{-vo$v<%$5P2?(1$ z1?zWk2I8<75QGTd(5s(ngp{7U0ln;ru!cf5U=S3~HDEm=2ysKg_-tGW5ETUsJ^fAq z=8S$!-QM2b+PWltI_z@A^szNE*|&PAx3Rb1ZOWQ!wriUv<@VFREU+$N$g=u@7MhayXcl zEDwK0!^c0?Wb+oXm&2oyFgn0#nehNpk~?|rOBdB zM9;|aJ3|^2(FDP|DnXUIQ?O+h+bgMkv4b8|oL}hK9ar}WgT^CU8wW-*cdh6y5;l$e zkgT-RS^02v>|))xw$+g!CMV;erNB)7oa;bN7KCYDS&6b}?)CgDD+y%n0kJ&u9dL(g4R^C&q68-H3)R1JI-4Y1b8`&Y@Hv zgQ#&Ml0M$oaz$y(JphfYzGO zE-lu$0iaDa@PWt$MpUYIH86M>5hklD#~dN&2k>aKVHW$Y=Am)JV60bfjnIoO5Kq}+ z=a+ySfy{`e$#((;_7LfFsNh?mGR3nILt1_SP?rLp`TJ8)3O)$f1p$*f5@B?}Fos=l zUdO{=^q2EwY@`By)Qo-K#K_1t+N7QOi<&oJ#_Ry|0sXmc!?V8}qSe3Nnm$G{SJK8AoMMAmYjNwu^6q51^sz=LO&QcY zoerxbolJUSV4{@RdAAVoL@kbxRYqJy7KCqG=ET+Yf3$Fz&vVwIK z(hxu5`&!5dvKbWNX7JvsRL}wxax9v`%hf|SLJV1`8bT<3rmS&8S8XenqkRg!>Sv#5 zp86AcYG-!aZGbkPKPjU;(zW-Bg;Ssz3(b)Y(@Z89t6y8AQ&CYpYr95y=4{&_Wp(7U zZ)1M_wpsmC#aAt{4m?vCP!Ztc_Dx{xn*|Lw&iqNng*+|AVKIZK7>7M$Se0n_=ldnK zz-xve?_R6YA}jiUT_~`2Ca50oYi@?|3tz;aiTvqHY zEylpq0YZ96ebBR_b;7%Ed%ol6SBWufDKmnj)snMeeLe7^`HUOX6pUC}PbL_A$s9^W z4detjgKEnBt>BJ83=ZHosQ)bOBbkx4gYKZ4#iHT2rwS^w>W3-^9+`KsqRw=|yFpA? z4nUzSr232NIh%$wGWJJKSbJmF@uL#f4U6OH5h;ykRny7TFHeKZ6@h`qsc5t7_()$3 zjiBnGmE`0HW>rJ}Z^W-(T(g% z|A!iH2*rHoprMqBdo!pbRjyNmhk!e35Kj<9Sh|qu&TQnW4q0U7#byPQ521YqJ~WIT z?M)w`AHznmAkZ-H?)p=XVIM#Zc%oY&8ugKi*vUWycD}g^tg!E6v-lXaD7$P`=VR;! z|CLJ>qH2yyOK~n0Aw<}%14W}gR!;;#_P?g55K!?E1LpUf3yTXs)Qsy&e}nVGqM3!# zN5U5aMAuKW<}WG~vLEQ^ysA8g-NUwzwdf8g?t8t*DGM(pgsaq%s|ow3Zm-V9w%88a zXs^Wj#39Hj_6$yUR4nx8@LtTS5@U| ztJ|S~SBXV4tllZE`#!j!kDH2lbBQp1m)y*Drr z6}~{PCgY`a$=K{FG^zfRfDyb@PY@akQuiA(;yo|tMbMB;L(E^t{KD^CX;^es(OkWB zGUnub88hfa$AYrjCHE2(MDFb1=J<`&{0=K3>`hBLI6|YKNRs^Kwcm4VF7&H7Pr=X*P$xImh$>_T-~k7CTmWC zq?*Jk(ocBnOpD8`2LDb9KI^imZN!i8`UGh$3-;0gOo~oOw<1qbdN4!RDL0?IaAVvy zgSjT*x-=0cgYMf}p-?Qn$+Kd}n3uuLwjzeF4;E#4J9SA@#=B?c#I9;Y-ceH-!`eth zg9+R!;QEg1%H*2siyGjd6!A@m2?;>uWrNt%=BC=q>Ys{38JN!g{(gbgT&8;35++oT z2wQq`HB^g0Pn-Qg$xcpQ@ow8*xG}J>HZ_APgso{S;Dbhdz*;hr(OI9fm42A!OLG~hO1oXK-qgoz0()EEl5Ji9dfJYCU0H@@1uYZwx$ZscD<^t7%1rio3zCX1?6Z+5hd{wNgdeEbrKy`ZFps zLUGOY*`;wha9VKOnht-;yGrlM?9ZP4!%ORJJ^S}{vv*8iuYDs(4gzZ^_tL3!5bKAX zJ*wK+akSmWN!o4$v3j8B^EHR*i*@vjyy@OPw(XZo8&t2^C>p;8x^MWVt~0P7QUWu( zlEprR;RVJU!Aw0WLsKF>dvPv|1yeDfp1nXy#+R%VXB|2u)0P9yom#F~`b|UT_OrGl z-kwM@5QH896O+CzJA!HH7@EbYOD!q$Yv^`SCi;6gm!0;0Z8G*zjh}WkhM`l{Z8&vs zk#6qpv6kp#(dgz7yi~IQlFm0iGsLsv>bCXp=gI6vZU=(|uzp#cBBngQ@s{z~20}4! zo5A)by7^+#*u~fMO00^Tr;jBVZIa_%ZQtd#AqtkA2&ddyBj-+eugx6#G0)FmGqU8g z&{KD5?xim}QL%-FH5@oT&6I?pORL#W@rjON-LZM(LiHoUB6y?WOnYzCrYq8|N^kWq zoW)xB&M8^Pr&6uc5L+&S)=#VLG#S9mG+fx0C|!qX-=6R843`?iI#AKVmUPv3-`A(j zJJ{GLN>5@H@D_bBG-}!gr*~#Std7*}glzjl`@}X&)si;vM8gj;Bi_<=r^OyixS(+{ zmyoj7L|K^EuDfmwz2XbJPQT{B3y!;b#zzIhqBs=yA)J;3j6^>gfy@NkId89g{vt;H z!>0ogH>^_cfEAC4v#p|=l!-r!tLM=VNo(A%1{>tBR8)@TJ`If@%=F0*7VJ0vi06pd zkT>z4DjH%ly507QBj${)!-1!IbMj^v2=wtza)4oisI}QV=@RKcQKZRU%S^#i^*v8q zY)MYOFJ^(jmfdLy#>umN+zS@S;y^p5&{k&5|7`m#|F&A_4$Mm%5%L0(V24k+6_dl{Pt>5}kgLUD1 zh;_`5M>Cv5j~676e z?eg<}{1}+Rjv0Xu(U{)i$&NQPHynd9PyL*T3J&jVWx}aaK;1d2O?Hi1b(!ws^S@Nk z9t{^rfX~?JMr3~iyRcF5?zET35?U7m=0*%cAJoNjdk|jf?~R0G`HRlw zfnD#Yc6VC-hSODlL_pohlFfH<(M+pM!woCn7DPyAtT;mtgxmklRLV)Sn`j2>j+r^+ znPUll)oxc)hBLRsn)#+TDyDOecA(}-?6jdU4eQlAW}zS`I#kafzj^keHp#22q$SJT zftie$*6HCi46H#o4H33Nq6vpJ&`^Io-bD@X37mo*!FmRwZ88Qk(n8T+OKOr<%=Y#} z5r!1k{U8>uOlOQ_ln;CH`-@&rkq@u+h3SU&#+;$d`aN>9T~;515{mjn3Q{QI>t(~L z@PabpLlNt;^8|XS!ZKAihq_t*vh>i-%p<&^t9tLaiWLDxNpBHqXap7uio^vatJJGmV|I~{F`h4$@x>G?^+_Xl?>_i}St)Uh64d<+Y=M}W3OHrl% z8U>g~&!6}lgC;RkpszYBrz;}(>b6O!oSW+|wq@5aXrm^2@oOs}WJf@fd8Ox_lgQdP zh1cd_m>0UDk!g3W|I;8|j={QO%V?{AtIfGQNENaalg5$D#O_AWd;e8m+( zSy7mAvwb(Tfw+g0a3uR(;&AS@~68Ok|x*xPWNQMRGnb z5%%`ZN6UStU|nNbw!5ye1xmk_)r#zlqX~bGoKH`LX}11u5UNhKL9=>t0Ho+;u)oV! zbBLWtnKoCznX%cbXn3UF6s;lLg`0a&0GWE+53)QZ8VBz!r?+r)g`Y$9zptqh8iS4+ zz9#?Hk#pwIH$EB zv|1B@+Mw6=`2K0^^?*lqUGfDa__}cV97T{}R2o2&_Ce92r zJ6@)=o;&)<@ZwBBmD5(S*7u7C;KC)3ls;LZL|tqdKuq%$lvPSM=M1K_n=^`1T#t;2 z)3-rKg9wVz^yseP+ zQ=0AqR!^9Gjc4%@Brhj_Qt>v)F%98#pa?SVj6T*`A}-yjZZVyYyvJ2Cb*%h+R1BDp zOmxj5&mKcrT}q475UoHg=69wRx=}H(qIgkjx+Z%iiJS%^6v32K03pf=9_qD1DdOXE zJQz6ASxOT7pm6J>ptjuZNw?57i6PaTiUZQJOf_kL6)QmdvjMBI%1R#UH9*DS_AnzL z{ns5Oa_$pdZI5ccrYG>p+X>#|GAD>%CeRavlI#Q}xy;G=Cqb|kbxd05vH}^}f$UOU zTNb*jfVlVtfOLv%C+N=&kj@Q2v-Z7sr#&lGfTopXBlc2zK^QlnuoonTh+~_!t7HMX zOby7e@%SWbURg1bR0r%Uo(gBpdkvHYcWw|@ha!obcAy`)Sb~TL*y8;5*WfeZuASLuG#gtIZv^#X0b z0=o*x-#X+oL<50wewUudSXTfFj6@hD%&}Kn%F3b(hfoz^e)~A}O(Oh|IUxjWXTYk` zMGzrTfL=Y=28~pm`lv010OMBQ6hbT<)g>pbfLZ%Eiy*&}#L5y`VdexY!`JlJ?lC%% z{{3ojeX}<7Gq@_mU9B_%k}A^h@nv@ug>$8(SH1q0D3T{i^bleNQzTHuAkr=#*Dp8Z zWhVD7`1SlkMJc2RnP0pv?eYB}jdfxHB*@Io*A1&z2~R0xGw`3rm97gZj6oyWt)AGn z;sY-;EY=SX8i7T52OsC&Yo&cJ`q?R9K?rUfYU`gFMdKU*k0oc1*f2%p5Ln7{zUgy! zc5fsnb9h8#-v`w7CiQj>&-K(_9BI!K8Z3YIF2>_W7*cNl8VPRwdiJSn5Z-504S61U zSUVdn>HefLgdH%i1)!-vvltH7A`t7rQ3`C$*OsqN+RBiwNp!V6{^@0dMh_oOgxhf` zUp82vEVF5ws~8MA;epF;iMTMc}t{3EmwljJ=PZoO&)q`SF@ZJsWqCI%eA7fcDflg6B*@iZMO__L*GZS`R>==Ae zFmnwcN=lldj_gPDe3&sKl2B7CX8q1qgcHu{TqGODpL2+HlOewm@!aG2D=>fO;BHBAT8${k%!|e~ZUtk|X!V{!VXe#^0|#HG{HYbv=(1{t z|Ngl;Mdt!4#i|bzIK~jKe5O)#=KKy1qguV=*^cVkt;vZVy31eD2TEdWzJ3d~6-Jgw zfAkd80?E0mXWlwjA>k9@_+xA^AtJ)ZJ~(YLJZai-s}iF)d{Lve^GqYrji=(e&z_fV zGX=pnI-V>pjr~&XuW3=d;S{{y+o7kYQZca|XAyc&{z!)Ch36_e3&u=E-JMr;6G46& zKE;{%*N6M|r&_58B`2hPOTui7q-_@3zC-vXV8oDJVYvDV{^zAHr))6o$RG*UyM)`@n>V@^=<~vEAluzs7)1uQ@ zzle-4)k>y2Idj`N8o&H(Bc4a~^j!|$>*aH+A2V0ENREo)Ca5F*zqo9TwJt8R9nD^|M)(>mK*h>A`$zU=I#dIl)|>8Y)D-nsQ#%d)R~u6 zqI~#>)h$uM$>Y`&eP1TO8cBb{sZwjHkC8jBKj{0XaftqUUU%m)|H6osZF8R!SAznT z9>XQ9lYzDd67I*(CRs4RznPVa)SFQ?>715Z&E8qY6slX_G5e8^$Kz9`Bk4kEonkSC zD%Ps|n1xMMtu^;<6QHVf#b&Ya5@;km-~Ckl#3<4HP*IaOM5-{2M>Wl4g>?U_gjNB4 z);`#!cG;trpK;1~02TKg(o-cgLPLkp!sc{FICyu*BlIN{{xCod(zET6-^mTU8~pI| z`$1mI+)NZ|;bX-bqro?BLgF51Xd%F9U=+K&DjlEX2@c!0&D;atEI>ajX+&XkQnu*1 zjjM_gEQsqc_lNjE`MUau8_z1|wkO$?>Z0G1GrtXzA-!R^BgQKv!H-X8&aZb;_0M;o zA&>Gfs@v9c^k`k6vrbQz8#jEyMVipp2pn@l7fx5|Us=DE zy!Ci%LeWDIZD7H_qTx=y^Z57v8R(72d135M`oRx8;a_xG&bd2uV{aE*_egGl-vn&nhlHF zTl$9k5A*$+3%CDKO5N$|dR26$%*WR)BUly59#F>rc}ilwr#=Idj|v};fCoOQVo1-M zGs-TeYfLF&Xs47$C>~?)%gLHv@~&;;eASZBwtp_+^b<0v-NuIO-gqJ7=&QIJaWQEH zZm&a8XJg=Jp?sq~WIQ9C`+Kfg4IiKH@A!4oPaCYSuBNPzCi+LKg*Ons8glFIn{sJa zv=b>+>YIm@bO{G>qq>ySKfGDt!kL(I$_fQR<+uE3X{CE{CUv9-0T1-=DqV1vN#^s+ zzmub=38%WL?Md9FwpsicO0i8u zVSY_>^4v<2gof`r+aCTa{q;AsL2$g*M;X$Bqgc(Cy>M9#`EE1jj0*6Ors@>0)kI)j zqB3i!1yCn)&to^|dm|L>VX=8SASr#M&m>M;iM!kn03lfll>9`G7)fqgtgFOf%--Q* zsNcW%D)TYzjcjZ~#MGF69-_rggXC68{m#teqnW|e5zDFG&aoqF+I`#st@gs^S`6ad zi1+&<)M5%G=A~%XNP7nk;~94Xrldq=e*~SbtE(HP93D+;t>FZG2u0ic;eib4n3%eT z3XxXaX-)r;_D7+|G{`~EeYx{e2h6wcFK%z=28*l<4yA4f{wQ9|XsGL|t-|;dw+bu! z=^M({1_tVjzK`-g52N7K!z3T!!3}c>RQx@4Vfx0&)&!}opY@8Ww%^_-R6AJ=5bc_dM_(ag#c1!xn<#@1@S%A8OptP>hcs>SinYb(+vAWNqL^5vO*o z?W81ouI~)>M|_^FFYmj0Y+NZCFmdgrQl>M^0wd}pNIe=ol;v59p*qiL^5pqV!^abY zE9E{r`2+35q`*xXp2n6&ZUZl*DK=Ny)R2{e`vv*a# z^OxoTF^II27t_j?rXAtw_QSnn4`O>!AXakl82@?e(v-G^*p=RH6I@?7=k8_klR0@x z)ZETSllZHFR?Nfi5pJkOaK4dQ&p2Xr7^UPsXD5s7S~%NZ%KLLr-SKr_rlMuJBtr`M zZ!{{OY3%%L5_fKDz4Y*Vh4-N0?&M!(ThUSXx)qzI%2#DG@rTv~hcfL*qsejG z3DSI*ojtB+Uo52Z@~7f$Hv;cjFkGg5FALe2=3Uh=m;co_ncze%Eq_QRjc(zJDitaU zV=!gK`&EUf$dl9I7|Ir|@)1OK#{O8uHpwTiS=`w~Szxl6r>bZF8n^KGh-*8U&JsP&wek)nS8+-D0o}_KOkONr( zfw{>bVXJrysQ#Kslei(OR}&3GZ-unDL7ANoqk|X-AQ$G4FR#s5?~XXSn!vc0us@>~ z33>Z$5_gon&*|{s9*BZ=9XJo6%d)(#2|o6NbRFCSz8-@(KRgU>Jn1z{nFcCWD0(I0VMhmT^-zr~v`=zrUkHOwfF-=T(!qQl-b4EYv0%qLhW)SpKLmW$H8l7W zF09wnHuls6_@6I#f71nYQss*}^~;PugRgo2(FbVDU#QDZp^^Wm98)!!cu1tPeM#AA%DEHKM<{v zT5SCNa#B#8g2hHi^-$-8R*>6x2jYu-!NN)E*Xfa9Gj<1b!=m2+Lv_>-K<`?6H@R80 znpn@i9sIqJDx1JoC=jbyz-$NC$oz)x%_hX+Z3I>+xFSspTs2Q zkTx+Wj|)k)4W)tDo%?nvDX*R{1C=MQPVkFp<=l4AyxZf<_}|R`-}>glFqj&@o%kc* z-5q3IA`dk4frI7^o$#qxhAS#!z%3_^4x=Egx@bCbe=V4^hpgE)A<3_WM zDikLKqLmZ$arJI*ztq-z241A1H-hN5Lk%DKF$sY@1E3|5Z&DviHS5v_FRn+c`hhwC zf$1*4$aL*pNcLT!`UVWPlQX~j*N=QeFfq|ZLgUEBb<-qJ95>Pi+Q}n@g?uG_^GXm+ z3sx$wvUQ1&vdoFKfV?GrGhFlExw^Q8Yb(5@K4#z4F6#U(f#0f@4szRltE~=G{d1bB zI!=KSJ#m{v)!i`f-=v2ypZv6e2xDb}1|lc8AlELm+NZ2n^fs4vWuo>%R9ehdNNK~o z5Y4@aHU#BoA^1L|ex&dm6!N1EV~$bpIHMyWDe>Q|_`m538rj*iZQOVDKi_-ozOkUa zy8y2#WaYqvoWI_D(BALR`f0CU{`tbkclXVNzKgvtc+Eq51Z4b=mplHnM?e4ad`9+P zDVzn5q5Qmj1WwwaTztF)WTj*kq?80sI{W#dPf3FpNTAcn2ZeTavGbMk@^KOq+vQ{A z<7Mw}2Od}b@6rDGL`n9a@1XwD;4I1yKhI;``^7(m#=cd>kC1G6J&lP?>-H z2q-Bk$}0*u3jFSql~+^-s?mM~Jpc5`$jGXI2P}W7Cf-| zTc6y2p{J|@o<;p_oU#&7^KX4B|1l28tN^_AZ{uX;mE`}kPgx1PDf?@jpAX8#-N9!! z225N+9B6i-y}B{>@&dm8{a%XJ*7I}(KK=dT3)s@v59Q;xYe`vo6y8N>`(3YjNMsVtSUWSMzFMw=~LWl2I3Av9u)K_+5`RNqPy(ni^3 z-`A2oQnHSn24iP5@8frV?|aUB{^$L_=e)-_&vQTbv)|9X_jB*(d0=+ISV%xh00x5z znGg&uU@$lUgYnv;;1J}~d*4$q7>}0Oxr;_zE*Hk|!{Km)gM*Haj+K>_SSCCx334=gzkjGPv9-E|<)$s^TW!<|aFFlkr?Ko=k?3ow#JDs_j?2 z7+DPZ!K`ko(arGv$^D72E*6EweFIA-1Ix?5i|z}S7v7v7e>n9dXoKv_HG<8P$=pyg zuGk}3YX=tw{c)pOXJD`v7>vAU&uthi1_tA*#@NWk$icXwFc+u=w7e>Z&s{ruEa8nbk3J-&^SXGcb2)_M*B_}6Cs0u8J?KsH`BkTZSI{_4#2f=fM48V3^ z5LGY=Y&RYQgQIwzVE?_-OE5t)HpUl59nA#%t>jUYW0`=_fE;SF9+XCafW_af7nKY$ zKv(6S)cm^-MMx8&bf9sRh`;eXvqv|EP3W5sO2R0fIcij{P8X9#tgF0+G33n^6dmbV z;XX)0EnF?zPY|CWrFta;mm4S%T=?8RqO^R#oV#>o>{2b-^4A{>V;&iQ@#j=C6Nqg- zHM?t4V5C&?rs-PW)d)(`Gm%VP%N!E1>y~_`Bm$Hj?x6gWh*@}G%NhQA<;G|mcr`&n z;e^YMn7!N!jT4zCj zm;PLUdcoQMH8W7N$HT0iqGxXX!J4Dvo%YI4>;wj3T{m}H=RR#%x>PgC<~;_XU+3LW z^|<75*XM>ks8rZquRFHW3sfmMN>n>&R7)6J>d|%0V<$lP=AO$;3o3B;#eHmmQR<>q zwSVPJTMgPv8^-b=9*Ft}Bt_E5+lUeFk6!D|lncz9MpF!<^mV=k={hDUqC_PH!W-@8 zx&b~4S!F{?7QRCdgHZ19Xav{jW_cTAz&!b4+)5SNJP06&y)U@tSe|g*I0g_&zD#?A zt)UpBDJbI{S!tZYbQsWS{;@GVNU0;o#)GYbpyceGW>fG?q%*0hAjfC|ZtmA%l*(KL*Hl8n^a5tJN!jKtZh z{}6=+r_(=7Pl1dRR>Y%i2osPtc;TvWoFDkL5HqG=^z5k0keMY1ZLc8(g%ZfKd#)?ve9ZaWHShoB`!9dqn`_m&o5LL14Nu~0@_3} z>xmD2m)J&*{@m$EtnBH|W)B(Nxdhec`};P+fV~1+wt}gr(BsaLr(R=kWi`#m#?5!g z5q2NfXwxY%uV!&WJveweW^$RCC22sDMTzph3ECM_#{GLx0INbl;B&BDy~(mT%&1S| zW=cA#Tds+jHWanr05fWyuhG_bL&DnoUoDf1n>qi4L<0JY-(wu~E?s9%sH502(Od2X zkQQwzeKZq#71BqTqH82Y0;pVF|J0wf07mNW&oZ)jx|KBXxw2`qn6u^`oft5OJ zi(;eFw``B&F{AOi-6{R(79Ir1dGrvAVxMDtZT}O+7LcIaH7>`dL#QJ!?m$+&V}FnE z#p#PH7eJ_Gse|Z0&@(vUAvOq(qTE$bgsdpm?N0tp(?H@db=ZSQdNOiBfU+jH-Sb7) zKe?Y)OnGyX3{cgUHscLpAjZ@k4SZM&14!*G<3r-ghy?>)l7c86Q^&`@o}vdK@^!M| zAfr`1;QMCB{!ujA&t3N9xKOm16%m zQpAI6 zApxw{dSuTnJRON*i~k{cN&`36!+^q$IpZ40?ZvSxu!I1lAjnF1f}KKQ5XBFCalbxc zQMo(7Zy0=jI8aDvK_h;Wfzye67$_1hr;Qm$Zxcp(;eNe@@ZjK13|>`f!P4oOpputCJaZqGAT@ye z^7YZ#R0PhVwh(Ck+B1<*Tk&Z(HSfE%%v$_iWb?E}dH)V#xv+0V(tC{9rq6Gf$t5}u zwA5sos`Y?Ehh4!_PphO6^j%4U39C$CJ#i_MwtA1@Y%YV3+rne8{~0iL8#x26Yw^Ih z`rOn0u%EantRBia))jK48$eYFf;Yuu5p)qmj^YU%CXEkN?vTL?GwNs%#4rboz`#+? z+w&glim8iB<0Ml~gag*@LVz0yXhKkWle=lm9VdBQc%auGXx=7COdanh9t5r;9FR~L z;}icVj$-dE=NX-}<82Z^tjAq~#14Y=EF~e}9mGKtu?|h6os?jMTd<${T^V73G-g!# z(w5e5ml&(hkhd?Vu86 zU!%hfGyp39W1hNEKKRsSOkKojT&u@zZ-9@wNC#?^e<^bc93?Iwj9z{7JJ?z~t4a&9IshkMW8j z-J9w5-{M0Tr_jLK)d!1VHSCb9t62xyrZon;zm5(DukpyY@(IMSC_#5dmZt}URNCg1 z1cA`dwFPwuo{{XqR{#a&ft>P>i`_qj}m zY0p-9-%2xqTpimCkDTMqd3nmM$zGJe&oecix-Nd+LK$aXqPHvF79jl|ukjo>=Dn%k z*B2u&HV|U|#}jKD$0|ttTHQ4@F7;_TKcdGk-kL*V&#-`Lman$LQ@$kq@BPZe5g_$z)mJGK2DqwTd4xooTV<$YmI zCz-TvG5mb(2~?^adFOI^=BCHJ=PJM9;6((zVk`_Ri^mA@xX_s@{G_5C;Bf>!KxXf6 zMh|}>h!N*;Np~EMx`{EQVV9hU01>={pnr&$Es@#CmQZNBb0P2t(^xY`^A5Ol6(>0} z=G3f}p^G&#g6-K#MvjCv=MlFi)|J+HPR|x@`AFau_*3_hAD%~Z=4&2{N215p6`pgT zMXuKv6MM5$9D(7?$835Hlsq?$R=pqoI%{kl@PHPpF~JE+qSa6CjJCcq_|1M#DC`-u zG~%xG+vWUy;?!YpPK#2YN-Esp_hRpnszW_Vh?B zYUDvif0<)op!bt6)^pi;HfO%x+H+4x5GPx;{uy@|EE!%6X^;24?6D9o@Wfd4CA!lg z7HIVL>2VB3^B5J?cA^mrI^S%;)Z-aI*-=Uq4(aQLy$ zwT8X;Cx0Ltsf#}JGA!Shs5Hew@D7yRFnn+c4%RL04#Wd5B_Drx`1l^)LD|#OioO%FKJWTjK05XC?iQZbSS$uMD2-^F4%Y3# z-)G|HnKx$p0RO%Uhq=?#Ro#$qtW#Sg6aC9eUJDt!ghH;r*iI%*ymOsBPkbDRsiw6X znZ%1B`JZnozMY!A#%Db_l1A{Q81BOF*Zw=sQ`2sACVtdf?DTNsskkS3eBlX~5b;=i zk8l2($de@m|3B-2!aeOH#%7uRi;O=%SjSCY5rlz8!EX!l$~}~amm7@q&6K)n_T?VD zI9_J(;`Aml$j}9*{rGh~ZT^bhUx4)``V{9M`y*dEMxtP|npv&08d(*u-r>>^^y(-z zNX$N-O^1)iGM<0^V%N@t&8YF(|xNljJj{FbYWOZ5igAI*6ztFkc&T*U2y7n%c~jIsPc zz=?y?HVgDte2OD3^Yf8!Ls7eD<*8!_3n7y5ELiTY+Ui^(4&C6`EzA#6(*1{S2Ny54za%cvW-;`jvXnIBiM(88dOHjvFQ;Xcd9kT9;F*2tI@Ebv>rT?=hN29U0v9t8U(lvv zLU((J8I5mn6%S*D9S8} z7A`h7^XmuueOr~4@)S>+4IGy?5-jttjNGWR8c94}I>VXaw-?ZASl4Rmi*wz5AaNgc z3DzlmKSC9rG!!f|N@94$&+$c_`OnckBi0HbqKM^FrMesxS~iQ_BA6tTNAdr!x;{ zIqwC({*Cxh674Qy=nDgBglwz{woo?KUBNyR8$ok7DU)Sf2SUmE;>fF) zYWxzQ@S=tdcm$9O&o;#84dA~J@w@Qt2FMM>!ar+B6!whSerh5_bII%gwK89gEbGy( zNTox7p!WUL!n8jrL2vl2C{z@%mLjl9`EF&|-4*_W9d@08ldp%C30RbtWZK8Ce;rL1 zzV-zdO25W|x{99Y;q4>bzpoLN7^B?~XlgLKr$^znB8RIjd5fgtzVX^q&oLZT<3M8JZx zi^Fa1 ze5)7O9*1=EeJB!HEfF;sVkpp6eAj{OHo{EeRL1~2Vlp>Qx zh)Kd-BvFHGR(>YY9h$&*DO99Se6-?ZymF2J+Xc2SSm9h6Wr$DLAfNKvv{rP9@LUT? zR3&5dMGu?-=iv;w)1d4!GbIVH-FtOnjYlU*uyZycq&#zs_$1p6eD6@r^B4IcFNVlSi@{=GbCj~^Ek*c` zgyp%YA4~gALTgUxEqYmRPu8zNIjSQLg}qIAQ&z43Ih_RPXL`)*V&J_2yyZAvzJ z-_y7TAM59KjTV&M4ein#7-(}Lu3nR@jqDP)`dOKM;q^VYo}F!KY9GHYOuzkh-_7Pm zi;-PL((0FQia*Uqdnzj`am|88yP#W`L$>96sxM7E0=k{qm#in|E-!{SKA9`cm)7%W z5S&@fu4DVJ;Q?QkyvzLTAt8fUT&LKKpH%>Y&wT`Ojyjl zOa>V2D-Kh?Xxd)q`pdk!3a^*#(J?5|e^WMdD}7ocD=)BU^4Hh?kYB6D1qNz8-`nme zHRVcXn0lxXH2f(lelOsL8vY?dv_bxmGtNqzftWf=M)Tp98zM#hIWC;yjj);x`$07} zE9togsY-WU(|p8=eGQq~9Ys@fDqV3j}hi;)MKr?l4Jj&&ElvvR>WN;2t0o zK+-B)CqS(K+#`0upi5{qNr1VKhssr?0+kzKd-=XIZw<6~{8Z`R`(t_~<{oaw^PSD6 z+r5B1^9nAh=7PT;6%UdY;W}}4o))=|B|0s>vlHtVW;}6Ot$7u>W1e4nF3{asZ<)WMqG@XP0%uH$EiV=s{;9c+W2Z@ww5&$kXFu1|Ek zZZ&vlhu!^B8FP2TWD7V?Zzk2PrQp8EQI2>0AZlu-7Q+cgDV&K~Y6%!d_n(}AqXGu;`LS5GSSvA>rR$m(sm79XBRAZN( zxgm>(_yT%6Y2+m)0I6W9W*zv&?d6||?zZ+X0amL0WmsA9x?kwoq8is5MRR?yj{gcxt*A8m84r{mzp&bN^Cl{rgH{w5iOU2~*GC>c|B8 ze7mKe3GA70Q?Z%V-&l;OD+V&IU8I?IJLvjD=Kg-WXa7RQ9QU=s9XIeT)BD{`iL*XY zH}hC`?@v|56YDdHF=F=6&i=k{eRnt0vS4Dcy%jxs8&p@xx)mqf@e=o!^~OtL^G?f# zSpRCltH_zH?}KB}{Vumvo7}8gi923tr%2^Q%!h&+D_URvavhnq*}?Lq5w4F{TB16ds?+ zm;AY}X$+x)up?RY4ku3SZ+!vj{zD{9P_>&=G2A@y%jW3!XQc-+{?zm-pyZmy)%XLn zTs^|Sw{+|)d7(kb(YFp1mgp}aT3UBK`WZ0%$@k;P%PftiB=vxzcT+*c(k|a^hW~bk zKIL2EZ)>}wdc_VoF9b+IIS%y`3n{=O*^8gbviI`YVN>tv3U`XdbWbiGEzC-3_BbI- z0n~pWzid{jFKXA8l75@6RSDB+eqc$J!}(DEmF}WKk$!bi5HW2)5Mpt1V4?+i&K*hERs2r)(B+*6@Ti!ftf930(;#_B!f7JC2S@f zOW#ajXAmu=L$zfQ9}fFHG4&~u(n471(?8!GyT#~EiP`raxTR?~eFKcX(%`~P^;*pc zGXTcmH)dEd5g&eRIjZm@iH{ze3U&(E3CZ8sUqMB?;F}J}y#ca*ahlDQn}#oC;+vF= z+GfN{5OAV-lxoe3W1H+|kZw#jkB&PjXio`yD`K`DVw}S<2j3ze%CY!BBN2;e>&{;q zA|h#%H8;{lA~c%f$}kp{)Aad_Jx||A>x2oJRik!$9Wn*3pUD0P?V?O2SMG_VJDxWsqj-}M|D%`w_kwo3ao#$2 zNmM(-dYJ;Sv&9vj=JJKRVhn}PKU==XcqMquzN0eHjaTAn1|K-u?$WQJTwp?9h-SU9 zXqT{w1~J;l{rkKnZem=@Mz#_Xg<|u=0bloqPr!?D&)3OQ(S--Io=f%Y0;n=?TOtZY zT&&0P`fp2~N%J)LB#nr#JwQ1}GtJuVP($a#L#*txJ{1ccb2wUYTuSabpUGp0?7n_F z`=*-pI`PJl(Je!j>+h18BCD zTUCQT5D3`*7T`TR4$_lzh+sXgc3(&@3|O$Vwpy2qSmXKcRsH|lWchzT78nZc%=0X= T(qM%`8o*49E*NH>a-#ePJ<. -# -import numpy as np - -import espressomd -espressomd.assert_features(["ELECTROSTATICS", "LENNARD_JONES"]) -from espressomd import code_info -from espressomd import analyze -from espressomd import integrate -from espressomd.interactions import * -from espressomd import reaction_ensemble - - -# System parameters -# -box_l = 60 - -# Integration parameters -# -system = espressomd.System(box_l=[box_l] * 3) -system.set_random_state_PRNG() -# system.seed = system.cell_system.get_state()['n_nodes'] * [1234] -np.random.seed(seed=12) - -system.time_step = 0.02 -system.cell_system.skin = 0.4 -system.cell_system.max_num_cells = 2744 - - -# -# Setup System # -# -mode = "reaction_ensemble" -# mode="constant_pH_ensemble" - -# Particle setup -# -# type 0 = HA -# type 1 = A- -# type 2 = H+ - -N0 = 50 # number of titratable units -K_diss = 0.000830 - -for i in range(N0): - system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=1) -for i in range(N0, 2 * N0): - system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=2) - -RE = None -if(mode == "reaction_ensemble"): - RE = reaction_ensemble.ReactionEnsemble( - temperature=1, - exclusion_radius=1, - seed=2) -elif(mode == "constant_pH_ensemble"): - RE = reaction_ensemble.ConstantpHEnsemble( - temperature=1, exclusion_radius=1, seed=2) - RE.constant_pH = 3 -RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], - product_types=[1, 2], product_coefficients=[1, 1], - default_charges={0: 0, 1: -1, 2: +1}) -print(RE.get_status()) -system.setup_type_map([0, 1, 2]) - - -alpha = [] - - -for i in range(100): - RE.reaction(100) - print("HA", system.number_of_particles(type=0), "A-", - system.number_of_particles(type=1), "H+", - system.number_of_particles(type=2)) - alpha.append(system.number_of_particles(type=1) / N0) - -alpha_av = np.mean(alpha) -alpha_err = np.std(alpha) / np.sqrt(len(alpha)) - -print("\n = {} (err = {})".format(alpha_av, alpha_err)) diff --git a/doc/tutorials/10-reaction_ensemble/scripts/RE_vs_cpH_poly.py b/doc/tutorials/10-reaction_ensemble/scripts/RE_vs_cpH_poly.py deleted file mode 100644 index a80b04d91d4..00000000000 --- a/doc/tutorials/10-reaction_ensemble/scripts/RE_vs_cpH_poly.py +++ /dev/null @@ -1,229 +0,0 @@ -# -# Copyright (C) 2013-2019 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 numpy as np - -import espressomd -espressomd.assert_features(["ELECTROSTATICS", "LENNARD_JONES"]) -from espressomd import code_info -from espressomd import analyze -from espressomd import integrate -from espressomd.interactions import * -from espressomd import reaction_ensemble -from espressomd import polymer -from espressomd import interactions -from espressomd import electrostatics -import sys - -# System parameters -# -box_l = 56.134 -l_bjerrum = 2.0 -temperature = 1. - -system = espressomd.System(box_l=[box_l] * 3) - -# Integration parameters -# - -system.set_random_state_PRNG() -np.random.seed(seed=system.seed) -system.time_step = 0.01 -system.cell_system.skin = 10. # only for tutorial purposes -system.cell_system.max_num_cells = 2744 - -system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42) - - -# -# Setup System # -# - -# reaction method -mode = "reaction_ensemble" -# mode="constant_pH_ensemble" - - -# Particle setup -# -N_P = 1 # number of chains -MPC = 50 # monomers per chain -N0 = N_P * MPC # total number of monomers -nNaOH = 0 # number of initial Na+OH- -nHCl = 0 # number of initial H+Cl- (additional H+'s) - - -type_HA = 0 # type 0 = HA -type_A = 1 # type 1 = A- -type_H = 2 # type 2 = H+ -type_OH = 3 # type 3 = OH- -type_Na = 4 # type 4 = Na+ -type_Cl = 5 # type 5 = Cl- - -charges = {} -charges[type_HA] = 0 -charges[type_A] = -1 -charges[type_H] = 1 -charges[type_OH] = -1 -charges[type_Na] = 1 -charges[type_Cl] = -1 - - -# bonding interaction parameter -bond_l = 1.2 # bond length -kbond = 100 # force constant for harmonic bond -harmonic_bond = interactions.HarmonicBond(k=kbond, r_0=bond_l) -system.bonded_inter.add(harmonic_bond) - -# non-bonding interactions (LJ) -lj_eps = 1.0 -lj_sig = 1.0 -lj_cut = 1.12246 -lj_shift = 0.0 - - -# setting up the polymer -positions = polymer.positions( - n_polymers=N_P, beads_per_chain=MPC, bond_length=bond_l, seed=13) -for polymer in positions: - for i, pos in enumerate(polymer): - id = len(system.part) - system.part.add(id=id, pos=pos, type=type_A, q=charges[type_A]) - if i > 0: - system.part[id].add_bond((harmonic_bond, id - 1)) -# setting up counterions -for i in range(N0): - system.part.add(pos=np.random.random(3) * - system.box_l, type=type_H, q=charges[type_H]) - -# setting up other ions -# - Na+ and OH- -for i in range(nNaOH): - system.part.add(pos=np.random.random(3) * - system.box_l, type=type_OH, q=charges[type_OH]) -for i in range(nNaOH): - system.part.add(pos=np.random.random(3) * - system.box_l, type=type_Na, q=charges[type_Na]) -# - (additional) H+ and Cl- -for i in range(nHCl): - system.part.add(pos=np.random.random(3) * - system.box_l, type=type_H, q=charges[type_H]) -for i in range(nHCl): - system.part.add(pos=np.random.random(3) * - system.box_l, type=type_Cl, q=charges[type_Cl]) - - -# setting up LJ-interactions -for i in range(1, 5): - for j in range(i + 1, 6): - system.non_bonded_inter[i, j].lennard_jones.set_params( - epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift=lj_shift) - - -# Setting up electrostatics -# -p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=1e-3) -system.actors.add(p3m) - - -K_diss = 0.002694 -K_w = 10.**(-14) * 0.02694**2 -RE = None -if(mode == "reaction_ensemble"): - RE = reaction_ensemble.ReactionEnsemble( - temperature=temperature, exclusion_radius=1, seed=3) -elif(mode == "constant_pH_ensemble"): - RE = reaction_ensemble.ConstantpHEnsemble( - temperature=temperature, exclusion_radius=1, seed=3) - RE.constant_pH = 0 - -# HA <--> A- + H+ -RE.add_reaction( - gamma=K_diss, reactant_types=[type_HA], reactant_coefficients=[1], - product_types=[type_A, type_H], product_coefficients=[1, 1], - default_charges={type_HA: charges[type_HA], type_A: charges[type_A], - type_H: charges[type_H]}) - -# H2O autoprotolysis -RE.add_reaction( - gamma=(1 / K_w), reactant_types=[type_H, type_OH], - reactant_coefficients=[1, 1], product_types=[], product_coefficients=[], - default_charges={type_H: charges[type_H], type_OH: charges[type_OH]}) - - -print(RE.get_status()) -system.setup_type_map([type_HA, type_A, type_H, type_OH, type_Na, type_Cl]) - - -alpha = [] -nHA = [] -nA = [] -nH = [] -nOH = [] -qdist = np.zeros(N0) -c = 0 -n_iterations = 500 # this is for tutorial only, too few integration steps -n_steps_production = 10000 -n_steps_thermalization = 2000 - -for i in range(n_steps_thermalization + n_steps_production): - RE.reaction() - system.integrator.run(n_iterations) - print(i, ") HA", system.number_of_particles(type=type_HA), "A-", - system.number_of_particles(type=type_A), "H+", - system.number_of_particles(type=type_H), 'OH-', - system.number_of_particles(type=type_OH), 'Cl-', - system.number_of_particles(type=type_Cl), 'NA+', - system.number_of_particles(type=type_Na)) - if (i > n_steps_thermalization): # just a bit of thermalization before starting to gain informations about the properties of the system - alpha.append(system.number_of_particles(type=type_A) / N0) - nHA.append(system.number_of_particles(type=type_HA)) - nA.append(system.number_of_particles(type=type_A)) - nH.append(system.number_of_particles(type=type_H)) - nOH.append(system.number_of_particles(type=type_OH)) - - c = c + 1 - - for n in range(N0): - qn = system.part[n].q - qdist[n] = qdist[n] + qn - print(qdist) - - -alpha_av = np.mean(alpha) -alpha_err = np.std(alpha) / np.sqrt(len(alpha)) - -nHA_av = np.mean(nHA) -nA_av = np.mean(nA) -nH_av = np.mean(nH) -nOH_av = np.mean(nOH) - - -print("\n = {} (err = {})".format(alpha_av, alpha_err)) -print("\n") -print("\n = {} ".format(nHA_av)) -print("\n = {} ".format(nA_av)) -print("\n = {} ".format(nH_av)) -print("\n = {} ".format(nOH_av)) - - -qdist = qdist / c -print('*******************************************') -print('*** charge distribution along the chain ***') -for i in range(N0): - print(i + 1, qdist[i]) diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt index 8419837e9f6..8c0af9d594e 100644 --- a/doc/tutorials/CMakeLists.txt +++ b/doc/tutorials/CMakeLists.txt @@ -76,7 +76,6 @@ add_subdirectory(05-raspberry_electrophoresis) add_subdirectory(06-active_matter) add_subdirectory(07-electrokinetics) add_subdirectory(08-visualization) -add_subdirectory(10-reaction_ensemble) add_subdirectory(11-ferrofluid) add_subdirectory(12-constant_pH) @@ -90,7 +89,6 @@ add_custom_target(tutorials DEPENDS tutorials_06 tutorials_07 tutorials_08 - tutorials_10 tutorials_11 tutorials_12 ) diff --git a/doc/tutorials/Readme.rst b/doc/tutorials/Readme.rst index 2a5f3e86715..900a67b2c8a 100644 --- a/doc/tutorials/Readme.rst +++ b/doc/tutorials/Readme.rst @@ -15,7 +15,6 @@ physical systems. Currently, the following tutorials are available: * :file:`06-active_matter`: Modelling of self-propelling particles * :file:`07-electrokinetics`: Modelling electrokinetics together with hydrodynamic interactions * :file:`08-visualization`: Using the online visualizers of ESPResSo -* :file:`10-reaction_ensemble`: Modelling chemical reactions by means of the reaction ensemble * :file:`11-ferrofluid`: Modelling a monolayer ferrofluid system * :file:`12-constant_pH`: Modelling an acid dissociation curve using the constant pH method From d2153982c65b7651dde29eabb0a77d362801ebc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 16:23:45 +0200 Subject: [PATCH 04/11] Remove unused LaTeX preambles They were listed as a dependencies of tutorial 06, but weren't actually included anywhere. --- doc/misc/common.tex | 124 ------------------ doc/tutorials/06-active_matter/CMakeLists.txt | 3 +- doc/tutorials/common/common.tex | 98 -------------- 3 files changed, 1 insertion(+), 224 deletions(-) delete mode 100644 doc/misc/common.tex delete mode 100644 doc/tutorials/common/common.tex diff --git a/doc/misc/common.tex b/doc/misc/common.tex deleted file mode 100644 index 6f3e2a35da5..00000000000 --- a/doc/misc/common.tex +++ /dev/null @@ -1,124 +0,0 @@ -% Do magic that \bfseries for keywords actually works -\usepackage{lmodern} % use the Latin Modern fonts - -% Custom colors -\definecolor{stringblue}{rgb}{0.09,0.211,0.57} -\definecolor{codered}{rgb}{0.655,0.1137,0.3747} -\definecolor{deepgreen}{rgb}{0,0.5,0} -\definecolor{commentgray}{rgb}{0.4,0.4,0.4} -\definecolor{framegray}{rgb}{0.8,0.8,0.8} -\lstdefinestyle{pypressostyle}{ - language=Python, - belowcaptionskip=1\baselineskip, - breaklines=true, - xleftmargin=\parindent, - showstringspaces=false, - stringstyle=\color{orange}, - belowskip=\bigskipamount, - aboveskip=\bigskipamount, - basicstyle=\ttfamily\small, - otherkeywords={self}, % Add keywords here - keywordstyle=\bfseries\color{codered}, - moredelim=*[s][]{.}{\ }, - moredelim=*[s][]{\ }{.}, - moredelim=*[s][]{.}{\[}, - emph={accuracy, - actors, - Actors, - Actor, - add, - add_bond, - Analysis, - analyze, - AngleCosine, - AngleCossquare, - AngleHarmonic, - bjerrum_length, - BondedInteraction, - BondedInteractionNotDefined, - BondedInteractions, - bonded_inter, - bonded_interaction_classes, - box_l, - cellsystem, - CellSystem, - COORDS_ALL_FIXED, - COORDS_FIX_MASK, - COORD_FIXED, - cuda_init, - CudaInitHandle, - cutoff, - Dihedral, - electrostatics, - epsilon, - espressomd, - features, - FeneBond, - galilei, - GalileiTransform, - harmonic, - HarmonicBond, - HarmonicDumbbellBond, - highlander, - id, - integrate, - integrator, - Integrator, - interactions, - lennard_jones, - minimize_energy, - MinimizeEnergy, - non_bonded_inter, - NonBondedInteraction, - NonBondedInteractionHandle, - NonBondedInteractions, - OifGlobalForces, - OifLocalForces, - OrderedDict, - Overlapped, - part, - particle_data, - PARTICLE_EXT_FORCE, - PARTICLE_EXT_TORQUE, - ParticleHandle, - ParticleList, - ParticleSlice, - polymer, - Polymer, - pos, - P3M, - q, - RigidBond, - run, - set_params, - set_steepest_descent, - shift, - sigma, - steps, - System, - Tabulated, - ThereCanOnlyBeOne, - thermostat, - Thermostat, - time_step, - update_wrapper, - utils, - v, - Virtual}, % Custom highlighting - emphstyle=\bfseries\color{deepgreen}, % Custom highlighting style - commentstyle=\color{commentgray}\ttfamily, - stringstyle=\color{stringblue}, - frame=tb, % Any extra options here - rulecolor=\color{framegray}, - showstringspaces=false % -} - -% Python environment -\lstnewenvironment{pypresso}{\lstset{style=pypressostyle}}{} - -% Python for external files -\newcommand\pypressoexternal[2][]{{ -\lstset{style=pypressostyle}\lstinputlisting[#1]{#2}}} - -% Python for inline -\newcommand\pypressoinline[1]{{\lstset{style=pypressostyle}\lstinline!#1!}} diff --git a/doc/tutorials/06-active_matter/CMakeLists.txt b/doc/tutorials/06-active_matter/CMakeLists.txt index ce94837c8fb..1f24e54a82d 100644 --- a/doc/tutorials/06-active_matter/CMakeLists.txt +++ b/doc/tutorials/06-active_matter/CMakeLists.txt @@ -6,11 +6,10 @@ get_filename_component(BASENAME ${CMAKE_CURRENT_SOURCE_DIR} NAME) add_custom_command( OUTPUT ${BASENAME}.pdf COMMAND sh ../../latexit.sh - ${CMAKE_CURRENT_SOURCE_DIR}:${CMAKE_CURRENT_SOURCE_DIR}/../common + ${CMAKE_CURRENT_SOURCE_DIR} ${BASENAME} DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${BASENAME}.tex ${CMAKE_CURRENT_SOURCE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR}/../common ) add_custom_target(tutorials_06 DEPENDS ${BASENAME}.pdf) diff --git a/doc/tutorials/common/common.tex b/doc/tutorials/common/common.tex deleted file mode 100644 index 852257252a0..00000000000 --- a/doc/tutorials/common/common.tex +++ /dev/null @@ -1,98 +0,0 @@ -% Copyright (C) 2010-2019 The ESPResSo project -% Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 -% Max-Planck-Institute for Polymer Research, Theory Group -% -% 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 . -% -\usepackage[draft]{varioref} % defines \vref -\usepackage{hyperref} % automatically creates links when - % using pdflatex, defines \url -\usepackage{ifpdf} % defines \ifpdf -\usepackage{graphicx} % handles graphics -\usepackage{color} % use colors - -\usepackage{amsmath} - -\usepackage{verbatim} % required for \verbatim and \endverbatim -\usepackage{fancyvrb} -\usepackage{calc} % compute length -\usepackage{ifthen} % provide ifthen -\usepackage{xspace} -\usepackage{units} -\usepackage[numbers]{natbib} - -\usepackage{listings} - -% For building the distribution docs, disable todo boxes. -%\usepackage[disable]{todonotes} -\usepackage{todonotes} - -\newcommand{\es}{\mbox{\textsf{ESPResSo}}\xspace} -\newcommand{\ie}{\textit{i.e.}\xspace} -\newcommand{\eg}{\textit{e.g.}\xspace} -\newcommand{\etal}{\textit{et al.}\xspace} - -\newcommand{\codebox}[1]% -{\texttt{#1}} - -\DefineVerbatimEnvironment{code}{Verbatim}% -{commandchars=\\\{\}} -\makeatletter -\newenvironment{tclcode} -{% - \addtolength{\linewidth}{-2em}% set the line length - \@minipagetrue%%%DPC%%% - \@tempswatrue%%%DPC%%% - \hsize=\linewidth% - \setbox0=\vbox\bgroup\verbatim -}{\endverbatim - \unskip\setbox0=\lastbox%%%DPC%%% - \egroup - \par% - \noindent\hspace{1em}% - \codebox{\box0}% - \par\noindent% -} -\makeatother - -% \newcommand{\todo}[1]{ -% \marginpar{% -% \setlength{\fboxrule}{1pt} -% \fcolorbox{red}{yellow}{% -% \parbox{\marginparwidth-2\fboxrule-2\fboxsep}{% -% \bf\raggedright\scriptsize #1% -% }% -% }% -% }% -% } - -\makeatletter -\renewcommand{\minisec}[1]{\@afterindentfalse \vskip 1.5ex - {\parindent \z@ - \raggedsection\normalfont\sffamily\itshape\nobreak#1\par\nobreak}% - \@afterheading} -\makeatother - -\newcommand{\esptitlehead}{ - \titlehead{ - \begin{center} - \includegraphics[width=5cm]{logo/logo.pdf} - \end{center} - } -} - -\input{misc/common} - From eabc23db466132facdd6fdd1f773aabdd9c081f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 17:38:07 +0200 Subject: [PATCH 05/11] Remove tests --- testsuite/scripts/tutorials/CMakeLists.txt | 2 -- ...0-reaction_ensemble__scripts__RE_vs_cpH.py | 31 ------------------ ...ction_ensemble__scripts__RE_vs_cpH_poly.py | 32 ------------------- 3 files changed, 65 deletions(-) delete mode 100644 testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH.py delete mode 100644 testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH_poly.py diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index b4cc6d1b030..665ce37ce62 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -36,8 +36,6 @@ tutorial_test(FILE test_06-active_matter__enhanced_diffusion.py) tutorial_test(FILE test_06-active_matter__rectification_simulation.py) tutorial_test(FILE test_07-electrokinetics.py LABELS "gpu") tutorial_test(FILE test_08-visualization.py) -tutorial_test(FILE test_10-reaction_ensemble__scripts__RE_vs_cpH.py) -tutorial_test(FILE test_10-reaction_ensemble__scripts__RE_vs_cpH_poly.py) tutorial_test(FILE test_11-ferrofluid_1.py) tutorial_test(FILE test_11-ferrofluid_2.py) tutorial_test(FILE test_11-ferrofluid_3.py) diff --git a/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH.py b/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH.py deleted file mode 100644 index fd16578bce1..00000000000 --- a/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH.py +++ /dev/null @@ -1,31 +0,0 @@ -# Copyright (C) 2019 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 - -tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/10-reaction_ensemble/scripts/RE_vs_cpH.py") - - -@skipIfMissingFeatures -class Tutorial(ut.TestCase): - system = tutorial.system - - -if __name__ == "__main__": - ut.main() diff --git a/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH_poly.py b/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH_poly.py deleted file mode 100644 index d3238529eb6..00000000000 --- a/testsuite/scripts/tutorials/test_10-reaction_ensemble__scripts__RE_vs_cpH_poly.py +++ /dev/null @@ -1,32 +0,0 @@ -# Copyright (C) 2019 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 - -tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/10-reaction_ensemble/scripts/RE_vs_cpH_poly.py", - n_iterations=50, n_steps_production=500, n_steps_thermalization=100) - - -@skipIfMissingFeatures -class Tutorial(ut.TestCase): - system = tutorial.system - - -if __name__ == "__main__": - ut.main() From 7397e6ab397deda4d4e370d1a662ce3fc2c38187 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 15 Oct 2019 15:36:30 +0200 Subject: [PATCH 06/11] core: particle_data: Removed redundant Cell::resize --- src/core/Cell.hpp | 4 ---- src/core/ParticleList.hpp | 16 +++++++--------- src/core/collision.cpp | 2 +- src/core/forces_inline.hpp | 2 +- src/core/global.cpp | 2 +- .../lb_particle_coupling.cpp | 2 +- src/core/integrators/steepest_descent.cpp | 2 +- src/core/integrators/velocity_verlet_npt.cpp | 2 +- src/core/particle_data.cpp | 5 ++--- 9 files changed, 15 insertions(+), 22 deletions(-) diff --git a/src/core/Cell.hpp b/src/core/Cell.hpp index fe1c68bf15d..0d6b65871ce 100644 --- a/src/core/Cell.hpp +++ b/src/core/Cell.hpp @@ -107,10 +107,6 @@ class Cell : public ParticleList { * @brief All neighbors of the cell. */ neighbors_type &neighbors() { return m_neighbors; } - - void resize(size_t size) { - realloc_particlelist(static_cast(this), this->n = size); - } }; #endif diff --git a/src/core/ParticleList.hpp b/src/core/ParticleList.hpp index 8693249ce57..16dc8a65bee 100644 --- a/src/core/ParticleList.hpp +++ b/src/core/ParticleList.hpp @@ -21,10 +21,10 @@ struct ParticleList { Utils::Span particles() { return {part, static_cast(n)}; } - int resize(int size) { -/** granularity of the particle buffers in particles */ - constexpr int PART_INCREMENT = 8; + /** granularity of the particle buffers in particles */ + static constexpr int INCREMENT = 8; + int resize(int size) { assert(size >= 0); int old_max = max; Particle *old_start = part; @@ -35,12 +35,10 @@ struct ParticleList { max = 0; else /* shrink not as fast, just lose half, rounded up */ - max = - PART_INCREMENT * - (((max + size + 1) / 2 + PART_INCREMENT - 1) / PART_INCREMENT); + max = INCREMENT * (((max + size + 1) / 2 + INCREMENT - 1) / INCREMENT); } else /* round up */ - max = PART_INCREMENT * ((size + PART_INCREMENT - 1) / PART_INCREMENT); + max = INCREMENT * ((size + INCREMENT - 1) / INCREMENT); if (max != old_max) part = Utils::realloc(part, sizeof(Particle) * max); return part != old_start; @@ -49,9 +47,9 @@ struct ParticleList { /** Allocate storage for local particles and ghosts. This version does \em not care for the bond information to be freed if necessary. - \param plist the list on which to operate + \param l the list on which to operate \param size the size to provide at least. It is rounded - up to multiples of \ref PART_INCREMENT. + up to multiples of \ref ParticleList::INCREMENT. \return true iff particle addresses have changed */ inline int realloc_particlelist(ParticleList *l, int size) { return assert(l), l->resize(size); diff --git a/src/core/collision.cpp b/src/core/collision.cpp index 90666eaf540..9c9757b1bd3 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -27,9 +27,9 @@ #include "event.hpp" #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "particle_data.hpp" #include "rotation.hpp" #include "virtual_sites/VirtualSitesRelative.hpp" -#include "particle_data.hpp" #include #include diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 06ebfc8e297..b13633791aa 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -60,8 +60,8 @@ #include "object-in-fluid/oif_global_forces.hpp" #include "object-in-fluid/oif_local_forces.hpp" #include "object-in-fluid/out_direction.hpp" -#include "thermostat.hpp" #include "particle_data.hpp" +#include "thermostat.hpp" #ifdef DIPOLES #include "electrostatics_magnetostatics/dipole_inline.hpp" diff --git a/src/core/global.cpp b/src/core/global.cpp index 194dd7bb992..7fa8d8dbab7 100644 --- a/src/core/global.cpp +++ b/src/core/global.cpp @@ -34,10 +34,10 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" #include "object-in-fluid/oif_global_forces.hpp" +#include "particle_data.hpp" #include "rattle.hpp" #include "thermostat.hpp" #include "tuning.hpp" -#include "particle_data.hpp" #include diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index d399187e4b0..9a935f0096f 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -27,8 +27,8 @@ #include "lb_interface.hpp" #include "lb_interpolation.hpp" #include "lbgpu.hpp" -#include "random.hpp" #include "particle_data.hpp" +#include "random.hpp" #include #include diff --git a/src/core/integrators/steepest_descent.cpp b/src/core/integrators/steepest_descent.cpp index c6d1f4d2b97..d83fb616c4c 100644 --- a/src/core/integrators/steepest_descent.cpp +++ b/src/core/integrators/steepest_descent.cpp @@ -24,8 +24,8 @@ #include "communication.hpp" #include "event.hpp" #include "integrate.hpp" -#include "rotation.hpp" #include "particle_data.hpp" +#include "rotation.hpp" #include diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 1302c8be4a8..7142e31e478 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -27,8 +27,8 @@ #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" -#include "thermostat.hpp" #include "particle_data.hpp" +#include "thermostat.hpp" #include diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index dba226c32b1..17c6019746b 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -25,8 +25,8 @@ */ #include "particle_data.hpp" -#include "Particle.hpp" #include "PartCfg.hpp" +#include "Particle.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "cells.hpp" #include "communication.hpp" @@ -543,8 +543,7 @@ void realloc_local_particles(int part) { if (part >= max_local_particles) { /* round up part + 1 in granularity INCREMENT */ - max_local_particles = - INCREMENT * ((part + INCREMENT) / INCREMENT); + max_local_particles = INCREMENT * ((part + INCREMENT) / INCREMENT); local_particles = Utils::realloc(local_particles, sizeof(Particle *) * max_local_particles); From ec8d9949c0cfa4bc5e0a5a0f1408e3aa11799e62 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 15 Oct 2019 19:38:12 +0200 Subject: [PATCH 07/11] maintainer: Added missing license header --- src/core/ParticleList.hpp | 40 +++++++++++++++---- .../p3m_gpu_cuda.cu | 1 + .../scafacos.cpp | 1 + src/core/rotation.cpp | 1 + 4 files changed, 36 insertions(+), 7 deletions(-) diff --git a/src/core/ParticleList.hpp b/src/core/ParticleList.hpp index 16dc8a65bee..00d31cced34 100644 --- a/src/core/ParticleList.hpp +++ b/src/core/ParticleList.hpp @@ -1,3 +1,23 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 . + */ #ifndef ESPRESSO_CORE_PARTICLE_LIST_HPP #define ESPRESSO_CORE_PARTICLE_LIST_HPP @@ -24,7 +44,16 @@ struct ParticleList { /** granularity of the particle buffers in particles */ static constexpr int INCREMENT = 8; - int resize(int size) { + /** + * @brief Allocate storage for local particles and ghosts. + * + * This version does \em not care for the bond information to be freed if + * necessary. + * @param size the size to provide at least. It is rounded + * up to multiples of @ref ParticleList::INCREMENT. + * @return true iff particle addresses have changed + */ + int realloc(int size) { assert(size >= 0); int old_max = max; Particle *old_start = part; @@ -43,14 +72,11 @@ struct ParticleList { part = Utils::realloc(part, sizeof(Particle) * max); return part != old_start; } + + int resize(int size) { return realloc(this->n = size); } }; -/** Allocate storage for local particles and ghosts. This version - does \em not care for the bond information to be freed if necessary. - \param l the list on which to operate - \param size the size to provide at least. It is rounded - up to multiples of \ref ParticleList::INCREMENT. - \return true iff particle addresses have changed */ +/** @copydoc ParticleList::resize */ inline int realloc_particlelist(ParticleList *l, int size) { return assert(l), l->resize(size); } diff --git a/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu b/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu index 5a6c1ea8606..e83fd3dc669 100644 --- a/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu +++ b/src/core/electrostatics_magnetostatics/p3m_gpu_cuda.cu @@ -65,6 +65,7 @@ #include "EspressoSystemInterface.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" #include "global.hpp" +#include "particle_data.hpp" #include using Utils::int_pow; diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index a4e0deaba5c..be47ea402df 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -41,6 +41,7 @@ #include "global.hpp" #include "grid.hpp" #include "integrate.hpp" +#include "particle_data.hpp" #include "tuning.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp index 45fb6bcfe58..abe1bb4ffa8 100644 --- a/src/core/rotation.cpp +++ b/src/core/rotation.cpp @@ -43,6 +43,7 @@ #include "global.hpp" #include "grid_based_algorithms/lb_interface.hpp" #include "integrate.hpp" +#include "particle_data.hpp" #include "thermostat.hpp" #include From f28097b192391426bf399c407b70b02355b9a544 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 16 Oct 2019 16:32:40 +0200 Subject: [PATCH 08/11] core: Replaced some uses of realloc_particlelist by new ParticleList::clear --- src/core/ParticleList.hpp | 45 ++++++++++++++----------- src/core/cells.cpp | 2 +- src/core/domain_decomposition.cpp | 8 ++--- src/core/ghosts.cpp | 12 ++++--- src/core/layered.cpp | 4 +-- src/core/nsquare.cpp | 2 +- src/core/particle_data.cpp | 16 ++++----- src/core/serialization/ParticleList.hpp | 2 +- src/core/unit_tests/link_cell_test.cpp | 8 ++--- src/core/unit_tests/verlet_ia_test.cpp | 8 ++--- 10 files changed, 53 insertions(+), 54 deletions(-) diff --git a/src/core/ParticleList.hpp b/src/core/ParticleList.hpp index 00d31cced34..bf55302576e 100644 --- a/src/core/ParticleList.hpp +++ b/src/core/ParticleList.hpp @@ -28,7 +28,6 @@ /** List of particles. The particle array is resized using a sophisticated * (we hope) algorithm to avoid unnecessary resizes. - * Access using \ref realloc_particlelist, ... */ struct ParticleList { ParticleList() : part{nullptr}, n{0}, max{0} {} @@ -36,23 +35,11 @@ struct ParticleList { Particle *part; /** Number of particles contained */ int n; + +private: /** Number of particles that fit in until a resize is needed */ int max; - Utils::Span particles() { return {part, static_cast(n)}; } - - /** granularity of the particle buffers in particles */ - static constexpr int INCREMENT = 8; - - /** - * @brief Allocate storage for local particles and ghosts. - * - * This version does \em not care for the bond information to be freed if - * necessary. - * @param size the size to provide at least. It is rounded - * up to multiples of @ref ParticleList::INCREMENT. - * @return true iff particle addresses have changed - */ int realloc(int size) { assert(size >= 0); int old_max = max; @@ -73,12 +60,30 @@ struct ParticleList { return part != old_start; } +public: + /** Current allocation size. */ + auto capacity() const { return max; } + + Utils::Span particles() { return {part, static_cast(n)}; } + + /** granularity of the particle buffers in particles */ + static constexpr int INCREMENT = 8; + + /** + * @brief Resize storage for local particles and ghosts. + * + * This version does \em not care for the bond information to be freed if + * necessary. + * @param size the size to provide at least. It is rounded + * up to multiples of @ref ParticleList::INCREMENT. + * @return true iff particle addresses have changed + */ int resize(int size) { return realloc(this->n = size); } -}; -/** @copydoc ParticleList::resize */ -inline int realloc_particlelist(ParticleList *l, int size) { - return assert(l), l->resize(size); -} + /** + * @brief Resize the List to zero. + */ + void clear() { resize(0); } +}; #endif diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 4fb7e100003..712f80e5ffb 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -410,7 +410,7 @@ void cells_resort_particles(int global_flag) { resort_particles = Cells::RESORT_NONE; rebuild_verletlist = true; - realloc_particlelist(&displaced_parts, 0); + displaced_parts.clear(); on_resort_particles(local_cells.particles()); } diff --git a/src/core/domain_decomposition.cpp b/src/core/domain_decomposition.cpp index f8ce2a35a5b..6409a6e3b31 100644 --- a/src/core/domain_decomposition.cpp +++ b/src/core/domain_decomposition.cpp @@ -700,7 +700,7 @@ void move_if_local(ParticleList &src, ParticleList &rest) { } } - realloc_particlelist(&src, src.n = 0); + src.resize(0); } /** @@ -761,7 +761,7 @@ void exchange_neighbors(ParticleList *pl, const Utils::Vector3i &grid) { Utils::Mpi::sendrecv(comm_cart, node_neighbors[2 * dir], 0, send_buf, node_neighbors[2 * dir], 0, recv_buf); - realloc_particlelist(&send_buf, 0); + send_buf.clear(); move_if_local(recv_buf, *pl); } else { @@ -784,8 +784,8 @@ void exchange_neighbors(ParticleList *pl, const Utils::Vector3i &grid) { move_if_local(recv_buf_l, *pl); move_if_local(recv_buf_r, *pl); - realloc_particlelist(&send_buf_l, 0); - realloc_particlelist(&send_buf_r, 0); + send_buf_l.clear(); + send_buf_r.clear(); } } } diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 93c44266cbb..f11f1e2ca68 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -193,11 +193,12 @@ void prepare_send_buffer(GhostCommunication *gc, int data_parts) { static void prepare_ghost_cell(Cell *cell, int size) { using Utils::make_span; - auto const old_cap = cell->max; + auto const old_cap = cell->capacity(); /* reset excess particles */ - if (size < cell->max) { - for (auto &p : make_span(cell->part + size, cell->max - size)) { + if (size < cell->capacity()) { + for (auto &p : + make_span(cell->part + size, cell->capacity() - size)) { p = Particle{}; p.l.ghost = true; } @@ -207,8 +208,9 @@ static void prepare_ghost_cell(Cell *cell, int size) { cell->resize(size); /* initialize new particles */ - if (old_cap < cell->max) { - auto new_parts = make_span(cell->part + old_cap, cell->max - old_cap); + if (old_cap < cell->capacity()) { + auto new_parts = + make_span(cell->part + old_cap, cell->capacity() - old_cap); std::uninitialized_fill(new_parts.begin(), new_parts.end(), Particle{}); for (auto &p : new_parts) { p.l.ghost = true; diff --git a/src/core/layered.cpp b/src/core/layered.cpp index fbea03314ef..f6775b8b003 100644 --- a/src/core/layered.cpp +++ b/src/core/layered.cpp @@ -482,6 +482,6 @@ void layered_exchange_and_sort_particles(int global_flag, } } - realloc_particlelist(&recv_buf_up, 0); - realloc_particlelist(&recv_buf_dn, 0); + recv_buf_up.clear(); + recv_buf_dn.clear(); } diff --git a/src/core/nsquare.cpp b/src/core/nsquare.cpp index e798db37b5e..08db52c0b09 100644 --- a/src/core/nsquare.cpp +++ b/src/core/nsquare.cpp @@ -168,7 +168,7 @@ void nsq_exchange_particles(int global_flag, ParticleList *displaced_parts) { auto const target_node = (p.identity() % n_nodes); send_buf.at(target_node).emplace_back(std::move(p)); } - realloc_particlelist(displaced_parts, displaced_parts->n = 0); + displaced_parts->resize(0); /* Exchange particles */ std::vector> recv_buf(n_nodes); diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 17c6019746b..b2b25d8a1a4 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -561,12 +561,12 @@ void update_local_particles(ParticleList *pl) { } void append_unindexed_particle(ParticleList *l, Particle &&part) { - realloc_particlelist(l, ++l->n); + l->resize(l->n + 1); new (&(l->part[l->n - 1])) Particle(std::move(part)); } Particle *append_indexed_particle(ParticleList *l, Particle &&part) { - auto const re = realloc_particlelist(l, ++l->n); + auto const re = l->resize(l->n + 1); auto p = new (&(l->part[l->n - 1])) Particle(std::move(part)); assert(p->p.identity <= max_seen_particle); @@ -582,7 +582,7 @@ Particle *move_unindexed_particle(ParticleList *dl, ParticleList *sl, int i) { assert(sl->n > 0); assert(i < sl->n); - realloc_particlelist(dl, ++dl->n); + dl->resize(dl->n + 1); auto dst = &dl->part[dl->n - 1]; auto src = &sl->part[i]; auto end = &sl->part[sl->n - 1]; @@ -592,14 +592,14 @@ Particle *move_unindexed_particle(ParticleList *dl, ParticleList *sl, int i) { new (src) Particle(std::move(*end)); } - realloc_particlelist(sl, --sl->n); + sl->resize(sl->n - 1); return dst; } Particle *move_indexed_particle(ParticleList *dl, ParticleList *sl, int i) { assert(sl->n > 0); assert(i < sl->n); - int re = realloc_particlelist(dl, ++dl->n); + int re = dl->resize(dl->n + 1); Particle *dst = &dl->part[dl->n - 1]; Particle *src = &sl->part[i]; Particle *end = &sl->part[sl->n - 1]; @@ -617,7 +617,7 @@ Particle *move_indexed_particle(ParticleList *dl, ParticleList *sl, int i) { new (src) Particle(std::move(*end)); } - if (realloc_particlelist(sl, --sl->n)) { + if (sl->resize(sl->n - 1)) { update_local_particles(sl); } else if (src != end) { local_particles[src->p.identity] = src; @@ -652,7 +652,7 @@ Particle extract_indexed_particle(ParticleList *sl, int i) { new (src) Particle(std::move(*end)); } - if (realloc_particlelist(sl, --sl->n)) { + if (sl->resize(sl->n - 1)) { update_local_particles(sl); } else if (src != end) { local_particles[src->p.identity] = src; @@ -1325,7 +1325,7 @@ void send_particles(ParticleList *particles, int node) { free_particle(&particles->part[pc]); } - realloc_particlelist(particles, particles->n = 0); + particles->clear(); } void recv_particles(ParticleList *particles, int node) { diff --git a/src/core/serialization/ParticleList.hpp b/src/core/serialization/ParticleList.hpp index 938784c897d..de15cb668dd 100644 --- a/src/core/serialization/ParticleList.hpp +++ b/src/core/serialization/ParticleList.hpp @@ -28,7 +28,7 @@ void load(Archive &ar, ParticleList &pl, const unsigned int /* version */) { int size; ar >> size; - realloc_particlelist(&pl, pl.n = size); + pl.resize(size); for (int i = 0; i < size; i++) { ar >> pl.part[i]; } diff --git a/src/core/unit_tests/link_cell_test.cpp b/src/core/unit_tests/link_cell_test.cpp index d732730c402..0c26dd08b5e 100644 --- a/src/core/unit_tests/link_cell_test.cpp +++ b/src/core/unit_tests/link_cell_test.cpp @@ -46,8 +46,8 @@ BOOST_AUTO_TEST_CASE(link_cell) { c.m_neighbors = Neighbors(neighbors, {}); - c.part = new Particle[n_part_per_cell]; - c.n = c.max = n_part_per_cell; + c.resize(n_part_per_cell); + std::uninitialized_fill(c.part, c.part + c.n, Particle()); for (unsigned i = 0; i < n_part_per_cell; ++i) { c.part[i].p.identity = id++; @@ -86,8 +86,4 @@ BOOST_AUTO_TEST_CASE(link_cell) { BOOST_CHECK((it->first == i) && (it->second == j)); ++it; } - - for (auto &c : cells) { - delete[] c.part; - } } diff --git a/src/core/unit_tests/verlet_ia_test.cpp b/src/core/unit_tests/verlet_ia_test.cpp index bf7117c6b23..d0cdadecc42 100644 --- a/src/core/unit_tests/verlet_ia_test.cpp +++ b/src/core/unit_tests/verlet_ia_test.cpp @@ -70,8 +70,8 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { c.m_neighbors = Neighbors(neighbors, {}); - c.part = new Particle[n_part_per_cell]; - c.n = c.max = n_part_per_cell; + c.resize(n_part_per_cell); + std::uninitialized_fill(c.part, c.part + c.n, Particle()); for (unsigned i = 0; i < n_part_per_cell; ++i) { c.part[i].p.identity = id++; @@ -146,8 +146,4 @@ BOOST_AUTO_TEST_CASE(verlet_ia) { [](int count) { return count == 1; })); check_pairs(n_part, pairs); - - for (auto &c : cells) { - delete[] c.part; - } } From e70998e647880d6d911ecceaa292f806541e1a9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 20:06:11 +0200 Subject: [PATCH 09/11] Remove pdbreader --- doc/sphinx/io.rst | 9 +- src/core/io/CMakeLists.txt | 1 - src/core/io/reader/CMakeLists.txt | 4 - src/core/io/reader/readpdb.cpp | 193 ---------------------------- src/core/io/reader/readpdb.hpp | 55 -------- src/script_interface/CMakeLists.txt | 2 +- 6 files changed, 2 insertions(+), 262 deletions(-) delete mode 100644 src/core/io/reader/CMakeLists.txt delete mode 100644 src/core/io/reader/readpdb.cpp delete mode 100644 src/core/io/reader/readpdb.hpp diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst index d2f137d9954..9706b8d2455 100644 --- a/doc/sphinx/io.rst +++ b/doc/sphinx/io.rst @@ -369,7 +369,7 @@ Writing various formats using MDAnalysis If the MDAnalysis package (https://mdanalysis.org) is installed, it is possible to use it to convert frames to any of the supported configuration/trajectory formats, including PDB, GROMACS, GROMOS, -CHARMM/NAMD, AMBER, LAMMPS, ...) +CHARMM/NAMD, AMBER, LAMMPS, ... To use MDAnalysis to write in any of these formats, one has first to prepare a stream from the |es| particle data using the class :class:`espressomd.MDA_ESP`, and then read from it @@ -397,10 +397,3 @@ using MDAnalysis. A simple example is the following: W.write_next_timestep(u.trajectory.ts) # append it to the trajectory For other examples, see :file:`/samples/MDAnalysisIntegration.py` - -.. _Parsing PDB Files: - -Parsing PDB Files ------------------ - -The feature allows the user to parse simple PDB files, a file format introduced by the protein database to encode molecular structures. Together with a topology file (here ) the structure gets interpolated to the grid. For the input you will need to prepare a PDB file with a force field to generate the topology file. Normally the PDB file extension is :file:`.pdb`, the topology file extension is :file:`.itp`. Obviously the PDB file is placed instead of and the topology file instead of . diff --git a/src/core/io/CMakeLists.txt b/src/core/io/CMakeLists.txt index cb6d8f64550..605ed7cf90c 100644 --- a/src/core/io/CMakeLists.txt +++ b/src/core/io/CMakeLists.txt @@ -1,3 +1,2 @@ add_subdirectory(writer) -add_subdirectory(reader) add_subdirectory(mpiio) diff --git a/src/core/io/reader/CMakeLists.txt b/src/core/io/reader/CMakeLists.txt deleted file mode 100644 index 6c2de627ac0..00000000000 --- a/src/core/io/reader/CMakeLists.txt +++ /dev/null @@ -1,4 +0,0 @@ -add_library(pdbreader SHARED readpdb.cpp) -target_link_libraries(pdbreader PUBLIC EspressoConfig EspressoCore pdbparser) -target_include_directories(pdbreader PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) -install(TARGETS pdbreader LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd) diff --git a/src/core/io/reader/readpdb.cpp b/src/core/io/reader/readpdb.cpp deleted file mode 100644 index f50e7a7fa16..00000000000 --- a/src/core/io/reader/readpdb.cpp +++ /dev/null @@ -1,193 +0,0 @@ -/* - * Copyright (C) 2010-2019 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 . - */ -#include "readpdb.hpp" -#include "grid.hpp" -#include "nonbonded_interactions/lj.hpp" -#include "particle_data.hpp" - -#include -#include -#include - -namespace Reader { -namespace PDB { - -#ifdef LENNARD_JONES -/* Add user requested Lennard-Jones interactions */ -static void add_lj_interaction( - std::set &types, - std::vector interactions, const double rel_cutoff) { - for (std::vector::const_iterator it = interactions.begin(); - it != interactions.end(); ++it) { - for (auto type : types) { - const double epsilon_ij = sqrt(it->epsilon * type.epsilon); - const double sigma_ij = 0.5 * (it->sigma + 10. * type.sigma); - const double cutoff_ij = rel_cutoff * sigma_ij; - const double shift_ij = - -(pow(sigma_ij / cutoff_ij, 12) - pow(sigma_ij / cutoff_ij, 6)); - if ((epsilon_ij <= 0) || (sigma_ij <= 0)) { - continue; - } - lennard_jones_set_params(it->other_type, type.espresso_id, epsilon_ij, - sigma_ij, cutoff_ij, shift_ij, 0.0, 0.0); - } - } -} - -/* Add Lennard-Jones interactions between particles added from pdb/itp file */ -static void add_lj_internal( - std::set &types, - const double rel_cutoff, bool only_diagonal) { - for (auto it = types.begin(); it != types.end(); ++it) { - for (auto type : types) { - if (it->espresso_id > type.espresso_id) - continue; - if (only_diagonal && (it->espresso_id != type.espresso_id)) - continue; - const double epsilon_ij = sqrtf(it->epsilon * type.epsilon); - const double sigma_ij = 0.5 * (10. * it->sigma + 10. * type.sigma); - const double cutoff_ij = rel_cutoff * sigma_ij; - const double shift_ij = - -pow(sigma_ij / cutoff_ij, 12) - pow(sigma_ij / cutoff_ij, 6); - if ((epsilon_ij <= 0) || (sigma_ij <= 0)) { - continue; - } - lennard_jones_set_params(it->espresso_id, type.espresso_id, epsilon_ij, - sigma_ij, cutoff_ij, shift_ij, 0.0, 0.0); - } - } -} -#endif /* LENNARD_JONES */ - -static int -add_particles(PdbParser::PdbParser &parser, int first_id, int default_type, - std::set - &seen_types, - int first_type = 0, bool fit = false) { - double pos[3]; - int id = first_id; - int stat; - int type; -#ifdef ELECTROSTATICS - double q; -#endif - PdbParser::BoundingBox bb; - bb.llx = bb.lly = bb.llz = 0.0; - double bb_l[3] = {box_geo.length()[0], box_geo.length()[1], - box_geo.length()[2]}; - - bb = parser.calc_bounding_box(); - - if (fit) { - bb_l[0] = (bb.urx - bb.llx); - bb_l[1] = (bb.ury - bb.lly); - bb_l[2] = (bb.urz - bb.llz); - - for (int i = 0; i < 3; i++) { - if (bb_l[i] > box_geo.length()[i]) { - rescale_boxl(i, bb_l[i]); - } - } - } - - int last_type = first_type; - for (std::vector::const_iterator it = - parser.pdb_atoms.begin(); - it != parser.pdb_atoms.end(); ++it) { - pos[0] = (it->x - bb.llx); - pos[1] = (it->y - bb.lly); - pos[2] = (it->z - bb.llz); - - stat = place_particle(id, pos); - - const std::map::const_iterator entry = - parser.itp_atoms.find(it->i); - switch (stat) { - case ES_PART_OK: - std::cerr << "Warning: position and type of particle " << id - << " was overwritten by value from pdb file." << std::endl; - /* Fall through */ - case ES_PART_CREATED: - /* See if we have a type from itp file, otherwise set default type */ - if (entry != parser.itp_atoms.end()) { - PdbParser::itp_atomtype itp_atomtype = - parser.itp_atomtypes[entry->second.type]; - /* See if we have seen that type before */ - auto type_iterator = seen_types.find(itp_atomtype); - if (type_iterator == seen_types.end()) { - itp_atomtype.espresso_id = last_type++; - type_iterator = seen_types.insert(itp_atomtype).first; - } - itp_atomtype = *type_iterator; -#ifdef ELECTROSTATICS - q = entry->second.charge; -#endif - type = itp_atomtype.espresso_id; - } else { - type = default_type; -#ifdef ELECTROSTATICS - q = 0; -#endif - } - set_particle_type(id, type); -#ifdef ELECTROSTATICS - set_particle_q(id, q); -#endif - id++; - break; - case ES_PART_ERROR: - std::cerr << "Warning: Illegal particle id " << id << std::endl; - return id - first_id; - break; - } - } - return id - first_id; -} - -int pdb_add_particles_from_file(char *pdb_file, int first_id, int type, - std::vector &ljInteractions, - double lj_rel_cutoff, char *itp_file, - int first_type, bool fit, bool lj_internal, - bool lj_diagonal) { - int n_part; - PdbParser::PdbParser parser; - if (!parser.parse_pdb_file(pdb_file)) - return 0; - - if (itp_file) { - if (!parser.parse_itp_file(itp_file)) - return 0; - } - - /* Unique set of types that actually have particles */ - std::set seen_types; - - n_part = add_particles(parser, first_id, type, seen_types, first_type, fit); - -#ifdef LENNARD_JONES - add_lj_interaction(seen_types, ljInteractions, lj_rel_cutoff); - if (lj_internal || lj_diagonal) - add_lj_internal(seen_types, lj_rel_cutoff, lj_diagonal); -#endif - - return n_part; -} - -} // namespace PDB -} // namespace Reader diff --git a/src/core/io/reader/readpdb.hpp b/src/core/io/reader/readpdb.hpp deleted file mode 100644 index 567e2a8cc8c..00000000000 --- a/src/core/io/reader/readpdb.hpp +++ /dev/null @@ -1,55 +0,0 @@ -/* - * Copyright (C) 2014-2019 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 . - */ - -#ifndef __READ_PDB_HPP -#define __READ_PDB_HPP - -#include "Particle.hpp" -#include "PdbParser.hpp" - -namespace Reader { -namespace PDB { - -struct PdbLJInteraction { - int other_type; - double epsilon, sigma; -}; - -/** Call only on the master node: Parse pdb file and add contained particles. - @param pdb_file Filename of the pdb file. - @param first_id Id of the first particle to add. - @param type Type for the particles. - @param ljInteractions LJ interactions vector. - @param lj_rel_cutoff LJ cutoff, as a multiple of the read LJ sigma value. - @param itp_file Optional ITP file. - @param first_type Number of atom types to skip in the ITP file. - @param fit Should the box be rescaled to hold the particles. - @param lj_internal Should LJ interactions within the molecule be added. - @param lj_diagonal Just the diagonal interaction terms of @p lj_internal. - @return Number of particles that were added. - */ -int pdb_add_particles_from_file(char *pdb_file, int first_id, int type, - std::vector &ljInteractions, - double lj_rel_cutoff = 2.5, - char *itp_file = nullptr, int first_type = 0, - bool fit = false, bool lj_internal = false, - bool lj_diagonal = false); -} // namespace PDB -} // namespace Reader -#endif diff --git a/src/script_interface/CMakeLists.txt b/src/script_interface/CMakeLists.txt index 848eab8dae2..e31ce8a7ba3 100644 --- a/src/script_interface/CMakeLists.txt +++ b/src/script_interface/CMakeLists.txt @@ -33,7 +33,7 @@ endif() add_library(EspressoScriptInterface SHARED ${EspressoScriptInterface_SRC}) install(TARGETS EspressoScriptInterface LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd) -target_link_libraries(EspressoScriptInterface PRIVATE EspressoConfig PUBLIC EspressoCore Shapes core_cluster_analysis mpiio pdbreader) +target_link_libraries(EspressoScriptInterface PRIVATE EspressoConfig PUBLIC EspressoCore Shapes core_cluster_analysis mpiio) if(HDF5_FOUND) target_link_libraries(EspressoScriptInterface PRIVATE H5mdCore) endif() From cf8bfc6244ec48a615f658084cea69f53686d401 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 20:10:09 +0200 Subject: [PATCH 10/11] Remove pdbparser --- src/CMakeLists.txt | 1 - src/core/CMakeLists.txt | 1 - src/pdbparser/CMakeLists.txt | 6 - src/pdbparser/include/PdbParser.hpp | 72 -------- src/pdbparser/src/PdbParser.cpp | 188 -------------------- src/pdbparser/unit_tests/CMakeLists.txt | 3 - src/pdbparser/unit_tests/PdbParser_test.cpp | 93 ---------- 7 files changed, 364 deletions(-) delete mode 100644 src/pdbparser/CMakeLists.txt delete mode 100644 src/pdbparser/include/PdbParser.hpp delete mode 100644 src/pdbparser/src/PdbParser.cpp delete mode 100644 src/pdbparser/unit_tests/CMakeLists.txt delete mode 100644 src/pdbparser/unit_tests/PdbParser_test.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7281e58f17a..23d7e01b1e4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,7 +31,6 @@ add_subdirectory(config) add_subdirectory(utils) add_subdirectory(profiler) -add_subdirectory(pdbparser) if(SCAFACOS) add_subdirectory(scafacos) diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index f692586062b..bb6cc82575a 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -173,7 +173,6 @@ endif() target_link_libraries(EspressoCore PUBLIC EspressoConfig utils) target_link_libraries(EspressoCore PUBLIC Boost::serialization Boost::mpi MPI::MPI_CXX Random123) target_link_libraries(EspressoCore PRIVATE Profiler) -target_link_libraries(EspressoCore PRIVATE pdbparser) if(FFTW3_FOUND) target_include_directories(EspressoCore PRIVATE SYSTEM ${FFTW3_INCLUDE_DIR}) diff --git a/src/pdbparser/CMakeLists.txt b/src/pdbparser/CMakeLists.txt deleted file mode 100644 index 3c4a66f1dc1..00000000000 --- a/src/pdbparser/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -add_library(pdbparser SHARED ${CMAKE_CURRENT_SOURCE_DIR}/src/PdbParser.cpp) -target_include_directories(pdbparser PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include) -install(TARGETS pdbparser LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd) -if(WITH_UNIT_TESTS) - add_subdirectory(unit_tests) -endif(WITH_UNIT_TESTS) diff --git a/src/pdbparser/include/PdbParser.hpp b/src/pdbparser/include/PdbParser.hpp deleted file mode 100644 index fa842aa5baa..00000000000 --- a/src/pdbparser/include/PdbParser.hpp +++ /dev/null @@ -1,72 +0,0 @@ -/* - * Copyright (C) 2014-2019 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 . - */ - -#ifndef __PDBPARSER_HPP -#define __PDBPARSER_HPP - -#include -#include -#include -#include -#include -#include - -namespace PdbParser { - -struct BoundingBox { - float llx, lly, llz; - float urx, ury, urz; -}; - -typedef struct { - int i; // index - float x, y, z; -} pdb_atom; - -typedef struct { - int i; - std::string type; - float charge; -} itp_atom; - -typedef struct { - int id, espresso_id; - float sigma, epsilon; -} itp_atomtype; - -struct itp_atomtype_compare { - bool operator()(const itp_atomtype &a, const itp_atomtype &b) const { - return a.id < b.id; - } -}; - -class PdbParser { -public: - bool parse_pdb_file(const std::string &filename); - bool parse_itp_file(const std::string &filename); - bool parse_file(const std::string &pdb_filename, - const std::string &itp_filename); - BoundingBox calc_bounding_box() const; - std::vector pdb_atoms; - std::map itp_atoms; - std::map itp_atomtypes; -}; -} // namespace PdbParser - -#endif diff --git a/src/pdbparser/src/PdbParser.cpp b/src/pdbparser/src/PdbParser.cpp deleted file mode 100644 index aa6ea3d3e4a..00000000000 --- a/src/pdbparser/src/PdbParser.cpp +++ /dev/null @@ -1,188 +0,0 @@ -/* - * Copyright (C) 2014-2019 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 . - */ - -#include "PdbParser.hpp" - -#include -#include - -using namespace std; - -namespace PdbParser { - -BoundingBox PdbParser::calc_bounding_box() const { - BoundingBox bb{}; - - bb.llx = std::numeric_limits::max(); - bb.lly = std::numeric_limits::max(); - bb.llz = std::numeric_limits::max(); - bb.urx = -std::numeric_limits::max(); - bb.ury = -std::numeric_limits::max(); - bb.urz = -std::numeric_limits::max(); - - for (auto pdb_atom : pdb_atoms) { - bb.llx = std::min(pdb_atom.x, bb.llx); - bb.lly = std::min(pdb_atom.y, bb.lly); - bb.llz = std::min(pdb_atom.z, bb.llz); - bb.urx = std::max(pdb_atom.x, bb.urx); - bb.ury = std::max(pdb_atom.y, bb.ury); - bb.urz = std::max(pdb_atom.z, bb.urz); - } - return bb; -} - -bool PdbParser::parse_pdb_file(const string &filename) { - ifstream file; - string tmp; - pdb_atom a; - - pdb_atoms.clear(); - - try { - file.open(filename.c_str()); - while (file.good()) { - - file >> tmp; - if (tmp == "ATOM") { - file.ignore(246, ' '); - file >> a.i; - file >> tmp >> tmp >> tmp >> tmp; - file >> a.x >> a.y >> a.z; - pdb_atoms.push_back(a); - } - } - } catch (ifstream::failure &e) { - return false; - } - - return true; -} - -bool PdbParser::parse_itp_file(const string &filename) { - ifstream file(filename.c_str()); - string tmp, buf; - itp_atom atom; - std::size_t pos; - - itp_atoms.clear(); - itp_atomtypes.clear(); - - while (file.good()) { - try { - buf = char(file.get()); - /* Skip leading whitespace */ - if (std::isspace(buf[0])) - continue; - - /* Comment, ignore rest of line */ - if (buf[0] == ';') { - std::getline(file, buf); - continue; - } - - /* Section statement */ - if (buf == "[") { - std::getline(file, buf); - pos = buf.find_first_not_of(" \t["); - if (pos == std::string::npos) - continue; - - std::string section = buf.substr(pos, std::string::npos); - pos = section.find_first_of(']'); - section = section.substr(0, pos); - pos = section.find_last_not_of(" \t"); - section = section.substr(0, pos + 1); - - if (section == "atoms") { - while (file.good()) { - buf = char(file.get()); - - /* Ignore leading whitespace, check for end of file */ - if (std::isspace(buf[0]) || file.eof()) { - continue; - } - /* End of atoms section */ - if (buf[0] == '[') { - file.unget(); - break; - } - /* Comment, ignore line */ - if (buf[0] == ';') { - std::getline(file, buf); - continue; - } - /* Push back first char */ - file.unget(); - /* Parse line */ - std::getline(file, buf); - std::istringstream line(buf); - line >> atom.i >> atom.type >> tmp >> tmp >> tmp >> tmp >> - atom.charge; - itp_atoms.insert(std::pair(atom.i, atom)); - } - } - if (section == "atomtypes") { - itp_atomtype type; - std::string type_name; - while (file.good()) { - buf = char(file.get()); - - /* Ignore leading whitespace */ - if (std::isspace(buf[0])) { - continue; - } - /* End of atoms section */ - if (buf[0] == '[') { - file.unget(); - break; - } - - /* Ignore leading whitespace, check for end of file */ - if (std::isspace(buf[0]) || file.eof()) { - continue; - } - - /* Push back first char */ - file.unget(); - /* Parse line */ - std::getline(file, buf); - std::istringstream line(buf); - line >> type_name >> tmp >> tmp >> tmp >> tmp >> type.sigma >> - type.epsilon; - /* Id is sequential number starting from zero */ - type.id = itp_atomtypes.size(); - itp_atomtypes.insert( - std::pair(type_name, type)); - } - } - } - } catch (...) { - return false; - } - } - - return true; -} - -bool PdbParser::parse_file(const string &pdb_filename, - const string &itp_filename) { - return parse_pdb_file(pdb_filename) && parse_itp_file(itp_filename); -} - -} // namespace PdbParser diff --git a/src/pdbparser/unit_tests/CMakeLists.txt b/src/pdbparser/unit_tests/CMakeLists.txt deleted file mode 100644 index c1ae079be42..00000000000 --- a/src/pdbparser/unit_tests/CMakeLists.txt +++ /dev/null @@ -1,3 +0,0 @@ -include(unit_test) - -unit_test(NAME PdbParser_test SRC PdbParser_test.cpp DEPENDS pdbparser) diff --git a/src/pdbparser/unit_tests/PdbParser_test.cpp b/src/pdbparser/unit_tests/PdbParser_test.cpp deleted file mode 100644 index 9cd0051fb83..00000000000 --- a/src/pdbparser/unit_tests/PdbParser_test.cpp +++ /dev/null @@ -1,93 +0,0 @@ -/* - * Copyright (C) 2019 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 . - */ - -/** \file - * Unit tests for the PdbParser class. - */ - -#define BOOST_TEST_MODULE PdbParser test -#define BOOST_TEST_DYN_LINK -#include - -#include - -#include "PdbParser.hpp" - -BOOST_AUTO_TEST_CASE(parser) { - // create input files - const std::string pdb_filepath("PdbParser.pdb"); - const std::string itp_filepath("PdbParser.itp"); - std::ofstream file; - file.open(pdb_filepath); - file << "ATOM 1 O1 LIG A 1 1.000 1.000 1.000 1.00\n"; - file << "ATOM 2 C1 LIG A 1 2.000 -0.500 1.000 1.00"; - file.close(); - file.open(itp_filepath); - file << "[ atomtypes ]\n"; - file << ";name bond_type mass charge ptype sigma epsilon\n"; - file << " c c 12.010 0.7700 A 1.000e-01 2.000e-01\n"; - file << "\n"; - file << "[ atoms ]\n"; - file << "; nr type resnr residue atom cgnr charge mass\n"; - file << " 1 O 1 LIG O1 1 -0.7700 16\n"; - file << " 2 c 1 LIG C1 2 0.7700 12.01"; - file.close(); - - PdbParser::PdbParser topology{}; - // check file parser - const bool success = topology.parse_file(pdb_filepath, itp_filepath); - BOOST_CHECK(success); - // check bounding box - const auto box = topology.calc_bounding_box(); - BOOST_CHECK(box.llx == 1.0); - BOOST_CHECK(box.lly == -0.5); - BOOST_CHECK(box.llz == 1.0); - BOOST_CHECK(box.urx == 2.0); - BOOST_CHECK(box.ury == 1.0); - BOOST_CHECK(box.urz == 1.0); - // check itp atomtypes - const auto O_atomtype = topology.itp_atomtypes["O"]; - BOOST_CHECK(O_atomtype.id == 0); - BOOST_CHECK(O_atomtype.sigma == 0.0f); - BOOST_CHECK(O_atomtype.epsilon == 0.0f); - const auto c_atomtype = topology.itp_atomtypes["c"]; - BOOST_CHECK(c_atomtype.id == 1); - BOOST_CHECK(c_atomtype.sigma == 0.1f); - BOOST_CHECK(c_atomtype.epsilon == 0.2f); - // check itp atoms - const auto O_itp_atom = topology.itp_atoms[1]; - BOOST_CHECK(O_itp_atom.i == 1); - BOOST_CHECK(O_itp_atom.type == "O"); - BOOST_CHECK(O_itp_atom.charge == -0.77f); - const auto c_itp_atom = topology.itp_atoms[0]; - BOOST_CHECK(c_itp_atom.i == 0); - BOOST_CHECK(c_itp_atom.type.empty()); - BOOST_CHECK(c_itp_atom.charge == 0.0f); - // check pdb atoms - const auto O_pdbatom = topology.pdb_atoms[0]; - BOOST_CHECK(O_pdbatom.i == 1); - BOOST_CHECK(O_pdbatom.x == 1.0f); - BOOST_CHECK(O_pdbatom.y == 1.0f); - BOOST_CHECK(O_pdbatom.z == 1.0f); - const auto c_pdbatom = topology.pdb_atoms[1]; - BOOST_CHECK(c_pdbatom.i == 2); - BOOST_CHECK(c_pdbatom.x == 2.0f); - BOOST_CHECK(c_pdbatom.y == -0.5f); - BOOST_CHECK(c_pdbatom.z == 1.0f); -} From f2fbd3ef0df32d0a69727c8838a9cb18e45f2e09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 21:26:51 +0200 Subject: [PATCH 11/11] Document parsing PDB with MDAnalysis --- doc/sphinx/io.rst | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst index 9706b8d2455..daced001ace 100644 --- a/doc/sphinx/io.rst +++ b/doc/sphinx/io.rst @@ -397,3 +397,39 @@ using MDAnalysis. A simple example is the following: W.write_next_timestep(u.trajectory.ts) # append it to the trajectory For other examples, see :file:`/samples/MDAnalysisIntegration.py` + +.. _Reading various formats using MDAnalysis: + +Reading various formats using MDAnalysis +---------------------------------------- + +MDAnalysis can read various formats, including MD topologies and trajectories. +To read a PDB file containing a single frame:: + + import MDAnalysis + import numpy as np + import espressomd + from espressomd.interactions import HarmonicBond + + # parse protein structure + universe = MDAnalysis.Universe("protein.pdb") + # extract only the C-alpha atoms of chain A + chainA = universe.select_atoms("name CA and segid A") + # use the unit cell as box + box_l = np.ceil(universe.dimensions[0:3]) + # setup system + system = espressomd.System(box_l=box_l) + system.time_step = 0.001 + system.cell_system.skin = 0.4 + # configure sphere size sigma and create a harmonic bond + system.non_bonded_inter[0, 0].lennard_jones.set_params( + epsilon=1, sigma=1.5, cutoff=2, shift="auto") + system.bonded_inter[0] = HarmonicBond(k=0.5, r_0=1.5) + # create particles and add bonds between them + system.part.add(pos=np.array(chainA.positions, dtype=float)) + for i in range(0, len(chainA) - 1): + system.part[i].add_bond((system.bonded_inter[0], system.part[i + 1].id)) + # visualize protein in 3D + from espressomd import visualization + visualizer = visualization.openGLLive(system, bond_type_radius=[0.2]) + visualizer.run(0)