Welcome to the raspberry electrophoresis **ESPResSo** tutorial! This tutorial assumes some basic knowledge of **ESPResSo**.

Electrophoresis denotes the process when charged colloids or macromolecules move through a fluid in an electric field. Charged colloids with a total surface charge $q$ are surrounded by a diffuse layer of ions whith equal but opposite charge $\approx -q$, as shown in Fig. 1. The colloid's overall charge including the ion shell is therefore close to zero. An external electric field of magnitude $E$ pulls the colloid and its surrounding ion shell in opposite directions, so that the overall force on the colloid is much smaller than one would intuitively expect, i.e. $F \ll |qE|$. Instead, the electric field causes the colloid and the ion shell to move relative to each other, allowing the colloid to *slip* through the surrounding ion shell. Hydrodynamic friction resists the relative motion of the colloid and the ion shell. Modelling electrophoresis therefore requires not only electrostatic interactions but also hydrodynamic interactions.

The first step is compiling **ESPResSo** with the appropriate flags, as listed in Sec. 2.
The tutorial starts by discussing how to build a colloid out of several MD beads. These particles typically
resemble a raspberry as can be seen in Fig. 1. After covering the construction of a raspberry colloid, we then
briefly discuss the inclusion of hydrodynamic interactions via a lattice-Boltzmann fluid. Finally we will cover
including ions via the restrictive primitive model (hard sphere ions) and the addition of an electric field
to measure the electrokinetic properties. This script will run a raspberry electrophoresis simulation and write the time and position of the colloid out to a file named `posVsTime.dat` in the same directory.
A sample set of data is included in the file `posVsTime_sample.dat`.

To run this tutorial, you will need to enable the following features in the `myconfig.hpp`

file when compiling **ESPResSo**:

```
#define ELECTROSTATICS
#define ROTATION
#define ROTATIONAL_INERTIA
#define EXTERNAL_FORCES
#define MASS
#define VIRTUAL_SITES_RELATIVE
#define LENNARD_JONES
```

These features are enabled in the default configuration, so if you have not created your own `myconfig.hpp`

, all of these options should already be active. To be able to use the GPU accelerated lattice-Boltzmann algorithm, CUDA should be activated, as explained in the user guide.

The first thing to do in any **ESPResSo** simulation is to import our espressomd features and set a few global simulation parameters:

In [1]:

```
import espressomd
import espressomd.interactions
import espressomd.electrostatics
import espressomd.lb
import sys
import tqdm
import logging
logging.basicConfig(level=logging.INFO, stream=sys.stdout)
espressomd.assert_features(["ELECTROSTATICS", "ROTATION", "ROTATIONAL_INERTIA", "EXTERNAL_FORCES",
"MASS", "VIRTUAL_SITES_RELATIVE", "LENNARD_JONES", "WALBERLA"])
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
# System parameters
#############################################################
box_l = 40. # size of the simulation box
skin = 0.3 # Skin parameter for the Verlet lists
time_step = 0.01
integ_steps = 150
temperature = 1.0
# Interaction parameters (Lennard-Jones for raspberry)
#############################################################
radius_col = 3.
harmonic_radius = 3.0
# the subscript c is for colloid and s is for salt (also used for the surface beads)
eps_ss = 1. # LJ epsilon between the colloid's surface particles.
sig_ss = 1. # LJ sigma between the colloid's surface particles.
eps_cs = 48. # LJ epsilon between the colloid's central particle and surface particles.
sig_cs = radius_col # LJ sigma between the colloid's central particle and surface particles (colloid's radius).
a_eff = 0.32 # effective hydrodynamic radius of a bead due to the discreteness of LB.
# Ion types
#############################################################
TYPE_CENTRAL = 0
TYPE_SURFACE = 1
TYPE_CATIONS = 2
TYPE_ANIONS = 3
# System setup
#############################################################
system = espressomd.System(box_l=[box_l] * 3)
system.time_step = time_step
```

