Molecular Dynamics Simulations¶
Overview¶
Questions¶
What is molecular dynamics?
How do I set up a molecular dynamics simulation in HOOMD-blue?
Objectives¶
Describe the equations of motion of the system and HOOMD-blue solves them with time steps.
Define forces, potential energy and explain how HOOMD-blue evaluates pair potentials within a cutoff.
Explain how the MD Integrator and integration methods solve the equations of motion and allow for different thermodynamic ensembles.
Boilerplate Code¶
[1]:
import hoomd
import matplotlib
import numpy
%matplotlib inline
matplotlib.style.use("ggplot")
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats("svg")
Equations of Motion¶
Molecular dynamics simulations model the movement of particles over time by solving the equations of motion numerically, advancing the state of the system forward by time dt
on each time step. You can use molecular dynamics to model dynamic, time dependent processes (like fluid flow) or thermodynamic equilibrium states (like crystals). This tutorial models a system of Lennard-Jones particles crystallizing.
The Integrator class in the md package implements molecular dynamics integration in HOOMD-blue:
[2]:
integrator = hoomd.md.Integrator(dt=0.005)
The dt
property sets the step size \(\Delta t\):
[3]:
integrator.dt
[3]:
0.005
Forces¶
The force term in the equations of motion determines how particles interact. The forces
Integrator property is the list of forces applied to the system. The default is no forces:
[4]:
integrator.forces[:]
[4]:
[]
You can add any number of Force objects to this list. Each will compute forces on all particles in the system state. Integrator will add all of the forces together into the net force used in the equations of motion.
HOOMD-blue provides a number of forces that are derived from potential energy. Pair potentials define the functional form of the potential energy between a single pair of particles given their separation distance r. The Lennard-Jones potential is:
[5]:
sigma = 1
epsilon = 1
r = numpy.linspace(0.95, 3, 500)
V_lj = 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
fig = matplotlib.figure.Figure(figsize=(5, 3.09))
ax = fig.add_subplot()
ax.plot(r, V_lj)
ax.set_xlabel("r")
ax.set_ylabel("V")
fig
[5]:
While pair potentials are nominally defined between all pairs of particles, Molecular dynamics simulations evaluate short ranged pair potentials only for \(r \lt r_{\mathrm{cut}}\) to make the computation fast through the use of a neighbor list. By default, HOOMD-blue introduces a discontinuity in V at \(r=r_{\mathrm{cut}}\), though there are options to shift or smooth the potential at the cutoff.
HOOMD-blue provides several neighbor list options to choose from. Cell performs well in most situations. The buffer
parameter sets an extra distance to include in the neighbor list, so it can be used for multiple timesteps without recomputing neighbors. Choose a value that minimizes the time to execute your simulation.
[6]:
cell = hoomd.md.nlist.Cell(buffer=0.4)
Construct the LJ pair force object to compute the Lennard-Jones interactions:
[7]:
lj = hoomd.md.pair.LJ(nlist=cell)
Pair potentials accept parameters for every pair of particle types in the simulation. Define a single A-A interaction for this tutorial:
[8]:
lj.params[("A", "A")] = dict(epsilon=1, sigma=1)
lj.r_cut[("A", "A")] = 2.5
Add the force to the Integrator:
[9]:
integrator.forces.append(lj)
Now, the Integrator will compute the net force term using the Lennard-Jones pair force.
Integration Methods¶
HOOMD-blue provides a number of integration methods, which define the equations of motion that apply to a subset of the particles in the system. The ConstantVolume method implements Newton’s laws while the thermostat scales the velocities to sample the canonical ensemble.
Lennard-Jones particles crystallize readily at constant temperature and volume. One of the points in the solid area of the phase diagram is \(kT=1.5\) at a number density \(\rho=1.2\) Apply the ConstantVolume method with the Bussi thermostat to all particles (later sections in this tutorial will set the density):
[10]:
nvt = hoomd.md.methods.ConstantVolume(
filter=hoomd.filter.All(), thermostat=hoomd.md.methods.thermostats.Bussi(kT=1.5)
)
kT
is the temperature multiplied by Boltzmann’s constant and has units of energy. tau
is a time constant that controls the amount of coupling between the thermostat and particle’s degrees of freedom. filter
is a particle filter object that selects which particles this integration method applies to. You can apply multiple integration methods to different sets of particles or leave some particles fixed.
Add the method to the Integrator:
[11]:
integrator.methods.append(nvt)
Now, you have defined a molecular dynamics integrator with a step size dt
that will use the ConstantVolume integration method with the Bussi thermostat on all particles interacting via the Lennard-Jones potential.
The remaining sections in this tutorial will show you how to initialize a random low density state, compress the system to a target density, and run the simulation long enough to self-assemble an ordered state.