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.
import espressomd
import espressomd.lb
import espressomd.shapes
import espressomd.observables
import espressomd.accumulators
espressomd.assert_features(["LENNARD_JONES", "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:
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)# 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)
In this scenario, we will sample the sedimentation dynamics using the Langevin thermostat.
Exercise:
# 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, 381.95it/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()
In this scenario, we want to sample the same system coupled to a lattice-Boltzmann fluid.
Exercise:
spacing
and LB time step equal to the MD time step and add this object to the system
(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:
# 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:
acc_lb_vel
(user guide)Hints:
n_height
by 2 units (parameters min_y
and max_y
are also affected)n_z_bins=1
to average the velocity along the z-directionsampling_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)
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]
0%| | 0/700 [00:00<?, ?it/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%|██████████| 700/700 [00:16<00:00, 41.85it/s]
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)