`box_l` sets the size of the simulation box. In general, one should check for finite
size effects which can be surprisingly large in simulations using hydrodynamic interactions. They
also generally scale as `box_l`$^{-1}$ or `box_l`$^{-3}$ depending on the transport mechanism
which sometimes allows for the infinite box limit to be extrapolated to, instead of using an
excessively large simulation box. As a rule of thumb, the box size should be five times greater than the characteristic
length scale of the object. Note that this example uses a small box
to provide a shorter simulation time.

In [2]:

```
system.cell_system.skin = skin
```

`skin`$\ \approx 0.3$.

In [3]:

```
system.periodicity = [True, True, True]
```

`periodicity` parameter indicates that the system is periodic in all three
dimensions. Note that the lattice-Boltzmann algorithm requires periodicity in all three directions (although
this can be modified using boundaries, a topic not covered in this tutorial).

Setting up the raspberry is a non-trivial task. The main problem lies in creating a relatively uniform distribution of beads on the surface of the colloid. In general one should take about 1 bead per lattice-Boltzmann grid point on the surface to ensure that there are no holes in the surface. The behavior of the colloid can be further improved by placing beads inside the colloid, though this is not done in this example script (for further reading see [1]). In our example we first define a harmonic interaction causing the surface beads to be attracted to the center, and a Lennard-Jones interaction preventing the beads from entering the colloid. There is also a Lennard-Jones potential between the surface beads to get them to distribute evenly on the surface.

In [4]:

```
# the LJ potential with the central bead keeps all the beads from simply collapsing into the center
system.non_bonded_inter[TYPE_SURFACE, TYPE_CENTRAL].wca.set_params(epsilon=eps_cs, sigma=sig_cs)
# the LJ potential (WCA potential) between surface beads causes them to be roughly equidistant on the
# colloid surface
system.non_bonded_inter[TYPE_SURFACE, TYPE_SURFACE].wca.set_params(epsilon=eps_ss, sigma=sig_ss)
# the harmonic potential pulls surface beads towards the central colloid bead
col_center_surface_bond = espressomd.interactions.HarmonicBond(k=3000., r_0=harmonic_radius)
system.bonded_inter.add(col_center_surface_bond)
```

Out[4]:

0

In [5]:

```
center = system.box_l / 2
colPos = center
# Charge of the colloid
q_col = -40
# Number of particles making up the raspberry (surface particles + the central particle).
n_col_part = int(4 * np.pi * np.power(radius_col, 2) + 1)
logging.info(f"Number of colloid beads = {n_col_part}")
# Place the central particle
central_part = system.part.add(pos=colPos, type=TYPE_CENTRAL, q=q_col,
fix=(True, True, True), rotation=(True, True, True))
```

INFO:root:Number of colloid beads = 114

Add `n_col_part-1`

particles of type `TYPE_SURFACE`

to the system and store the returned particle slice in `surface_parts`

(see the user guide section on how to add several particles at once). The particles shall be at random positions with a distance of exactly `radius_col`

from the colloid's center at `col_pos`

.

In [6]:

```
# Create surface beads uniformly distributed over the surface of the central particle
colSurfPos = np.random.uniform(low=-1, high=1, size=(n_col_part - 1, 3))
colSurfPos = colSurfPos / np.linalg.norm(colSurfPos, axis=1)[:, np.newaxis] * radius_col + colPos
colSurfTypes = np.full(n_col_part - 1, TYPE_SURFACE)
surface_parts = system.part.add(pos=colSurfPos, type=colSurfTypes)
```

In [7]:

```
# Relax bead positions. The LJ potential with the central bead combined with the
# harmonic bond keep the monomers roughly radius_col away from the central bead. The LJ
# between the surface beads cause them to distribute more or less evenly on the surface.
# We use a gradient descent algorithm with an additional constraint for surface particles.
system.integrator.set_steepest_descent(f_max=0, gamma=30, max_displacement=0.01 * sig_ss)
def constrain_surface_particles():
# This loop moves the surface beads such that they are once again exactly radius_col
# away from the center. For the scalar distance, we use system.distance() which
# considers periodic boundaries and the minimum image convention
colPos = central_part.pos
for p in surface_parts:
p.pos = (p.pos - colPos) / np.linalg.norm(system.distance(p, central_part)) * radius_col + colPos
p.pos = (p.pos - colPos) / np.linalg.norm(p.pos - colPos) * radius_col + colPos
logging.info("Relaxation of the raspberry surface particles")
for _ in range(100):
system.integrator.run(50)
constrain_surface_particles()
force_max = np.max(np.linalg.norm(system.part.all().f, axis=1))
logging.info(f"maximal force: {force_max:.1f}")
if force_max < 10.:
break
system.integrator.set_vv()
```

