New feature: per-particle propagation flags

[Jean-Noël Grad]

Dear ESPResSo users,

with this mail, we'd like to inform people using the development version 
of ESPResSo, i.e., the python branch, of API changes.

Background:
One important feature in the next ESPResSo release will be to give users 
much more control over the equations of motion being used for particles. 
Examples include:
* using VirtualSitesRelative (rigid body mechanics) and 
VirtualSitesInertialessTracers (particles following a Lattice-Boltzmann 
fluid as tracers) together in a simulation
* coupling the translational degrees of freedom of a particle to a 
lattice-Boltzmann fluid while using the Langevin equation for the 
rotational degrees of freedom
* Deciding on a per-particle level, whether relative virtual sites' 
rotational degree of freedom is controlled by the central particle (e.g. 
dipole swimmers) or is propagated using Euler or Langevin equations 
(e.g., magnetic moments of a nanoparticle re-aligning without a rotation 
of the entire particle, Neel mechanism)

Summary:
* Due to the complexity of the changes and interactions with other 
ongoing development work, we have decided to merge the new functionality 
in steps.
   * Step 1 (now): change how relative virtual sites and LB inertialess 
tracers are set up and introduces the infrastructure to control the 
equation of motion on a per-particle level
   * Step 2 (Q1/2024): change how the system-wide defaults for 
equations of motions (integrator, thermostat) are set up
* Only a part of the new functionality described above will be 
introduced with step 1.

## Step 1 - being merged now.

Introduces the new infrastructure but limits the selection of equations 
of motions to mainly those previously present in ESPResSo.

Users not using virtual sites (relative or inertialess tracers) are not 
affected by the API changes in step 1 and can jump to step 2.

### Background

Particles have gained a new `propagation` attribute that allows setting 
per-particle integrators, which works similarly to per-particle 
thermostat friction coefficients. The per-particle propagation flag is a 
bitmask constructed with the enum values documented in the user guide in 
chapter espressomd.propagation [1]. The behavior of the 
`system.integrator.set_*()` methods remain unchanged, and they set up 
the "main" integrator. Depending on which main integrator is selected, 
different "secondary" integrators become available. For now, only the 
velocity Verlet integrator is available as a secondary integrator, using 
flags `Propagation.TRANS_NEWTON` for translation following Newton's 
equations of motion and `Propagation.ROT_EULER` for rotation following 
Euler's equations of rotation; in this way, selected particles can be 
decoupled from a thermostat.

### API changes

* The new `ParticleHanlde.propagation` attribute on each particle 
controls which equation of motion is applied to that particle. By 
default `Propagation.SYSTEM_DEFAULT` is set, which means that the 
system-wide thermostat and integrator settings are used.
* The `System.virtual_sites` attribute and the 
`espressomd.virtual_sites` module have been removed, as this is now 
handled by setting the appropriate propagation for the particle.
* Relative virtual sites are still set up by 
`ParticleHanlde.vs_auto_relate_to()`, however, there are additional 
keyword arguments for controlling whether the virtual site should be 
coupled to a lattice-Boltzmann fluid (`couple_to_lb`) or receive noise 
and friction terms from the Langevin equation (`couple_to_langevin`). 
This will set the appropriate propagation flags on the virtual site.
* The `act_on_virtual` keyword argument to `system.thermostat.set_lb()` 
has been removed, as this is now controlled by the propagation attribute 
of the individual virtual site. Warning: until now, relative virtual 
sites would couple to the LB fluid by default, but now, they aren't 
coupled by default! This is not considered a silent change, because an 
old Python script won't silently run the simulation with the new default 
settings (it will fail early due to the missing VirtualSitesRelative 
class). When updating a LB script that contains virtual sites, be sure 
to add `couple_to_lb=True` where appropriate.
* The `act_on_virtual` keyword argument of all the other thermostats has 
been removed.
* To set up an LB inertialess tracers, set the `propagation` attribute 
of the particle to `propagation.TRANS_LB_TRACER`. This is one of the 
rare use cases where you will need to use the 
`espressomd.propagation.Propagation` enum.

### Example script

```python
from espressomd.propagation import Propagation
import espressomd
from espressomd.lb import LBFluidWalberla

# System setup
system = espressomd.System(box_l=(10, 10, 10))
system.time_step = 0.01
system.cell_system.skin = 0.4

# Add lattice-Boltzmann fluid
lbf = LBFluidWalberla(agrid=1, density=1,
                      kinematic_viscosity=1, tau=system.time_step)
system.lb = lbf
system.thermostat.set_lb(gamma=1,LB_fluid=lbf)

# Add virtual sites relative
base_particle = system.part.add(pos=(0, 0, 0))
virtual_site = system.part.add(pos=(1, 0, 0))
virtual_site.vs_auto_relate_to(base_particle, couple_to_lb=True)
print(virtual_site.propagation) # yields 
Propagation.ROT_VS_RELATIVE|TRANS_LB_MOMENTUM_EXCHANGE|TRANS_VS_RELATIVE

# Add LB tracer
lb_tracer = system.part.add(pos=(3, 3, 3), 
propagation=Propagation.TRANS_LB_TRACER)

system.integrator.run(1)
```

### Outlook on step 1

It is now possible to set up per-particle propagation flags such that:
* translations are integrated by velocity Verlet using the LB 
thermostat, rotations are integrated using either Euler or Langevin 
equations of rotation;
* translations are integrated by velocity Verlet without thermalization 
while all other particles in the system are coupled to a Langevin 
thermostat;
* rotations are integrated by the Euler integrator without 
thermalization while all other particles in the system are coupled to a 
Langevin thermostat;
* for relative virtual sites, it can be controlled on a per-particle 
level, whether they are coupled to a lattice-Boltzmann fluid via 
momentum exchange (e.g., as in the virtual sties making up a raspberry) 
or receive Langevin noise and friction terms.

To find out how to set up such flags, have a look at 
`test_propagation()` in `propagation_langevin.py` [2]. There is a test 
case for each integrator, whose names follow the pattern 
`testsuite/python/propagation_*.py`, where you can see which 
combinations of flags are currently allowed.

## Step 2 - due next year

In a second step, we will change the way in which defaults for the 
equations of motion are set for the system, i.e., setting integrators 
and thermostats.

We also plan on enable more options for a particle's propagation in the 
future, starting with relative virtual sites whose rotational degrees of 
freedom can be integrated according to the Euler or Langevin equation 
rather than being derived from the base particle. Finally, we will 
rewrite the thermostats as python code, at which point it will be 
necessary to use a fresh build folder (or delete build/src/python to 
force a rebuild) to avoid error messages about python modules not being 
found.

Best regards,
Jean-Noël Grad and Rudolf Weeber

Links:
[1] 
https://espressomd.github.io/doc/espressomd.html#module-espressomd.propagation
[2] 
https://github.com/espressomd/espresso/blob/cc3a319/testsuite/python/propagation_langevin.py#L310

--
Dr. Jean-Noël Grad, Research Software Engineer
Institute for Computational Physics
University of Stuttgart
Allmandring 3
70569 Stuttgart
Germany
Phone: +49(0)711/685-67715
Email: jgrad@icp.uni-stuttgart.de