Running a Simulation¶
This page contains a sample python script to run a simulation using the FSD plugin. If the script is called run.py
, the simulation can be run with the command python3 run.py
.
A sample of the script run.py
that shears a simple cubic lattice is:
# Import modules
import hoomd
from hoomd import _hoomd
from hoomd.md import _md
import hoomd.PSEv3
import os
import math
import numpy as np
# Initialize HOOMD
hoomd.context.initialize()
# Time step
dt = 1e-3
# Number of time steps for simulation
nsteps = 1000
# Temperature and size
radius = 1.0
diameter = 2.0 * radius
T = 1.0
phi = 0.2
# Random seed for stochastic calculations
seed = np.random.randint( 1, 10000 )
# Simple cubic crystal
# n = number of particles in each direction
# N = total number of particles
# L = box size
# a = spacing between particles
n = 10;
N = int(n*n*n)
L = ( 4.0*math.pi/3.0 * N * radius**3.0 / phi )**(1.0/3.0)
a = L / n
# Create the box and particles
system = hoomd.init.create_lattice(
unitcell=hoomd.lattice.sc(a=a),
n=int(n)
)
# Set up sinusoidal shearing function
shear_function = hoomd.PSEv3.shear_function.sine(
dt = dt,
shear_rate = 1.0,
shear_freq = 0.2
)
# To shear, need to resize the box for Lees-Edwards BCs
# !! It is very important to turn "scale_particles" off, or else HOOMD
# will shift particle coordinates by deformation, resulting in wrong
# dynamics
box_resize = hoomd.update.box_resize(
xy = hoomd.PSEv3.variant.shear_variant(
function_form = shear_function,
total_timestep = nsteps
),
scale_particles = False
)
# Set up pse integrator
#
# The integrator takes the following parameters, default values are given in parentheses.
#
# group = group of particles for integrator to work on
# seed = seed for stochastic sampling (0)
# xi = Ewald splitting parameter (0.5)
# error = global relative error tolerance (0.001)
# function_form = function form for shearing (None)
# max_strain = maximum total strain for shearing before box reset (0.5)
# fileprefix = name for stresslet output file (stresslet)
# period = frequency to output stresslet file. 0 = no output. (0)
#
# The only parameter that needs to be set is "group". For the rest, the default values
# can be assumed if there's no need to change them.
itg = hoomd.md.integrate.mode_standard(dt=dt)
pse = hoomd.PSEv3.integrate.PSEv3(
group = hoomd.group.all(),
seed = seed,
T = T,
function_form = shear_function
)
# Run the simulation
hoomd.run( nsteps )
Sheared Simulations¶
- Simple shear is implemented in the FSD plugin using the Lees-Edwards boundary conditions, which translates the periodic images above and below the sample relative to one another. The flow, gradient, and vorticity directions are the cartesian directions x, y, and z, respectively. This is accomplished in HOOMD-blue using the built-in box resizer
box_resize = hoomd.update.box_resize( xy = hoomd.PSEv3.variant.shear_variant( function_form = shear_function, total_timestep = nsteps ), scale_particles = False )
In this command, the first argument sets a variant resizer for the box deformation using a specific shear function and the second argument turns off HOOMD’s automatic particle rescaling. If rescaling is not turned off, the particle coordinates are updated proportional to the box deformation automatically by HOOMD, which produces incorrect particle trajectories, because particle advection by shear is included in the hydrodynamic integration performed by FSD.
- There are three options for the type of time-dependent shear signal used in a simulation (a sample call to set up the shear function for each type is shown as well):
- Steady, constant strain rate –
shear_function = hoomd.PSEv3.shear_function.steady( dt, shear_rate )
- Sinusoidal, oscillatory strain rate with constant frequency –
shear_function = hoomd.PSEv3.shear_function.sine( dt = dt, shear_rate, shear_freq )
- Chirp, oscillatory strain rate with exponentially increasing frequency –
shear_function = hoomd.PSEv3.shear_function.chirp( dt, amp, omega_0, omega_f, periodT )
- Steady, constant strain rate –
- The chirp signal (which has an exponentially increasing frequency), is most efficiently analyzed using a Tukey window to restrict the amplitude. The plugin also supports a windowed shear function, which is implemented as
shear_function = hoomd.PSEv3.shear_function.chirp( dt, amp, omega_0, omega_f, periodT ) shear_window = hoomd.PSEv3.shear_function.tukey_window( dt, periotT, tukey_param ) windowed_shear_function = hoomd.PSEv3.shear_function.windowed( function_form = shear_function, window = shear_window )
where windowed_shear_function
can be used as shear_function
in e.g. the box resizer command.
In most cases, when the box resizer resets the box to orthorhombic (without rescaling particles) after a total relative strain of 0.5. This ensures that the deformation doesn’t get so large so as to introduce uncontrolled numerical errors.