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 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):
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", "w") 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.ConstantPressure(
filter=hoomd.filter.All(),
S=p,
tauS=tauP,
couple="xyz",
thermostat=hoomd.md.methods.thermostats.MTTK(kT=kT, tau=tau),
)
],
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]:
with warnings.catch_warnings():
# Hide deprecation warnings from hoomd
warnings.simplefilter("ignore")
project.run()
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"))
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.
[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