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)
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);
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')