Sedimentation in a fluid

The purpose of this tutorial is to demonstrate how hydrodynamic interactions can have a dramatic impact on the overall dynamics of a molecular system. We will set up a simple semi-two-dimensional system of sedimenting particles and simulate their interaction with a solvent.

In the first scenario, Langevin dynamics (LD) will model the friction from a "simple" implicit solvent. Although LD in general doesn't fully model an implicit solvent due to the lack of contribution for electrostatic screening and the hydrophobic effect, these limitations won't affect this particular simulation, since we represent the sediment particles as electrically neutral and fully solvated beads.

LD adds two extra terms to Newton's second law, which account for the friction and random collisions from a solvent:

\begin{equation} m_i \dot{v}_i(t) = f_i(t) - \gamma v_i(t) + \sqrt{2\gamma k_{\mathrm{B}} T} \cdot \eta_i(t) \end{equation}

with $m_i$, $v_i(t)$, $f_i(t)$ the mass, velocity and forces acting on particle $i$ at time $t$, $\gamma$ the dampening constant, $k_{\mathrm{B}}$ the Boltzmann constant, $T$ the temperature, $\eta_i(t)$ a random uniform number drawn for particle $i$ at time $t$. The second term in the right-hand side is responsible for drag due to friction with the solvent, while the third term is responsible for heat exchange due to collisions with the solvent. This method doesn't capture hydrodynamic effects, since the particles only interact with each other via their inter-atomic potentials; when two particles are further away than their potential cutoff, they cannot interact.

In the second scenario, we will use the lattice-Boltzmann (LB) method to model the solvent using a grid-based solver for the Navier-Stokes equation in the limit of a low Reynolds number. Particle–fluid interaction is achieved through momentum exchange via frictional coupling:

\begin{equation} m_i \dot{v}_i(t) = f_i(t) - \gamma \left( v_i(t) - u(x_i(t), t) \right) + \sqrt{2\gamma k_{\mathrm{B}} T} \cdot \eta_i(t) \end{equation}

with $u(x_i(t), t)$ the interpolated fluid velocity at the position of the particle. The frictional term is also back-propagated to the nearest nodes of the LB fluid. While computationally more demanding, this method recovers hydrodynamic effects, i.e. particles in close proximity influence each other due to coupling with the fluid inbetween them, even when their separation distance exceeds the inter-atomic potential cutoff. In addition, the fluid can develop a macroscopic flow field, in which case all particles will influence each other over large distances.

We will compare both scenarios by generating a video of the sedimentation dynamics. To make the effect more clearly visible, we will not couple the particles and fluid to a thermal bath.

System setup

import espressomd
import espressomd.lb
import espressomd.shapes
import espressomd.observables
import espressomd.accumulators

espressomd.assert_features(["LENNARD_JONES", "EXTERNAL_FORCES", "WALBERLA"])

# imports for data handling, plotting, and progress bar
import tqdm
import numpy as np
import matplotlib.pyplot as plt

The initial particles positions will be chosen to form a two-dimensional hexagonal Bravais lattice structure. The spacing we use is large enough for the particles not to interact initially. If particles were positioned closer, they would rearrange to minimize the interaction energy, breaking up the initial structure.

We setup the simulation box such that the particle structure is periodically repeated in $x$-direction.

def hexagonal_lattice(n_rows, parts_per_row, spacing):
    positions = np.array([[x + 0.12 + (0.5 if y % 2 == 1 else 0),
                           y * np.sqrt(3) / 2, 0]
                          for x in range(parts_per_row)
                          for y in range(n_rows)]) * spacing
    return positions

Now, we are ready to define the system parameters and initialize the simulation system.

# WCA potential parameters
lj_sigma = 1.
lj_epsilon = 1.
lj_cutoff = 2**(1. / 6.) * lj_sigma

# lattice spacing and number of lattice rows to use
spacing = lj_cutoff
n_rows = 10

# system size in units of lattice spacing
n_height = 50
n_width = 20
n_depth = 2

# resulting box geometry
box_height = n_height * spacing
box_width = n_width * spacing
box_depth = n_depth * spacing

# system setup
system = espressomd.System(box_l=[box_width, box_height, box_depth])
system.time_step = 0.01
system.cell_system.skin = 0.4

# add non-bonded WCA interaction
system.non_bonded_inter[0, 0].lennard_jones.set_params(
    epsilon=lj_epsilon, sigma=lj_sigma, cutoff=lj_cutoff, shift="auto")

We add a wall constraint on bottom and top of the simulation box, respectively.

Exercise:

• set up two wall constraints wall_shape_b (bottom) and wall_shape_t (top) in the xz-plane at a distance of 1 unit from the top and bottom sides of the box (user guide)
• add these two shapes to the system constraints using particle type 0 (user guide)

