Example Scripts

Examples:

• Example 1: Using dump files
  The simulation trajectory is written to a file which can be viewed in visualization software such as VMD
• Example 2: Using IMD
  A running simulation is connected to VMD which displays the current system state updated in real-time
• Example 3: Using the polymer generator
  The polymer generator can generate complicated initial conditions for bead-spring polymer systems
• Example 4: Using arbitrary input files
  None of the built-in generators works for you? Write an input file describing the initial condition and load it in.
• Example 5: Using a python loop to run many simulations at different parameters
  One simulation in a script not enough? Use a simple python loop to run many in a row!

  ---

Example 1: Using dump files

This simple example is a simple adaptation of the quick start script. It performs a simulation of a Lennard-Jones liquid, dumping snapshots of the system every 100 time steps.

from hoomd_script import *

# create 1000 random particles of name A
init.create_random(N=1000, 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, alpha=1.0)

# integrate at constant temperature
integrate.nvt(dt=0.005, T=1.2, tau=0.5)

# dump a .mol2 file for the structure information
mol2 = dump.mol2()
mol2.write(filename='example1.mol2')

# dump a .dcd file for the trajectory
dump.dcd(filename='example1.dcd', period=100)

# run 10,000 time steps
run(10e3)

Running this quick simulation should result in two output files being generated in the current working directory: example.mol2 and example.dcd. The .mol2 file generated by dump.mol2 contains the particle names and coordinates at time step 0. If there were any bonds specified, they would be included too. VMD or other applications can read in the .mol2 to obtain this information.

dump.dcd includes snapshots of the system state (particle position coordinates only) written every 100 time steps. This file can be loaded into visualization software such as VMD and played as a movie or read for analysis purposes.

If you have VMD installed, you can load up the entire simulation trajectory by running

vmd example1.mol2 example1.dcd

on the command line or by loading these files using VMD's GUI. For the best visualization, open VMD's Graphical Representation menu and set the Drawing Method to VDW. The default of lines will draw a seemingly random line through the simulation box. This is actually a dummy bond between particles 0 and 1, as VMD refuses to load a MOL2 file without any bonds specified in it.

---

Example 2: Using IMD

Here is the same simulation as Example 1, this time configured for real-time display in VMD using the IMD interface.

from hoomd_script import *

# create 1000 random particles of name A
init.create_random(N=1000, 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, alpha=1.0)

# integrate at constant temperature
integrate.nvt(dt=0.005, T=1.2, tau=0.5)

# dump a .mol2 file for the structure information
mol2 = dump.mol2()
mol2.write(filename='example2.mol2')

# setup the IMD server
analyze.imd(port=54321, period=500)

# run a very long time so the simulation can be watched in VMD
run(1e9)

Start the simulation running in HOOMD, then load up VMD. Inside VMD, create a new molecule and load the file example2.mol2 generated at the beginning of the simulation. Then go to the VMD command window and run the command

imd connect localhost 54321

The particles in the display window should begin moving. The display is of the current state of the simulation, updated in real-time. Again, the best visualization is obtained by setting the Drawing Method to VDW in VMD's Graphical Representation menu.

Switch back to the terminal where HOOMD is running and press CTRL-C to kill the simulation. It is set to run for an extremely long time on purpose to allow ample time to launch VMD and issue the imd command.

---

Example 3: Using the polymer generator

Here is a more complicated script that generates a system of bead-spring polymers which self-assemble into a hex phase when run for a few million time steps. The polymers are A6B7A6 block copolymers in an implicit solvent. The script also shows a few examples of how writing python code in the script can be handy: here the concentration phi_P is a parameter and math operations are performed to calculate the length of the box.

For more information on the model in this script, see
"Micellar crystals in solution from molecular dynamics simulations"
J. Chem. Phys. 128, 184906 (2008); DOI:10.1063/1.2913522
http://link.aip.org/link/?JCPSA6/128/184906/1

Any of the polymer systems in the paper could be easily run just by changing a few parameters in this script.

from hoomd_script import *
import math

