Skip to content
Permalink
Browse files

Merge #3116

3116: Constant pH ipython notebook r=jngrad a=jonaslandsgesell

partially Fixes #3087 via simplifying the tutorial. We now do not perform reaction ensemble simulations or simulations with p3m in the tutorial.
This makes the tutorial easy and super fast.

What needs to be done is the integration with tests. @jngrad could you do this?
Also the old tutorial needs to be deleted and the folder renamed to 10-constant_pH_method

Co-authored-by: Jonas Landsgesell <jonaslandsgesell@users.noreply.github.com>
Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>
  • Loading branch information...
3 people committed Sep 10, 2019
2 parents e679e45 + 12cb87c commit 87905ae562f4df82a8a58b3279e3e61bf551a26e
@@ -299,6 +299,7 @@ Currently, the following tutorials are available:
* :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

.. _Sample scripts:

@@ -529,16 +529,13 @@ Particle number counting feature



Knowing the number of particles of a certain type in simulations in the grand canonical ensemble,
or other purposes, when particles of certain types are created and
deleted frequently is often of interest. Particle ids can be stored in a map for each
individual type and so random ids of particles of a certain type can be
drawn. ::
Knowing the number of particles of a certain type in simulations where particle numbers can fluctuate is of interest.
Particle ids can be stored in a map for each
individual type::

import espressomd
system = espressomd.System()
system.setup_type_map([_type])
system.find_particle(_type)
system.number_of_particles(_type)

If you want to keep track of particle ids of a certain type you have to
@@ -547,9 +544,7 @@ initialize the method by calling ::
system.setup_type_map([_type])

After that will keep track of particle ids of that type. Keeping track of particles of a given type is not enabled by default since it requires memory.
When using the
keyword ``find_particle`` and a particle type, the command will return a randomly
chosen particle id, for a particle of the given type. The keyword
The keyword
``number_of_particles`` as argument will return the number of
particles which have the given type. For counting the number of particles of a given type you could also use :meth:`espressomd.particle_data.ParticleList.select` ::

@@ -558,7 +553,7 @@ particles which have the given type. For counting the number of particles of a g
...
number_of_particles = len(system.part.select(type=type))

However calling ``select(type=type)`` results in looping over all particles. Therefore calling ``select()`` is slow compared to using :meth:`espressomd.system.System.number_of_particles` which directly can return the number of particles with that type.
However calling ``select(type=type)`` results in looping over all particles which is slow. In contrast, :meth:`espressomd.system.System.number_of_particles` directly can return the number of particles with that type.

.. _Self-propelled swimmers:

@@ -522,13 +522,13 @@
"_uncorrelated_ samples:\n",
"\n",
"\\begin{equation}\n",
" SE = \\sqrt{\\frac{\\sigma}{N}},\n",
" SE = \\sqrt{\\frac{\\sigma^2}{N}},\n",
"\\end{equation}\n",
"\n",
"where $\\sigma$ is the standard deviation\n",
"where $\\sigma^2$ is the variance\n",
"\n",
"\\begin{equation}\n",
" \\sigma = \\left\\langle x^2 - \\langle x\\rangle^2 \\right\\rangle\n",
" \\sigma^2 = \\left\\langle x^2 - \\langle x\\rangle^2 \\right\\rangle\n",
"\\end{equation}"
]
},
@@ -92,7 +92,7 @@ \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{,}
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}

Large diffs are not rendered by default.

@@ -0,0 +1,6 @@
configure_file(12-constant_pH.ipynb ${CMAKE_CURRENT_BINARY_DIR}/12-constant_pH.ipynb COPYONLY)

add_custom_target(tutorials_12)

nb_export(TUTORIAL tutorials_12 FILE "12-constant_pH.ipynb" HTML_RUN)

@@ -78,6 +78,7 @@ add_subdirectory(07-electrokinetics)
add_subdirectory(08-visualization)
add_subdirectory(10-reaction_ensemble)
add_subdirectory(11-ferrofluid)
add_subdirectory(12-constant_pH)


### here: add the appropriate tutorial target after DEPENDS
@@ -91,6 +92,7 @@ add_custom_target(tutorials DEPENDS
tutorials_08
tutorials_10
tutorials_11
tutorials_12
)

