Skip to content

Commit

Permalink
Make lb swimmer force center a particle (#4745)
Browse files Browse the repository at this point in the history 
Description of changes:
- Replace magic force center by a virtual site relative particle that applies the force to the fluid
- Add helper to attach said particle to a swimmer and restore the old behavior
- Remove v_swim feature
  • Loading branch information
@kodiakhq
kodiakhq[bot] committed Jun 20, 2023
2 parents a78400b + b8a6442 commit 52923c7
Show file tree
Hide file tree
Showing 19 changed files with 509 additions and 373 deletions.
96 changes: 55 additions & 41 deletions doc/sphinx/particles.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -629,60 +629,74 @@ Self-propelled swimmers
Langevin swimmers
~~~~~~~~~~~~~~~~~

.. note::

Requires feature ``ENGINE``.

::

import espressomd
system = espressomd.System(box_l=[1, 1, 1])
system.part.add(pos=[1, 0, 0], swimming={'f_swim': 0.03})

This enables the particle to be self-propelled in the direction determined by
its quaternion. For setting the particle's quaternion see
:attr:`~espressomd.particle_data.ParticleHandle.quat`. The self-propulsion
speed will relax to a constant velocity, that is specified by ``v_swim``.
Alternatively it is possible to achieve a constant velocity by imposing a
constant force term ``f_swim`` that is balanced by friction of a (Langevin)
thermostat. The way the velocity of the particle decays to the constant
terminal velocity in either of these methods is completely determined by the
friction coefficient. You may only set one of the possibilities ``v_swim`` *or*
``f_swim`` as you cannot relax to constant force *and* constant velocity at the
same time. Note that there is no real difference between ``v_swim`` and
``f_swim``, since the latter may always be chosen such that the same terminal
velocity is achieved for a given friction coefficient.
This enables the particle to be self-propelled along its director.
The terminal propulsion speed is determined by the friction of a (Langevin)
thermostat (``v_swim = f_swim / gamma``).

.. _Lattice-Boltzmann swimmers:

Lattice-Boltzmann swimmers
~~~~~~~~~~~~~~~~~~~~~~~~~~

.. note::

Requires features ``ENGINE`` and ``VIRTUAL_SITES_RELATIVE``.

::

import espressomd
system = espressomd.System(box_l=[1, 1, 1])
system.part.add(pos=[2, 0, 0], rotation=[True, True, True], swimming={
'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0})

For an explanation of the parameters ``v_swim`` and ``f_swim`` see the previous
item. In lattice-Boltzmann self-propulsion is less trivial than for regular MD,
because the self-propulsion is achieved by a force-free mechanism, which has
strong implications for the far-field hydrodynamic flow field induced by the
self-propelled particle. In |es| only the dipolar component of the flow field
of an active particle is taken into account. This flow field can be generated
by a *pushing* or a *pulling* mechanism, leading to change in the sign of the
dipolar flow field with respect to the direction of motion. You can specify the
nature of the particle's flow field by using one of the modes: ``pusher`` or
``puller``. You will also need to specify a ``dipole_length`` which determines
the distance of the source of propulsion from the particle's center. Note that
you should not put this distance to zero; |es| (currently) does not support
mathematical dipole flow fields.

You may ask: "Why are there two methods ``v_swim`` and ``f_swim`` for the
self-propulsion using the lattice-Boltzmann algorithm?" The answer is
straightforward. When a particle is accelerating, it has a monopolar flow-field
contribution which vanishes when it reaches its terminal velocity (for which
there will only be a dipolar flow field). The major difference between the
above two methods is that with ``v_swim`` the flow field *only* has a monopolar
moment and *only* while the particle is accelerating. As soon as the particle
reaches a constant speed (given by ``v_swim``) this monopolar moment is gone
and the flow field is zero! In contrast, ``f_swim`` always, i.e., while
accelerating *and* while swimming at constant force possesses a dipolar flow
field.
swimmer = system.part.add(pos=[2, 0, 0],
rotation=[True, True, True],
swimming={'f_swim': 0.01})

In lattice-Boltzmann, self-propulsion is less trivial than for regular MD, because for hydrodynamic
interactions it is important that the propulsion force is generated by the swimmer itself and
not by an external force. Since the swimmer can only push itself forward by pushing fluid backward,
the total system is net force-free at all times. The resulting flow-field can thus not contain
monopolar contributions in the far field and the slowest-decaying flow-field mode is a dipole.
In |es|, the propulsion mechanism can be mimicked by applying a force to the fluid that is equal in
magnitude but opposite in sign to the forward force on the swimming particle.
For this, particles can be marked as force appliers as shown in the following example:

::

dipole_partcl = system.part.add(pos=[1, 0, 0],
virtual = True,
swimming={'f_swim': swimmer.swimming['f_swim'],
'is_engine_force_on_fluid': True})

This makes the particle not experience any friction or stochastic forces, but apply ``f_swim``
along its internal orientation (``director``).
To make the force applier follow the swimmer and also update its orientation accordingly,
the :class:`espressomd.virtual_sites.VirtualSitesRelative` mechanism is used.
|es| provides a helper function :func:`~espressomd.swimmer_helpers.add_dipole_particle` to set
up the virtual particle with the correct distance, relative position and orientation::

import espressomd.swimmer_helpers.add_dipole_particle as add_dip
dipole_partcl = add_dip(system, swimmer, 2., 0, mode="pusher")

It creates pushers with the propulsion behind the swimmer and pullers with the
propulsion in front of the swimmer.

Notes:

* The terminal velocity is **not** ``v_swim = gamma_lb * f_swim``. The flow-field generated
by the dipole particle interacts in a non-trivial way with the swimmer particle.
You will have to calibrate your propulsion force to the desired swim velocity.
* For the same reason, do not place the dipole particle too close to the swimmer
(at least one grid spacing is recommended). Lattice-Boltzmann cannot create an exact,
mathematical dipole, which would require zero distance but diverging forces.
* Since the swimmer is a point particle, it cannot experience shear or rotational flow
components, only random rotational noise. If shear or rotational flow are important in your
system, consider creating an extended, raspberry particle as described in :ref:`Rigid arrangements of particles`.
77 changes: 45 additions & 32 deletions doc/tutorials/active_matter/active_matter.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@
"import espressomd.accumulators\n",
"\n",
"espressomd.assert_features(\n",
" [\"ENGINE\", \"ROTATION\", \"MASS\", \"ROTATIONAL_INERTIA\", \"WALBERLA\"])"
" [\"ENGINE\", \"ROTATION\", \"MASS\", \"ROTATIONAL_INERTIA\", \"WALBERLA\", \"VIRTUAL_SITES_RELATIVE\"])"
]
},
{
Expand All @@ -136,7 +136,7 @@
"ED_PARAMS = {'time_step': 0.01,\n",
" 'box_l': 3*[10.],\n",
" 'skin': 0.4,\n",
" 'active_velocity': 5,\n",
" 'f_swim': 5,\n",
" 'kT': 1,\n",
" 'gamma': 1,\n",
" 'gamma_rotation': 1,\n",
Expand Down Expand Up @@ -218,7 +218,7 @@
},
"source": [
"```python\n",
"part_act = system.part.add(pos=[5.0, 5.0, 5.0], swimming={'v_swim': ED_PARAMS['active_velocity']},\n",
"part_act = system.part.add(pos=[5.0, 5.0, 5.0], swimming={'f_swim': ED_PARAMS['f_swim']},\n",
" mass=ED_PARAMS['mass'], rotation=3 * [True], rinertia=ED_PARAMS['rinertia'])\n",
"part_pass = system.part.add(pos=[5.0, 5.0, 5.0],\n",
" mass=ED_PARAMS['mass'], rotation=3 * [True], rinertia=ED_PARAMS['rinertia'])\n",
Expand Down Expand Up @@ -506,7 +506,7 @@
" 'funnel_angle': np.pi / 4.0,\n",
" 'funnel_thickness': 0.1,\n",
" 'n_particles': 500,\n",
" 'active_velocity': 7,\n",
" 'f_swim': 7,\n",
" 'time_step': 0.01,\n",
" 'wca_sigma': 0.5,\n",
" 'wca_epsilon': 0.1,\n",
Expand Down Expand Up @@ -636,7 +636,7 @@
"source": [
"### Exercise\n",
"* Place an equal number of swimming particles (the total number should be `RECT_PARAMS['n_particles']`) in the left and the right part of the box such that the center of mass is exactly in the middle. (Hint: Particles do not interact so you can put multiple in the same position)\n",
"* Particles must be created with a random [orientation](https://espressomd.github.io/doc/particles.html#rotational-degrees-of-freedom-and-particle-anisotropy)"
"* Particles should be created with a random [orientation](https://espressomd.github.io/doc/particles.html#rotational-degrees-of-freedom-and-particle-anisotropy)"
]
},
{
Expand All @@ -657,7 +657,7 @@
" np.sin(theta) * np.cos(phi),\n",
" np.cos(theta)]\n",
"\n",
" system.part.add(pos=pos, swimming={'v_swim': RECT_PARAMS['active_velocity']},\n",
" system.part.add(pos=pos, swimming={'f_swim': RECT_PARAMS['f_swim']},\n",
" director=director, rotation=3*[True])\n",
"```"
]
Expand Down Expand Up @@ -708,9 +708,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"metadata": {},
"outputs": [],
"source": []
},
Expand Down Expand Up @@ -769,7 +767,7 @@
"that active systems are net force free. \n",
"This is because the particles have to generate their propulsion on their own\n",
"as opposed to being dragged by an external force.\n",
"As a consequence, the flow field around a swimmer does not contain a monopolar (Stokeslet) contribution. \n",
"As a consequence, the far field flow around a swimmer does not contain a monopolar (Stokeslet) contribution. \n",
"In the case of an E. Coli cell, see <a href='#fig:pusher-puller'>Fig.&nbsp;3</a>(a), the reasoning is as follows:\n",
"The corkscrew tail pushes against the fluid and the fluid pushes against the\n",
"tail, which generates the forward motion.\n",
Expand Down Expand Up @@ -806,14 +804,18 @@
"since it is much faster. Moreover, the current implementation of the CPU\n",
"self-propulsion with hydrodynamics is limited to one core.\n",
"\n",
"**ESPResSo** implements swimming coupled with hydrodynamics by applying a force to the particle along its director as in the 'dry' case above, but now also applying a counterforce to the fluid that implicitly models the propulsion mechanism (flagella, cilia, etc).\n",
"The friction between the particle and the fluid then results in net force-free swimming.\n",
"For the setup of swimming particles with hydrodynamics we cannot use the `v_swim` argument anymore because it is not trivial to determine the friction acting on the particle. Instead, we have to provide the keys `f_swim` and `dipole_length`. Together they determine the dipole strength and the terminal velocity of the swimmer.\n",
"**ESPResSo** implements swimming coupled with hydrodynamics by applying a force to the particle along its director as in the 'dry' case above, just with a different thermostat that provides the friction.\n",
"Additionally, it also gives you the tools to apply a counterforce to the fluid that implicitly models the propulsion mechanism (flagella, cilia, etc).\n",
"This is done by creating a virtual particle that follows the swimmer particle in a fixed relative configuration, i.e., it moves at constant distance with the swimmer and also adapts to the swimmer rotation.\n",
"The orientation of the virtual particle must be pointing in the opposite direction as the swimmer.\n",
"The virtual particle then marked as ``\"is_engine_force_on_fluid\"``, which results in it not experiencing any friction but instead applying the propulsion force ``f_swim`` to the fluid.\n",
"``espressomd.swimmer_helper.add_dipole_particle`` is a helper function that performs these steps for you and creates a pusher or puller type swimmer for you.\n",
"\n",
"The keyword ``mode`` (possible values: ``\"pusher\"``, ``\"puller\"``) determines where the counterforce on the fluid is applied, see again\n",
"<a href='#fig:pusher-puller'>Fig.&nbsp;3</a>(c, d).\n",
"\n",
"Caveats of the algorithm are:\n",
"- One should make sure that the `dipole_length` is at least one\n",
"- One should make sure that the dipole length is at least one\n",
"lattice Boltzmann grid spacing, since velocities are interpolated from that grid. \n",
"If the length is less than one grid spacing, you can easily run into discretization artifacts\n",
"or cause the particle not to move.\n",
Expand All @@ -835,7 +837,12 @@
"metadata": {},
"outputs": [],
"source": [
"import espressomd.lb"
"import espressomd.lb\n",
"import espressomd.virtual_sites\n",
"import espressomd.swimmer_helpers\n",
"system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(\n",
" have_quaternion=True\n",
")"
]
},
{
Expand All @@ -862,8 +869,9 @@
" 'gamma': 2,\n",
" 'mass': 1,\n",
" 'dipole_length': 4,\n",
" 'active_force': 0.01,\n",
" 'mode': 'pusher'}\n",
" 'f_swim': 0.01,\n",
" 'dipole_particle_type': 1\n",
" }\n",
"\n",
"HYDRO_N_STEPS = 2000\n",
"\n",
Expand Down Expand Up @@ -926,7 +934,8 @@
},
"source": [
"### Exercise\n",
"* Using `HYDRO_PARAMS`, place particle at `pos` that swims in `z`-direction. The particle handle should be called `particle`."
"* Using `HYDRO_PARAMS`, place particle at `pos` that swims in `z`-direction. The particle handle should be called `particle`.\n",
"* Using `espressomd.swimmer_helpers.add_dipole_particle` and `HYDRO_PARAMS`, create a dipole particle that applies the propulsion force to the fluid."
]
},
{
Expand All @@ -936,11 +945,13 @@
},
"source": [
"```python\n",
"director = np.array([0,0,1])\n",
"particle = system.part.add(\n",
" pos=pos, mass=HYDRO_PARAMS['mass'], rotation=3*[False],\n",
" swimming={'f_swim': HYDRO_PARAMS['active_force'],\n",
" 'mode': HYDRO_PARAMS['mode'],\n",
" 'dipole_length': HYDRO_PARAMS['dipole_length']})\n",
" pos=pos, \n",
" director=director,\n",
" mass=HYDRO_PARAMS['mass'], rotation=3*[False],\n",
" swimming={'f_swim': HYDRO_PARAMS['f_swim']})\n",
"espressomd.swimmer_helpers.add_dipole_particle(system, particle, HYDRO_PARAMS['dipole_length'], HYDRO_PARAMS['dipole_particle_type'])\n",
"```"
]
},
Expand Down Expand Up @@ -999,18 +1010,20 @@
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"vtk_base_dir = os.path.join('vtk_out', 'RESULTS_FLOW_FIELD')\n",
"vtk_identifier = f'T_{HYDRO_PARAMS[\"mode\"]}_P_{pos[2]}'\n",
"vtk_outdir = os.path.join(vtk_base_dir, vtk_identifier)\n",
"lb_vtk = espressomd.lb.VTKOutput(identifier=vtk_identifier,\n",
"import pathlib\n",
"vtk_base_dir = pathlib.Path('./vtk_out').resolve()\n",
"vtk_identifier = 'RESULTS_FLOW_FIELD'\n",
"vtk_outdir = vtk_base_dir / vtk_identifier\n",
"vtk_outdir.mkdir(parents=True, exist_ok=True)\n",
"lb_vtk = espressomd.lb.VTKOutput(lb_fluid=lbf, \n",
" identifier=vtk_identifier,\n",
" observables=[\"velocity_vector\"],\n",
" base_folder=vtk_base_dir,\n",
" prefix=\"lb_velocity\")\n",
" base_folder=vtk_base_dir.as_posix(),\n",
" prefix='lb_velocity')\n",
"lbf.add_vtk_writer(vtk=lb_vtk)\n",
"for i in range(HYDRO_N_STEPS // 100):\n",
" system.integrator.run(100)\n",
" system.part.writevtk(os.path.join(vtk_outdir, f'position_{i}.vtk'), types=[0])\n",
" system.part.writevtk(vtk_outdir / f'position_{i}.vtk', types=[0])\n",
" lb_vtk.write()"
]
},
Expand Down Expand Up @@ -1070,7 +1083,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -1084,7 +1097,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
"version": "3.10.11"
}
},
"nbformat": 4,
Expand Down
27 changes: 11 additions & 16 deletions src/core/Particle.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,30 +45,21 @@ inline bool get_nth_bit(uint8_t const bitfield, unsigned int const bit_idx) {
}
} // namespace detail

#ifdef ENGINE
/** Properties of a self-propelled particle. */
struct ParticleParametersSwimming {
/** Is the particle a swimmer. */
bool swimming = false;
/** Imposed constant force. */
double f_swim = 0.;
/** Constant velocity to relax to. */
double v_swim = 0.;
/** Flag for the swimming mode in a LB fluid.
* Values:
* - -1: pusher
* - +1: puller
* - 0: no swimming
*/
int push_pull = 0;
/** Distance of the source of propulsion from the particle
* center in a LB fluid.
*/
double dipole_length = 0.;
/** Is the particle a swimmer. */
bool swimming = false;
/** Whether f_swim is applied to the particle or to the fluid. */
bool is_engine_force_on_fluid = false;

template <class Archive> void serialize(Archive &ar, long int /* version */) {
ar &swimming &f_swim &v_swim &push_pull &dipole_length;
ar &f_swim &swimming &is_engine_force_on_fluid;
}
};
#endif

/** Properties of a particle which are not supposed to
* change during the integration, but have to be known
Expand Down Expand Up @@ -603,7 +594,9 @@ struct Particle { // NOLINT(bugprone-exception-escape)
};

BOOST_CLASS_IMPLEMENTATION(Particle, object_serializable)
#ifdef ENGINE
BOOST_CLASS_IMPLEMENTATION(ParticleParametersSwimming, object_serializable)
#endif
BOOST_CLASS_IMPLEMENTATION(ParticleProperties, object_serializable)
BOOST_CLASS_IMPLEMENTATION(ParticlePosition, object_serializable)
BOOST_CLASS_IMPLEMENTATION(ParticleMomentum, object_serializable)
Expand All @@ -617,7 +610,9 @@ BOOST_CLASS_IMPLEMENTATION(decltype(ParticleProperties::vs_relative),
object_serializable)
#endif

#ifdef ENGINE
BOOST_IS_BITWISE_SERIALIZABLE(ParticleParametersSwimming)
#endif
BOOST_IS_BITWISE_SERIALIZABLE(ParticleProperties)
BOOST_IS_BITWISE_SERIALIZABLE(ParticlePosition)
BOOST_IS_BITWISE_SERIALIZABLE(ParticleMomentum)
Expand Down
2 changes: 1 addition & 1 deletion src/core/forces.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ inline ParticleForce external_force(Particle const &p) {
#ifdef ENGINE
// apply a swimming force in the direction of
// the particle's orientation axis
if (p.swimming().swimming) {
if (p.swimming().swimming and !p.swimming().is_engine_force_on_fluid) {
f.f += p.swimming().f_swim * p.calc_director();
}
#endif
Expand Down
Loading

0 comments on commit 52923c7

Please sign in to comment.