1.7 signac-flow minimal example (with HOOMD-blue)

About

Note: This example is in the process of being updated for the latest versions of signac and signac-flow.

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

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"
project = signac.init_project(name="FlowTutorialHOOMDProject", root=project_root)
HOOMD-blue 2.9.3 DOUBLE HPMC_MIXED TBB SSE SSE2 SSE3
Compiled: 10/09/2020
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 bundle all operations modifying our data space in a Operations class. This will make it easier later to execute these operations by name.

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. We define all operations as static methods of an Operations class in order to keep them in a separate namespace.

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


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()
[4]:
for job in project:
    print(job.document)

Now that we have defined all operations it’s time to embed them into a general workflow. For this purpose we specialize a flow.FlowProject.

The workflow is controlled by two core functions: classify() and next_operation(): - The classify() function allows us to label our jobs and get a good overview of the project status. This is especially important once the data space becomes larger and more complex and operations more expensive. - The next_operation() functions helps to automate the workflow by identifying the next required operation for each job.

[5]:
class MyProject(flow.FlowProject):
    def classify(self, job):
        yield "init"
        if "V" in job.document:
            yield "estimated"
        if job.document.get("sample_step", 0) >= 5000:
            yield "sampled"

    def next_operation(self, job):
        labels = set(self.classify(job))
        if "V" not in job.document:
            return "estimate"
        if "sampled" not in labels:
            return "sample"

We need to use the get_project() class method to get a project handle for this special project class.

[6]:
project = MyProject.get_project(root=project_root)

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

[7]:
# 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)
    with project.open_job(sp):
        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:

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



Overview:

Total # of jobs: 10

label
-------


operation
-----------


Detailed View:

job id                            operation      p  labels
--------------------------------  -----------  ---  --------
18c063206c37ab03d06de130c6ae2b70  [ ]          1
a89bc9f5cceaf268f18efed936a70333  [ ]          0.5
a5d4427e8d285c78239005f740a2625b  [ ]          2
99f5397cf2f6250afe34927606aa2f37  [ ]          4.5
5b3c8f0fbd9458a304222d83bb82bc30  [ ]          3.5
2a05a487f254db90b08f576d199a83b7  [ ]          5
6ff6b952b42f0ad3e1315e137cf67cdb  [ ]          4
4ef567c185064a4b27b73f62124991c0  [ ]          3
c342fff540bfa0a2128efac700fc5aef  [ ]          2.5
6681f083e2b03f5b9816f35502e411f7  [ ]          1.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 operations by cycling through jobs and operations until no next operations are defined anymore.

We limit the max. number of cycles to prevent accidental infinite loops, the number of cycles is arbitrary.

[9]:
for i in range(3):
    for job in project:
        for j in range(5):
            next_op = project.next_operation(job)
            if next_op is None:
                break
            print("execute", job, next_op)
            globals()[next_op](job)
            assert next_op != project.next_operation(job)
        else:
            raise RuntimeError("Reached max. # cycle limit!")
execute 5b3c8f0fbd9458a304222d83bb82bc30 estimate
execute 5b3c8f0fbd9458a304222d83bb82bc30 sample
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 1504.51 | ETA 00:00:00
Average TPS: 1503.67
---------
-- 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 **
execute a89bc9f5cceaf268f18efed936a70333 estimate
execute a89bc9f5cceaf268f18efed936a70333 sample
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:02 | Step 5000 / 5000 | TPS 1839.88 | ETA 00:00:00
Average TPS: 1839.17
---------
-- 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 **
execute 6681f083e2b03f5b9816f35502e411f7 estimate
execute 6681f083e2b03f5b9816f35502e411f7 sample
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 1663.52 | ETA 00:00:00
Average TPS: 1663.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 **
execute a5d4427e8d285c78239005f740a2625b estimate
execute a5d4427e8d285c78239005f740a2625b sample
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 1595.41 | ETA 00:00:00
Average TPS: 1594.94
---------
-- 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 **
execute 6ff6b952b42f0ad3e1315e137cf67cdb estimate
execute 6ff6b952b42f0ad3e1315e137cf67cdb sample
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 1489.42 | ETA 00:00:00
Average TPS: 1489.08
---------
-- 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 **
execute 18c063206c37ab03d06de130c6ae2b70 estimate
execute 18c063206c37ab03d06de130c6ae2b70 sample
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:02 | Step 5000 / 5000 | TPS 1733.42 | ETA 00:00:00
Average TPS: 1732.88
---------
-- 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 **
execute 99f5397cf2f6250afe34927606aa2f37 estimate
execute 99f5397cf2f6250afe34927606aa2f37 sample
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 1459.47 | ETA 00:00:00
Average TPS: 1459.14
---------
-- 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 **
execute 2a05a487f254db90b08f576d199a83b7 estimate
execute 2a05a487f254db90b08f576d199a83b7 sample
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 1466.86 | ETA 00:00:00
Average TPS: 1466.45
---------
-- 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 **
execute 4ef567c185064a4b27b73f62124991c0 estimate
execute 4ef567c185064a4b27b73f62124991c0 sample
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 1534.76 | ETA 00:00:00
Average TPS: 1533.67
---------
-- 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 **
execute c342fff540bfa0a2128efac700fc5aef estimate
execute c342fff540bfa0a2128efac700fc5aef sample
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 1568.46 | ETA 00:00:00
Average TPS: 1568.09
---------
-- 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.

[10]:
project.print_status();



Overview:

Total # of jobs: 10

label
-------


operation
-----------



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

[11]:
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.

[12]:
%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()
../../_images/examples_notebooks_signac_107_HOOMD-blue_22_0.png

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

[13]:
# %rm -r projects/tutorial-signac-flow-hoomd-blue/workspace
[ ]: