In [1]:

    

import hoomd
import hoomd.md
import gsd
import gsd.hoomd

import scipy.constants as const


    

In [10]:

    

#Calculate sigma and epsilon based on TraPPE FF
epsilon = 98.0 # units of K
# Convert to units of
epsilon = epsilon / const.physical_constants["electron volt-kelvin relationship"][0]
print(epsilon)


    

    




0.008444983701655409

In [4]:

    

hoomd.context.initialize("");
ethane_snap = hoomd.data.make_snapshot(N=2, box=hoomd.data.boxdim(L=5), particle_types=['C3', 'C3'],  bond_types=['alkanes'])


    

In [5]:

    

# set particle positions
ethane_snap.particles.position[:] = [[0,0,0], [1.597649,0.000000, 0.000000]]

# set particle masses
ethane_snap.particles.mass[:] = [15.03452, 15.03452]

# create bonds
ethane_snap.bonds.resize(1)
ethane_snap.bonds.group[:] = [[0,1]]


    

In [7]:

    

hoomd.init.read_snapshot(ethane_snap);

# set LJ and bond parameters - Actually don't need LJ for ethane. 
nl = hoomd.md.nlist.cell();
lj = hoomd.md.pair.lj(r_cut=99.0, nlist=nl);



sigma = 3.750 # units of angstrom

lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=2.0);


    

    




**ERROR**: Cannot initialize more than once






    




---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-7-7187270d3223> in <module>()
----> 1 hoomd.init.read_snapshot(ethane_snap);
      2 
      3 # set LJ and bond parameters - Actually don't need LJ for ethane.
      4 nl = hoomd.md.nlist.cell();
      5 lj = hoomd.md.pair.lj(r_cut=99.0, nlist=nl);

/home/jessica/anaconda3/lib/python3.6/site-packages/hoomd/init.py in read_snapshot(snapshot)
    237     if is_initialized():
    238         hoomd.context.msg.error("Cannot initialize more than once\n");
--> 239         raise RuntimeError("Error initializing");
    240 
    241     # broadcast snapshot metadata so that all ranks have _global_box (the user may have set box only on rank 0)

RuntimeError: Error initializing

In [ ]:

    

print(const.physical_constants["electron volt-kelvin relationship"])


    

In [ ]:

    

hoomd.context.initialize("");

# Initialize snapshot
snap = hoomd.data.make_snapshot(N=4, box=hoomd.data.boxdim(L=10), particle_types=['A', 'B'], bond_types=['polymer'])

# Assign coordinates
snap.particles.position[0] = [1,2,3];
snap.particles.position[1] = [-1,-2,-3];
snap.particles.position[2] = [3,2,1];
snap.particles.position[3] = [-3,-2,-1];

# Define bonds
snap.bonds.resize(2);
snap.bonds.group[:] = [[0,1], [1, 2]]

# Initialize hoomd blue (can't set forces without doing this)
hoomd.init.read_snapshot(snap);

# Set LJ parameters
nl = hoomd.md.nlist.cell();
lj = hoomd.md.pair.lj(r_cut=3.0, nlist=nl);
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=2.0);
lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0);
lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0);
all = hoomd.group.all();

# Define bond types
bond1 = hoomd.md.bond.harmonic(name="polymer")
bond1.bond_coeff.set('polymer', k=330.0, r0=0.84)

# Set up dynamics
# hoomd.md.integrate.mode_standard(dt=0.001);
# hoomd.md.integrate.langevin(group=all, kT=1.0, seed=987);


    

In [ ]:

    

s = gsd.hoomd.Snapshot()  
s.particles.N = 4  
s.particles.types = ['A', 'B']  
s.particles.typeid = [0,0,1,1]  
s.particles.position = [[0,0,0],[1,1,1], [-1,-1,-1], [1,-1,-1]]  
s.configuration.box = [3, 3, 3, 0, 0, 0] 
s.bonds.group = [[0,1], [1, 2], [2,3]];
s.bonds.typeid = [1,2,1]
#s.bonds.types 
s.bonds.N = 3

print(s.bonds.N, s.bonds.group, s.bonds.typeid, s.bonds.types)
traj = gsd.hoomd.open(name='test.gsd', mode='wb')
traj.append(s)


    

In [ ]:

    

# Write snapshot
hoomd.dump.gsd("init.gsd", period=None, group=all, overwrite=True);


    

In [ ]:

    

# Read snapshot
snap = gsd.hoomd.open('init.gsd','rb')