Source code for arbdmodel.ssdna_two_bead

# -*- coding: utf-8 -*-
import numpy as np
import sys

## Local imports
from .logger import logger,  get_resource_path
from . import ParticleType, PointParticle, Group
from .polymer import PolymerBeads, PolymerModel
from .interactions import TabulatedNonbonded, HarmonicBond, HarmonicAngle, HarmonicDihedral

"""Define particle types"""

## units "295 k K/(160 amu * 1.24/ps)" "AA**2/ns"
## units "295 k K/(180 amu * 1.24/ps)" "AA**2/ns"
_P = ParticleType("P",
                 diffusivity = 1621,
                 mass = 121,
                 radius = 5,
                 nts = 0.5      # made compatible with nbPot
)

_B = ParticleType("B",
                 diffusivity = 1093,
                 mass = 181,    # thymine
                 radius = 3,                 
                 nts = 0.5      # made compatible with nbPot
)

[docs] class DnaStrandBeads(PolymerBeads): p = PointParticle(_P, (0,0,0), "P") b = PointParticle(_B, (3,0,1), "B") nt = Group( name = "nt", children = [p,b]) nt.add_bond( i=p, j=b, bond = get_resource_path('two_bead_model/BPB.dat'), exclude = True ) def __init__(self, polymer, sequence=None, **kwargs): self.polymer = polymer if sequence is not None and sequence != 'T'*len(sequence): logger.warn('Sequence given to DNA model does not appear to be poly(dT)') PolymerBeads.__init__(self, polymer, sequence, **kwargs) assert(self.monomers_per_bead_group == 1) def _generate_ith_bead_group(self, i, r, o): new = DnaStrandBeads.nt.duplicate() new.orientation = o new.position = r new.resid=i+1 return new def _join_adjacent_bead_groups(self, ids): nts = self.children ## Two consecutive groups if len(ids) == 2: p1,b1 = nts[ids[0]].children p2,b2 = nts[ids[1]].children self.add_bond( i=b1, j=p2, bond = get_resource_path('two_bead_model/BBP.dat'), exclude=True ) self.add_bond( i=p1, j=p2, bond = get_resource_path('two_bead_model/BPP.dat'), exclude=True ) self.add_angle( i=p1, j=p2, k=b2, angle = get_resource_path('two_bead_model/p1p2b2.dat') ) self.add_angle( i=b1, j=p2, k=b2, angle = get_resource_path('two_bead_model/b1p2b2.dat') ) self.add_dihedral( i=b1, j=p1, k=p2, l=b2, dihedral = get_resource_path('two_bead_model/b1p1p2b2.dat') ) self.add_exclusion( i=b1, j=b2 ) self.add_exclusion( i=p1, j=b2 ) elif len(ids) == 3: p1,b1 = nts[ids[0]].children p2,b2 = nts[ids[1]].children p3,b3 = nts[ids[2]].children self.add_angle( i=p1, j=p2, k=p3, angle = get_resource_path('two_bead_model/p1p2p3.dat') ) self.add_angle( i=b1, j=p2, k=p3, angle = get_resource_path('two_bead_model/b1p2p3.dat') ) self.add_dihedral( i=b1, j=p2, k=p3, l=b3, dihedral = get_resource_path('two_bead_model/b1p2p3b3.dat') ) self.add_exclusion( i=p1, j=p3 ) self.add_exclusion( i=b1, j=p3 ) elif len(ids) == 4: p1,b1 = nts[ids[0]].children p2,b2 = nts[ids[1]].children p3,b3 = nts[ids[2]].children p4,b4 = nts[ids[3]].children self.add_dihedral( i=p1, j=p2, k=p3, l=p4, dihedral = get_resource_path('two_bead_model/p0p1p2p3.dat') ) else: raise Exception("Programming error; should not be here") def _apply_grids(self): for p in self.polymer: ...
# def hybridize(strand1, strand2, parent=None, num_bp=None, start1=None, end2=None): # """ hybridize num_bp basepairs between strand1 and strand2, # starting with nt at start1 and ending with nucleotide and end2 """ # num_nt1 = len(strand1.children) # num_nt2 = len(strand2.children) # if parent is None: # assert(strand1.parent == strand2.parent) # parent = strand1.parent # if start1 is None: # start1 = 0 # if end2 is None: # end2 = num_nt2 # assert(start1 >= 0) # assert(end2 > 0) # assert(start1 < num_nt1) # assert(end2 <= num_nt2) # if num_bp is None: # num_bp = min(num_nt1-start1, end2) # if num_bp > num_nt1-start1: # raise ValueError("Attempted to hybridize too many basepairs ({}) for strand1 ({})".format(num_bp,num_nt1-start1)) # if num_bp > end2: # raise ValueError("Attempted to hybridize too many basepairs ({}) for strand1 ({})".format(num_bp,end2)) # nts1 = strand1.children[start1:start1+num_bp] # nts2 = strand2.children[end2-num_bp:end2][::-1] # assert( len(nts1) == len(nts2) ) # kAngle = 0.0274155 # 90 kcal_mol per radian**2 # ## every bp # for i in range(num_bp): # p11,b11 = nts1[i].children # p21,b21 = nts2[i].children # parent.add_bond( i=b11, j=b21, bond=HarmonicBond(10,7.835) ) # parent.add_angle( i=p11, j=b11, k=b21, angle=HarmonicAngle(kAngle,162.0) ) # parent.add_angle( i=p21, j=b21, k=b11, angle=HarmonicAngle(kAngle,162.0) ) # ## every 2 bp # for i in range(num_bp-1): # p11,b11 = nts1[i].children # p12,b12 = nts1[i+1].children # p21,b21 = nts2[i].children # p22,b22 = nts2[i+1].children # parent.add_bond( i=b11, j=b22, bond=HarmonicBond(1,8.41) ) # parent.add_bond( i=b12, j=b21, bond=HarmonicBond(1,7.96) ) # parent.add_angle( i=p11, j=p12, k=b12, angle=HarmonicAngle(kAngle,87.0) ) # parent.add_angle( i=p22, j=p21, k=b11, angle=HarmonicAngle(kAngle,87.0) ) # ## every 3 bp # for i in range(num_bp-2): # p11,b11 = nts1[i].children # p12,b12 = nts1[i+1].children # p13,b13 = nts1[i+2].children # p21,b21 = nts2[i].children # p22,b22 = nts2[i+1].children # p23,b23 = nts2[i+2].children # parent.add_angle( i=p11, j=p12, k=p13, angle=HarmonicAngle(kAngle,150.0) ) # parent.add_angle( i=p21, j=p22, k=p23, angle=HarmonicAngle(kAngle,150.0) )
[docs] class DnaModel(PolymerModel): def __init__(self, polymers, sequences = None, DEBUG=False, **kwargs): if 'timestep' not in kwargs: kwargs['timestep'] = 20e-6 if 'cutoff' not in kwargs: kwargs['cutoff'] = 35 PolymerModel.__init__(self, polymers, sequences, monomers_per_bead_group=1, **kwargs) self.strands = self.children # make a nice alias self.add_nonbonded_interaction( TabulatedNonbonded(get_resource_path('two_bead_model/NBBB.dat')), typeA=_B, typeB=_B ) self.add_nonbonded_interaction( TabulatedNonbonded(get_resource_path('two_bead_model/NBPB.dat')), typeA=_P, typeB=_B ) self.add_nonbonded_interaction( TabulatedNonbonded(get_resource_path('two_bead_model/NBPP.dat')), typeA=_P, typeB=_P ) self.generate_beads() def _generate_polymer_beads(self, polymer, sequence, polymer_index=None): return DnaStrandBeads( polymer, sequence, polymer_index = polymer_index )
if __name__ == "__main__": """ ## old code strands = [] for i in range(5,60,5): strands.extend( [strand.duplicate() for j in range(int(round(600/i**1.2)))] ) ## Randomly place strands through system model = ArbdModel( strands, dimensions=(200, 200, 200) ) old_schemes = model.nbSchemes model.nbSchemes = [] model.useNonbondedScheme( TabulatedPotential(get_resource_path('two_bead_model/NBBB.dat')), typeA=B, typeB=B ) model.useNonbondedScheme( TabulatedPotential(get_resource_path('two_bead_model/NBPB.dat')), typeA=P, typeB=B ) model.useNonbondedScheme( TabulatedPotential(get_resource_path('two_bead_model/NBPP.dat')), typeA=P, typeB=P ) model.nbSchemes.extend(old_schemes) for s in strands: s.orientation = randomOrientation() s.position = np.array( [(a-0.5)*b for a,b in zip( np.random.uniform(size=3), model.dimensions )] ) model.simulate( output_name = 'many-strands', output_period=1e4, num_steps=1e6 ) """