`VirtualSitesRelative`

scheme. This converts the raspberry to a rigid body
in which the surface particles follow the translation and rotation of the central particle,
and transfer the interpolated LB momentum back to the central particle.
Newton's equations of motion are only integrated for the central particle.
It is given an appropriate mass and moment of inertia tensor (note that the inertia tensor
is given in the frame in which it is diagonal.)

In [8]:

```
# Setting min_global_cut is necessary when there is no interaction defined with a range larger than
# the colloid such that the virtual particles are able to communicate their forces to the real particle
# at the center of the colloid
system.min_global_cut = radius_col
```

- Compute the center of mass of all particles in
`surface_parts`

and store its position in the variable`com`

. - Compute the moment of inertia of the particles in
`surface_parts`

using the previously computed`com`

, and store it in the variable`momI`

.*Hint:*we assume that the colloid is spherically symmetric. Therefore it suffices to compute a scalar value. Each particle has a mass of 1.

In [9]:

```
# Calculate the center of mass position (com) and the moment of inertia (momI) of the colloid
com = np.average(surface_parts.pos, 0) # surface_parts.pos returns an n-by-3 array
momI = 0
for p in surface_parts:
momI += np.power(np.linalg.norm(com - p.pos), 2)
```

In [10]:

```
# note that the real particle must be at the center of mass of the colloid because of the integrator
logging.info(f"Moving central particle from {central_part.pos} to {com}")
central_part.fix = [False, False, False]
central_part.pos = com
central_part.mass = n_col_part
central_part.rinertia = np.ones(3) * momI
# Convert the surface particles to virtual sites related to the central particle
# The id of the central particles is 0, the ids of the surface particles start at 1.
for p in surface_parts:
p.vs_auto_relate_to(central_part, couple_to_lb=True)
```

INFO:root:Moving central particle from [20. 20. 20.] to [19.98401984 19.88027911 19.97011589]

Next we insert enough ions at random positions (outside the radius of the colloid) with opposite charge to the colloid such that the system is electro-neutral. In addition, ions of both signs are added to represent the salt in the solution.

In [11]:

```
salt_rho = 0.001 # Number density of ions
volume = system.volume()
N_counter_ions = int(round(volume * salt_rho + abs(q_col)))
i = 0
while i < N_counter_ions:
pos = np.random.random(3) * system.box_l
# make sure the ion is placed outside of the colloid
if np.linalg.norm(pos - center) > radius_col + 1:
system.part.add(pos=pos, type=TYPE_CATIONS, q=1)
i += 1
logging.info(f"Added {N_counter_ions} positive ions")
N_co_ions = N_counter_ions - abs(q_col)
i = 0
while i < N_co_ions:
pos = np.random.random(3) * system.box_l
# make sure the ion is placed outside of the colloid
if np.linalg.norm(pos - center) > radius_col + 1:
system.part.add(pos=pos, type=TYPE_ANIONS, q=-1)
i += 1
logging.info(f"Added {N_co_ions} negative ions")
```

INFO:root:Added 104 positive ions INFO:root:Added 64 negative ions

We then check that charge neutrality is maintained

In [12]:

```
# Check charge neutrality
assert np.abs(np.sum(system.part.all().q)) < 1E-10
```

In [13]:

