#coding: utf-8
import copy
import random
import os
from math import *
try:
import DNA3SPN.make_bp_params
import DNA3SPN.pdb2cg_dna
except ImportError:
print('Error while importing DNA3SPN subpackages, LinearChain will be unable to build this DNA model')
import numpy as np
from hoobas import CGBead
from hoobas.Units import SimulationUnits as SimUnits
from hoobas import Util
import warnings
import operator
from hoobas import Composite
from hoobas.Quaternion import Quat
import networkx
[docs]class LinearChain(Composite.CompositeObject):
"""
This is a class for linear chains in the model. Subclasses are expected to give constants for angles and other funct
ions. initial chain is expected to run along Z, while DNA is added on x-y plane
Properties :
zero_pos : position of the first bead in the chain
center_position : position of the center of the chain, calculated by arithmetic average
bond_types : returns the bond types in the chain in the form [['name', k, r0], ... ]
angle_types : returns the angle types in the chain in the form [['name', k, t0], ...]
dihedral_types : returns the dihedral types in the form [['name', params0, params1, ... ], ...], parameter number
depends on the potential type, OPLS has 4
pnum_offset : shifts all bead numbers in potentials by this constant
Methods :
add_dna(#, n_ss, n_ds, sticky_end) : adds DNA defined by the parameters to the linear chain. DNA are grafted to
remaining sites
randomize_dirs : randomizes the chain directions
change_remaining_att_sites(key_search, max_num_bindings = 1) : changes the remaining attachment sites in the chain
by using a key search in the beads. the key is a dictionary (e.g. {'beadtype':'P'}). max_num_bindings is the maximum
number of grafts a site can have, including all previous bindings
graft_N_ext_obj(N, obj, connecting_bond = None, add_attachment_sites = False) : grafts N copies (N<remaining sites)
of obj to the chain, connecting the newly added chain by connecting_bond[0]. If add_attachment_sites is true, the
grafted object attachment sites will be appended to the remaining sites
random_roration():
Rotates the chain by a random matrix
"""
def __init__(self, n_monomer=None, kuhn_length=None, units=None):
super(LinearChain, self).__init__(units)
self.number_monomers = n_monomer
self.kuhn_length = kuhn_length
self.dir = np.array([0., 0., 1.0], dtype=np.float)
self.chain_name = None
@Composite.deferred
def random_rotation(self):
_r_mat = Util.gen_random_mat()
self.rotate(_r_mat)
def __str__(self):
if self.chain_name is None:
return str(self.__class__.__name__)
else:
return str(self.chain_name)
def build(self):
if hasattr(self.number_monomers, '__call__'):
self.number_monomers = self.number_monomers()
super(LinearChain, self).build()
[docs] def XY_rotate(self, angle):
"""
rotate around the polymer backbone
:param angle: angle to rotate
:return:
"""
op = Quat(np.array([sin(angle/2.0) * self.dir[0],
sin(angle/2.0) * self.dir[1],
sin(angle/2.0) * self.dir[2],
cos(angle/2.0)]))
self.rotate(op)
@Composite.deferred
def rotate(self, operation):
q = Quat(operation)
tr = q.transform
self.dir = tr.dot(self.dir)
super(LinearChain, self).rotate(operation)
@property
def zero_pos(self):
return self.positions[0, :]
@zero_pos.setter
def zero_pos(self, val):
diff = val - self.positions[0, :]
self.positions += diff
@property
def thickness(self): # approximation to the thickness for multilayers; layers are by definition along z
maxz = -np.infty
minz = np.infty
for bead in self.beads:
maxz = np.maximum(bead.position[2], maxz)
minz = np.minimum(bead.position[2], minz)
return maxz - minz
@property
def center_position(self):
maxv = np.array([-np.infty, -np.infty, -np.infty], dtype=np.single)
minv = np.array([np.infty, np.infty, np.infty], dtype=np.single)
for bead in self.beads:
maxv = np.maximum(bead.position, maxv)
minv = np.minimum(bead.position, minv)
return 0.5*(maxv+minv)
@center_position.setter
def center_position(self, value):
cacp = self.center_position
diff = value - cacp
self.positions += diff
@Composite.deferred
def randomize_dirs(self, tol=1e-1):
for i in range(1, self.positions.shape[0]):
old_vec = np.array(self.positions[i, :] - self.positions[i - 1, :])
mat = Util.gen_random_mat()
new_vec = mat.dot(old_vec)
self.positions[i:, :] += (tol * new_vec - tol * old_vec)
def do_on_copy(self): # virtual method for special cases where copying objects needs to do something (increment molecule index or whatever)
pass
def set_residues(self, res):
for bead in self.beads:
bead.residue = res
def type_by_idx(self, idx):
return str(self.beads[idx])
[docs]class Polyaniline(LinearChain):
def __init__(self, n_monomer=None):
_units = SimUnits()
_units.set_energy('kJ/mol')
_units.set_length('R')
_units.set_mass('amu')
LinearChain.__init__(self, n_monomer=n_monomer, kuhn_length=[5.56 / 20, 5.52 / 20, 5.47 / 20], units=_units)
# set types, this is a list initializer since bonds are harmonics
self.b_types += Composite.BondType('', ['paA-paA', 71.84 * 4.18 * (20 ** 2), 5.56 / 20])
self.b_types += Composite.BondType('', ['paA-paB', 111.28 * 4.18 * 20 ** 2, 5.52 / 20])
self.b_types += Composite.BondType('', ['paB-paB', 143.03 * 4.18 * 20 ** 2, 5.47 / 20])
self.a_types += Composite.AngleType('', ['paA-paA-paB', 19.46 * 4.18 * sin(1.94) ** 2, 1.94])
self.a_types += Composite.AngleType('', ['paA-paB-paB', 38.80 * 4.18 * sin(2.45) ** 2, 2.45])
self.d_types += Composite.DihedralType('paA-paA-paB-paB',
{'phi': 0.5 * 0.8 * 4.18, '2phi': 0.5 * 0.73 * 4.18,
'3phi': 0.5 * -0.03 * 4.18, '4phi': 0.5 * 0.22 * 4.18},
{'phi': 'E', '2phi': 'E', '3phi': 'E', '4phi': 'E'})
self.d_types += Composite.DihedralType('paA-paB-paB-paA',
{'phi': 0.5 * 0.66 * 4.18, '2phi': 0.5 * -0.54 * 4.18,
'3phi': 0.5 * -0.1 * 4.18, '4phi': 0.5 * 0.02 * 4.18},
{'phi': 'E', '2phi': 'E', '3phi': 'E', '4phi': 'E'})
self.d_types += Composite.DihedralType('paB-paA-paA-paB',
{'phi': 0.5 * 0.45 * 4.18, '2phi': 0.5 * -1.33 * 4.18,
'3phi': 0.5 * 0.13 * 4.18, '4phi': 0.5 * 0.67 * 4.18},
{'phi': 'E', '2phi': 'E', '3phi': 'E', '4phi': 'E'})
self.wdv = [['paA', 'paA', 0.32 * 4.18, 5.14 / 20], ['paB', 'paB', 0.34 * 4.18, 5.14 / 20],
['paA', 'paB', (0.32 * 0.34) ** 0.5 * 4.18, 5.14 / 20]]
self.mass = [163.0, 162.0]
def build(self):
self.__build_chain()
self.randomize_dirs()
self.build_finalize()
def __build_chain(self):
# base moomer ABBA
self.beads += CGBead.Bead(beadtype='A', position=[0., 0., 0.], mass=163.0)
self.beads += CGBead.Bead(beadtype='B', position=[0., 0., self.kuhn_length[1]], mass=162.0)
self.beads += CGBead.Bead(beadtype='B', position=[0., 0., self.kuhn_length[1] + self.kuhn_length[2]], mass=162.0)
self.beads += CGBead.Bead(beadtype='A', position=[0., 0., self.kuhn_length[1] + 2 * self.kuhn_length[2]], mass=162.0)
self.bonds += Composite.Bond(['paA-paB', 0, 1])
self.bonds += Composite.Bond(['paB-paB', 1, 2])
self.bonds += Composite.Bond(['paA-paB', 2, 3])
self.angles += Composite.Angle(['paA-paB-paB', 0, 1, 2])
self.angles += Composite.Angle(['paA-paB-paB', 1, 2, 3])
self.dihedrals += Composite.Dihedral(['paA-paB-paB-paA', 0, 1, 2, 3])
self.p_num = [0, 1, 2, 3]
# struct ia ABBA type repeated
for i in range(self.number_monomers - 1):
self.beads += CGBead.Bead(position=[0.0, 0.0, self.positions[-1, 2] + self.kuhn_length[0]], beadtype='A',
mass=162.0)
self.p_num.append(self.p_num[-1] + 1)
self.bonds += Composite.Bond(['paA-paA', self.p_num[-1], self.p_num[-2]])
self.angles += Composite.Angle(['paA-paA-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3]])
self.dihedrals += Composite.Dihedral(
['paA-paA-paB-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3], self.p_num[-4]])
self.beads += CGBead.Bead(position=[0.0, 0.0, self.positions[-1, 2] + self.kuhn_length[1]], beadtype='B',
mass=162.0)
self.p_num.append(self.p_num[-1] + 1)
self.bonds += Composite.Bond(['paA-paB', self.p_num[-1], self.p_num[-2]])
self.angles += Composite.Angle(['paA-paA-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3]])
self.dihedrals += Composite.Dihedral(
['paB-paA-paA-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3], self.p_num[-4]])
self.beads += CGBead.Bead(position=[0.0, 0.0, self.positions[-1, 2] + self.kuhn_length[2]], beadtype='B',
mass=162.0)
self.p_num.append(self.p_num[-1] + 1)
self.bonds += Composite.Bond(['paB-paB', self.p_num[-1], self.p_num[-2]])
self.angles += Composite.Angle(['paA-paB-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3]])
self.dihedrals += Composite.Dihedral(
['paA-paA-paB-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3], self.p_num[-4]])
self.beads += CGBead.Bead(position=[0.0, 0.0, self.positions[-1, 2] + self.kuhn_length[1]], beadtype='A',
mass=163.0)
self.p_num.append(self.p_num[-1] + 1)
self.bonds += Composite.Bond(['paA-paB', self.p_num[-1], self.p_num[-2]])
self.angles += Composite.Angle(['paA-paB-paB', self.p_num[-1], self.p_num[-2], self.p_num[-3]])
self.dihedrals += Composite.Dihedral(
['paA-paB-paB-paA', self.p_num[-1], self.p_num[-2], self.p_num[-3], self.p_num[-4]])
[docs]class GenericPolymer(LinearChain):
def __init__(self, n_mono=100, kuhn_length=1.0, beadname=None, rigid=False, units=None):
if units is None:
units = SimUnits()
LinearChain.__init__(self, n_monomer=n_mono, kuhn_length=kuhn_length, units=units)
if beadname is None:
self.BeadGenericName = 'GenericPolymer'
else:
self.BeadGenericName = str(beadname)
self.wdv = [[self.BeadGenericName, self.BeadGenericName, 1.0, kuhn_length]]
self.b_types += [self.BeadGenericName + 'Backbone', 1.0, self.kuhn_length]
self.is_rigid = rigid
def build(self):
super(GenericPolymer, self).build()
self.__build_chain()
if self.is_rigid:
self.a_types += [self.BeadGenericName + 'BackboneAngle', 1.0, pi]
self.__build_angles()
self.randomize_dirs()
def __build_chain(self):
self.beads += CGBead.Bead(position=[0.0, 0.0, 0.0], beadtype=self.BeadGenericName)
self.p_num = [0]
self.mass = [1.0]
for i in range(self.number_monomers - 1):
self.beads += CGBead.Bead(position=[0.0, 0.0, self.positions[-1, 2] + self.kuhn_length], beadtype=self.BeadGenericName)
self.p_num.append(self.p_num[-1]+1)
self.bonds += [self.BeadGenericName+'Backbone', self.p_num[-1], self.p_num[-2]]
def __build_angles(self):
for i in range(self.number_monomers - 2):
self.angles += [self.BeadGenericName + 'BackboneAngle', self.p_num[i], self.p_num[i+1], self.p_num[i+2]]
[docs]class GenericRingPolymer(LinearChain):
def __init__(self, n_mono=100, kuhn_length=1.0, beadname=None, rigid=False, units=None):
if units is None:
units = SimUnits()
LinearChain.__init__(self, n_monomer=n_mono, kuhn_length=kuhn_length, units=units)
if beadname is None:
self.BeadGenericName = beadname
else:
self.BeadGenericName = str(beadname)
self.wdv = [[self.BeadGenericName, self.BeadGenericName, 1.0, kuhn_length]]
self.b_types += [self.BeadGenericName + 'Backbone', 1.0, self.kuhn_length]
self.is_rigid = rigid
self.R = None
def build(self):
super(GenericRingPolymer, self).build()
self.R = self.kuhn_length * self.number_monomers / (2 * pi)
self.__build_chain()
if self.is_rigid:
self.__build_angles()
self.randomize_dirs()
def __build_chain(self):
self.beads += CGBead.Bead(position=[self.R, 0.0, 0.0], beadtype=self.BeadGenericName)
self.p_num = [0]
self.mass = [1.0]
for i in range(self.number_monomers - 1):
self.beads += CGBead.Bead(position=[self.R * cos(2 * pi * (i + 1) / (self.number_monomers + 1)),
self.R * sin(2 * pi * (i + 1) / (self.number_monomers + 1)), 0.0],
beadtype=self.BeadGenericName)
self.p_num.append(self.p_num[-1] + 1)
self.bonds.append([self.BeadGenericName + 'Backbone', self.p_num[-1], self.p_num[-2]])
self.bonds += [self.BeadGenericName + 'Backbone', self.p_num[-1], self.p_num[0]]
def __build_angles(self):
self.a_types += [self.BeadGenericName + 'BackboneAngle', 1.0, pi]
for i in range(self.number_monomers):
# pnum lists should range as [i, i+1, i+2, ..., i+self.number_monomers-1] so we contain overflow by % self.number_monomers
self.angles += [self.BeadGenericName + 'BackboneAngle',
self.p_num[i],
(self.p_num[i + 1] % self.number_monomers) + self.p_num[0],
(self.p_num[i + 2] % self.number_monomers) + self.p_num[0]]
[docs]class GenericBlockRingPolymer(LinearChain):
def __init__(self, n_mono=tuple([100]), kuhn_length=1.0, rigid=None, units=None, randomize_blocks=False):
super(GenericBlockRingPolymer, self).__init__(n_monomer=sum(n_mono), kuhn_length=kuhn_length, units=units)
self.blocks = n_mono
self.nblocks = len(n_mono)
self.beadtypes = [chr(ord('A') + i) for i in range(self.nblocks)]
if rigid is None:
self.rigid = [False for i in range(self.blocks)]
else:
self.rigid = rigid
self.random_blocks = randomize_blocks
self.R = self.kuhn_length * self.number_monomers / (2.0 * pi)
def build(self):
if self.random_blocks:
self.__build_chain_random()
else:
self.__build_chain()
self.__build_bonds()
self.__build_angles()
self.randomize_dirs()
def __build_chain(self):
for beadt in range(self.nblocks):
for bead in range(self.blocks[beadt]):
self.beads += CGBead.Bead(position=[self.R * cos(2*pi*len(self.beads)/self.number_monomers),
self.R * sin(2*pi*len(self.beads)/self.number_monomers), 0.],
beadtype=chr(ord('A')+beadt))
try:
self.p_num .append(self.p_num[-1] + 1)
except IndexError:
self.p_num.append(0)
def __build_chain_random(self):
remains = [self.beadtypes[j] for j in range(len(self.blocks)) for i in range(self.blocks[j])]
for beadi in range(self.number_monomers):
self.beads += CGBead.Bead(position=[self.R * cos(2*pi*len(self.beads)/self.number_monomers),
self.R * sin(2*pi*len(self.beads)/self.number_monomers), 0.],
beadtype=remains.pop(np.random.randint(0, len(remains))))
try:
self.p_num.append(self.p_num[-1] + 1)
except IndexError:
self.p_num.append(0)
def __build_bonds(self):
for beadi in range(len(self.beads)):
nexti = (beadi + 1) % self.number_monomers
self.bonds += [self.beads[beadi].beadtype + '-' + self.beads[nexti].beadtype, beadi, nexti]
def __build_angles(self):
for beadi in range(len(self.beads)):
nexti = (beadi + 1) % self.number_monomers
previ = (beadi - 1 + self.number_monomers) % self.number_monomers
if self.rigid[self.beadtypes.index(self.beads[beadi].beadtype)]:
self.angles += [self.beads[previ].beadtype+'-'+self.beads[beadi].beadtype+'-'+self.beads[nexti].beadtype,
previ, beadi, nexti]
[docs]class PolystyreneSulfonateMonomer(LinearChain):
def __init__(self, tact):
_units = SimUnits(length='nm', energy='kJ/mol', mass='amu')
super(PolystyreneSulfonateMonomer, self).__init__(n_monomer=3, kuhn_length=1.0, units=_units)
self.beads += CGBead.Bead(beadtype='PSSA', mass=27.0, charge=0.0, position=[0.0, 0.0, 0.0])
self.beads += CGBead.Bead(beadtype='PSSB', mass=72.0, charge=-0.505, position=[1.0, 0.0, 0.0])
self.beads += CGBead.Bead(beadtype='PSSC', mass=80.0, charge=-0.495, position=[2.0, 0.0, 0.0])
self.chain_side = tact
self.bonds += ['PSSAB', 0, 1]
self.bonds += ['PSSBC', 1, 2]
self.att_list += Composite.AttachmentSite(index=0,
orientation=np.array([0.0, 0.0, 1.0]),
qualifier='forward',
restriction='backward')
self.att_list += Composite.AttachmentSite(index=1,
orientation=np.array([0.0, 0.0, -1.0]),
qualifier='backward',
restriction='forward')
self.p_num = [0, 1, 2]
self.chain_name = 'PSS-monomer' + str(self.chain_side)
[docs]class DNAChain(LinearChain): # this in reduced units
"""
DNA chain based on works published:
C. Knorowski, S. Burleigh, A. Travesset, "Dynamics and statics of DNA-programmable nanoparticle self-assembly
and crystallization", Phys. Rev. Lett. 106, 2011
T. I. N. G. Li, R. Sknepnek, R. J. MacFarlane, C. A. Mirkin, M. Olvera de la Cruz, "Modeling the crystallization
of spherical nucleic acid nanoparticle conjugates with molecular dynamics simulations", J. Am. Chem. Soc. 135, 2013
"""
def __init__(self, n_ss, n_ds, sticky_end, flexor=None, **kwargs):
LinearChain.__init__(self, n_monomer=n_ss + n_ds + len(sticky_end), kuhn_length=0.84)
self.num_ssDNA_bases = n_ss
self.num_dsDNA_bases = n_ds
self.sticky_end = sticky_end
""" Bead types:
S: ssDNA
A: dsDNA
B: sticky-end
NP: nanoparticle
Fl: flanking bead
C:
"""
# Bond types
self.b_types += [['S-NP', 330.0, 0.84],
['S-S', 330.0, 0.84 * 0.5],
['S-A', 330.0, 0.84 * 0.75],
['backbone', 330.0, 0.84],
['A-B', 330.0, 0.84 * 0.75],
['B-B', 330.0, 0.84 * 0.5],
['C-FL', 330.0, 0.84 * 0.5],
['B-FL', 330.0, 0.84 * 0.5 * 1.41],
['C-C', 330.0, 0.84 * 0.5]]
# Angle types
self.a_types += [['flexor', 2.0, pi],
['dsDNA', 30.0, pi],
['C-C-C', 10.0, pi],
['B-B-B', 10.0, pi],
['FL-C-FL', 100.0, pi],
['A-B-C', 120.0, 0.5 * pi],
['A-A-B', 2.0, pi],
['A-B-B', 2.0, pi],
['C-C-endFL', 50.0, pi]]
self.rem_sites = [] # we dont attach anything to the DNA chain
self.inner_types = []
if flexor is not None:
self.flexors = flexor - self.num_ssDNA_bases
else:
self.flexors = flexor
if self.num_dsDNA_bases < 2:
self.flexors = []
if self.number_monomers < 1:
raise StandardError('DNA chain does not have a single monomer, unable to build')
self.mass = []
self.__build_chain()
self.__build_bonds()
self.__build_angles()
self.__build_beads()
#############################
# to be deprecated
# overloading older methods
@property
def harmonic_bonds_in_one_dna_chain(self):
return self.bonds
@property
def angles_in_one_dna_chain(self):
return self.angles
@property
def beads_in_one_dna_chain(self):
return self.beads
########################
def __build_chain(self):
num_ss_left = self.num_ssDNA_bases
num_ds_left = self.num_dsDNA_bases
num_stick_left = len(self.sticky_end)
monomers_left = self.number_monomers
sticky_left = copy.deepcopy(self.sticky_end)
# Build the chain one bead at a time, in the order of
# 1) ssDNA
# 2) dsDNA
# 3) sticky-end
while monomers_left > 0:
# Add ssDNA
if num_ss_left > 0:
num_ss_left -= 1
self.inner_types.append('S')
self.mass.append(1.0)
try:
self.p_num.append(self.p_num[-1] + 1)
except IndexError:
self.p_num.append(0)
# Add dsDNA
elif num_ds_left > 0:
num_ds_left -= 1
self.inner_types.append('A')
self.mass.append(4.0)
try:
self.p_num.append(self.p_num[-1] + 1)
except IndexError:
warnings.warn('DNA chain starts with dsDNA', UserWarning)
self.p_num.append(0)
# Add the sticky end
else:
num_stick_left -= 1
self.inner_types.append('B')
self.mass.append(1.0)
try:
self.p_num.append(self.p_num[-1] + 1)
except IndexError:
warnings.warn('DNA chain starts with sticky ends, there is no linker / duplexor, check for errors', UserWarning)
self.p_num.append(0)
self.inner_types.append(sticky_left.pop(0))
self.mass.append(1.0)
self.p_num.append(self.p_num[-1] + 1)
# Add flanker beads to make sure the chain doesn't bind to multiple other chains
self.inner_types.append('FL')
self.mass.append(1.0)
self.p_num.append(self.p_num[-1] + 1)
self.inner_types.append('FL')
self.mass.append(1.0)
self.p_num.append(self.p_num[-1] + 1)
monomers_left -= 1
# If there is a sticky end, the ends of the sticky ends could also bind
# Add another flanking bead to protect against this
if len(self.sticky_end) > 1:
self.inner_types.append('FL')
self.mass.append(1.0)
self.p_num.append(self.p_num[-1] + 1)
def __build_bonds(self):
for idx in range(len(self.p_num) - 1):
# ssDNA - ssDNA
if self.inner_types[idx] == 'S' and self.inner_types[idx + 1] == 'S':
self.bonds += ['S-S', idx, idx + 1]
# ssDNA - dsDNA
if self.inner_types[idx] == 'S' and self.inner_types[idx + 1] == 'A':
self.bonds += ['S-A', idx, idx + 1]
# dsDNA - dsDNA
if self.inner_types[idx] == 'A' and self.inner_types[idx + 1] == 'A':
self.bonds += ['backbone', idx, idx + 1]
# dsDNA - sticky-end
if self.inner_types[idx] == 'A' and self.inner_types[idx + 1] == 'B':
self.bonds += ['A-B', idx, idx + 1]
# ssDNA - sticky-end
if self.inner_types[idx] == 'S' and self.inner_types[idx + 1] == 'B': # special case with no dsDNA
self.bonds += ['S-B', idx, idx + 1]
###
# Stuff get messier afterwards. for each B, there is 1 sticky + 2 flanking
# plus the end flanking to check.
###
# sticky-end - sticky-end
if self.inner_types[idx] == 'B' and self.inner_types[idx + 4] == 'B':
self.bonds += ['B-B', idx, idx + 4]
# sticky-end - flanker1
if self.inner_types[idx] == 'B' and self.inner_types[idx + 2] == 'FL':
self.bonds += ['B-FL', idx, idx + 2]
# sticky-end - flanker2
if self.inner_types[idx] == 'B' and self.inner_types[idx + 3] == 'FL':
self.bonds += ['B-FL', idx, idx + 3]
# backbone - sticky-end
if self.inner_types[idx] == 'B' and self.inner_types[idx + 1] in self.sticky_end:
self.bonds += ['B-C', idx, idx + 1]
if self.inner_types[idx] == 'B' and self.inner_types[idx + 4] == 'FL':
self.bonds += ['B-FL', idx, idx + 4]
if self.inner_types[idx] in self.sticky_end and self.inner_types[idx + 1] == 'FL':
self.bonds += ['C-FL', idx, idx + 1]
if self.inner_types[idx] in self.sticky_end and self.inner_types[idx + 2] == 'FL':
self.bonds += ['C-FL', idx, idx + 2]
try:
if self.inner_types[idx] in self.sticky_end and self.inner_types[idx + 4] in self.sticky_end:
self.bonds += ['C-C', idx, idx + 4]
except IndexError:
pass
if self.inner_types[idx] in self.sticky_end and self.inner_types[idx + 3] == 'FL':
self.bonds += ['C-FL', idx, idx + 3]
def __build_angles(self):
for idx in range(len(self.p_num) - 2):
if self.inner_types[idx] == 'S':
continue
if self.inner_types[idx] == 'A':
if self.inner_types[idx + 1] == 'A' and self.inner_types[idx + 2] == 'A':
if self.flexors is not None and idx + 1 in self.flexors:
self.angles += ['flexor', idx, idx + 1, idx + 2]
else:
self.angles += ['dsDNA', idx, idx + 1, idx + 2]
if self.inner_types[idx + 1] == 'A' and self.inner_types[idx + 2] == 'B':
self.angles += ['A-A-B', idx, idx + 1, idx + 2]
try:
# offset by 4 see generator above
if self.inner_types[idx + 1] == 'B' and self.inner_types[idx + 5] == 'B':
self.angles += ['A-B-B', idx, idx + 1, idx + 5]
except IndexError:
pass
try:
if self.inner_types[idx + 1] == 'B' and self.inner_types[idx + 2] in self.sticky_end:
self.angles += ['A-B-C', idx, idx + 1, idx + 2]
except IndexError:
pass
try:
if self.inner_types[idx] == 'B' \
and self.inner_types[idx + 4] == 'B'\
and self.inner_types[idx + 8] == 'B':
self.angles += ['B-B-B', idx, idx + 4, idx + 8]
except IndexError:
pass
# have to check angles FL - C - FL, C - C - C and C - C - FL
if self.inner_types[idx] in self.sticky_end:
# should be, sanity check
if self.inner_types[idx + 1] == 'FL' and self.inner_types[idx + 2] == 'FL':
self.angles += ['FL-C-FL', idx + 1, idx, idx + 2]
try:
if self.inner_types[idx + 4] in self.sticky_end and self.inner_types[idx + 8] in self.sticky_end:
self.angles += ['C-C-C', idx, idx + 4, idx + 8]
except IndexError:
pass
try:
if self.inner_types[idx + 4] in self.sticky_end and self.inner_types[idx + 7] == 'FL':
self.angles += ['C-C-endFL', idx, idx+4, idx+7]
except IndexError:
pass
def __build_beads(self):
base_length = self.kuhn_length
for idx in range(len(self.p_num)):
if idx == 0:
self.beads += CGBead.Bead(position=[0., 0., 0.], mass=self.mass[idx], beadtype=self.inner_types[idx])
continue
# S - S bond
if self.inner_types[idx - 1] == 'S' and self.inner_types[idx] == 'S':
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.0, 0.50 * base_length], mass=self.mass[idx], beadtype=self.inner_types[idx])
# S - A bond
elif self.inner_types[idx - 1] == 'S' and self.inner_types[idx] == 'A':
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.0, 0.50 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
# A - A bond
elif self.inner_types[idx - 1] == 'A' and self.inner_types[idx] == 'A':
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.0, 1.0 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
elif self.inner_types[idx] == 'B':
if self.inner_types[idx - 1] == 'A':
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.0, 0.50 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
# special case where sticky ends are bound to ssDNA
elif self.inner_types[idx - 1] == 'S':
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.0, 0.25 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
# s
elif self.inner_types[idx - 4] == 'B':
self.beads += CGBead.Bead(position=self.positions[-4, :] + [0.0, 0.0, 0.50 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
# 3 scenarios, 1st flanking for C, 2nd flanking for C or chain end FL
elif self.inner_types[idx] == 'FL':
if self.inner_types[idx - 2] == 'B':
self.beads += CGBead.Bead(position=self.positions[-2, :] + [0.3 * base_length, 0.4 * base_length, 0.0],
mass=self.mass[idx], beadtype=self.inner_types[idx])
elif self.inner_types[idx - 3] == 'B':
self.beads += CGBead.Bead(
position=self.positions[-3, :] + [-0.3 * base_length, 0.4 * base_length, 0.0],
mass=self.mass[idx], beadtype=self.inner_types[idx])
# chain end flanking
elif self.inner_types[idx - 4] == 'B':
self.beads += CGBead.Bead(position=self.positions[-4, :] + [0.0, 0.0, 0.30 * base_length],
mass=self.mass[idx], beadtype=self.inner_types[idx])
elif self.inner_types[idx] in self.sticky_end:
self.beads += CGBead.Bead(position=self.positions[-1, :] + [0.0, 0.40 * base_length, 0.0],
mass=self.mass[idx], beadtype=self.inner_types[idx])
X3DNA_LICENSE = """
3DNA is a suite of software programs for the analysis,
rebuilding and visualization of 3-Dimensional Nucleic Acid
structures. Permission to use, copy, modify, and distribute
this suite for any purpose, with or without fee, is hereby
granted, and subject to the following conditions:
At least one of the 3DNA papers must be cited, including the
following two primary ones:
1. Lu, X. J., & Olson, W. K. (2003). "3DNA: a software
package for the analysis, rebuilding and visualization
of three-dimensional nucleic acid structures." Nucleic
Acids Research, 31(17), 5108-5121.
2. Lu, X. J., & Olson, W. K. (2008). "3DNA: a versatile,
integrated software system for the analysis,
rebuilding and visualization of three-dimensional
nucleic-acid structures." Nature Protocols, 3(7),
1213-1227.
THE 3DNA SOFTWARE IS PROVIDED "AS IS", WITHOUT EXPRESSED OR
IMPLIED WARRANTY OF ANY KIND.
Any 3DNA-related questions, comments, and suggestions are
welcome and should be directed to the open 3DNA Forum
(http://forum.x3dna.org/).
"""
DNA3SPN_LICENSE = """
### Oct. 12, 2015 ###
The GCGI/ directory was added, which contains necessary files for using
the general coarse- grained ions model presented by Hinckley and de Pablo,
JCTC (2015) (DOI:10.1021/acs.jctc.5b00341).
The following website will have information regarding on-going improvements to
the 3SPN.2 model, as well as recent publications by the de Pablo group and
others.
http://ime.uchicago.edu/de_pablo_lab/research/dna_folding_and_hybridization/3spn.2/
"""
[docs]class DNA3SPNChain(LinearChain):
_LICENSE_DISPLAY = True
_MOL_INDEX = 1 # this has chain index that need to be incremented for each molecule. Check do_on_copy() for usage
_MASS_TABLE_AMU = {'Phos': 94.97,
'Sug': 83.11,
'Ade': 134.10,
'Thy': 125.10,
'Cyt': 110.10,
'Gua': 150.10}
_ALL_TYPES = ['Phos', 'Sug', 'Ade', 'Thy', 'Gua', 'Cyt']
_BASE_TYPES = ['Ade', 'Thy', 'Gua', 'Cyt']
_EXCL_DIAM_LIST = {'Phos': 4.5,
'Sug': 6.2,
'Ade': 5.4,
'Thy': 7.1,
'Gua': 4.9,
'Cyt': 6.4}
INTERMEDIATE_FILES_DIR = '.'
try:
X3DNA_REL_PATH = os.path.abspath(os.environ['PBS_O_WORKDIR'])
except KeyError:
X3DNA_REL_PATH = os.path.abspath('./')
print('Could not find PBS_O_WORKDIR variable, setting X3 DNA path to ' + X3DNA_REL_PATH + '/DNA3SPN/X3DNA')
print('X3DNA is only required for the DNA 3SPN chain model '
'and path can manually be set using LinearChain.Dna3spnChain.X3DNA_REL_PATH')
def __init__(self, ss_dna_sequence=None, ds_dna_sequence=None, sticky_sequence=None, flexor=None, mpi_rank=None,
mpi_sync_method=None, intermediate_files_dir=None, **kwargs):
# Set the units
_units = SimUnits()
_units.set_energy('kJ/mol')
_units.set_length('Ang')
_units.set_mass('amu')
# Initialize the data structure for a linear chain
LinearChain.__init__(self,
n_monomer=(len(ss_dna_sequence) + len(ds_dna_sequence) + len(sticky_sequence)) * 3,
kuhn_length=0.00,
units=_units)
# set where intermediate files will be dumped, if passed in
if intermediate_files_dir is not None:
self.intermediate_files_dir = intermediate_files_dir
else:
self.intermediate_files_dir = DNA3SPNChain.INTERMEDIATE_FILES_DIR
# Check that the files necessary to build a DNA chain are present
if not os.path.isfile(DNA3SPNChain.X3DNA_REL_PATH + '/DNA3SPN/X3DNA/bin/x3dna_utils'):
print('X3 DNA package is required to build 3SPN DNA model')
raise IOError('Could not find X3 DNA package at '
+ DNA3SPNChain.X3DNA_REL_PATH + '/DNA3SPN/X3DNA/bin/x3dna_utils')
if not os.access(DNA3SPNChain.X3DNA_REL_PATH + '/DNA3SPN/X3DNA/bin/x3dna_utils', os.X_OK):
print('Packages requires execute permission on /bin/x3dna_utils of X3DNA')
raise IOError('x3dna_utils reported not having os.X_OK bit set')
if not os.access(DNA3SPNChain.X3DNA_REL_PATH + '/DNA3SPN/X3DNA/bin/rebuild', os.X_OK):
print('Packages requires execute permissions on /bin/rebuild of X3DNA')
raise IOError('rebuild reported not having os.X_OK bit set')
# Display the license information the first time the code is run
if DNA3SPNChain._LICENSE_DISPLAY:
print('This model uses X3 DNA to build the representation')
print(X3DNA_LICENSE)
print('This model makes use of parts of the USER-3SPN2 package')
print(DNA3SPN_LICENSE)
print('LICENSE information is displayed only once')
DNA3SPNChain._LICENSE_DISPLAY = False
# check that sequences are composed of the characters ATCG
try:
for b in ss_dna_sequence:
if b not in ['A', 'T', 'C', 'G']:
raise SyntaxError('Sequences of bases require iterables of A, T, C, G')
except TypeError:
raise TypeError('squences must be of iterable types')
try:
for b in ds_dna_sequence:
if b not in ['A', 'T', 'C', 'G']:
raise SyntaxError('Sequences of bases require iterables of A, T, C, G')
except TypeError:
raise TypeError('squences must be of iterable types')
try:
for b in sticky_sequence:
if b not in ['A', 'T', 'C', 'G']:
raise SyntaxError('Sequences of bases require iterables of A, T, C, G')
except TypeError:
raise TypeError('squences must be of iterable types')
if flexor is not None:
raise NotImplementedError('flexors will be implemented in a future version')
# Copy the sequences
self.ss_sequence = [nucleotide for nucleotide in ss_dna_sequence]
self.ds_sequence = [nucleotide for nucleotide in ds_dna_sequence]
self.sticky_sequence = [self.get_complement(nucleotide) for nucleotide in sticky_sequence]
self.three_prime_end = 0
self.five_prime_end = 0
# If this is a serial simulation or the first rank of an MPI simulation,
# create an atomistic structure of the DNA
if mpi_rank is None or mpi_rank == 0:
# becaues write_dna_structure call other people's files that
# assume that intermediate files will be saved in current dir,
# hack a little to be able to dump intermediate files elsewhere
cur_dir = os.getcwd()
os.chdir(self.intermediate_files_dir)
self.write_dna_structure()
os.chdir(cur_dir)
if mpi_sync_method is not None:
mpi_sync_method()
# Set force-field information
self.b_types += Composite.BondType('Phos-Sug', {'K2': 0.6, 'K4': 60.0, 'r0': 3.899},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.b_types += Composite.BondType('Sug-Phos', {'K2': 0.6, 'K4': 60.0, 'r0': 3.559},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.b_types += Composite.BondType('Sug-Ade', {'K2': 0.6, 'K4': 60.0, 'r0': 4.670},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.b_types += Composite.BondType('Sug-Thy', {'K2': 0.6, 'K4': 60.0, 'r0': 4.189},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.b_types += Composite.BondType('Sug-Gua', {'K2': 0.6, 'K4': 60.0, 'r0': 4.829},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.b_types += Composite.BondType('Sug-Cyt', {'K2': 0.6, 'K4': 60.0, 'r0': 3.844},
{'K2': 'E/LL', 'K4': 'E/LLLL', 'r0': 'L'})
self.a_types += Composite.AngleType('Sug-Phos-Sug', {'k': 200.0, 't0': 94.49}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Phos-Sug-Phos', {'k': 200.0, 't0': 120.15}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Phos-Sug-Ade', {'k': 200.0, 't0': 103.53}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Phos-Sug-Thy', {'k': 200.0, 't0': 92.06}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Phos-Sug-Gua', {'k': 200.0, 't0': 107.40}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Phos-Sug-Cyt', {'k': 200.0, 't0': 103.79}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Ade-Sug-Phos', {'k': 200.0, 't0': 112.07}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Thy-Sug-Phos', {'k': 200.0, 't0': 116.68}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Gua-Sug-Phos', {'k': 200.0, 't0': 110.12}, {'k': 'E', 't0': '1'})
self.a_types += Composite.AngleType('Cyt-Sug-Phos', {'k': 200.0, 't0': 110.33}, {'k': 'E', 't0': '1'})
e_tab = np.array([[14.39, 14.34, 13.25, 14.51], [10.37, 13.36, 10.34, 12.89], [14.81, 15.57, 14.93, 15.39],
[11.42, 12.79, 10.52, 13.24]])
r0_tab = np.array([[3.716, 3.675, 3.827, 3.744], [4.238, 3.984, 4.416, 4.141], [3.576, 3.598, 3.664, 3.635],
[4.038, 3.798, 4.208, 3.935]])
t0_tab = np.array(
[[101.15, 85.94, 105.26, 89.00], [101.59, 89.50, 104.31, 91.28], [100.89, 84.83, 105.48, 88.28],
[106.49, 93.31, 109.54, 95.46]])
t0_tab *= pi / 180.0
for i in range(4):
for j in range(4):
self.a_types += Composite.AngleType('Sug-' + str(self._BASE_TYPES[i]) + '-' + str(self._BASE_TYPES[j]),
{'epsilon': e_tab[i][j], 'alpha': 3.0, 'r0': r0_tab[i][j],
't0': t0_tab[i][j]},
{'epsilon': 'E', 'alpha': '1/L', 'r0': 'L', 't0': '1'})
self.d_types += Composite.DihedralType('Phos-Sug-Phos-Sug', {'kphi': 6.0, 'phi0': -154.79, 'sigma_phi': 0.30},
{'kphi': 'E', 'phi0': '1', 'sigma_phi': 0.30})
self.d_types += Composite.DihedralType('Sug-Phos-Sug-Phos', {'kphi': 6.0, 'phi0': -179.17, 'sigma_phi': 0.30},
{'kphi': 'E', 'phi0': '1', 'sigma_phi': 0.30})
# Use 3SPN to coarse-grain the atomistic structure and generate a structure and topology
self.subsystem = DNA3SPN.pdb2cg_dna.System()
DNA3SPN.pdb2cg_dna.read_pdb(os.path.join(os.path.abspath(self.intermediate_files_dir),
'atomistic.pdb'),
self.subsystem)
self.subsystem.do_coarse_graining()
self.subsystem.generate_topology()
# Convert this 3SPN topology to a hoobas topology
self.generate_from_3spn_topologies()
# the hoomd code swaps A-T and C-G parameters depending on quat.x
self.set_quaternion_inversion()
# set charges
self.set_charge_phosphates()
# move the chain upwards in Z
self.set_position_array()
self.generate_3base_bonds()
# the 3' end is bound to the NP and not the 5' end, we need to reverse the chain order, this is just semantics
# self.swap_ends()
self.zero_pos = np.array([0.0, 0.0, 5.0])
self.strict_topology = True
self.set_masses()
def write_dna_structure(self):
# Write the configuration file
DNA3SPN.make_bp_params.write_bp_parameters(self.ss_sequence + self.ds_sequence + self.sticky_sequence)
# Modify the parameters to use B-DNA
# Make sure the X3DNA environmental variable is set so it can find itself
try:
os.environ['X3DNA']
except KeyError:
os.environ['X3DNA'] = os.path.join(DNA3SPNChain.X3DNA_REL_PATH, 'DNA3SPN/X3DNA/')
parameter_modifier = os.path.join(DNA3SPNChain.X3DNA_REL_PATH, 'DNA3SPN/X3DNA/bin/x3dna_utils')
use_b_dna_command = '{} {} {}'.format(parameter_modifier, 'cp_std', 'BDNA')
os.system(use_b_dna_command)
# Build the DNA structure and write it a file
dna_builder = os.path.join(DNA3SPNChain.X3DNA_REL_PATH, 'DNA3SPN/X3DNA/bin/rebuild')
parameter_path = os.path.join(os.path.abspath('.'), 'bp_step.par')
dna_structure_path = os.path.join(os.path.abspath('.'), 'atomistic.pdb')
build_dna_command = '{} -atomic {} {}'.format(dna_builder, parameter_path, dna_structure_path)
os.system(build_dna_command)
return {'dna_parameter_path': parameter_path,
'dna_structure_path': dna_structure_path}
def generate_from_3spn_topologies(self):
# convert subsystem topology to local topology
# note that subsystem topology always run 5' to 3', with beads Sugar - Base - Phosphate - Sugar - ....
# The system we want looks something like this:
# ssDNA | dsDNA | sticky-end
# Chain 5: 5' ACTGACTGACTG ACTGACTGACTGACTGACTGACTGACTG 3'
# Chain 3: 3' TGACTGACTGACTGACTGACTGACTGAC ACTGACTGACTGACTG 5'
#
# However, it is easiest to initialize the strands as only double-stranded, ie
# Chain 5: 5' ACTGACTGACTG ACTGACTGACTGACTGACTGACTGACTG ACTGACTGACTGACTG 3'
# Chain 3: 3' TGACTGACTGAC TGACTGACTGACTGACTGACTGACTGAC TGACTGACTGACTGAC 5'
#
# We need to excise the unwanted parts of the DNA
beadtype_dict = {'S': 'Sug',
'P': 'Phos',
'A': 'Ade',
'T': 'Thy',
'C': 'Cyt',
'G': 'Gua'}
# Each particle has an index associated with it.
# track the indices for the 5' --> 3' chain, and 3' --> 5' chain separately
# These two indices contain the LAMMPS indices, not the the hoobas indices
chain5to3_indices = []
chain3to5_indices = []
chain5to3 = self.subsystem.molecules[0]
chain3to5 = self.subsystem.molecules[1]
num_bases_chain5to3 = len(self.ss_sequence) + len(self.ds_sequence)
num_bases_chain3to5 = len(self.sticky_sequence) + len(self.ds_sequence)
chain3to5_offset = len(self.subsystem.molecules[0].beads)
# Chain 5' --> 3'
# Copy the beads from the cg model, excluding the sticky end
first_sticky_bead_idx = min(3 * num_bases_chain5to3, # End of ss-DNA / ds-DNA
len(chain5to3.beads)) # End of chain
for local_bead in chain5to3.beads[0:first_sticky_bead_idx]:
# The chain_index represents the chain for which a base belongs
# If the bead is a dna base, it is +1 for the 5' --> 3' chain, and -1 for the 3' --> 5' chain
# else: it is 0 (not a dna base)
# If the dna chain is copied, these numbers will be incremented by +1 and -1 for the 5' --> 3'
# and 3' --> 5' chains respectively
if beadtype_dict[local_bead.beadtype] in ['Ade', 'Thy', 'Cyt', 'Gua']:
chain_index = 1.0
else:
chain_index = 0.0
self.beads += CGBead.Bead(position=[local_bead.x, local_bead.y, local_bead.z],
charge=local_bead.charge,
beadtype=beadtype_dict[local_bead.beadtype],
diameter=chain_index)
chain5to3_indices.append(local_bead.beadnumber)
num_beads_chain5to3 = len(self.beads)
self.three_prime_end = len(self.beads) - 1
# Chain 3' --> 5'
# Copy the second strand of DNA, excluding the ssDNA and including the sticky end
# The chain starts at the 5' end, so the sticky end bp come first, then the dsDNA
# this cuts the unused section of the 5' - 3' sense strand (removes the ssDNA sequence)
# Exclude the last bead of the 3' --> 5' strand
for idx in range(3 * len(self.sticky_sequence) + 3 * len(self.ds_sequence)):
lb = self.subsystem.molecules[1].beads
if beadtype_dict[lb[idx].beadtype] in ['Ade', 'Thy', 'Cyt', 'Gua']:
chain_index = -1.0
else:
chain_index = 0.0
self.beads += CGBead.Bead(position=[lb[idx].x, lb[idx].y, lb[idx].z],
charge=lb[idx].charge,
beadtype=beadtype_dict[lb[idx].beadtype],
diameter=chain_index)
chain3to5_indices.append(lb[idx].beadnumber + chain3to5_offset) # who the hell offset half his lists
# Set the ID for each bead
self.p_num = [particle_id for particle_id in range(len(self.beads))]
# Translate the bonds from lammps
# Chain 5' --> 3'
for local_bond in chain5to3.bonds:
# Add the bond if
# 1) we kept both beads
# 2) the type is not -1
if local_bond.stea in chain5to3_indices and local_bond.steb in chain5to3_indices and local_bond.type != -1:
bead1_idx = local_bond.stea - 1
bead2_idx = local_bond.steb - 1
bond_name = "{}-{}".format(self.beads[bead1_idx].beadtype, self.beads[bead2_idx].beadtype)
new_bond = Composite.Bond([bond_name, bead1_idx, bead2_idx])
self.bonds += new_bond
# Chain 3' --> 5'
for local_bond in chain3to5.bonds:
if local_bond.stea in chain3to5_indices and local_bond.steb in chain3to5_indices and local_bond.type != -1:
bead1_idx = chain3to5_indices.index(local_bond.stea) + num_beads_chain5to3
bead2_idx = chain3to5_indices.index(local_bond.steb) + num_beads_chain5to3
bond_name = "{}-{}".format(self.beads[bead1_idx].beadtype, self.beads[bead2_idx].beadtype)
new_bond = Composite.Bond([bond_name, bead1_idx, bead2_idx])
self.bonds += new_bond
# Angles
# Chain 5' --> 3'
for local_angle in chain5to3.bends:
if local_angle.stea in chain5to3_indices and local_angle.steb in chain5to3_indices \
and local_angle.stec in chain5to3_indices and local_angle.type != -1:
bead1_idx = local_angle.stea - 1
bead2_idx = local_angle.steb - 1
bead3_idx = local_angle.stec - 1
angle_name = '{}-{}-{}'.format(self.beads[bead1_idx].beadtype,
self.beads[bead2_idx].beadtype,
self.beads[bead3_idx].beadtype)
self.angles += [angle_name, bead1_idx, bead2_idx, bead3_idx]
# Chain 3' --> 5'
for local_angle in chain3to5.bends:
if local_angle.stea in chain3to5_indices and local_angle.steb in chain3to5_indices\
and local_angle.stec in chain3to5_indices and local_angle.type != -1:
bead1_idx = chain3to5_indices.index(local_angle.stea) + num_beads_chain5to3
bead2_idx = chain3to5_indices.index(local_angle.steb) + num_beads_chain5to3
bead3_idx = chain3to5_indices.index(local_angle.stec) + num_beads_chain5to3
angle_name = '{}-{}-{}'.format(self.beads[bead1_idx].beadtype,
self.beads[bead2_idx].beadtype,
self.beads[bead3_idx].beadtype)
self.angles += [angle_name, bead1_idx, bead2_idx, bead3_idx]
# Dihedrals
# Chain 5' --> 3'
for local_diehedral in chain5to3.torsions:
if local_diehedral.stea in chain5to3_indices and \
local_diehedral.steb in chain5to3_indices and \
local_diehedral.stec in chain5to3_indices and \
local_diehedral.sted in chain5to3_indices and \
local_diehedral.type != -1:
bead1_idx = local_diehedral.stea - 1
bead2_idx = local_diehedral.steb - 1
bead3_idx = local_diehedral.stec - 1
bead4_idx = local_diehedral.sted - 1
dihedral_name = '{}-{}-{}-{}'.format(self.beads[bead1_idx], self.beads[bead2_idx],
self.beads[bead3_idx], self.beads[bead4_idx])
self.dihedrals += [dihedral_name, bead1_idx, bead2_idx, bead3_idx, bead4_idx]
# Chain 3' --> 5'
for local_diehedral in chain3to5.torsions:
if local_diehedral.stea in chain3to5_indices and \
local_diehedral.steb in chain3to5_indices and \
local_diehedral.stec in chain3to5_indices and \
local_diehedral.sted in chain3to5_indices and \
local_diehedral.type != -1:
bead1_idx = chain3to5_indices.index(local_diehedral.stea) + num_beads_chain5to3
bead2_idx = chain3to5_indices.index(local_diehedral.steb) + num_beads_chain5to3
bead3_idx = chain3to5_indices.index(local_diehedral.stec) + num_beads_chain5to3
bead4_idx = chain3to5_indices.index(local_diehedral.sted) + num_beads_chain5to3
dihedral_name = '{}-{}-{}-{}'.format(self.beads[bead1_idx], self.beads[bead2_idx],
self.beads[bead3_idx], self.beads[bead4_idx])
self.dihedrals += [dihedral_name, bead1_idx, bead2_idx, bead3_idx, bead4_idx]
def generate_3base_bonds(self):
# Using the diameter attribute to represent the chain each bead belongs to is confusing
# Create a helper function to make this explicit
chain = operator.attrgetter('diameter')
# generate 'dummy' bonds between bases less than 4 away to avoid intra-strand base pairing
num_beads = len(self.beads)
for idx, bead in enumerate(self.beads):
# Check that it is a nucleotide and not a sugar/phosphate
if str(bead) in ['Ade', 'Thy', 'Cyt', 'Gua']:
# Check the next three neighbors
for offset in range(1, 4):
next_bead_idx = idx + 3 * offset
# Make sure we are still on the same chain
if next_bead_idx < num_beads and chain(bead) == chain(self.beads[next_bead_idx]):
self.bonds += ['dummy', idx, next_bead_idx]
def set_quaternion_inversion(self):
for bead in self.beads:
if str(bead) == 'Thy':
bead.orientation = np.array([-1.0, 0.0, 0.0, 0.0], dtype=np.single)
if str(bead) == 'Gua':
bead.orientation = np.array([0.5, 0.0, 0.0, 0.0], dtype=np.single)
if str(bead) == 'Ade':
bead.orientation = np.array([1.0, 0.0, 0.0, 0.0], dtype=np.single)
if str(bead) == 'Cyt':
bead.orientation = np.array([-0.50, 0.0, 0.0, 0.0], dtype=np.single)
# Using the diameter attribute to represent the chain each bead belongs to is confusing
# Create a helper function to make this explicit
chain = operator.attrgetter('diameter')
# find terminal bases
for idx, bead in enumerate(self.beads):
if str(bead) in DNA3SPNChain._BASE_TYPES:
if 6 <= idx <= (len(self.beads) - 1 - 6):
this_chain = chain(bead)
if this_chain == chain(self.beads[idx - 6]) and this_chain == chain(self.beads[idx + 6]):
bead.orientation += np.array([0.0, 1.0, 1.0, -1.0], dtype=np.single)
elif this_chain == chain(self.beads[idx - 6]) and this_chain == chain(self.beads[idx + 3]):
bead.orientation += np.array([0.0, -1.0, 1.0, -1.0])
elif this_chain == chain(self.beads[idx - 6]) and this_chain != chain(self.beads[idx + 3]):
bead.orientation += np.array([0.0, -1.0, 1.0, 1.0])
elif this_chain == chain(self.beads[idx + 6]) and this_chain == chain(self.beads[idx - 3]):
bead.orientation += np.array([0.0, 1.0, -1.0, -1.0])
elif this_chain == chain(self.beads[idx + 6]) and this_chain != chain(self.beads[idx - 3]):
bead.orientation += np.array([0.0, 1.0, -1.0, 1.0])
# First base
elif 0 <= idx < 3:
bead.orientation += np.array([0.0, 1.0, -1.0, 1.0])
# Second base
elif 3 <= idx < 6:
bead.orientation += np.array([0.0, 1.0, -1.0, -1.0])
# Last base
elif (len(self.beads) - 1 - 3) < idx <= (len(self.beads) - 1):
bead.orientation += np.array([0.0, -1.0, 1.0, 1.0])
# Second to last base
elif (len(self.beads) - 1 - 6) < idx <= (len(self.beads) - 1 - 3):
bead.orientation += np.array([0.0, -1.0, 1.0, -1.0])
else:
raise ArithmeticError('Bad number of DNA bases somewhere')
def set_charge_phosphates(self):
for bead in self.beads:
if str(bead) == 'Phos':
bead.charge = -1.0
def set_mass_particles(self):
for bead in self.beads:
bead.mass = DNA3SPNChain._MASS_TABLE_AMU[str(bead)]
def do_on_copy(self):
for bead in self.beads:
bead.diameter *= DNA3SPNChain._MOL_INDEX
DNA3SPNChain._MOL_INDEX += 1
def set_position_array(self):
self.positions = np.zeros((len(self.beads), 3), dtype=np.single)
for idx in range(len(self.beads)):
self.positions[idx, :] = self.beads[idx].position
@staticmethod
def force_field_types():
return DNA3SPNChain._ALL_TYPES
@staticmethod
def force_field_excluded_volumes():
return DNA3SPNChain._EXCL_DIAM_LIST
@staticmethod
def force_field_amus():
return DNA3SPNChain._MASS_TABLE_AMU
@staticmethod
def get_excl_sig():
# 0 1 2 3 4 5
# _ALL_TYPES = ['Phos', 'Sug', 'Ade', 'Thy', 'Gua', 'Cyt']
for ite1 in range(len(DNA3SPNChain._ALL_TYPES)):
for ite2 in range(ite1, len(DNA3SPNChain._ALL_TYPES)):
if ite1 == 2 and ite2 == 3 or ite1 == 4 and ite2 == 5:
yield (DNA3SPNChain._ALL_TYPES[ite1],
DNA3SPNChain._ALL_TYPES[ite2],
0.0,
0.0)
else:
sig1 = DNA3SPNChain._EXCL_DIAM_LIST[DNA3SPNChain._ALL_TYPES[ite1]]
sig2 = DNA3SPNChain._EXCL_DIAM_LIST[DNA3SPNChain._ALL_TYPES[ite2]]
yield (DNA3SPNChain._ALL_TYPES[ite1],
DNA3SPNChain._ALL_TYPES[ite2],
(sig1 + sig2) * 0.5,
(sig1 + sig2) * 0.5 * 2.0 ** (1.0 / 6.0))
def randomize_dirs(self, tol=1e-1): # can't randomize this thing or it'll have insane topology
pass
def set_masses(self):
for bead in self.beads:
bead.mass = DNA3SPNChain._MASS_TABLE_AMU[str(bead)]
[docs] def swap_ends(self):
"""
swaps ends of the chain to bind to the three prime end instead of 5 prime
:return:
"""
self.rotate(Util.get_v1_v2_rotmat([0., 0., 1.], [0., 0., -1.]))
# Get the index of the 3' bead
old_3p_idx = self.three_prime_end
new_5p_idx = old_3p_idx
# Get the index of the 5' bead
old_5p_idx = self.five_prime_end
new_3p_idx = old_5p_idx
# Swap the beads at these locations
self.beads[old_3p_idx], self.beads[old_5p_idx] = self.beads[old_5p_idx], self.beads[old_3p_idx]
#
new_chain_5to3, new_chain_3to5 = (copy.deepcopy(self.positions[old_5p_idx, :]),
copy.deepcopy(self.positions[old_3p_idx, :]))
self.positions[new_5p_idx, :] = new_chain_5to3
self.positions[new_3p_idx, :] = new_chain_3to5
# self.p_num[old_3p_idx], self.p_num[old_5p_idx] = self.p_num[old_5p_idx], self.p_num[old_3p_idx]
for topo_item in self.bonds + self.angles + self.dihedrals + self.impropers:
for index in range(1, len(topo_item)):
if topo_item[index] == old_3p_idx:
topo_item[index] = new_3p_idx
elif topo_item[index] == old_5p_idx:
topo_item[index] = new_5p_idx
@staticmethod
def get_complement(base):
base_pairs = {'A': 'T',
'T': 'A',
'C': 'G',
'G': 'C'}
return base_pairs[base]
[docs]class BeadMonomer(LinearChain):
def __init__(self, Bead, units=None):
super(BeadMonomer, self).__init__(n_monomer=1, kuhn_length=1.0, units=units)
self.beads += copy.deepcopy(Bead)
self.p_num.append(0)
[docs]class PolymerBySequence(LinearChain):
def __init__(self, sequence, units=None):
super(PolymerBySequence, self).__init__(n_monomer=len(sequence), kuhn_length=1.0, units=units)
for index, item in enumerate(sequence):
self.beads += CGBead.Bead(beadtype=item, position=np.array([0.0, 0.0, 1.0 * index]), mass=1.0)
self.p_num.append(index)
if index > 0:
self.bonds += ['PolymerBySequenceBond:'+str(item)+'-'+str(self.beads[index-1]), index-1, index]
[docs]class RandomPolymer(LinearChain):
"""
General random polymers comprised of arbitrary :py:class:`hoobas.Composite.CompositeObject` monomers
Args:
n_mono: Number of monomers to assemble. This parameter can be an integer or function. The function should not take any arguments
monomers (dict): list of monomers that will be used for building the chain. If any of the monomers are beads, they will be converted to BeadMonomer objects
distribution (function): function that determines the next monomer to be added. Extra and superseding behavior
-None: Entirely random distribution
-foo(\**kwargs): queries kwargs from self and passes them in foo
-foo() returns a CompositeObject: use the composite object as the next monomer
-foo() returns None: terminate building for this chain
"""
def __init__(self, n_mono, monomers=None, distribution=None, units=None):
super(RandomPolymer, self).__init__(n_monomer=0, units=units)
self.monomer_list = monomers
if self.monomer_list is not None:
for key, value in self.monomer_list.items():
if value.__class__ == CGBead.Bead:
self.monomer_list[key] = BeadMonomer(copy.deepcopy(value), units)
self.dist = distribution
self.max_monomers = n_mono
def build(self):
self.is_built = True
if hasattr(self.max_monomers, '__call__'):
self.max_monomers = self.max_monomers()
prev_monomer = None
while self.number_monomers < self.max_monomers:
# parse the next monomer to be added
if self.dist is None:
next_monomer = self.monomer_list[random.choice(list(self.monomer_list.keys()))]
else:
# parse arguments for distribution
dist_return = self.dist(**self.resolve_args(self.dist))
if dist_return is None:
break
elif issubclass(dist_return.__class__, Composite.CompositeObject):
next_monomer = dist_return
else:
next_monomer = self.monomer_list[dist_return]
# check whether the next monomer is a function (constructor was passed in)
if hasattr(next_monomer, '__call__'):
ncp = copy.deepcopy(next_monomer())
else:
ncp = copy.deepcopy(next_monomer)
# if this is the first monomer, just merge it
if self.number_monomers == 0:
self.merge(ncp)
prev_monomer = next_monomer
self.number_monomers += 1
self.att_list += ncp.att_list
continue
# make a copy of the next monomer and graft it. Note that since the object is in built state, this is done
# immediately so that representations of build() are finished before any deferred call is done
n = self.graft_external_objects(ncp,
1,
connecting_topology='Placeholder',
merge_attachment_sites=True,
strict_compliance=True)
if n == 0: # we failed to graft anything due to constraints
continue
# remap the bond to something very generic
bond_indexes = self.bonds[-1].topology_tuple
beadnameA = str(self.beads[bond_indexes[0]])
beadnameB = str(self.beads[bond_indexes[1]])
self.remap_bond('Placeholder', 'RandomPolymer:' + str(prev_monomer)+'-'+str(next_monomer)+':'+beadnameA+'-'+beadnameB)
self.number_monomers += 1
[docs]class PolyStyreneSulfonateChain(RandomPolymer):
def __init__(self, n_monomers, tacticity_distribution=None):
monomer_dictionary = {'+': PolystyreneSulfonateMonomer('+'), '-': PolystyreneSulfonateMonomer('-')}
if tacticity_distribution == 'isotactic':
tacticity_distribution = lambda: '+'
elif tacticity_distribution == 'syntactic':
tacticity_distribution = lambda subunits: '+' if len(subunits) == 0 or str(subunits[-1]) == '-' else '-'
super(PolyStyreneSulfonateChain, self).__init__(n_mono=n_monomers,
monomers=monomer_dictionary,
distribution=tacticity_distribution)