Propagators (#4820)

Fixes #4603

Description of changes:
- add a per-particle `propagation` flag to control which equations of motion are propagated during integration
- all virtual sites implementations can now be used simultaneously (`system.virtual_sites` was removed)
- virtual sites relative can now be coupled only to the Langevin and LB thermostats (possibly both: Langevin for rotation, LB for translation), using the new `particle.vs_relate_to()` keyword arguments `couple_to_lb=True` and `couple_to_langevin=True` (both are `False` by default)
- finer control over the propagation mode combinations can be obtained by setting the `particle.propagation` bitmask using combinations of `espressomd.propagation.Propagation` enum values; only a small number of combinations are currently allowed, but the list will grow as new propagation schemes are developed

doc/sphinx/CMakeLists.txt

@@ -87,6 +87,8 @@ if(SPHINX_FOUND)
     ${Python_EXECUTABLE} "${CMAKE_CURRENT_BINARY_DIR}/samples.py"
     COMMAND ${SPHINX_API_DOC_EXE} -f -o ${CMAKE_CURRENT_BINARY_DIR}
             ${SPHINX_PYTHON_DIR} ${EXCLUDE}
+    COMMAND ${Python_EXECUTABLE}
+            "${CMAKE_CURRENT_SOURCE_DIR}/autodoc_overrides.py"
     COMMAND
       ${SPHINX_EXECUTABLE} -q -W -b html -c "${CMAKE_CURRENT_BINARY_DIR}" -d
       "${SPHINX_CACHE_DIR}" "${CMAKE_CURRENT_BINARY_DIR}" "${SPHINX_HTML_DIR}"

doc/sphinx/advanced_methods.rst

@@ -68,7 +68,6 @@ Several modes are available for different types of binding.
   the point of contact or you can use :class:`espressomd.interactions.Virtual` which acts as a marker, only.
   The method is setup as follows::
 
-      system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
       system.collision_detection.set_params(
           mode="bind_at_point_of_collision",
           distance=0.1,
@@ -116,7 +115,6 @@ Several modes are available for different types of binding.
 
   The method is used as follows::
 
-      system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
       system.collision_detection.set_params(
             mode="glue_to_surface",
             distance=0.1,

src/script_interface/virtual_sites/CMakeLists.txt → doc/sphinx/autodoc_overrides.py

@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2020-2022 The ESPResSo project
+# Copyright (C) 2023 The ESPResSo project
 #
 # This file is part of ESPResSo.
 #
@@ -17,5 +17,21 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
-target_sources(espresso_script_interface
-               PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/initialize.cpp)
+
+def update_file(filepath, function):
+    with open(filepath, "r+") as f:
+        content = function(f.read())
+        f.seek(0)
+        f.truncate()
+        f.write(content)
+
+
+def update_espressomd(content):
+    token = ".. automodule:: espressomd.propagation"
+    assert token in content
+    assert content.count(token) == 1
+    content = content.replace(token, token + "\n   :member-order: bysource", 1)
+    return content
+
+
+update_file("espressomd.rst", update_espressomd)

doc/sphinx/particles.rst

@@ -251,26 +251,8 @@ according to their respective particle type. Before the next integration
 step, the forces accumulated on a virtual site are distributed back to
 those particles, from which the virtual site was derived.
 
-
-There are different schemes for virtual sites, described in the
-following sections. To switch the active scheme, the system
-:attr:`~espressomd.system.System.virtual_sites` property can be used::
-
-    import espressomd
-    import espressomd.virtual_sites
-
-    system = espressomd.System(box_l=[1, 1, 1])
-    system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(have_quaternion=False)
-    # or
-    system.virtual_sites = espressomd.virtual_sites.VirtualSitesOff()
-
-By default, :class:`espressomd.virtual_sites.VirtualSitesOff` is selected.
-This means that virtual particles are not touched during integration.
-The ``have_quaternion`` parameter determines whether the quaternion of
-the virtual particle is updated (useful in combination with the
-:attr:`~espressomd.particle_data.ParticleHandle.vs_quat` property of the
-virtual particle which defines the orientation of the virtual particle
-in the body fixed frame of the related real particle).
+There are different schemes for virtual sites, described in the following sections.
+The scheme acting on a virtual site is dependent on its propagation mode.
 
 .. _Rigid arrangements of particles:
 
@@ -298,19 +280,17 @@ the non-virtual particle rotates, the virtual sites rotates on an orbit
 around the non-virtual particles center.
 
 To use this implementation of virtual sites, activate the feature
-``VIRTUAL_SITES_RELATIVE``. Furthermore, an instance of
-:class:`~espressomd.virtual_sites.VirtualSitesRelative` has to be set as the
-active virtual sites scheme (see above). To set up a virtual site:
+``VIRTUAL_SITES_RELATIVE``. Furthermore, particles have to be set up with the
+propagation modes :attr:`~espressomd.propagation.Propagation.TRANS_VS_RELATIVE`
+and :attr:`~espressomd.propagation.Propagation.ROT_VS_RELATIVE`.
 
 #. Place the particle to which the virtual site should be related.
    It needs to be in the center of mass of the rigid arrangement of
    particles you create::
 
        import espressomd
-       import espressomd.virtual_sites
 
        system = espressomd.System(box_l=[10., 10., 10.])
-       system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
        p1 = system.part.add(pos=[1., 2., 3.])
 
 #. Place a particle at the desired relative position, make it virtual
@@ -320,8 +300,8 @@ active virtual sites scheme (see above). To set up a virtual site:
        p2 = system.part.add(pos=p1.pos + rel_offset)
        p2.vs_auto_relate_to(p1)
 
-   This will set the :attr:`~espressomd.particle_data.ParticleHandle.virtual`
-   attribute on particle ``p2`` to ``True``.
+   The :meth:`~espressomd.particle_data.ParticleHandle.is_virtual`
+   method on particle ``p2`` will now return ``True``.
 
 #. Repeat the previous step with more virtual sites, if desired.
 
@@ -371,22 +351,16 @@ Please note:
 Inertialess lattice-Boltzmann tracers
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-:class:`espressomd.virtual_sites.VirtualSitesInertialessTracers`
-
-When this implementation is selected, the virtual sites follow the motion of a
-lattice-Boltzmann fluid (both, CPU and GPU). This is achieved by integrating
+Using the propagation mode :attr:`~espressomd.propagation.Propagation.TRANS_LB_TRACER`,
+the virtual sites follow the motion of a LB fluid. This is achieved by integrating
 their position using the fluid velocity at the virtual sites' position.
 Forces acting on the virtual sites are directly transferred as force density
 onto the lattice-Boltzmann fluid, making the coupling free of inertia.
+Please note that the velocity attribute of the virtual particles
+does not carry valid information for this virtual sites scheme.
 The feature stems from the implementation of the
 :ref:`Immersed Boundary Method for soft elastic objects`, but can be used independently.
 
-For correct results, the LB thermostat has to be deactivated for virtual sites::
-
-   system.thermostat.set_lb(kT=0, act_on_virtual=False)
-
-Please note that the velocity attribute of the virtual particles does not carry valid information for this virtual sites scheme.
-
 .. _Interacting with groups of particles:
 
 Interacting with groups of particles
@@ -667,24 +641,22 @@ the total system is net force-free at all times. The resulting flow-field can th
 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:
-
-::
+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})
+    dipole = system.part.add(pos=[1, 0, 0],
+                             swimming={"f_swim": swimmer.swimming["f_swim"],
+                                       "is_engine_force_on_fluid": True})
+    dipole.vs_auto_relate_to(swimmer)
+    dipole.vs_quat = calc_quaternions_from_angles(np.pi, 0.)
 
-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.
+This makes the dipole particle not experience any friction or stochastic forces,
+but apply ``f_swim`` along its internal orientation (``director``)
+and follow the swimmer while updating its orientation accordingly.
 |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")
+    dipole = 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.

doc/tutorials/active_matter/active_matter.ipynb

@@ -874,11 +874,7 @@
    "outputs": [],
    "source": [
     "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",
-    ")"
+    "import espressomd.swimmer_helpers"
    ]
   },
   {

doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb

@@ -71,7 +71,6 @@
     "import espressomd.interactions\n",
     "import espressomd.electrostatics\n",
     "import espressomd.lb\n",
-    "import espressomd.virtual_sites\n",
     "\n",
     "import sys\n",
     "import tqdm\n",
@@ -335,43 +334,15 @@
    "id": "478a2301",
    "metadata": {},
    "source": [
-    "Now that the beads are arranged in the shape of a raspberry, the surface beads are made virtual particles\n",
+    "Now that the beads are arranged in the shape of a raspberry, the surface beads need to be set up as virtual particles\n",
     "using the `VirtualSitesRelative` scheme. This converts the raspberry to a rigid body\n",
-    "in which the surface particles follow the translation and rotation of the central particle.\n",
+    "in which the surface particles follow the translation and rotation of the central particle,\n",
+    "and transfer the interpolated LB momentum back to the central particle.\n",
     "Newton's equations of motion are only integrated for the central particle.\n",
     "It is given an appropriate mass and moment of inertia tensor (note that the inertia tensor\n",
     "is given in the frame in which it is diagonal.)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "cb7fc3a7",
-   "metadata": {},
-   "source": [
-    "#### Exercise\n",
-    "Activate the `VirtualSitesRelative` implementation in **ESPResSo**. See the [documentation](https://espressomd.github.io/doc/particles.html#virtual-sites) for help."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "fae9dbe4",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# SOLUTION CELL\n",
-    "# Select the desired implementation for virtual sites\n",
-    "system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "9f256d8e",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -435,7 +406,7 @@
     "# Convert the surface particles to virtual sites related to the central particle\n",
     "# The id of the central particles is 0, the ids of the surface particles start at 1.\n",
     "for p in surface_parts:\n",
-    "    p.vs_auto_relate_to(central_part)"
+    "    p.vs_auto_relate_to(central_part, couple_to_lb=True)"
    ]
   },
   {

samples/drude_bmimpf6.py

@@ -34,7 +34,6 @@
 import espressomd.electrostatics
 import espressomd.interactions
 import espressomd.drude_helpers
-import espressomd.virtual_sites
 import espressomd.visualization
 
 required_features = ["LENNARD_JONES", "P3M", "MASS", "ROTATION",
@@ -80,7 +79,6 @@
 print("\n-->Ion pairs:", n_ionpairs, "Box size:", box_l)
 
 system = espressomd.System(box_l=[box_l, box_l, box_l])
-system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
 
 if args.visu:
     d_scale = 0.988 * 0.5

samples/immersed_boundary/addSoft.py

@@ -31,7 +31,7 @@ def AddSoft(system, comX, comY, comZ, k1, k2):
             X = float(line[0]) + comX
             Y = float(line[1]) + comY
             Z = float(line[2]) + comZ
-            system.part.add(id=i, pos=[X, Y, Z], virtual=True)
+            system.part.add(id=i, pos=[X, Y, Z])
 
     # triangles
     import espressomd.interactions

samples/immersed_boundary/sampleImmersedBoundary.py

@@ -26,7 +26,6 @@
 import espressomd
 import espressomd.lb
 import espressomd.shapes
-import espressomd.virtual_sites
 
 required_features = ["VIRTUAL_SITES_INERTIALESS_TRACERS", "WALBERLA"]
 espressomd.assert_features(required_features)
@@ -48,7 +47,6 @@
 system = espressomd.System(box_l=(20, 20, boxZ))
 system.time_step = 1 / 6.
 system.cell_system.skin = 0.1
-system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers()
 print(f"Parallelization: {system.cell_system.node_grid}")
 
 force = 0.001
@@ -79,7 +77,7 @@
     agrid=1, density=1, kinematic_viscosity=1, tau=system.time_step,
     ext_force_density=[force, 0, 0])
 system.lb = lbf
-system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0, act_on_virtual=False)
+system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0)
 
 # Setup boundaries
 wall_shapes = [None] * 2