# parameters
phi_P = 0.25
n_poly = 600
T = 1.2
polymer1 = dict(bond_len=1.2, type=['A']*6 + ['B']*7 + ['A']*6, 
                bond="linear", count=n_poly)

# perform some simple math to find the length of the box
N = len(polymer1['type']) * polymer1['count']
L = math.pow(math.pi * N / (6.0 * phi_P), 1.0/3.0)

# generate the polymer system
init.create_random_polymers(box=hoomd.BoxDim(L), polymers=[polymer1], 
                            separation=dict(A=0.35, B=0.35), seed=12)

# force field setup
harmonic = bond.harmonic()
harmonic.set_coeff('polymer', k=330.0, r0=0.84)
lj = pair.lj(r_cut=3.0)
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0, alpha=1.0)

# dump every 100,000 steps
mol2 = dump.mol2()
mol2.write(filename="example3.mol2")
dump.dcd(filename="example3.dcd", period=100000)

# integrate NVT for a bunch of time steps
integrate.nvt(dt=0.005, T=T, tau=0.5)
run(2000)

# uncomment the next run() command if you have a few hours to spare 
# running this on a GPU the resulting dump files should show the
# polymers self-assembling into the hex phase
# run(10e9)

---

Example 4: Using arbitrary input files

Of course, HOOMD is not limited by the built-in random initial condition generators used in the previous example. You can load in an arbitrary initial condition from a formatted xml file. Here is a simple example demonstrating most of the types of data that can be input (see XML File Format for full documentation of this format):

<?xml version="1.0" encoding="UTF-8"?>
<hoomd_xml>
<configuration time_step="0">
<box units="sigma"  Lx="10" Ly="10" Lz="10"/>
<!-- Setup the initial condition to place all particles in a line -->
<position units="sigma">
-3 0 0
-2 0 0
-1 0 0
0 0 0
1 0 0
2 0 0
3 0 0
</position>
<!-- Name the first 3 particles A and the rest B -->
<type>
A
A
A
B
B
B
B
</type>
<!-- Bond the particles together into a polymer chain -->
<bond>
polymer 0 1
polymer 1 2
polymer 2 3
polymer 3 4
polymer 4 5
polymer 5 6
</bond>
<!-- Give the particles a little kick with an initial velocity -->
<velocity units="sigma/tau">
1 2 3
3 2 1
1 0 0
0 1 0
0 0 1
-1 -2 -3
-3 -2 -1
</velocity>
</configuration>
</hoomd_xml>

Copy and paste this data to a file example4.xml. The initial conditions can be read into a simulation using the command init.read_xml as shown in the example script below.

from hoomd_script import *
import math

# read in the file
init.read_xml(filename="example4.xml")

# example4.xml defines a single polymer: use the same force field as in example 3
# force field setup
harmonic = bond.harmonic()
harmonic.set_coeff('polymer', k=330.0, r0=0.84)
lj = pair.lj(r_cut=3.0)
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0, alpha=1.0)

# dump every few steps
mol2 = dump.mol2()
mol2.write(filename="example4.mol2")
dump.dcd(filename="example4.dcd", period=10)

# integrate NVT for a bunch of time steps
integrate.nvt(dt=0.005, T=1.2, tau=0.5)
run(2000)

Example 5: Using a python loop to run many simulations at different parameters

In a loop, this example performs a simulation of a Lennard-Jones liquid at different system sizes one after the other. It shows how the power of python code in a hoomd script can be utilized to accomplish complicated tasks. For more information on python and a tutorial on the python language, check out http://www.python.org.

from hoomd_script import *

# create a list ranging from 1000 to 2000 in steps of 100
size_array = range(1000,2000,100)

# loop over the system sizes in the array
for size in size_array:
        # create 1000 random particles of name A
        init.create_random(N=size, 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, alpha=1.0)
        
        # integrate at constant temperature
        integrate.nvt(dt=0.005, T=1.2, tau=0.5)
        
        # run 10,000 time steps
        run(10e3)
        
        # we need to clear all saved variables before resetting and
        # starting the next simulation so there is no memory leak
        del lj
        init.reset()