Skip to content

Commit

Permalink
Merge branch 'python' into patch-6
Browse files Browse the repository at this point in the history
  • Loading branch information
fweik authored Jun 29, 2019
2 parents 9ab98b7 + bb1776b commit 160d172
Show file tree
Hide file tree
Showing 173 changed files with 2,078 additions and 6,494 deletions.
3 changes: 3 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,9 @@ cmake_minimum_required(VERSION 3.4)
include(FeatureSummary)
project(ESPResSo)
include(cmake/FindPythonModule.cmake)
if(NOT ${CMAKE_VERSION} VERSION_LESS "3.12")
cmake_policy(SET CMP0074 NEW)
endif()

enable_language(CXX)

Expand Down
3 changes: 2 additions & 1 deletion bors.toml
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
status = [
"ICP GitLab CI"
]
timeout-sec = 14400
timeout_sec = 14400
required_approvals = 1
106 changes: 1 addition & 105 deletions doc/sphinx/advanced_methods.rst
Original file line number Diff line number Diff line change
Expand Up @@ -119,109 +119,6 @@ The following limitations currently apply for the collision detection:
* The ``"bind at point of collision"`` approach cannot handle collisions
between virtual sites

.. _Swimmer Reactions:

Swimmer Reactions
-----------------


With the help of the feature ``SWIMMER_REACTIONS``, one can define three particle types to act as reactant (e.g. :math:`\mathrm{H_2 O_2}`), catalyzer (e.g. platinum), and product (e.g. :math:`\mathrm{O_2}` and :math:`\mathrm{H_2 O}`). The current setup allows one to simulate active swimmers and their chemical propulsion.

For a Janus swimmer consisting of platinum on one hemisphere and gold on the other hemisphere, both surfaces catalytically induce a reaction. We assume an initial abundance of hydrogen peroxide and absence of products, so that back (recombination) reactions seldom occur at the surface. A typical model for the propulsion of such a particle assumes

.. math::
\begin{aligned}
\mathrm{H_2 O_2} &\xrightarrow{\text{Pt}} \mathrm{2 H^{+} + 2 e^{-} + O_2} \\
\mathrm{2 H^{+} + 2 e^{-} + H_2 O_2} &\xrightarrow{\text{Au}} \mathrm{2 H_2 O}
\end{aligned}
That is, catalytic surfaces induce a reactions that produce charged species by consuming hydrogen peroxide. It is the change in distribution of charged species that leads to motion of the swimmer, a process referred to as self-electrophoresis. A minimal model for this would be

.. math::
\begin{aligned}
A &\xrightarrow{C^{+}} B \\
B &\xrightarrow{C^{-}} A
\end{aligned}
where on the upper half of the catalyst :math:`C^{+}` a species :math:`A` is converted into :math:`B`, and on the lower half :math:`C^{-}` the opposite reaction takes place. Note that when :math:`A` and :math:`B` are charged, this reaction conserves charge, provided the rates are equal. Note that this feature uses the word catalyst in a meaning which cannot be brought into agreement with the definition of a catalyst. If the catalyst :math:`C^{+}` catalyzes (on average) the reaction, where :math:`A` is converted to :math:`B`, then it is impossible that a catalyst :math:`C^{-}` performs (on average) the reverse reaction. For the example with hydrogen peroxide this would mean that hydrogen peroxide is created spontaneously using a catalyst (under the same environment where another catalyst wants to split hydrogen peroxide). This is chemically impossible. What is meant to be modeled is that hydrogen peroxide is constantly flowing into the system from the bulk and therefore it is not depleted. This behaviour cannot be modeled using a catalyst (in the defined meaning of the word catalyst).
In |es| the orientation of a catalyzer particle is used to define hemispheres; half spaces going through the particle's center. The reaction region is bounded by the *reaction range*: :math:`r`. Inside the reaction range, we react only reactant-product pairs. The particles in a pair are swapped from hemisphere to another with a rate prescribed by

