1.6 signac-flow HOOMD-blue Example
About
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.
Author
Carl Simon Adorf, Bradley Dice
Before you start
This example requires signac, signac-flow, HOOMD-blue and numpy. You can install these package for example via conda:
conda install -c conda-forge hoomd numpy signac signac-flow
[1]:
import flow
import hoomd
import numpy as np
import signac
# Initialize the HOOMD-blue execution context
hoomd.context.initialize("")
project_root = "projects/tutorial-signac-flow-hoomd"
class MyProject(flow.FlowProject):
pass
HOOMD-blue 2.9.6 DOUBLE HPMC_MIXED TBB SSE SSE2 SSE3
Compiled: 03/17/2021
Copyright (c) 2009-2019 The Regents of the University of Michigan.
-----
You are using HOOMD-blue. Please cite the following:
* J A Anderson, J Glaser, and S C Glotzer. "HOOMD-blue: A Python package for
high-performance molecular dynamics and hard particle Monte Carlo
simulations", Computational Materials Science 173 (2020) 109363
-----
HOOMD-blue is running on the CPU
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.
[2]:
from math import ceil
def init(N):
with hoomd.context.SimulationContext():
n = ceil(pow(N, 1 / 3))
assert n ** 3 == N
hoomd.init.create_lattice(unitcell=hoomd.lattice.sc(a=1.0), n=n)
hoomd.dump.gsd("init.gsd", period=None, group=hoomd.group.all())
def sample_lj(N, sigma, seed, kT, tau, p, tauP, steps, r_cut):
from hoomd import md
hoomd.init.read_gsd("init.gsd", restart="restart.gsd")
group = hoomd.group.all()
hoomd.dump.gsd("restart.gsd", truncate=True, period=100, phase=0, group=group)
lj = md.pair.lj(r_cut=r_cut, nlist=md.nlist.cell())
lj.pair_coeff.set("A", "A", epsilon=1.0, sigma=sigma)
md.integrate.mode_standard(dt=0.005)
md.integrate.npt(group=group, kT=kT, tau=tau, P=p, tauP=tauP)
hoomd.analyze.log("dump.log", ["volume"], 100, phase=0)
hoomd.run_upto(steps)
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.
[3]:
@MyProject.label
def estimated(job):
return "V" in job.document
@MyProject.label
def sampled(job):
return job.document.get("sample_step", 0) >= 5000
@MyProject.operation
@MyProject.post(estimated)
def estimate(job):
sp = job.statepoint()
job.document["V"] = sp["N"] * sp["kT"] / sp["p"]
@MyProject.operation
@MyProject.post(sampled)
def sample(job):
import hoomd
with job:
with hoomd.context.SimulationContext():
try:
sample_lj(steps=5000, **job.statepoint())
finally:
job.document["sample_step"] = hoomd.get_step()
Now it’s time to actually generate some data! Let’s initialize the data space!
[4]:
project = MyProject.init_project(name="FlowTutorialHOOMDProject", root=project_root)
# 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=p, tau=1.0, tauP=1.0, r_cut=2.5)
job = project.open_job(sp).init()
with job:
init(N=sp["N"])
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
notice(2): Group "all" created containing 512 particles
The print_status()
function allows to get a quick overview of our project’s status:
[5]:
project.print_status(detailed=True, parameters=["p"])
Overview: 10 jobs, 10 jobs with eligible operations.
label
-------
operation number of eligible jobs submission status
----------- ------------------------- -------------------
estimate 10 [U]: 10
sample 10 [U]: 10
Detailed View:
job id operation 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
The next cell will attempt to execute all eligible operations.
[6]:
project.run()
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:04 | Step 5000 / 5000 | TPS 1227.09 | ETA 00:00:00
Average TPS: 1226.85
---------
-- Neighborlist stats:
316 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 104 / n_neigh_avg: 51.2246
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 60 / n_max: 68 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1523.76 | ETA 00:00:00
Average TPS: 1523.43
---------
-- Neighborlist stats:
287 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 30 / n_neigh_avg: 10.3984
shortest rebuild period: 3
-- Cell list stats:
Dimension: 4, 4, 4
n_min : 1 / n_max: 16 / n_avg: 8
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1366.9 | ETA 00:00:00
Average TPS: 1366.02
---------
-- Neighborlist stats:
305 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 93 / n_neigh_avg: 43.7891
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 60 / n_max: 67 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1318.05 | ETA 00:00:00
Average TPS: 1317.6
---------
-- Neighborlist stats:
315 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 100 / n_neigh_avg: 45.3223
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 62 / n_max: 67 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1253.11 | ETA 00:00:00
Average TPS: 1252.47
---------
-- Neighborlist stats:
309 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 108 / n_neigh_avg: 52.6016
shortest rebuild period: 4
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 62 / n_max: 69 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1428.05 | ETA 00:00:00
Average TPS: 1427.45
---------
-- Neighborlist stats:
301 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 66 / n_neigh_avg: 30.8477
shortest rebuild period: 3
-- Cell list stats:
Dimension: 3, 3, 3
n_min : 11 / n_max: 25 / n_avg: 18.963
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:04 | Step 5000 / 5000 | TPS 1247.66 | ETA 00:00:00
Average TPS: 1247.02
---------
-- Neighborlist stats:
322 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 107 / n_neigh_avg: 52.8262
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 59 / n_max: 68 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1280.47 | ETA 00:00:00
Average TPS: 1280.24
---------
-- Neighborlist stats:
306 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 114 / n_neigh_avg: 53.2871
shortest rebuild period: 4
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 58 / n_max: 68 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1311.12 | ETA 00:00:00
Average TPS: 1310.35
---------
-- Neighborlist stats:
312 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 106 / n_neigh_avg: 49.9785
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 60 / n_max: 66 / n_avg: 64
** run complete **
notice(2): Group "all" created containing 512 particles
notice(2): -- Neighborlist exclusion statistics -- :
notice(2): Particles with 0 exclusions : 512
notice(2): Neighbors included by diameter : no
notice(2): Neighbors excluded when in the same body: no
** starting run **
Time 00:00:03 | Step 5000 / 5000 | TPS 1367.57 | ETA 00:00:00
Average TPS: 1367.29
---------
-- Neighborlist stats:
306 normal updates / 50 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 95 / n_neigh_avg: 47.293
shortest rebuild period: 3
-- Cell list stats:
Dimension: 2, 2, 2
n_min : 57 / n_max: 73 / n_avg: 64
** run complete **
Let’s double check the project status.
[7]:
project.print_status(detailed=True)
Overview: 10 jobs, 0 jobs with eligible operations.
label ratio
--------- ----------------------------------------------------------
estimated |████████████████████████████████████████| 10/10 (100.00%)
sampled |████████████████████████████████████████| 10/10 (100.00%)
operation
-----------
Detailed View:
job id operation 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
After running all operations we can make a brief examination of the collected data.
[8]:
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 629.75796286 146.28571428571428
0.5 1146.859093668 1024.0
1.5 789.514068364 341.3333333333333
2.0 723.460173872 256.0
4.0 614.63846208 128.0
1.0 905.471509388 512.0
4.5 599.651218704 113.77777777777777
5.0 589.6517064440001 102.4
3.0 657.398353356 170.66666666666666
2.5 692.0376509839999 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.
[9]:
%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$]")
plt.legend()
plt.show()

Uncomment and execute the following line to remove all data and start over.
[10]:
# %rm -r projects/tutorial-signac-flow-hoomd-blue/workspace