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:
import espressomd
import espressomd.interactions
import espressomd.electrostatics
import espressomd.lb
import espressomd.zn
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
The parameter 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.
system.cell_system.skin = skin
The skin is used for constructing the Verlet lists and is purely an optimization parameter. Whatever value provides the fastest integration speed should be used. For the type of simulations covered in this tutorial, this value turns out to be skin$\ \approx 0.3$.
system.periodicity = [True, True, True]
The 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.
# 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)
0
We set up the central bead and the other beads are initialized at random positions on the surface of the colloid. The beads are then allowed to relax using an integration loop where the forces between the beads are capped.
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.
# SOLUTION CELL
# 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)
# 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()
INFO:root:Relaxation of the raspberry surface particles INFO:root:maximal force: 1640.1 INFO:root:maximal force: 11.8 INFO:root:maximal force: 4.1
The best way to ensure a relatively uniform distribution of the beads on the surface is to simply take a look at a VMD snapshot of the system after this integration. Such a snapshot is shown in Fig. 1.
Now that the beads are arranged in the shape of a raspberry, the surface beads need to be set up as virtual particles
using the 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.)
# 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
surface_parts and store its position in the variable com.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.# SOLUTION CELL
# 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)
# 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 [20.17393276 19.99958468 19.87291289]
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.
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
# Check charge neutrality
assert np.abs(np.sum(system.part.all().q)) < 1E-10
A WCA potential acts between all of the ions. This potential represents a purely repulsive version of the Lennard-Jones potential, which approximates hard spheres of diameter $\sigma$. The ions also interact through a WCA potential with the central bead of the colloid, using an offset of around $\mathrm{radius\_col}-\sigma +a_{\mathrm{grid}}/2$. This makes the colloid appear as a hard sphere of radius roughly $\mathrm{radius\_col}+a_{\mathrm{grid}}/2$ to the ions, which is approximately equal to the hydrodynamic radius of the colloid
# 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)
After inserting the ions, again a short integration is performed with a force cap to prevent strong overlaps between the ions.
# 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:
# 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
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}$.
# SOLUTION CELL
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.
system.part.all().v = (0, 0, 0)
The important parameters for the LB fluid are the density, the viscosity, the time step, and the friction coefficient used to couple the particle motion to the fluid. The time step should generally be comparable to the MD time step. While large time steps are possible, a time step of $0.01$ turns out to provide more reasonable values for the root mean squared particle velocities. Both density and viscosity should be around $1$, while the friction should be set around $20$. The grid spacing should be comparable to the ions' size.
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.
system.thermostat.turn_off()
system.lb = lb
system.thermostat.set_lb(LB_fluid=lb, seed=123, gamma=20.0)
We use the ZnDraw visualizer and run our simulation for a few steps in order to get an idea what is happening in our system:
# Initializing ZnDraw visualizer
colors = {
TYPE_CENTRAL: "#fca000", #orange
TYPE_SURFACE: "#7fc454", #green
TYPE_CATIONS: "#ff0000", #red
TYPE_ANIONS: "#030ffc" #blue
}
radii = {
TYPE_CENTRAL: sig_ss/2,
TYPE_SURFACE: sig_ss/2,
TYPE_CATIONS: sig_ss/2,
TYPE_ANIONS: sig_ss/2
}
arrows_config = {'colormap': [[-0.5, 0.9, 0.5], [-0.0, 0.9, 0.5]],
'normalize': True,
'colorrange': [0, 1],
'scale_vector_thickness': True,
'opacity': 1.0}
lbfield = espressomd.zn.LBField(system, step_x=2, step_y=2, step_z=5, scale=0.3)
vis = espressomd.zn.Visualizer(system, colors=colors, radii=radii, folded=True, vector_field=lbfield)
for _ in range(50):
system.integrator.run(10)
vis.update()
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 0.66
10 6 4.39595e-01 5.83813e+00 9.970e-04 7.071e-04 7.028e-04 0.63
10 5 4.74224e-01 5.39622e+00 9.986e-04 7.071e-04 7.051e-04 0.65
10 4 4.92500e-01 5.18846e+00 1.593e-03 7.071e-04 1.427e-03 accuracy not achieved
14 7 3.10807e-01 8.36520e+00 9.860e-04 7.071e-04 6.872e-04 0.40
14 6 3.24545e-01 7.99829e+00 9.882e-04 7.071e-04 6.904e-04 0.39
14 5 3.50302e-01 7.38915e+00 9.985e-04 7.071e-04 7.050e-04 0.52
14 4 4.12979e-01 6.22911e+00 9.996e-04 7.071e-04 7.065e-04 0.58
14 3 4.39595e-01 5.83813e+00 2.567e-03 7.071e-04 2.467e-03 accuracy not achieved
18 7 2.46578e-01 1.06341e+01 9.970e-04 7.071e-04 7.029e-04 0.46
18 6 2.57988e-01 1.01470e+01 9.902e-04 7.071e-04 6.932e-04 0.44
18 5 2.79539e-01 9.33729e+00 9.958e-04 7.071e-04 7.011e-04 0.44
18 4 3.24545e-01 7.99829e+00 1.064e-03 7.071e-04 7.952e-04 accuracy not achieved
22 7 2.05080e-01 1.28713e+01 9.988e-04 7.071e-04 7.054e-04 0.65
22 6 2.14654e-01 1.22771e+01 9.933e-04 7.071e-04 6.977e-04 0.61
22 5 2.33297e-01 1.12620e+01 9.960e-04 7.071e-04 7.014e-04 0.60
22 4 2.57988e-01 1.01470e+01 1.292e-03 7.071e-04 1.081e-03 accuracy not achieved
26 7 1.75884e-01 1.50903e+01 9.993e-04 7.071e-04 7.061e-04 0.92
26 6 1.84086e-01 1.43946e+01 9.979e-04 7.071e-04 7.042e-04 0.86
26 5 2.00490e-01 1.31766e+01 1.000e-03 7.071e-04 7.071e-04 0.88
26 4 2.33297e-01 1.12620e+01 1.091e-03 7.071e-04 8.306e-04 accuracy not achieved
30 7 1.54604e-01 1.72458e+01 9.846e-04 7.071e-04 6.852e-04 1.34
30 6 1.61794e-01 1.64530e+01 9.894e-04 7.071e-04 6.920e-04 1.24
30 5 1.76895e-01 1.50010e+01 9.858e-04 7.071e-04 6.869e-04 1.21
30 4 1.84086e-01 1.43946e+01 1.728e-03 7.071e-04 1.577e-03 accuracy not achieved
34 7 1.37508e-01 1.94696e+01 9.951e-04 7.071e-04 7.001e-04 3.58
34 6 1.44418e-01 1.85063e+01 9.852e-04 7.071e-04 6.860e-04 3.53
34 5 1.57547e-01 1.69123e+01 9.973e-04 7.071e-04 7.033e-04 3.54
34 4 1.76895e-01 1.50010e+01 1.267e-03 7.071e-04 1.051e-03 accuracy not achieved
resulting parameters: mesh: (14, 14, 14), cao: 6, r_cut_iL: 3.2454e-01,
alpha_L: 7.9983e+00, accuracy: 9.8823e-04, time: 0.39
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
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!
# 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)
vis.update()
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:03, 5.90it/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:03<00:00, 5.62it/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:
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)