# SOLUTION CELL
# create wall shapes bottom (b) and top (t)
wall_shape_b = espressomd.shapes.Wall(normal=[0, 1, 0], dist=1)
wall_shape_t = espressomd.shapes.Wall(
    normal=[0, -1, 0], dist=-(box_height - 1))

# add wall constraints
for wall_shape in [wall_shape_b, wall_shape_t]:
    system.constraints.add(shape=wall_shape, particle_type=0)

We will now calculate the particle initial positions and introduce a small crystalline defect to help break the symmetry of the system. We will also configure an external force acting on all particles, which models the effect of gravity.

# gravitational force
f_gravity = [0., -3., 0.]

# number of frames in the output trajectory
sampling_steps = 700

# lattice positions
initial_positions = hexagonal_lattice(n_rows, n_width, spacing)

# introduce a small imperfection in the initial lattice structure
initial_positions[99, 0] += 0.15 * spacing
initial_positions[109, 0] -= 0.15 * spacing

# shift initial positions to the top
y_max = np.amax(initial_positions[:, 1])
initial_positions += (box_height - 1 - lj_cutoff) - y_max
initial_positions = np.remainder(initial_positions, np.copy(system.box_l))

# total number of particles
n_parts = n_rows * n_width
assert initial_positions.shape == (n_parts, 3)

Langevin dynamics

In this scenario, we will sample the sedimentation dynamics using the Langevin thermostat.

Exercise:

• set up an unthermalized Langevin thermostat to add a purely frictional term to the equations of motion with $\gamma=15$ (user guide)

# SOLUTION CELL
system.thermostat.set_langevin(kT=0., gamma=15., seed=12)

We can now sample the particle positions as a function of time. We will use a particle observable and an accumulator to record the trajectory.

parts = system.part.add(pos=initial_positions, ext_force=[f_gravity] * n_parts)
obs_particle_pos = espressomd.observables.ParticlePositions(ids=list(range(n_parts)))
acc_particle_pos = espressomd.accumulators.TimeSeries(obs=obs_particle_pos, delta_N=1)

system.integrator.run(0)

for step in tqdm.tqdm(range(sampling_steps)):
    acc_particle_pos.update()
    system.integrator.run(25)

data_ld = np.remainder(np.reshape(acc_particle_pos.time_series(), (-1, n_parts, 3)), system.box_l)

100%|██████████| 700/700 [00:01<00:00, 457.03it/s]

We will now disable the thermostat, reset the particles to their initial positions and zero out particle velocities.

parts.pos = initial_positions
parts.v = [0, 0, 0]
system.thermostat.turn_off()

Hydrodynamics

In this scenario, we want to sample the same system coupled to a lattice-Boltzmann fluid.

Exercise:

• create an unthermalized lattice-Boltzmann object with viscosity 1, density 1, grid size equal to spacing and LB time step equal to the MD time step and add this object to the system (user guide)
• activate particle coupling to the fluid by setting the LB thermostat with $\gamma=15$ (user guide)

# SOLUTION CELL
lbf = espressomd.lb.LBFluidWalberla(agrid=spacing,
                                    density=1.,
                                    kinematic_viscosity=1.,
                                    tau=system.time_step, kT=0.)
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, gamma=15., seed=0)

The wall constraints that were previously added now have to be registered as LB boundaries.

Exercise:

• convert the wall shapes to LB boundaries and add them to the system list of LB boundaries (user guide)

# SOLUTION CELL
# add LB boundaries
for wall_shape in [wall_shape_b, wall_shape_t]:
    lbf.add_boundary_from_shape(wall_shape)

We will plot the fluid flow field in the final video using 2D vectors. To this end, we need to record the fluid trajectory with the same frequency as the particle positions.

Exercise:

• create a LB velocity profile in Cartesian coordinates (user guide)
• register that observable in a time series accumulator named acc_lb_vel (user guide)

Hints:

• the velocity observable takes parameters in MD units, not LB units
• there is no fluid inside the top an bottom boundaries, therefore the number of bins for the y-axis is smaller than n_height by 2 units (parameters min_y and max_y are also affected)
• use n_z_bins=1 to average the velocity along the z-direction
• for sampling_delta_* and sampling_offset_*, use spacing and 0.5 * spacing respectively

