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