```
# WCA interactions for the ions, essentially giving them a finite volume
system.non_bonded_inter[TYPE_CENTRAL, TYPE_CATIONS].lennard_jones.set_params(
epsilon=eps_ss, sigma=sig_ss,
cutoff=sig_ss * pow(2., 1. / 6.), shift="auto", offset=sig_cs - 1 + a_eff)
system.non_bonded_inter[TYPE_CENTRAL, TYPE_ANIONS].lennard_jones.set_params(
epsilon=eps_ss, sigma=sig_ss,
cutoff=sig_ss * pow(2., 1. / 6.), shift="auto", offset=sig_cs - 1 + a_eff)
system.non_bonded_inter[TYPE_CATIONS, TYPE_CATIONS].wca.set_params(epsilon=eps_ss, sigma=sig_ss)
system.non_bonded_inter[TYPE_CATIONS, TYPE_ANIONS].wca.set_params(epsilon=eps_ss, sigma=sig_ss)
system.non_bonded_inter[TYPE_ANIONS, TYPE_ANIONS].wca.set_params(epsilon=eps_ss, sigma=sig_ss)
```

In [14]:

```
# Langevin thermostat for warmup before turning on the LB.
system.time_step = time_step / 10
system.thermostat.set_langevin(kT=temperature, gamma=1., seed=42)
logging.info("Removing overlap between ions")
ljcap = 100
for _ in range(100):
system.force_cap = ljcap
system.integrator.run(integ_steps)
ljcap += 5
system.force_cap = 0
system.time_step = time_step
```

INFO:root:Removing overlap between ions

Electrostatics are simulated using the Particle-Particle Particle-Mesh (P3M) algorithm. In **ESPResSo** this can be added to the simulation rather trivially:

In [15]:

```
# Turning on the electrostatics
# Note: Production runs would typically use a target accuracy of 10^-4
logging.info("Tuning P3M parameters...")
bjerrum = 2.
p3m = espressomd.electrostatics.P3M(prefactor=bjerrum * temperature, accuracy=0.001)
system.electrostatics.solver = p3m
logging.info("Tuning complete")
```

INFO:root:Tuning P3M parameters... INFO:root:Tuning complete CoulombP3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 2.00000e+00 System: box_l = 4.00000e+01 # charged part = 169 Sum[q_i^2] = 1.76800e+03 mesh cao r_cut_iL alpha_L err rs_err ks_err time [ms] 6 7 cao too large for this mesh 10 7 4.20356e-01 6.11568e+00 9.986e-04 7.071e-04 7.052e-04 1.28 10 6 4.39595e-01 5.83813e+00 9.970e-04 7.071e-04 7.028e-04 1.53 10 5 4.74224e-01 5.39622e+00 9.986e-04 7.071e-04 7.051e-04 1.48 10 4 4.92500e-01 5.18846e+00 1.593e-03 7.071e-04 1.427e-03 accuracy not achieved 14 7 3.10341e-01 8.37822e+00 9.938e-04 7.071e-04 6.983e-04 0.82 14 6 3.24298e-01 8.00459e+00 9.918e-04 7.071e-04 6.954e-04 0.79 14 5 3.50571e-01 7.38328e+00 9.955e-04 7.071e-04 7.008e-04 1.05 14 4 4.12967e-01 6.22929e+00 9.997e-04 7.071e-04 7.066e-04 1.21 14 3 4.20356e-01 6.11568e+00 3.051e-03 7.071e-04 2.968e-03 accuracy not achieved 18 7 2.47024e-01 1.06142e+01 9.873e-04 7.071e-04 6.891e-04 1.05 18 6 2.57792e-01 1.01550e+01 9.939e-04 7.071e-04 6.984e-04 0.99 18 5 2.79327e-01 9.34463e+00 9.987e-04 7.071e-04 7.053e-04 0.99 18 4 3.24298e-01 8.00459e+00 1.067e-03 7.071e-04 7.987e-04 accuracy not achieved 22 7 2.05131e-01 1.28680e+01 9.974e-04 7.071e-04 7.034e-04 1.40 22 6 2.14406e-01 1.22918e+01 9.989e-04 7.071e-04 7.056e-04 1.34 22 5 2.33500e-01 1.12519e+01 9.927e-04 7.071e-04 6.967e-04 1.33 22 4 2.78236e-01 9.38263e+00 9.964e-04 7.071e-04 7.020e-04 1.36 22 3 2.79327e-01 9.34463e+00 3.402e-03 7.071e-04 3.327e-03 accuracy not achieved 26 7 1.76037e-01 1.50767e+01 9.945e-04 7.071e-04 6.993e-04 1.94 26 6 1.84246e-01 1.43817e+01 9.937e-04 7.071e-04 6.982e-04 1.87 26 5 2.00664e-01 1.31648e+01 9.966e-04 7.071e-04 7.023e-04 1.85 26 4 2.33500e-01 1.12519e+01 1.088e-03 7.071e-04 8.266e-04 accuracy not achieved 30 7 1.54417e-01 1.72673e+01 9.911e-04 7.071e-04 6.945e-04 2.61 30 6 1.6147

WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:64) on node 0

2e-01 1.64870e+01 9.990e-04 7.071e-04 7.057e-04 3.17 30 5 1.76365e-01 1.50477e+01 9.971e-04 7.071e-04 7.031e-04 2.48 30 4 2.00664e-01 1.31648e+01 1.191e-03 7.071e-04 9.582e-04 accuracy not achieved 34 7 1.37785e-01 1.94291e+01 9.841e-04 7.071e-04 6.845e-04 6.24 34 6 1.43986e-01 1.85639e+01 9.996e-04 7.071e-04 7.065e-04 6.15 34 5 1.57764e-01 1.68883e+01 9.921e-04 7.071e-04 6.959e-04 6.12 34 4 1.76365e-01 1.50477e+01 1.282e-03 7.071e-04 1.069e-03 accuracy not achieved resulting parameters: mesh: (14, 14, 14), cao: 6, r_cut_iL: 3.2430e-01, alpha_L: 8.0046e+00, accuracy: 9.9179e-04, time: 0.79

Generally a Bjerrum length of $2$ is appropriate when using WCA interactions with $\sigma=1$, since a typical ion has a radius of $0.35\ \mathrm{nm}$, while the Bjerrum length in water is around $0.7\ \mathrm{nm}$.

The external electric field is simulated by simply adding a constant force equal to the simulated field times the particle charge. Generally the electric field is set to $0.1$ in MD units, which is the maximum field before the response becomes nonlinear. Smaller fields are also possible, but the required simulation time is considerably larger. Sometimes, Green-Kubo methods are also used, but these are generally only feasible in cases where there is either no salt or a very low salt concentration.

Add an external force $F=q\vec{E}$ to every particle according to its charge, where the electric field is $\vec{E}=\begin{pmatrix} 0.1 \\ 0 \\ 0 \end{pmatrix}$.

In [16]:

```
E = 0.1 # an electric field of 0.1 is the upper limit of the linear response regime for this model
Efield = np.array([E, 0, 0])
for p in system.part:
p.ext_force = p.q * Efield
```

Before creating the LB fluid it is a good idea to set all of the particle velocities to zero. This is necessary to set the total momentum of the system to zero. Failing to do so will lead to an unphysical drift of the system, which will change the values of the measured velocities.

In [17]:

```
system.part.all().v = (0, 0, 0)
```

In [18]:

```
lb = espressomd.lb.LBFluidWalberla(kT=temperature, seed=42,
density=1., kinematic_viscosity=3.,
agrid=1., tau=system.time_step)
```

A logical way of picking a specific set of parameters is to choose them such that the hydrodynamic radius of an ion roughly matches its physical radius determined by the WCA potential ($R=0.5\sigma$). Using the following equation:

\begin{equation} \frac{1}{\Gamma}=\frac{1}{6\pi \eta R_{\mathrm{H0}}}=\frac{1}{\Gamma_0} +\frac{1}{g\eta a} \end{equation}one can see that the set of parameters grid spacing $a=1\sigma$, fluid density $\rho=1$, a kinematic viscosity of $\nu=3 $ and a friction of $\Gamma_0=50$ leads to a hydrodynamic radius of approximately $0.5\sigma$.

The last step is to first turn off all other thermostats, followed by turning on the LB thermostat. The temperature is typically set to 1, which is equivalent to setting $k_\mathrm{B}T=1$ in molecular dynamics units.

In [19]:

```
system.thermostat.turn_off()
system.lb = lb
system.thermostat.set_lb(LB_fluid=lb, seed=123, gamma=20.0)
```

