PolyethylenimineΒΆ
Generates a polydisperse set of polyethylenimine using Martini force fields.
imports
import Hoobas
import numpy as np
import hoomd
import math
from hoomd import md
initialize the execution context
hoomd.context.initialize()
set up simulation units
Hoobas.Units.SimulationUnits.set_default('amu', 'nm', 'kJ/mol')
define random direction
def random_direction():
return np.random.uniform(-1., 1., 3)
set up primary branching site
primary = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='p'))
add attachment sites
primary.add_free_attachment_sites(key_search={}, properties={'orientation': random_direction, 'qualifier': 'uncharged'})
set up charged beads
primary_charged = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='pq', charge=1.0))
add attachment site to charged beads
primary_charged.add_free_attachment_sites(key_search={}, properties={'orientation': random_direction, 'restriction': 'uncharged'})
set up secondary branching sites, attachment sites, and charges on them
secondary = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='s'))
secondary.add_free_attachment_sites(key_search={}, max_binding_per_bead=2, properties={'orientation': random_direction, 'qualifier': 'uncharged'})
secondary_charged = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='sq', charge=1.0))
secondary_charged.add_free_attachment_sites(key_search={}, max_binding_per_bead=2, properties={'orientation': random_direction, 'restriction': 'uncharged'})
set up tertiary branching sites, attachment sites, and charges on them
ternary = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='t', charge=0.0))
ternary.add_free_attachment_sites(key_search={}, max_binding_per_bead=3, properties={'orientation': random_direction, 'qualifier': 'uncharged'})
ternary_charged = Hoobas.LinearChain.BeadMonomer(Hoobas.CGBead.Bead(beadtype='tq', charge=1.0))
ternary_charged.add_free_attachment_sites(key_search={}, max_binding_per_bead=3, properties={'orientation': random_direction, 'restriction': 'uncharged'})
generate an array containing all types of monomers and assign distribution weights to each
monomer_arrays = [primary, primary_charged, secondary, secondary_charged, ternary, ternary_charged]
probabilities = np.array([0.01, 0.02, 0.5, 0.3, 0.15, 0.02])
probabilities_short = np.array([0.00, 0.00, 0.53, 0.3, 0.15, 0.02])
def polymer_rules(number_monomers):
if number_monomers < 10:
return np.random.choice(monomer_arrays, p=probabilities_short)
else:
return np.random.choice(monomer_arrays, p=probabilities)
polyethylenimine = hoobas.LinearChain.RandomPolymer(n_mono=100, distribution=polymer_rules)
domain = hoobas.SimulationDomain.EmptyCube(12.5)
builder = hoobas.Build.HOOMDBuilder(domain)
# normally the whole force-field should be mapped, ideally in a separate class
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:t-t', topodict={'energy_constant': 20000.0, 'distance_constant': 0.335}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:t-pq', topodict={'energy_constant': 7000.0, 'distance_constant': 0.296}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:sq-s', topodict={'energy_constant': 5000.0, 'distance_constant': 0.357}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:s-s', topodict={'energy_constant': 10000.0, 'distance_constant': 0.364}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:t-sq', topodict={'energy_constant': 6000.0, 'distance_constant': 0.286}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:t-p', topodict={'energy_constant': 3000.0, 'distance_constant': 0.303}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:s-t', topodict={'energy_constant': 20000.0, 'distance_constant': 0.360}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:s-pq', topodict={'energy_constant': 5000.0, 'distance_constant': 0.360}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:t-s', topodict={'energy_constant': 20000.0, 'distance_constant': 0.330}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:sq-t', topodict={'energy_constant': 6000.0, 'distance_constant': 0.342}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:s-sq', topodict={'energy_constant': 5000.0, 'distance_constant': 0.357}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.bond_types += hoobas.Composite.BondType('RandomPolymer:BeadMonomer-BeadMonomer:s-p', topodict={'energy_constant': 12000.0, 'distance_constant': 0.376}, unitdict={'energy_constant': 'E/LL', 'distance_constant': 'L'})
polyethylenimine.add_angles()
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-t-t', topodict={'energy_constant': 100.0, 'angle_constant': math.radians(132.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-t-p', topodict={'energy_constant': 50.0, 'angle_constant': math.radians(145.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-s-s', topodict={'energy_constant': 500.0, 'angle_constant': math.radians(152.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-t-t', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(91.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-s-t', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(159.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-t-t', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(138.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-t-pq', topodict={'energy_constant': 50.0, 'angle_constant': math.radians(145.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-sq-s', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(168.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-s-s', topodict={'energy_constant': 500.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-t-s', topodict={'energy_constant': 50.0, 'angle_constant': math.radians(157.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-sq-s', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(134.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-s-pq', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-t-s', topodict={'energy_constant': 50.0, 'angle_constant': math.radians(125.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-s-sq', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(121.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-t-sq', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(133.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-t-p', topodict={'energy_constant': 30.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-s-t', topodict={'energy_constant': 500.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-s-pq', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(132.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-t-pq', topodict={'energy_constant': 100.0, 'angle_constant': math.radians(128.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-s-sq', topodict={'energy_constant': 70.0, 'angle_constant': math.radians(151.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:t-s-p', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-t-pq', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(120.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:sq-s-s', topodict={'energy_constant': 50.0, 'angle_constant': math.radians(141.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-t-s', topodict={'energy_constant': 100.0, 'angle_constant': math.radians(154.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-sq-t', topodict={'energy_constant': 20.0, 'angle_constant': math.radians(141.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('Angle:s-s-sq', topodict={'energy_constant': 200.0, 'angle_constant': math.radians(180.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:t-t-s', topodict={'energy_constant': 1000.0, 'angle_constant': math.radians(68.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:s-t-t', topodict={'energy_constant': 1000.0, 'angle_constant': math.radians(68.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:sq-t-sq', topodict={'energy_constant': 100.0, 'angle_constant': math.radians(87.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:t-t-t', topodict={'energy_constant': 1700.0, 'angle_constant': math.radians(79.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:s-t-s', topodict={'energy_constant': 700.0, 'angle_constant': math.radians(88.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:s-t-pq', topodict={'energy_constant': 700.0, 'angle_constant': math.radians(88.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:pq-t-s', topodict={'energy_constant': 700.0, 'angle_constant': math.radians(88.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:s-t-p', topodict={'energy_constant': 700.0, 'angle_constant': math.radians(88.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.angle_types += hoobas.Composite.AngleType('NAngle:p-t-s', topodict={'energy_constant': 700.0, 'angle_constant': math.radians(88.0)}, unitdict={'energy_constant': 'E', 'angle_constant': ''})
polyethylenimine.remap_beadtype('p', 'P2')
polyethylenimine.remap_beadtype('t', 'P2')
polyethylenimine.remap_beadtype('s', 'P2')
polyethylenimine.remap_beadtype('pq', 'Qd')
polyethylenimine.remap_beadtype('sq', 'Qd')
polyethylenimine.remap_beadtype('tq', 'Qd')
builder.add_N_ext_obj(polyethylenimine, 50)
builder.fix_remaining_charge(ntype='Q0', nion_mass=72.0)
builder.set_attribute_by_beadtype('P2', 'mass', 72.0)
builder.set_attribute_by_beadtype('Qd', 'mass', 72.0)
builder.add_rho_molar_ions(14.0 * 0.9 * 0.9, qtype='P4', ion_mass=72.0, q=0.0)
builder.add_rho_molar_ions(14.0 * 0.1 * 0.9, qtype='BP4', ion_mass=72.0, q=0.0)
# set N to N angles in the structure
for bond in builder.bonds:
if bond.topology_tuple[1] == min(bond.topology_tuple):
bond.topology_name = 'N' + bond.topology_name
sn = hoomd.data.make_snapshot(N=builder.num_beads,
particle_types=builder.bead_types,
box=hoomd.data.boxdim(*builder.current_box()))
builder.set_snapshot(sn)
system = hoomd.init.read_snapshot(sn)
# if the force-field has been clearly defined in hoobas it can be directly imported
bonded = hoomd.md.bond.harmonic()
for bond_t in builder.bond_types:
bonded.bond_coeff.set(bond_t.typename, k=bond_t['energy_constant'], r0=bond_t['distance_constant'])
angle = hoomd.md.angle.cosinesq()
for angle_t in builder.angle_types:
angle.angle_coeff.set(angle_t.typename, k=angle_t['energy_constant'], t0=angle_t['angle_constant'])
nlist = hoomd.md.nlist.cell()
charges = hoomd.md.charge.pppm(hoomd.group.all(), nlist)
charges.set_params(Nx=16, Ny=16, Nz=16, order=6, rcut=2.5*0.62)
nonbonded = hoomd.md.pair.lj(nlist=nlist, r_cut=2.5*0.62)
def set_lj(m):
hoomd.util.quiet_status()
for beadtypeA in builder.bead_types:
for beadtypeB in builder.bead_types:
params = hoobas.MartiniModels.MartiniForceMappings.get_pair(beadtypeA, beadtypeB)
nonbonded.pair_coeff.set(beadtypeA,
beadtypeB,
epsilon=params['eps'],
sigma=m*params['sigma'],
r_cut=2.5*params['sigma'])
hoomd.util.unquiet_status()
traj = hoomd.dump.gsd('PE.gsd', 10000, hoomd.group.all(), overwrite=True)
log = hoomd.analyze.log('PE.log', period=2500, quantities=['temperature', 'pressure', 'potential_energy', 'pair_lj_energy'], overwrite=True)
neq_steps = 10
neq_values = np.linspace(0.0, 1.0, neq_steps)
nvei = hoomd.md.integrate.nve(group=hoomd.group.all(), limit=.002)
hoomd.md.integrate.mode_standard(dt=.001)
for step in range(neq_steps):
set_lj(neq_values[step])
print 'setting LJ to :' + str(neq_values[step])
hoomd.run(1000)
nvei.disable()
set_lj(1.0)
..
hoomd.md.integrate.mode_standard(dt=0.0001)
nvt = hoomd.md.integrate.langevin(group=hoomd.group.all(), seed=1255345, kT=2.5)
nvt.set_gamma('P2', 2.0)
nvt.set_gamma('Qd', 5.0)
hoomd.run(10000)
hoomd.md.integrate.mode_standard(dt=0.0005)
hoomd.run(1e6)
hoomd.md.integrate.mode_standard(dt=0.003)
hoomd.run(1e6)