University of Michigan Logo

Here is a script that generates a system of bead-spring polymers that 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

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

To run this example, copy the examples directory from share/hoomd/examples your hoomd installation directory.

cp -R HOOMD-INSTALL-DIR/share/hoomd/examples ~/hoomd-examples
cd ~/hoomd-examples

Then run:

The results of this example may be visualized with VMD:

hoomd_script::init::create_random_polymers hoomd_script::dump::xml hoomd_script::dump::dcd hoomd_script::pair::lj hoomd_script::bond::harmonic hoomd_script::integrate::mode_standard hoomd_script::integrate::nvt
# Due to deficiencies in doxygen, the commands used in this example are listed explicitly here
# run this script with "python -x filename" to skip the first line, or remove this header
# ---- ----
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']
# generate the polymer system
init.create_random_polymers(box=data.boxdim(volume=math.pi * N / (6.0 * phi_P)), polymers=[polymer1],
separation=dict(A=0.35, B=0.35), seed=12)
# force field setup
harmonic = bond.harmonic()
harmonic.bond_coeff.set('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
dump.xml(filename="create_random_polymers.xml", vis=True)
dump.dcd(filename="create_random_polymers.dcd", period=100000)
# integrate NVT for a bunch of time steps
all = group.all()
integrate.nvt(group=all, T=1.2, tau=0.5)
# 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(10e6)