add_custom_target(tutorials_html DEPENDS
@@ -107,6 +109,7 @@ add_custom_target(tutorials_html DEPENDS
tutorials_11_1_html
tutorials_11_2_html
tutorials_11_3_html
tutorials_12_html
)

add_custom_target(tutorials_python_only DEPENDS
@@ -122,6 +125,7 @@ add_custom_target(tutorials_python_only DEPENDS
tutorials_11_1_python
tutorials_11_2_python
tutorials_11_3_python
tutorials_12_python
)

add_custom_target(tutorials_python DEPENDS
@@ -17,13 +17,16 @@ physical systems. Currently, the following tutorials are available:
* :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

Using the tutorials
-------------------
For using the tutorials, you need ESPResSo running. For installation
instructions, please see: http://espressomd.org/html/doc/installation.html

Tutorials 1, 2, 4, 5, and 8 are available as IPython notebooks. I.e., they consist of a `.ipynb` file which contains both, the source code and the corresponding explanations.
Tutorials 1, 2, 4, 5, 8 and 12 are available as IPython notebooks, i.e.
they consist of a `.ipynb` file which contains both the source code
and the corresponding explanations.
They can be viewed, changed and run interactively.


@@ -1491,11 +1491,11 @@ void remove_id_from_map(int part_id, int type) {
particle_type_map.at(type).erase(part_id);
}

int get_random_p_id(int type) {
if (particle_type_map.at(type).empty())
throw std::runtime_error("No particles of given type could be found");
int rand_index = i_random(particle_type_map.at(type).size());
return *std::next(particle_type_map[type].begin(), rand_index);
int get_random_p_id(int type, int random_index_in_type_map) {
if (random_index_in_type_map + 1 > particle_type_map.at(type).size())
throw std::runtime_error("The provided index exceeds the number of "
"particles listed in the type_map");
return *std::next(particle_type_map[type].begin(), random_index_in_type_map);
}

void add_id_to_type_map(int part_id, int type) {
@@ -930,7 +930,7 @@ void auto_exclusions(int distance);
void init_type_map(int type);

/* find a particle of given type and return its id */
int get_random_p_id(int type);
int get_random_p_id(int type, int random_index_in_type_map);
int number_of_particles_with_type(int type);

// The following functions are used by the python interface to obtain
@@ -251,7 +251,8 @@ bool ReactionAlgorithm::all_reactant_particles_exist(int reaction_id) {
*/
void ReactionAlgorithm::append_particle_property_of_random_particle(
int type, std::vector<StoredParticleProperty> &list_of_particles) {
int p_id = get_random_p_id(type);
int random_index_in_type_map = i_random(number_of_particles_with_type(type));
int p_id = get_random_p_id(type, random_index_in_type_map);
StoredParticleProperty property_of_part = {p_id, charges_of_types[type],
type};
list_of_particles.push_back(property_of_part);
@@ -791,12 +792,15 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
std::vector<int> p_id_s_changed_particles;

// save old_position
int p_id = get_random_p_id(type);
int random_index_in_type_map = i_random(number_of_particles_with_type(type));
int p_id = get_random_p_id(type, random_index_in_type_map);
for (int i = 0; i < particle_number_of_type_to_be_changed; i++) {
// determine a p_id you have not touched yet
while (is_in_list(p_id, p_id_s_changed_particles)) {
random_index_in_type_map = i_random(number_of_particles_with_type(type));
p_id = get_random_p_id(
type); // check whether you already touched this p_id, then reassign
type, random_index_in_type_map); // check whether you already touched
// this p_id, then reassign
}

auto part = get_particle_data(p_id);
@@ -182,6 +182,12 @@ class ReactionAlgorithm {
void restore_properties(std::vector<StoredParticleProperty> &property_list,
int number_of_saved_properties);

/**
* @brief draws a random integer from the uniform distribution in the range
* [0,maxint-1]
*
* @param maxint range.
*/
int i_random(int maxint) {
std::uniform_int_distribution<int> uniform_int_dist(0, maxint - 1);
return uniform_int_dist(m_generator);
@@ -47,5 +47,4 @@ cdef bool skin_set

cdef extern from "particle_data.hpp":
int init_type_map(int type) except +
int get_random_p_id(int type) except +
int number_of_particles_with_type(int type) except +
@@ -522,13 +522,3 @@ cdef class System:
self.check_valid_type(type)
number = number_of_particles_with_type(type)
return int(number)

def find_particle(self, type=None):
"""
The command will return a randomly chosen particle id, for a particle of
the given type.
"""
self.check_valid_type(type)
pid = get_random_p_id(type)
return int(pid)
@@ -84,7 +84,7 @@ class ReactionEnsembleTest(ut.TestCase):
WLRE.add_collective_variable_degree_of_association(
associated_type=0, min=0, max=1, corresponding_acid_types=[0, 1])
WLRE.set_wang_landau_parameters(
final_wang_landau_parameter=1e-2,
final_wang_landau_parameter=0.5 * 1e-2,
do_not_sample_reaction_partition_function=True,
full_path_to_output_filename="WL_potential_out.dat")

@@ -40,6 +40,7 @@ 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)
tutorial_test(FILE test_12-constant_pH.py)

add_custom_target(check_tutorials COMMAND ${CMAKE_CTEST_COMMAND} --timeout ${TEST_TIMEOUT} -j${TEST_NP} $(ARGS) --output-on-failure)

@@ -0,0 +1,44 @@
# 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 <http://www.gnu.org/licenses/>.

import unittest as ut
import importlib_wrapper
import numpy as np


tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
"@TUTORIALS_DIR@/12-constant_pH/12-constant_pH.py")


@skipIfMissingFeatures
class Tutorial(ut.TestCase):
system = tutorial.system

def test(self):
expected_values = tutorial.ideal_degree_of_dissociation(
tutorial.pHs, tutorial.pK)
simulated_values = tutorial.degrees_of_dissociation
simulated_values_error = tutorial.std_dev_degree_of_dissociation / \
np.sqrt(tutorial.num_samples)
# test alpha +/- 0.05 and standard error of alpha less than 0.10
np.testing.assert_allclose(expected_values, simulated_values, rtol=0,
atol=0.05)
self.assertLess(np.max(simulated_values_error), 0.10)


if __name__ == "__main__":
ut.main()

0 comments on commit 87905ae

Please sign in to comment.
You can’t perform that action at this time.