In this tutorial we explore the ways to simulate self-propulsion in the
simulation software package **ESPResSo**. We consider three examples that illustrate
the properties of these systems. First, we study the concept of enhanced
diffusion of a self-propelled particle. Second, we investigate rectification in
an asymmetric geometry. Finally, we determine the flow field around a
self-propelled particle using lattice-Boltzmann simulations (LB). These three
subsections should give insight into the basics of simulating active matter
with **ESPResSo**. The tutorial assumes basic knowledge of Python and **ESPResSo**,
as well as the use of lattice-Boltzmann within **ESPResSo**. It is therefore
recommended to go through the relevant tutorials first, before attempting this one.

Active matter is a term that describes a class of systems, in which agents (particles) constantly consume energy to perform work, for example in the form of directed motion. These systems are therefore highly out-of-equilibrium and thus defy description using the standard framework of statistical mechanics. Active systems are, however, ubiquitous. On our length scale, we encounter flocks of birds, schools of fish, and, of course, humans; on the mesoscopic level examples are found in bacteria, sperm, and algae; and on the nanoscopic level, transport along the cytoskeleton is achieved by myosin motors.

The `ENGINE` feature offers intuitive syntax for adding self-propulsion to
a particle. Propulsion will occur along the vector that defines the
orientation of the particle (henceforth referred to as ‘director’).
Within the `ENGINE` feature there are two ways of setting up a self-propelled
particle, with and without coupling to a fluid.

For this type of self-propulsion the Langevin thermostat can be used. The
Langevin thermostat imposes a velocity-dependent friction on a particle.
When a constant force is applied along the director,
the friction causes the particle to attain a terminal velocity, due to the balance
of driving and friction force, see Fig. 1. The timescale on
which the particle's velocity relaxes towards this value depends on the
strength of the friction and the mass of the particle. The `ENGINE`
feature implies that rotation of the particles (the `ROTATION` feature) is
compiled into **ESPResSo**. The particle can thus reorient due to external torques or
due to thermal fluctuations, whenever the rotational degrees of freedom are
thermalized. Note that the rotation of the particles has to be enabled
explicitly via their `rotation`

property. The ‘engine’ building block can
be connected to other particles, *e.g.*, via the virtual sites (rigid
body) to construct complex self-propelled objects.

First we import the necessary modules, define the parameters and set up the system.

In [1]:

```
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
import tqdm
import numpy as np
import espressomd.observables
import espressomd.accumulators
espressomd.assert_features(
["ENGINE", "ROTATION", "MASS", "ROTATIONAL_INERTIA", "WALBERLA", "VIRTUAL_SITES_RELATIVE"])
```

In [2]:

```
ED_PARAMS = {'time_step': 0.01,
'box_l': 3*[10.],
'skin': 0.4,
'f_swim': 5,
'kT': 1,
'gamma': 1,
'gamma_rotation': 1,
'mass': 0.1,
'rinertia': 3*[1.],
'corr_tmax': 100}
ED_N_SAMPLING_STEPS = 5000000
```

In [3]:

```
system = espressomd.System(box_l=ED_PARAMS['box_l'])
system.cell_system.skin = ED_PARAMS['skin']
system.time_step = ED_PARAMS['time_step']
```

- Use
`ED_PARAMS`

to set up a Langevin thermostat for translation and rotation of the particles.

In [4]:

```
system.thermostat.set_langevin(kT=ED_PARAMS['kT'],
gamma=ED_PARAMS['gamma'],
gamma_rotation=ED_PARAMS['gamma_rotation'],
seed=42)
```

`ParticleHandle` class of **ESPResSo**, which represents a particle in the
simulation. You can either set up the self-propulsion during the creation of a
particle or at a later stage.

In [5]:

```
part_act = system.part.add(pos=[5.0, 5.0, 5.0], swimming={'f_swim': ED_PARAMS['f_swim']},
mass=ED_PARAMS['mass'], rotation=3 * [True], rinertia=ED_PARAMS['rinertia'])
part_pass = system.part.add(pos=[5.0, 5.0, 5.0],
mass=ED_PARAMS['mass'], rotation=3 * [True], rinertia=ED_PARAMS['rinertia'])
```

