1.6 signac-flow HOOMD-blue Example


This notebook contains a minimal example for running a signac-flow project from scratch. The example demonstrates how to compare an ideal gas with a Lennard-Jones (LJ) fluid by calculating a p-V phase diagram.

This examples uses the general-purpose simulation toolkit HOOMD-blue for the execution of the molecular-dynamics (MD) simulations.


Carl Simon Adorf, Bradley Dice

Before you start

This example requires signac, signac-flow, HOOMD-blue, gsd, and numpy. You can install these package for example via conda:

conda install -c conda-forge gsd hoomd numpy signac signac-flow
import itertools
import os
import warnings

import flow
import gsd.hoomd
import hoomd
import numpy as np
import signac

project_path = "projects/tutorial-signac-flow-hoomd"

class MyProject(flow.FlowProject):

We want to generate a pressure-volume phase diagram for a Lennard-Jones fluid with molecular dynamics (MD) using the general-purpose simulation toolkit HOOMD-blue (http://glotzerlab.engin.umich.edu/hoomd-blue/).

We start by defining two functions, one for the initialization of our simulation and one for the actual execution.

from contextlib import contextmanager
from math import ceil

def init(N):
    frame = gsd.hoomd.Frame()
    frame.particles.N = N
    frame.particles.types = ["A"]
    frame.particles.typeid = np.zeros(N)

    n = ceil(pow(N, 1 / 3))
    assert n**3 == N
    spacing = 1.2
    L = n * spacing
    frame.configuration.box = [L, L, L, 0.0, 0.0, 0.0]
    x = np.linspace(-L / 2, L / 2, n, endpoint=False)
    position = list(itertools.product(x, repeat=3))
    frame.particles.position = position
    with gsd.hoomd.open("init.gsd", "w") as traj:

def get_sim(N, sigma, seed, kT, tau, p, tauP, r_cut):
    sim = hoomd.Simulation(hoomd.device.CPU())
    state_fn = "init.gsd"
    if os.path.exists("restart.gsd"):
        state_fn = "restart.gsd"
    sim.operations += hoomd.write.GSD(
        filename="restart.gsd", truncate=True, trigger=100, filter=hoomd.filter.All()
    lj = hoomd.md.pair.LJ(default_r_cut=r_cut, nlist=hoomd.md.nlist.Cell(0.4))
    lj.params[("A", "A")] = {"epsilon": 1.0, "sigma": sigma}
    integrator = hoomd.md.Integrator(
                thermostat=hoomd.md.methods.thermostats.MTTK(kT=kT, tau=tau),
    sim.operations.integrator = integrator
    logger = hoomd.logging.Logger(categories=("scalar", "string"))
    logger["Volume"] = (lambda: sim.state.box.volume, "scalar")

    # This context manager keeps the log file open while the simulation runs.
    with open("dump.log", "w+") as outfile:
        table = hoomd.write.Table(
            trigger=100, logger=logger, output=outfile, pretty=False
        sim.operations += table
        yield sim

We want to use signac to manage our simulation data and signac-flow to define a workflow acting on the data space.

Now that we have defined the core simulation logic above, it’s time to embed those functions into a general workflow. For this purpose we add labels and operations with preconditions/postconditions to MyProject, a subclass of flow.FlowProject.

The estimate operation stores an ideal gas estimation of the volume for the given system. The sample operation actually executes the MD simulation as defined in the previous cell.

def estimated(job):
    return "V" in job.document

def sampled(job):
    return job.document.get("sample_step", 0) >= 5000

def estimate(job):
    sp = job.statepoint()
    job.document["V"] = sp["N"] * sp["kT"] / sp["p"]

def sample(job):
    with job:
        with get_sim(**job.statepoint()) as sim:
                job.document["sample_step"] = sim.timestep

Now it’s time to actually generate some data! Let’s initialize the data space!

project = MyProject.init_project(project_path)

# Uncomment the following two lines if you want to start over!
# for job in project:
#     job.remove()

for p in np.linspace(0.5, 5.0, 10):
    sp = dict(
        N=512, sigma=1.0, seed=42, kT=1.0, p=float(p), tau=1.0, tauP=1.0, r_cut=2.5
    job = project.open_job(sp).init()
    with job:

The print_status() function allows to get a quick overview of our project’s status:

project.print_status(detailed=True, parameters=["p"])

Overview: 10 jobs/aggregates, 10 jobs/aggregates with eligible operations.


operation/group      number of eligible jobs  submission status
-----------------  -------------------------  -------------------
estimate                                  10  [U]: 10
sample                                    10  [U]: 10

Detailed View:

job id                            operation/group      p  labels
--------------------------------  -----------------  ---  --------
18c063206c37ab03d06de130c6ae2b70  estimate [U]       1
                                  sample [U]         1
2a05a487f254db90b08f576d199a83b7  estimate [U]       5
                                  sample [U]         5
4ef567c185064a4b27b73f62124991c0  estimate [U]       3
                                  sample [U]         3
5b3c8f0fbd9458a304222d83bb82bc30  estimate [U]       3.5
                                  sample [U]         3.5
6681f083e2b03f5b9816f35502e411f7  estimate [U]       1.5
                                  sample [U]         1.5
6ff6b952b42f0ad3e1315e137cf67cdb  estimate [U]       4
                                  sample [U]         4
99f5397cf2f6250afe34927606aa2f37  estimate [U]       4.5
                                  sample [U]         4.5
a5d4427e8d285c78239005f740a2625b  estimate [U]       2
                                  sample [U]         2
a89bc9f5cceaf268f18efed936a70333  estimate [U]       0.5
                                  sample [U]         0.5
c342fff540bfa0a2128efac700fc5aef  estimate [U]       2.5
                                  sample [U]         2.5

[U]:unknown [R]:registered [I]:inactive [S]:submitted [H]:held [Q]:queued [A]:active [E]:error [GR]:group_registered [GI]:group_inactive [GS]:group_submitted [GH]:group_held [GQ]:group_queued [GA]:group_active [GE]:group_error

The next cell will attempt to execute all eligible operations.

with warnings.catch_warnings():
    # Hide deprecation warnings from hoomd

Let’s double check the project status.


Overview: 10 jobs/aggregates, 0 jobs/aggregates with eligible operations.

label      ratio
---------  ----------------------------------------------------------
estimated  |████████████████████████████████████████| 10/10 (100.00%)
sampled    |████████████████████████████████████████| 10/10 (100.00%)


Detailed View:

job id                            operation/group    labels
--------------------------------  -----------------  ------------------
18c063206c37ab03d06de130c6ae2b70  [ ]                estimated, sampled
2a05a487f254db90b08f576d199a83b7  [ ]                estimated, sampled
4ef567c185064a4b27b73f62124991c0  [ ]                estimated, sampled
5b3c8f0fbd9458a304222d83bb82bc30  [ ]                estimated, sampled
6681f083e2b03f5b9816f35502e411f7  [ ]                estimated, sampled
6ff6b952b42f0ad3e1315e137cf67cdb  [ ]                estimated, sampled
99f5397cf2f6250afe34927606aa2f37  [ ]                estimated, sampled
a5d4427e8d285c78239005f740a2625b  [ ]                estimated, sampled
a89bc9f5cceaf268f18efed936a70333  [ ]                estimated, sampled
c342fff540bfa0a2128efac700fc5aef  [ ]                estimated, sampled

[U]:unknown [R]:registered [I]:inactive [S]:submitted [H]:held [Q]:queued [A]:active [E]:error [GR]:group_registered [GI]:group_inactive [GS]:group_submitted [GH]:group_held [GQ]:group_queued [GA]:group_active [GE]:group_error

After running all operations we can make a brief examination of the collected data.

def get_volume(job):
    log = np.genfromtxt(job.fn("dump.log"), names=True)
    N = len(log)
    return log[int(0.5 * N) :]["Volume"].mean(axis=0)

for job in project:
    print(job.statepoint["p"], get_volume(job), job.document.get("V"))
3.5 633.09071616 146.28571428571428
0.5 1238.75831864 1024.0
1.5 789.8835313200001 341.3333333333333
2.0 732.7477800800001 256.0
4.0 616.6703303999999 128.0
1.0 931.94641252 512.0
4.5 603.5544676 113.77777777777777
5.0 592.8023312 102.4
3.0 663.6705322 170.66666666666666
2.5 676.3909288000001 204.8

For a better presentation of the results we need to aggregate all results and sort them by pressure.

The following code requires matplotlib.

%matplotlib inline
# Display plots within the notebook

from matplotlib import pyplot as plt

V = dict()
V_idg = dict()

for job in project:
    V[job.statepoint()["p"]] = get_volume(job)
    V_idg[job.statepoint()["p"]] = job.document["V"]

p = sorted(V.keys())
V = [V[p_] for p_ in p]
V_idg = [V_idg[p_] for p_ in p]

plt.plot(p, V, label="LJ")
plt.plot(p, V_idg, label="idG")
plt.xlabel(r"pressure [$\epsilon / \sigma^3$]")
plt.ylabel(r"volume [$\sigma^3$]")

Uncomment and execute the following line to remove all data and start over.

# %rm -r projects/tutorial-signac-flow-hoomd-blue/workspace