samples/rigid_body.py

@@ -26,17 +26,14 @@
 import math
 
 import numpy as np
-
 import espressomd
-import espressomd.virtual_sites
 import espressomd.rotation
 
 required_features = ["VIRTUAL_SITES_RELATIVE", "MASS", "ROTATIONAL_INERTIA"]
 espressomd.assert_features(required_features)
 
 
 system = espressomd.System(box_l=[10.0] * 3)
-system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
 system.time_step = 0.01
 system.thermostat.set_langevin(kT=1.0, gamma=20.0, seed=42)
 

src/core/CMakeLists.txt

@@ -39,6 +39,7 @@ add_library(
   particle_node.cpp
   polymer.cpp
   pressure.cpp
+  propagation.cpp
   rattle.cpp
   rotation.cpp
   Observable_stat.cpp

src/core/virtual_sites/VirtualSitesOff.hpp → src/core/CellParticleIterator.hpp

@@ -16,15 +16,10 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef VIRTUAL_SITES_VIRTUAL_SITES_OFF_HPP
-#define VIRTUAL_SITES_VIRTUAL_SITES_OFF_HPP
 
-#include "config/config.hpp"
-#ifdef VIRTUAL_SITES
-#include "VirtualSites.hpp"
+#pragma once
 
-/** @brief Do-nothing virtual-sites implementation */
-class VirtualSitesOff : public VirtualSites {};
+#include "ParticleIterator.hpp"
+#include "cell_system/Cell.hpp"
 
-#endif
-#endif
+using CellParticleIterator = ParticleIterator<Cell *const *>;