**ESPResSo** correlators for the Mean Square Displacement (MSD), Velocity Autocorrelation Function (VACF) and the Angular Velocity Autocorrelation Function (AVACF).

In [6]:

```
pos_obs = espressomd.observables.ParticlePositions(
ids=[part_act.id, part_pass.id])
msd = espressomd.accumulators.Correlator(obs1=pos_obs,
corr_operation="square_distance_componentwise",
delta_N=1,
tau_max=ED_PARAMS['corr_tmax'],
tau_lin=16)
system.auto_update_accumulators.add(msd)
vel_obs = espressomd.observables.ParticleVelocities(
ids=[part_act.id, part_pass.id])
vacf = espressomd.accumulators.Correlator(obs1=vel_obs,
corr_operation="componentwise_product",
delta_N=1,
tau_max=ED_PARAMS['corr_tmax'],
tau_lin=16)
system.auto_update_accumulators.add(vacf)
ang_obs = espressomd.observables.ParticleAngularVelocities(
ids=[part_act.id, part_pass.id])
avacf = espressomd.accumulators.Correlator(obs1=ang_obs,
corr_operation="componentwise_product",
delta_N=1,
tau_max=ED_PARAMS['corr_tmax'],
tau_lin=16)
system.auto_update_accumulators.add(avacf)
```

No more setup needed! We can run the simulation and plot our observables.

In [7]:

```
for i in tqdm.tqdm(range(100)):
system.integrator.run(int(ED_N_SAMPLING_STEPS/100))
```

100%|██████████| 100/100 [02:15<00:00, 1.36s/it]

In [8]:

```
system.auto_update_accumulators.remove(msd)
msd.finalize()
system.auto_update_accumulators.remove(vacf)
vacf.finalize()
system.auto_update_accumulators.remove(avacf)
avacf.finalize()
```

In [9]:

```
taus_msd = msd.lag_times()
msd_result = msd.result()
msd_result = np.sum(msd_result, axis=2)
taus_vacf = vacf.lag_times()
vacf_result = np.sum(vacf.result(), axis=2)
taus_avacf = avacf.lag_times()
avacf_result = np.sum(avacf.result(), axis=2)
```

In [10]:

```
fig_msd = plt.figure(figsize=(10, 6))
plt.plot(taus_msd, msd_result[:, 0], label='active')
plt.plot(taus_msd, msd_result[:, 1], label='passive')
plt.xlim((taus_msd[1], None))
plt.loglog()
plt.xlabel('t')
plt.ylabel('MSD(t)')
plt.legend()
plt.show()
```

The Mean Square Displacement of an active particle is characterized by a longer ballistic regime and an increased diffusion coefficient for longer lag times. In the overdamped limit it is given by

$$ \langle r^{2}(t) \rangle = 6 D t + \frac{v^{2} \tau^{2}_{R}}{2} \left[ \frac{2 t}{\tau^{2}_{R}} + \exp\left( \frac{-2t}{\tau^{2}_{R}} \right) - 1 \right], $$where $\tau_{R} = \frac{8\pi\eta R^{3}}{k_{B} T}$ is the characteristic time scale for rotational diffusion and $ D = \frac{k_B T}{\gamma}$ is the translational diffusion coefficient.

For small times ($t \ll \tau_{R}$) the motion is ballistic

$$\langle r^{2}(t) \rangle = 6 D t + v^{2} t^{2},$$while for long times ($t \gg \tau_{R}$) the motion is diffusive

$$\langle r^{2}(t) \rangle = (6 D + v^{2}\tau_{R}) t.$$Note that no matter the strength of the activity, provided it is some finite value, the crossover between ballistic motion and enhanced diffusion is controlled by the rotational diffusion time.

The passive particle also displays a crossover from a ballistic to a diffusive motion. However, the crossover time $\tau_{C}=\frac{m}{\gamma}$ is not determined by the rotational motion but instead by the mass of the particles.