Now the main simulation can begin! The only important thing is to make sure the system has enough time to equilibrate. There are two separate equilibration times: 1) the time for the ion distribution to stabilize, and 2) the time needed for the fluid flow profile to equilibrate. In general, the ion distribution equilibrates fast, so the needed warmup time is largely determined by the fluid relaxation time, which can be calculated via $\tau_\mathrm{relax} = \mathrm{box\_length}^2/\nu$. This means for a box of size 40 with a kinematic viscosity of 3 as in our example script, the relaxation time is $\tau_\mathrm{relax} = 40^2/3 = 533 \tau_\mathrm{MD}$, or 53300 integration steps. In general it is a good idea to run for many relaxation times before starting to use the simulation results for averaging observables. To be on the safe side $10^6$ integration steps is a reasonable equilibration time. Please feel free to modify the provided script and try and get some interesting results!

In [20]:

```
# Reset the simulation clock
system.time = 0
initial_pos = central_part.pos
num_iterations = 20
num_steps_per_iteration = 20
with open('posVsTime.dat', 'w') as f: # file where the raspberry trajectory will be written to
for _ in tqdm.tqdm(range(num_iterations)):
system.integrator.run(num_steps_per_iteration)
pos = central_part.pos - initial_pos
f.write(f"{system.time:.2f} {pos[0]:.4f} {pos[1]:.4f} {pos[2]:.4f}\n")
logging.info("Finished")
```

5%|▌ | 1/20 [00:00<00:08, 2.22it/s]WARNING: Recalculating forces, so the LB coupling forces are not included in the particle force the first time step. This only matters if it happens frequently during sampling. in function void Thermostat::Thermostat::lb_coupling_deactivate() (/builds/espressomd/espresso/src/core/thermostat.cpp:91) on node 0 100%|██████████| 20/20 [00:08<00:00, 2.25it/s]

INFO:root:Finished

The above code cell saves the trajectory of the raspberry into the file `posVsTime.dat`

. For this purpose, the particle's `pos`

member should be used, as opposed to its `pos_folded`

member, which returns the particle's position folded back into the simulation box. In systems with periodic boundary conditions, particles can "leave" the simulation box. When a particle leaves the box, one of its periodic images enters the box from the opposite side, so it might appear that the particle never leaves the box. The truth, however, is that particles *can* leave the simulation box and therefore, their coordinates can end up outside of it. `pos`

returns these "true" coordinates. On the other hand, `pos_folded`

returns the position folded back into the simulation box, which is used to compute interactions between particles, and also for visualization of the simulation box.
Since the process of folding the particle position back into the simulation box destroys the information on how far it has actually travelled, `pos`

needs to be used to obtain a particle's trajectory.

Finally, we plot the raspberry trajectory with `matplotlib`:

In [21]:

```
from mpl_toolkits.mplot3d import Axes3D
trajectory_file = "posVsTime_sample.dat"
trajectory = np.loadtxt(trajectory_file)[:, 1:4]
# optional: trajectory smoothing with a running average
N = 6
trajectory = np.array(
[np.convolve(trajectory[:, i], np.ones((N,)) / N, mode='valid') for i in range(3)])
# calculate bounding box (cubic box to preserve scaling)
trajectory_range = np.max(trajectory, axis=1) - np.min(trajectory, axis=1)
mid_range = np.median(trajectory, axis=1)
max_range = 1.01 * np.max(np.abs(trajectory_range))
bbox = np.array([mid_range - max_range / 2, mid_range + max_range / 2])
# 3D plot
fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
ax.set_xlim(*bbox[:, 0])
ax.set_ylim(*bbox[:, 1])
ax.set_zlim(*bbox[:, 2])
ax.text(*trajectory[:, 0], '\u2190 start', 'y')
ax.scatter(*trajectory[:, 0])
ax.plot(*trajectory)
plt.tight_layout()
plt.rcParams.update({'font.size': 14})
```

[1] Lukas P. Fischer et al. The raspberry model for hydrodynamic interactions revisited. I. Periodic arrays of spheres and dumbbells, J. Chem. Phys. **143**, 084107 (2015)