# SOLUTION CELL
obs_lb_vel = espressomd.observables.LBVelocityProfile(
    n_x_bins=n_width,
    n_y_bins=n_height - 2,  # skip data inside the LB boundaries (top and bottom walls)
    n_z_bins=1,             # averaged velocity along the z-direction
    min_x=0.0,
    min_y=spacing,
    min_z=0.0,
    max_x=system.box_l[0],
    max_y=system.box_l[1] - spacing,
    max_z=system.box_l[2],
    sampling_delta_x=spacing,
    sampling_delta_y=spacing,
    sampling_delta_z=spacing,
    sampling_offset_x=0.5 * spacing,
    sampling_offset_y=0.5 * spacing,
    sampling_offset_z=0.5 * spacing,
    allow_empty_bins=True)
acc_lb_vel = espressomd.accumulators.TimeSeries(obs=obs_lb_vel, delta_N=1)

We can now sample the particle positions and fluid velocity as a function of time.

obs_particle_pos = espressomd.observables.ParticlePositions(ids=list(range(n_parts)))
acc_particle_pos = espressomd.accumulators.TimeSeries(obs=obs_particle_pos, delta_N=1)

system.integrator.run(0)

for step in tqdm.tqdm(range(sampling_steps)):
    acc_lb_vel.update()
    acc_particle_pos.update()
    system.integrator.run(25)

data_lb = np.remainder(np.reshape(acc_particle_pos.time_series(), (-1, n_parts, 3)), system.box_l)
data_flowfield = acc_lb_vel.time_series()[:, :, :, 0, 0:2]

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:94) on node 0
100%|██████████| 700/700 [00:24<00:00, 28.98it/s]

Visualization

Now let's visualize the data. First some imports and definitions for inline visualization.

import matplotlib.animation as animation
import matplotlib.quiver
import tempfile
import base64

VIDEO_TAG = """<video controls>
 <source src="data:video/x-m4v;base64,{0}" type="video/mp4">
 Your browser does not support the video tag.
</video>"""

# set ignore 'divide' and 'invalid' errors
# these occur when plotting the flowfield containing a zero velocity
np.seterr(divide='ignore', invalid='ignore')

def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with tempfile.NamedTemporaryFile(suffix='.mp4') as f:
            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
            with open(f.name, "rb") as g:
                video = g.read()
        anim._encoded_video = base64.b64encode(video).decode('ascii')
        plt.close(anim._fig)
    return VIDEO_TAG.format(anim._encoded_video)

animation.Animation._repr_html_ = anim_to_html

And now the actual visualization code.

Note: if Jupyter encapsulates the video in a scrollable area, click on the Out[]: text to the left of the output cell to toggle auto-scrolling off. This can also be achieved by hitting the key combination Shift+O after highlighting the output cell.

# setup figure and prepare axes
fig = plt.figure(figsize=(2 * 5, 5 / box_width * box_height))
gs = fig.add_gridspec(1, 2, wspace=0.1)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1], sharey=ax1)

ax1.set_title("Langevin")
ax1.set_xlim((0, box_width))
ax1.set_ylim((0, box_height))

ax2.set_title("LB Fluid")
ax2.set_xlim((0, box_width))
ax2.set_ylim((0, box_height))

# draw walls
for ax in [ax1, ax2]:
    ax.hlines((1, box_height-1), 0, box_width, color="gray")

# create meshgrid for quiver plot
xs = np.array([x for x in range(n_width)]) * spacing
ys = np.array([y for y in range(1, n_height-1)]) * spacing
X, Y = np.meshgrid(xs, ys)

# create a transposed flow field for quiver plot
data_flowfield_t = np.transpose(data_flowfield, axes=(0, 2, 1, 3))

# set quiver scale (fraction of the highest velocity in the XY plane)
lb_vel_max = np.sum(np.square(data_flowfield), axis=-1)
quiver_scale = np.sqrt(np.max(lb_vel_max))

def plot_lb_vel(ax, X, Y, flowfield, t, scale):
    return ax.quiver(X, Y,
                     flowfield[t, :, :, 0],
                     flowfield[t, :, :, 1],
                     scale_units="xy", scale=scale)

# initialize plot objects
lb_ff = plot_lb_vel(ax2, X, Y, data_flowfield_t, 0, quiver_scale)
lb_particles, = ax2.plot([], [], 'o')
ld_particles, = ax1.plot([], [], 'o')

def draw_frame(t):
    # manually remove Quivers from ax2
    for artist in ax2.get_children():
        if isinstance(artist, matplotlib.quiver.Quiver):
            artist.remove()

    # draw new quivers
    lb_ff = plot_lb_vel(ax2, X, Y, data_flowfield_t, t, quiver_scale)

    # draw particles
    ld_particles.set_data(data_ld[t, :, 0], data_ld[t, :, 1])
    lb_particles.set_data(data_lb[t, :, 0], data_lb[t, :, 1])

    return [ld_particles, lb_particles, lb_ff]

animation.FuncAnimation(fig, draw_frame, frames=sampling_steps, blit=True, interval=0, repeat=False)