From the longterm MSD of the active particles we can define an effective diffusion coefficient $D_{\mathrm{eff}} = D + v^{2}\tau_{R}/6$. One could be tempted to also connect this increased diffusion with an effective temperature. However, this apparent equivalence leads to problems when one then attempts to apply statistical mechanics to such systems at the effective temperature. That is, there is typically more to being out-of-equilibrium than can be captured by a simple remapping of equilibrium parameters, as we will see in the second part of the tutorial.

In [11]:

```
def acf_stable_regime(x, y):
"""
Remove the noisy tail in autocorrelation functions of finite time series.
"""
cut = np.argmax(y <= 0.) - 2
assert cut >= 1
return (x[1:cut], y[1:cut])
```

In [12]:

```
fig_vacf = plt.figure(figsize=(10, 6))
plt.plot(*acf_stable_regime(taus_vacf, vacf_result[:, 0]), label='active')
plt.plot(*acf_stable_regime(taus_vacf, vacf_result[:, 1]), label='passive')
plt.xlim((taus_vacf[1], None))
plt.loglog()
plt.xlabel('t')
plt.ylabel('VACF(t)')
plt.legend()
plt.show()
```

In [13]:

```
fig_avacf = plt.figure(figsize=(10, 6))
plt.plot(*acf_stable_regime(taus_avacf, avacf_result[:, 0]), label='active')
plt.plot(*acf_stable_regime(taus_avacf, avacf_result[:, 1]), label='passive')
plt.xlim((taus_avacf[1], None))
plt.loglog()
plt.xlabel('t')
plt.ylabel('AVACF(t)')
plt.legend()
plt.show()
```

Before we go to the second part, it is important to clear the state of the system.

In [14]:

```
def clear_system(system):
system.part.clear()
system.thermostat.turn_off()
system.constraints.clear()
system.auto_update_accumulators.clear()
system.time = 0.
```

In [15]:

```
clear_system(system)
```

In [16]:

```
import espressomd.shapes
import espressomd.math
```

In [17]:

```
RECT_PARAMS = {'length': 100,
'radius': 20,
'funnel_inner_radius': 3,
'funnel_angle': np.pi / 4.0,
'funnel_thickness': 0.1,
'n_particles': 500,
'f_swim': 7,
'time_step': 0.01,
'wca_sigma': 0.5,
'wca_epsilon': 0.1,
'skin': 0.4,
'kT': 0.1,
'gamma': 1.,
'gamma_rotation': 1}
RECT_STEPS_PER_SAMPLE = 100
RECT_N_SAMPLES = 500
TYPES = {'particles': 0,
'boundaries': 1}
box_l = np.array(
[RECT_PARAMS['length'], 2*RECT_PARAMS['radius'], 2*RECT_PARAMS['radius']])
system.box_l = box_l
system.cell_system.skin = RECT_PARAMS['skin']
system.time_step = RECT_PARAMS['time_step']
system.thermostat.set_langevin(
kT=RECT_PARAMS['kT'], gamma=RECT_PARAMS['gamma'], gamma_rotation=RECT_PARAMS['gamma_rotation'], seed=42)
```

In [18]:

```
cylinder = espressomd.shapes.Cylinder(
center=0.5 * box_l,
axis=[1, 0, 0], radius=RECT_PARAMS['radius'], length=RECT_PARAMS['length'], direction=-1)
system.constraints.add(shape=cylinder, particle_type=TYPES['boundaries'])
# Setup walls
wall = espressomd.shapes.Wall(dist=0, normal=[1, 0, 0])
system.constraints.add(shape=wall, particle_type=TYPES['boundaries'])
wall = espressomd.shapes.Wall(dist=-RECT_PARAMS['length'], normal=[-1, 0, 0])
system.constraints.add(shape=wall, particle_type=TYPES['boundaries'])
```

Out[18]:

<espressomd.constraints.ShapeBasedConstraint at 0x73b64e1c4fe0>

In [19]:

