Analyzing simulation data from HOOMD-blue at runtimeΒΆ

The following script shows how to use freud to compute the radial distribution function \(g(r)\) on data generated by the molecular dynamics simulation engine HOOMD-blue during a simulation run.

Generally, most users will want to run analyses as post-processing steps, on the saved frames of a particle trajectory file. However, it is possible to use analysis callbacks in HOOMD-blue to compute and log quantities at runtime, too. By using analysis methods at runtime, it is possible to stop a simulation early or change the simulation parameters dynamically according to the analysis results.

HOOMD-blue can be installed with conda install -c conda-forge hoomd.

The simulation script runs a Monte Carlo simulation of spheres, with outputs parsed with numpy.genfromtxt.

%matplotlib inline
import itertools

import freud
import h5py
import hoomd
import matplotlib.pyplot as plt
import numpy as np
from hoomd import hpmc
def make_simple_cubic_snapshot(a, n):
    """Make a snapshot with a simple cubic lattice.

        a (float): Lattice spacing
        n (int): Number of particles

        hoomd.Snapshot: The initial system snapshot.
    k = int(np.ceil(n ** (1 / 3)))
    L = k * a
    x = np.linspace(-L / 2, L / 2, k, endpoint=False)
    position = list(itertools.product(x, repeat=3))
    position = position[:n]

    snap = hoomd.Snapshot()
    snap.particles.N = n
    snap.particles.types = ["A"]
    snap.particles.typeid[:] = [0] * n
    snap.particles.position[:] = position = [L, L, L, 0, 0, 0]

    return snap

class FreudRDFCalc(hoomd.custom.Action):
    """Compute RDF periodically as the simulation progresses."""

    def __init__(self, sim, rdf):
        self._sim = sim
        self._rdf = rdf

    def act(self, timestep):
        snap = self._sim.state.get_snapshot()
        self._rdf.compute(system=snap, reset=False)

class FreudW6Calc(hoomd.custom.Action):
    """Compute w6 for logging during a simulation."""

    def __init__(self, sim, w6):
        self._sim = sim
        self._w6 = w6

    def w6(self):
        """Calculate w6 order parameter."""
        snap = self._sim.state.get_snapshot()
        self._w6.compute(snap, neighbors=dict(num_neighbors=12))
        return np.mean(w6.particle_order)

    def act(self, timestep):
# set up simulation
sim = hoomd.Simulation(device=hoomd.device.CPU())
sim.create_state_from_snapshot(make_simple_cubic_snapshot(a=1.0, n=1000))
sim.seed = 42
mc = hpmc.integrate.Sphere(default_d=0.1, default_a=0.1)
mc.shape["A"] = dict(diameter=0.5)
sim.operations.integrator = mc

# freud-related classes
rdf = freud.density.RDF(bins=50, r_max=4)
w6 = freud.order.Steinhardt(l=6, wl=True)
rdf_calc = FreudRDFCalc(sim, rdf)
w6_calc = FreudW6Calc(sim, w6)

# Equilibrate the system a bit before accumulating the RDF.

# add logger/w6/rdf calculations to the simulation
logger = hoomd.logging.Logger(hoomd.write.HDF5Log.accepted_categories)
logger.add(w6_calc, quantities=["w6"])
logger.add(sim, quantities=["timestep"])
hdf5_writer = hoomd.write.HDF5Log(
    trigger=100, filename="w6_data.h5", mode="w", logger=logger

rdf_action = hoomd.write.CustomWriter(action=rdf_calc, trigger=100)

# Store the computed RDF in a file
    "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)"
rdf_data = np.genfromtxt("rdf.csv", delimiter=",")
plt.plot(rdf_data[:, 0], rdf_data[:, 1])
plt.title("Radial Distribution Function")
w6_data = h5py.File("w6_data.h5", "r")
plt.title("$w_6$ Order Parameter")