.. math::
P_{\text{move}} = 1 - \mathrm{e}^{-k_{\mathrm{ct}}\,\Delta t} ,
with the reaction rate :math:`k_{\mathrm{ct}}` and the simulation time step :math:`\Delta t`. A pair may be swapped only once per MD time step, to avoid a no-net-effect situation. That is, we allow an exchange move only when the following conditions are met:

1. Both partners of the reactant-product pair have to reside within the reaction range.
2. The product has to reside in the upper half-space of the reaction range.
3. The reactant has to reside in the lower half-space of the reaction range.

Self-propulsion is achieved by imposing an interaction asymmetry between the partners of a swapped pair. That is, the heterogeneous distribution of chemical species induced by the swapping leads to a net force on the particle, counter balanced by friction.

To set up the system for catalytic reactions the class :class:`espressomd.swimmer_reaction.Reaction`
can be used. ::

from espressomd.swimmer_reaction import Reaction

system = espressomd.System()

# setting up particles etc

r = Reaction(product_type=1, reactant_type=2, catalyzer_type=0,
ct_range=2, ct_rate=0.2, eq_rate=0)
r.start()
r.stop()

print r

* the first invocation of ``Reaction``, in the above example, defines a
reaction with particles of type number 2 as reactant, type 0 as catalyzer and
type 1 as product [#1]_. The catalytic reaction rate constant is given by :math:`\mathrm{ct\_rate}`
[#2]_ and to override the default rate constant for the equilibrium reaction
( = 0), one can specify it by as ``eq_rata``. By default each reactant particle is checked
against each catalyst particle (``react_once=False``). However, when creating
smooth surfaces using many catalyst particles, it can be desirable to let the
reaction rate be independent of the surface density of these particles. That
is, each particle has a likelihood of reacting in the vicinity of the surface
(distance is less than :math:`r`) as specified by the rate constant, i.e.,
*not* according to :math:`P_{\text{cvt}} = 1 - \exp \left( - n k\Delta t
\right)`, with :math:`n` the number of local catalysts. To accomplish this,
each reactant is considered only once each time step by using the option
``react_once=True`` . The reaction command is set up such that the different
properties may be influenced individually.

* ``r.stop()`` disables the reaction. Note that at the moment, there can
only be one reaction in the simulation.

* ``print r`` returns the current reaction parameters.

In future versions of |es| the capabilities of the ``SWIMMER_REACTIONS`` feature may be generalized
to handle multiple reactant, catalyzer, and product types, as well as
more general reaction schemes. Other changes may involve merging the
current implementation with the ``COLLISION_DETECTION`` feature.

.. rubric:: Footnotes

.. [#1]
Only one type of particle can be assigned to each of these three
reaction species and no particle type may be assigned to multiple
species. That is, currently does not support particles of type 1 and
2 both to be reactants, nor can particles of type 1 be a reactant as
well as a catalyst. Moreover, only one of these reactions can be
implemented in a single Tcl script. If, for instance, there is a
reaction involving particle types 1, 2, and 4, there cannot be a
second reaction involving particles of type 5, 6, and 8. It is
however possible to modify the reaction properties for a given set of
types during the simulation.
.. [#2]
Currently only strictly positive values of the catalytic conversion
rate constant are allowed. Setting the value to zero is equivalent to
``r.stop()``.
..
.. _Lees-Edwards boundary conditions:

Expand Down Expand Up @@ -253,8 +150,7 @@ friction coefficient (this is different from the Object-in-Fluid method
below). We implement these marker points as virtual tracer
particles which are not integrated using the usual velocity-Verlet
scheme, but instead are propagated using a simple Euler algorithm with
the local fluid velocity (if the ``IMMERSED_BOUNDARY`` feature is turned
on).
the local fluid velocity.

The immersed boundary method consists of two components, which can be used independently:

Expand Down
2 changes: 1 addition & 1 deletion doc/sphinx/analysis.rst
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ For example, ::

Momentum of the System
~~~~~~~~~~~~~~~~~~~~~~
:meth:`espressomd.analyze.Analysis.analyze_linear_momentum`
:meth:`espressomd.analyze.Analysis.linear_momentum`

This command returns the total linear momentum of the particles and the
lattice Boltzmann (LB) fluid, if one exists. Giving the optional
Expand Down
7 changes: 5 additions & 2 deletions doc/sphinx/conf.py.in
Original file line number Diff line number Diff line change
Expand Up @@ -193,8 +193,8 @@ html_show_sourcelink = False
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
html_show_sphinx = False

# If true, "(C) Copyright (C) 2018 The ESPResSo project
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
# If true, "(C) Copyright (C) 2018 The ESPResSo project" is shown
# in the HTML footer. Default is True.
html_show_copyright = True

# If true, an OpenSearch description file will be output, and all pages will
Expand Down Expand Up @@ -241,6 +241,9 @@ rst_prolog = """
rst_epilog = """
.. |es| replace:: `ESPResSo`
"""
# Re-enable grouping of arguments in numpydoc
# https://github.com/sphinx-doc/sphinx/issues/2227
napoleon_use_param = False


# -- Other options --------------------------------------------------------
Expand Down
46 changes: 24 additions & 22 deletions doc/sphinx/constraints.rst
Original file line number Diff line number Diff line change
Expand Up @@ -78,19 +78,19 @@ The extra argument ``particle_type`` specifies the non-bonded interaction to be
that constraint.

There are two additional optional parameters
to fine tune the behavior of the constraint. If ``penetrable`` is
set to ``True`` then particles can move through the constraint. In this case the
other option ``only_positive`` controls whether the particle is subject to the interaction
potential of the wall. If set to ``False``, then the constraint will only act in
the direction of the normal vector.
to fine tune the behavior of the constraint. If ``penetrable`` is set to
``True`` then particles can move through the constraint. In this case the
other option ``only_positive`` controls whether the particle is subject to the
interaction potential of the wall. If set to ``False``, then the constraint
will only act in the direction of the normal vector.

If we wanted to add a non-penetrable pore constraint to our simulation,
we could do the following::

pore = espressomd.shapes.SimplePore(
axis=[1, 0, 0], length=2, pos=[15, 15, 15], radius=1, smoothing_radius=0.5)
pore_constraint = espressomd.constraints.ShapeBasedConstraint(
shape=pore, penetrable=0, particle_type=1)
shape=pore, penetrable=False, particle_type=1)
system.constraints.add(pore_constraint)

Interactions between the pore and other particles are then defined
Expand Down Expand Up @@ -197,11 +197,11 @@ Pictured is an example constraint with a ``Wall`` shape created with ::
wall = Wall(dist=20, normal=[0.1, 0.0, 1])
system.constraints.add(shape=wall, particle_type=0)

In variant (1) if the ``only_positive`` flag is set to 1, interactions are only calculated if
the particle is on the side of the wall in which the normal vector is
pointing.
In variant (1) if the ``only_positive`` flag is set to ``True``, interactions
are only calculated if the particle is on the side of the wall in which the
normal vector is pointing.
This has only an effect for penetrable walls. If the flag is
set to 1, then slip boundary interactions apply that are essential for
set to ``True``, then slip boundary interactions apply that are essential for
microchannel flows like the Plane Poiseuille or Plane Couette Flow.
You also need to use the tunable_slip interaction (see [sec:tunableSlip])
for this too work.
Expand All @@ -211,7 +211,8 @@ for this too work.
A sphere.

The resulting surface is a sphere with center ``center`` and radius ``radius``.
The direction ``direction`` determines the force direction, ``-1`` for inward and ``+1`` for outward.
The direction ``direction`` determines the force direction, ``-1`` for inward
and ``+1`` for outward.

.. _shape-sphere:

Expand Down Expand Up @@ -287,7 +288,7 @@ The direction ``direction`` determines the force direction, ``-1`` for inward an
b=[0.0, 0.0, 1.0],
c=[0.0, 1.0, 0.0],
direction=1)
system.constraints.add(shape=rhomboid, particle_type=0, penetrable=1)
system.constraints.add(shape=rhomboid, particle_type=0, penetrable=True)

creates a rhomboid defined by one corner located at ``[5.0, 5.0, 5.0]`` and three
adjacent edges, defined by the three vectors connecting the corner with its three neighboring corners, ``(1,1,0)`` , ``(0,0,1)`` and ``(0,1,0)``.
Expand All @@ -313,7 +314,7 @@ Pictured is an example constraint with a ``SimplePore`` shape created with ::
radius=12.5,
smoothing_radius=2,
center=[25, 25, 25])
system.constraints.add(shape=pore, particle_type=0, penetrable=1)
system.constraints.add(shape=pore, particle_type=0, penetrable=True)


:class:`espressomd.shapes.Stomatocyte`
Expand Down Expand Up @@ -352,7 +353,7 @@ Pictured is an example constraint with a ``Stomatocyte`` shape (with a closeup o
center=[25, 25, 25],
layer_width=3,
direction=1)
system.constraints.add(shape=stomatocyte, particle_type=0, penetrable=1)
system.constraints.add(shape=stomatocyte, particle_type=0, penetrable=True)



Expand Down Expand Up @@ -399,7 +400,7 @@ Pictured is an example constraint with a ``Slitpore`` shape created with ::
pore_width=5,
dividing_plane=40)

system.constraints.add(shape=slitpore, particle_type=0, penetrable=1)
system.constraints.add(shape=slitpore, particle_type=0, penetrable=True)


:class:`espressomd.shapes.SpheroCylinder`
Expand Down Expand Up @@ -457,11 +458,12 @@ Pictured is an example constraint with a ``Hollowcone`` shape created with ::
center=[25, 25, 25],
width=2,
direction=1)
system.constraints.add(shape=hollowcone, particle_type=0, penetrable=1)
system.constraints.add(shape=hollowcone, particle_type=0, penetrable=True)


For the shapes ``wall``; ``sphere``; ``cylinder``; ``rhomboid``; ``maze``; ``pore`` and ``stomatocyte``, constraints are able to be penetrated if ``penetrable`` is set to ``True``.
Otherwise, when the ``penetrable`` option is
For the shapes ``wall``, ``sphere``, ``cylinder``, ``rhomboid``, ``maze``,
``pore`` and ``stomatocyte``, constraints are able to be penetrated if
``penetrable`` is set to ``True``. Otherwise, when the ``penetrable`` option is
ignored or is set to ``False``, the constraint cannot be violated, i.e. no
particle can go through the constraint surface (|es| will exit if it does).

Expand Down Expand Up @@ -506,19 +508,19 @@ words, this option helps defines the sign of the normal surface vector.

For now, this may not sound useful but it can be practical when used with
together with constraint options such as ``penetrable`` or ``only_positive``.
In the former case, using non-penetrable surfaces with ``penetrable=0`` will
In the former case, using non-penetrable surfaces with ``penetrable=False`` will
cause |es| to throw an error is any distances between interacting particles and
constraints are found to be negative. This can be used to stop a simulation if
for one reason or another particles end up in an unwanted location.

The ``only_positive`` constraint option is used to define if a force should be
applied to a particle that has a negative distance. For example, consider the
same probe particle as in the previous case. The plot below shows the particle
force with ``only_positive=1``. Notice that when the distance is negative,
force with ``only_positive=True``. Notice that when the distance is negative,
forces are not applied at all to the particle. Thus the constraint surface is
either purely radially outwards (when ``direction=1``) or radially inwards
(when ``direction=-1``). Note that in both cases the constraint was set to be
penetrable with ``penetrable=1`` or else the simulation would crash whenever
penetrable with ``penetrable=True`` or else the simulation would crash whenever
the particle was found in any location that yields a negative distance.

.. figure:: figures/constraint-force.png
Expand All @@ -527,7 +529,7 @@ the particle was found in any location that yields a negative distance.
:height: 8.00000cm

The next figure shows what happens if we turn off the ``only_positive`` flag by
setting ``only_positive=0``. In this case the particle is pushed radially
setting ``only_positive=False``. In this case the particle is pushed radially
inward if it is inside the sphere and radially outward if it is outside. As
with the previous example, the constraint was set to be penetrable for this to
make sense.
Expand Down
Loading

0 comments on commit 160d172

Please sign in to comment.