```
funnel_length = (RECT_PARAMS['radius']-RECT_PARAMS['funnel_inner_radius']
)/np.tan(RECT_PARAMS['funnel_angle'])
```

- Using
`funnel_length`

and the geometric parameters in`RECT_PARAMS`

, set up the funnel cone (Hint: Use a Conical Frustum)

In [20]:

```
ctp = espressomd.math.CylindricalTransformationParameters(
axis=[1, 0, 0], center=box_l/2.)
hollow_cone = espressomd.shapes.HollowConicalFrustum(
cyl_transform_params=ctp,
r1=RECT_PARAMS['funnel_inner_radius'], r2=RECT_PARAMS['radius'],
thickness=RECT_PARAMS['funnel_thickness'],
length=funnel_length,
direction=1)
system.constraints.add(shape=hollow_cone, particle_type=TYPES['boundaries'])
```

Out[20]:

<espressomd.constraints.ShapeBasedConstraint at 0x73b64dec23e0>

- Set up a WCA potential between the walls and the particles using the parameters in
`RECT_PARAMS`

In [21]:

```
system.non_bonded_inter[TYPES['particles'], TYPES['boundaries']].wca.set_params(
epsilon=RECT_PARAMS['wca_epsilon'], sigma=RECT_PARAMS['wca_sigma'])
```

- 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) - Particles should be created with a random orientation

In [22]:

```
for i in range(RECT_PARAMS['n_particles']):
pos = box_l / 2
pos[0] += (-1)**i * 0.25 * RECT_PARAMS['length']
# https://mathworld.wolfram.com/SpherePointPicking.html
theta = np.arccos(2. * np.random.random() - 1)
phi = 2. * np.pi * np.random.random()
director = [np.sin(theta) * np.cos(phi),
np.sin(theta) * np.cos(phi),
np.cos(theta)]
system.part.add(pos=pos, swimming={'f_swim': RECT_PARAMS['f_swim']},
director=director, rotation=3*[True])
```

In [23]:

```
com_deviations = list()
times = list()
```

- Run the simulation using
`RECT_N_SAMPLES`

and`RECT_STEPS_PER_SAMPLE`

and calculate the deviation of the center of mass from the center of the box in each sample step. (Hint: Center of mass) - Save the result and the corresponding time of the system in the lists given above.

In [24]:

```
for _ in tqdm.tqdm(range(RECT_N_SAMPLES)):
system.integrator.run(RECT_STEPS_PER_SAMPLE)
com_deviations.append(system.galilei.system_CMS()[0] - 0.5 * box_l[0])
times.append(system.time)
```

100%|██████████| 500/500 [00:52<00:00, 9.60it/s]

In [25]:

```
def moving_average(data, window_size):
return np.convolve(data, np.ones(window_size), 'same') / window_size
```

In [26]:

```
smoothing_window = 10
com_smoothed = moving_average(com_deviations, smoothing_window)
fig_rect = plt.figure(figsize=(10, 6))
plt.plot(times[smoothing_window:-smoothing_window],
com_smoothed[smoothing_window:-smoothing_window])
plt.xlabel('t')
plt.ylabel('center of mass deviation')
plt.show()
```

In situations where hydrodynamic interactions between swimmers or swimmers and
objects are of importance, we use lattice-Boltzmann (LB) to simulate the fluid.
We recommend the GPU-based variant of LB in **ESPResSo**,
since it is much faster. Moreover, the current implementation of the CPU
self-propulsion with hydrodynamics is limited to one core.

