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, gsd, and numpy. You can install these package for example via conda:
conda install -c conda-forge gsd hoomd numpy signac signac-flow
[1]:
import itertools
import os
import flow
import gsd.hoomd
import hoomd
import numpy as np
import signac
project_path = "projects/tutorial-signac-flow-hoomd"
class MyProject(flow.FlowProject):
pass
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 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", "wb") as traj:
traj.append(frame)
@contextmanager
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.create_state_from_gsd(state_fn)
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(
dt=0.005,
methods=[
hoomd.md.methods.NPT(
filter=hoomd.filter.All(), kT=kT, tau=tau, S=p, tauS=tauP, couple="xyz"
)
],
forces=[lj],
)
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.
[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.post(estimated)
@MyProject.operation
def estimate(job):
sp = job.statepoint()
job.document["V"] = sp["N"] * sp["kT"] / sp["p"]
@MyProject.post(sampled)
@MyProject.operation
def sample(job):
with job:
with get_sim(**job.statepoint()) as sim:
try:
sim.run(5_000)
finally:
job.document["sample_step"] = sim.timestep
Now it’s time to actually generate some data! Let’s initialize the data space!
[4]:
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:
init(N=sp["N"])
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/aggregates, 10 jobs/aggregates with eligible operations.
label
-------
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.
[6]:
project.run()
/home/bdice/mambaforge/envs/signac-examples/lib/python3.10/site-packages/hoomd/md/methods/methods.py:365: FutureWarning: NPT is deprecated and wil be removed in hoomd 4.0. In version 4.0, use the ConstantPressure method with the desired thermostat from hoomd.md.methods.thermostats.
warnings.warn(
Let’s double check the project status.
[7]:
project.print_status(detailed=True)
Overview: 10 jobs/aggregates, 0 jobs/aggregates with eligible operations.
label ratio
--------- ----------------------------------------------------------
estimated |████████████████████████████████████████| 10/10 (100.00%)
sampled |████████████████████████████████████████| 10/10 (100.00%)
operation/group
-----------------
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.
[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"))
2.0 735.545045936 256.0
1.5 798.7716894079999 341.3333333333333
5.0 588.32751952 102.4
3.0 662.30139688 170.66666666666666
4.5 602.9928609199999 113.77777777777777
4.0 620.5540934320001 128.0
2.5 697.688737468 204.8
3.5 638.340131088 146.28571428571428
0.5 1168.8474468079999 1024.0
1.0 915.2272913159999 512.0
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