University of Michigan Logo
Quick Start Tutorial

Example script

So you have HOOMD-blue installed. Now what!?
Let's start with the classic MD simulation, the Lennard-Jones liquid. Place N particles randomly in a box and allow them to interact with the following potential between pairs of particles:

\[ V(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right] \]

To configure HOOMD-blue to perform this simulation, a simple Python script must be written.

from hoomd_script import *
# create 100 random particles of name A
init.create_random(N=100, phi_p=0.01, name='A')
# specify Lennard-Jones interactions between particle pairs
lj = pair.lj(r_cut=3.0)
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
# integrate at constant temperature
all = group.all();
integrate.mode_standard(dt=0.005)
integrate.nvt(group=all, T=1.2, tau=0.5)
# run 10,000 time steps
run(10e3)

If you don't know Python, don't worry. You can learn everything about it that you need to know for HOOMD-blue scripts here. Of course if you do know Python (which is a full-fledged object oriented programming language: tutorials at http://www.python.org) then you could make use of its capabilities in setting up complicated simulations.


Running the example

For now, copy and paste the above code into a file test.hoomd. Assuming you have installed HOOMD-blue, you can run the simulation script from the command line:

$ hoomd test.hoomd

And you should see output that looks something like this:

HOOMD-blue 0.9.0
Compiled: Wed Oct 28 06:58:46 EDT 2009
Copyright 2008, 2009 Ames Laboratory Iowa State University and the Regents of the University of Michigan
-----
http://codeblue.umich.edu/hoomd-blue/
This code is the implementation of the algorithms discussed in:
Joshua A. Anderson, Chris D. Lorenz, and Alex Travesset - 'General
Purpose Molecular Dynamics Fully Implemented on Graphics Processing
Units', Journal of Computational Physics 227 (2008) 5342-5359
-----
test.hoomd:004 | init.create_random(N=100, phi_p=0.01, name='A')
test.hoomd:007 | lj = pair.lj(r_cut=3.0)
test.hoomd:008 | lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
test.hoomd:011 | all = group.all();
Group "all" created containing 100 particles
test.hoomd:012 | integrate.mode_standard(dt=0.005)
test.hoomd:013 | integrate.nvt(group=all, T=1.2, tau=0.5)
test.hoomd:016 | run(10e3)
** starting run **
Time 00:00:00 | Step 10000 / 10000 | TPS 35417.9 | ETA 00:00:00
Average TPS: 35405
---------
-- Neighborlist stats:
370 normal updates / 100 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 10 / n_neigh_avg: 2.41
bins_min: 0 / bins_max: 6 / bins_avg: 1.5625
** run complete **

That's it! You've just run your first simulation with HOOMD-blue.


Understanding the output

The first few lines of output are just a header notifying you which version of HOOMD-blue you are running and when it was compiled along with a link to the website and the reference to the paper discussing the algorithms used in HOOMD-blue.

The simulation output starts at:

test.hoomd:004 | init.create_random(N=100, phi_p=0.01, name='A')

Each hoomd_script command prints the file and line where it was run along with the entire text of the command. This can be potentially very useful in debugging problems, as it will allow you to hone in on the command in the script that is producing the error message.

When a run() command is executed, hoomd steps the simulation forward that many time steps. While it is doing so, it periodically prints out status lines like the one seen above.

Time 00:00:00 | Step 10000 / 10000 | TPS 41990.9 | ETA 00:00:00
  • Time is the total time spent (so far) in the current simulation in HH:MM:SS (totaled over multiple run() commands)
  • Step current / final prints the current time step the simulation is at and the final time step of the run()
  • TPS is the current rate (in Time steps Per Second) at which the simulation is progressing.
  • ETA is the estimated time to completion of the current run() in HH:MM:SS

Since this run was so short, only one line was printed at the end. Modify the script to run a few million time steps and run it to see what happens in longer simulations. Or just see here...

Time 00:00:10 | Step 488501 / 10000000 | TPS 48850.1 | ETA 00:03:14
Time 00:00:20 | Step 976312 / 10000000 | TPS 48781.1 | ETA 00:03:04
Time 00:00:30 | Step 1462718 / 10000000 | TPS 48640.5 | ETA 00:02:55
Time 00:00:40 | Step 1950647 / 10000000 | TPS 48792.8 | ETA 00:02:44
Time 00:00:50 | Step 2436905 / 10000000 | TPS 48625.4 | ETA 00:02:35
Time 00:01:00 | Step 2924701 / 10000000 | TPS 48779.5 | ETA 00:02:25
Time 00:01:10 | Step 3410821 / 10000000 | TPS 48612 | ETA 00:02:15

The final bit of output at the end of the run prints statistics from various parts of the computation. In this example, only the neighbor list prints statistics.

---------
-- Neighborlist stats:
370 normal updates / 100 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 10 / n_neigh_avg: 2.41
bins_min: 0 / bins_max: 6 / bins_avg: 1.5625

Differently configured simulation scripts may print additional information here.


Line by line details

  1. The first line of every hoomd job script (except for comment lines starting with #) must be
    from hoomd_script import *
    This line takes the python code of hoomd_script, compiles it and loads it in so it can be used. hoomd_script contains the code for commands such as init.create_random() which is why this line must be first.

  2. After hoomd_script has been imported, the system must be initialized before any other command can be executed. In this example, we create a 100 random particles named A.
    init.create_random(N=100, phi_p=0.01, name='A')
    Here is a good point to call attention to one of HOOMD-blue's nifty features. You can name a particle type anything you want. If you want to name a particular particle type 'My ridiculously long particle type name', be my guest. hoomd doesn't care one way or another as it just stores the string you give it. Just be warned that some software packages only handle names up to a certain length, so some information may be lost when they load mol2 (dump.mol2) or xml (dump.xml) files written by HOOMD-blue.
    Documentation for init.create_random

  3. The next line specifies the pair force between particle pairs in the simulation. In the example, we create a Lennard-Jones pair force with a cutoff radius of 3.0.
    lj = pair.lj(r_cut=3.0)
    This line has the structure variable=command which saves the result of the command for later modification (see why on the next line of the script). In HOOMD-blue, any number of forces can be specified, even zero if that is what you need (use a separate line and variable name for each one). All specified forces are added together during the simulation.
    Documentation for pair.lj

  4. The parameters of the force must also be specified before the simulation is run() (you will get an error if you forget). The next line sets the parameters epsilon and sigma for the Lennard-Jones force between particles of types A and A.
    lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
    The example script only has a single particle type in the simulation so one line is sufficient. In more complicated simulations with more than one particle type every unique pair must be specified (i.e. 'A','A', 'A','B', and 'B',B'). Use one line like that above for each pair, or specify a list of types to assign the given parameters (see pair_coeff.set() for examples).

  5. After that, we choose the integrator to move particles forward in time. As with forces, more than one integration method can be applied, as long as they are set on different groups. However, only one integration mode can be active at a time (specifying another will overwrite the first). Also, there must be an integration mode set. If you try running the simulation without one hoomd will print an error message. Here, we specify a standard integration mode with a time step dt set to 0.005. The Nosé-Hoover thermostat is applied to a group containing all particles, with the temperature target T set at 1.2 and the parameter tau at 0.5. Like with the pair force above, the group is saved in a variable to be used a few lines later in the script.
    all = group.all();
    integrate.mode_standard(dt=0.005)
    integrate.nvt(group=all, T=1.2, tau=0.5)
    After the initialization and before the run() command, the order in which commands are called doesn't matter. The NVT integrator could just as easily been specified before the pair force.
    Documentation for group
    Documentation for integrate
    Documentation for integrate.mode_standard
    Documentation for integrate.nvt

  6. Finally, the run command actually takes the job settings previously specified and runs the simulation through time.
    run(10e3)
    This simple example only runs for 10,000 steps but real simulations might be run for 10's of millions, spending days of computation time in this single command. There is no limit that there be a single run command in a given script. If your simulation needs to turn off a certain force or change integrators and then continue, you can do that. Just execute the commands to make the changes after the first run and before the second.
    Documentation for run()

Components of a hoomd_script command

The init line is as good an example as any.

init.create_random(N=100, phi_p=0.01, name='A')

Parts of the command

  • init - The package. Every command is in its own package to keep it organized
  • . - Python syntax needed to access a member of the package
  • create_random - The command name to run
  • ( - Python required syntax to note the start of an argument list
  • N=100, phi_p=0.01, name='A' - Arguments (more on these below)
  • ) - Python required syntax to note the end of an argument list
  • <enter> - Python required syntax to execute the command

About the arguments
Multiple arguments of the form name=value are separated by commas. Whitespace is ignored so name = value works too. The order of arguments doesn't matter (as long as you specify them by name). I.e. all of the following are identical:

init.create_random(N=100, name='A', phi_p=0.01)
init.create_random(phi_p = 0.01, N = 100, name = 'A')
init.create_random(phi_p=0.01, name='A', N=100)

Check the documentation for a specific command to see what the arguments are and what they mean (i.e. see init.create_random). Here is a copy of the documentation for init.create_random:

init.create_random
( N,
phi_p,
name = "A",
min_dist = 1.0
)
Parameters:
N Number of particles to create
phi_p Packing fraction of particles in the simulation box
name Name of the particle type to create
min_dist Minimum distance particles will be separated by

First of all, notice in the header for the command that some arguments are listed with an = sign, like name="A". This means that argument has a default value associated with it. If you are happy with the default value, you don't need to specify that argument in the list. In our example, we have always been setting name='A' anyways so using

init.create_random(N=100, phi_p=0.01)

is identical.

Those arguments that are listed without the = sign (N, phi_P here) have no default value and must be specified. Python will give you an error if you don't.

init.create_random() doesn't have any optional arguments, but some commands do. Optional arguments will be labeled as such and have a default value of None. The documentation for any optional arguments will clearly indicate what occurs when you do or do not specify it in the argument list.


Where to go from here

This quick tutorial is only the tip of the iceberg. There are a lot more resources in the documentation.