**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.
Additionally, it also gives you the tools to apply a counterforce to the fluid that implicitly models the propulsion mechanism (flagella, cilia, etc).
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.
The orientation of the virtual particle must be pointing in the opposite direction as the swimmer.
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.
`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.

The keyword `mode`

(possible values: `"pusher"`

, `"puller"`

) determines where the counterforce on the fluid is applied, see again
Fig. 3(c, d).

Caveats of the algorithm are:

- One should make sure that the dipole length is at least one lattice Boltzmann grid spacing, since velocities are interpolated from that grid. If the length is less than one grid spacing, you can easily run into discretization artifacts or cause the particle not to move.
- The stability of the swimming algorithm is governed by the stability limitations of the LB method.
- Since the particle is
essentially a point particle, there is no rotation caused by shear
flow,
*e.g.*, a swimmer in a Poiseuille flow. If the thermostat is switched on, the rotational degrees of freedom will be thermalized, but there is no contribution of rotation due to ‘external’ flow fields. It is recommended to use an alternative means of obtaining rotations in your LB swimming simulations. For example, by constructing a raspberry particle.

In [27]:

```
import espressomd.lb
import espressomd.swimmer_helpers
```

In [28]:

```
clear_system(system)
```

In [29]:

```
HYDRO_PARAMS = {'box_l': 3*[35],
'time_step': 0.01,
'skin': 1,
'agrid': 1,
'dens': 0.1,
'visc': 5,
'gamma': 2,
'mass': 1,
'dipole_length': 4,
'f_swim': 0.01,
'dipole_particle_type': 1
}
HYDRO_N_STEPS = 2000
system.box_l = HYDRO_PARAMS['box_l']
system.cell_system.skin = HYDRO_PARAMS['skin']
system.time_step = HYDRO_PARAMS['time_step']
system.min_global_cut = HYDRO_PARAMS['dipole_length']
```

- Using
`HYDRO_PARAMS`

, set up a lattice-Boltzmann fluid and activate it as a thermostat (Hint: lattice-Boltzmann)

In [30]:

```
lbf = espressomd.lb.LBFluidWalberla(agrid=HYDRO_PARAMS['agrid'],
density=HYDRO_PARAMS['dens'],
kinematic_viscosity=HYDRO_PARAMS['visc'],
tau=HYDRO_PARAMS['time_step'])
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, gamma=HYDRO_PARAMS['gamma'], seed=42)
```

In [31]:

```
box_l = np.array(HYDRO_PARAMS['box_l'])
pos = box_l/2.
pos[2] = -15.
```

- Using
`HYDRO_PARAMS`

, place particle at`pos`

that swims in`z`

-direction. The particle handle should be called`particle`

. - Using
`espressomd.swimmer_helpers.add_dipole_particle`

and`HYDRO_PARAMS`

, create a dipole particle that applies the propulsion force to the fluid.

In [32]:

```
director = np.array([0,0,1])
particle = system.part.add(
pos=pos,
director=director,
mass=HYDRO_PARAMS['mass'], rotation=3*[False],
swimming={'f_swim': HYDRO_PARAMS['f_swim']})
espressomd.swimmer_helpers.add_dipole_particle(system, particle, HYDRO_PARAMS['dipole_length'], HYDRO_PARAMS['dipole_particle_type'])
```

Out[32]:

<espressomd.particle_data.ParticleHandle at 0x73b6841b5e40>

In [33]:

```
system.integrator.run(HYDRO_N_STEPS)
```

In [34]:

```
vels = np.squeeze(lbf[:, int(system.box_l[1]/2), :].velocity)
vel_abs = np.linalg.norm(vels, axis=2)
lb_shape = lbf.shape
xs, zs = np.meshgrid(np.linspace(0.5, box_l[0] - 0.5, num=lb_shape[0]),
np.linspace(0.5, box_l[2] - 0.5, num=lb_shape[2]))
fig_vels, ax_vels = plt.subplots(figsize=(10, 6))
im = plt.pcolormesh(vel_abs.T, cmap='YlOrRd')
vels = np.copy(vels)
ax_vels.streamplot(xs, zs, vels[:, :, 0].T, vels[:, :, 2].T, density=1.5)
circ = plt.Circle(particle.pos_folded[[0, 2]], 0.5, color='blue')
ax_vels.add_patch(circ)
ax_vels.set_aspect('equal')
plt.xlabel('x')
plt.ylabel('z')
cb = plt.colorbar(im, label=r'$|v_{\mathrm{fluid}}|$')
plt.show()
```