Added parser for maya files and fixed some issues with segmentmodel

arbdmodel.py

@@ -104,8 +104,9 @@ class Parent():
 
     def add_exclusion(self, i,j):
         ## TODO: how to handle duplicating and cloning bonds
-        beads = [b for b in self]
-        for b in (i,j): assert(b in beads)
+        ## TODO: perform following check elsewhere
+        # beads = [b for b in self]
+        # for b in (i,j): assert(b in beads)
         self.exclusions.append( (i,j) )
 
     def get_bonds(self):
@@ -262,7 +263,7 @@ class ParticleType():
     def _equal_check(a,b):
         if a.name == b.name:
             if a.__hash_key() != b.__hash_key():
-                raise Exception("Badly formed ParticleTypes have 'name' collision")
+                raise Exception("Two different ParticleTypes have same 'name' attribute")
 
     def __eq__(a,b):
         a._equal_check(b)

segmentmodel.py

@@ -64,7 +64,7 @@ class Segment(ConnectableElement, Group):
                                         radius = 1,
                                     )
     # orientation_bond = HarmonicBond(10,2)
-    orientation_bond = HarmonicBond(30,1.5)
+    orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )
 
     ssDNA_particle = ParticleType("S",
                                   diffusivity = 43.5,
@@ -148,12 +148,27 @@ class Segment(ConnectableElement, Group):
         # self._bead_model_max_nts_per_bead = max_nts_per_bead
 
         direction = self.end_position - self.start_position
+        
+
         num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
         nts_per_bead = float(self.num_nts)/num_beads
         twist_per_bead = nts_per_bead * self.twist_per_nt
 
         last = None
 
+        if num_beads <= 2:
+            ## not yet implemented for dsDNA
+            assert( isinstance(self, SingleStrandedSegment) )
+            b = self._generate_one_bead(self.start_position,
+                                        self.num_nts,
+                                        self.start_orientation)
+            self.children.append(b)
+            self.beads.append(b) # don't add orientation bead
+            self._assign_particles_to_locations()
+            return
+
+        assert( np.linalg.norm(direction) > 0 )
+
         for i in range(num_beads+1):
             nts = nts_per_bead
             if i == 0 or i == num_beads: 
@@ -177,7 +192,7 @@ class Segment(ConnectableElement, Group):
             if "orientation_bead" in b.__dict__:
                 o = b.orientation_bead
                 self.children.append(o)
-                self.add_bond(b,o, Segment.orientation_bond)
+                self.add_bond(b,o, Segment.orientation_bond, exclude=True)
                 
             # if last is not None:
             #     self.add_bond( i=last, j=b, bond="ssdna" )
@@ -228,26 +243,26 @@ class DoubleStrandedSegment(Segment):
         self.end3 = Location( self, address=-1, type_ = "end3" )
 
     ## Convenience methods
-    def connect_start5(self, end3, force_connection=False):
+    def connect_start5(self, end3, type_="intrahelical", force_connection=False):
         if isinstance(end3, SingleStrandedSegment):
             end3 = end3.end3
-        self._connect_ends( self.start5, end3, force_connection = force_connection )
-    def connect_start3(self, end5, force_connection=False):
+        self._connect_ends( self.start5, end3, type_, force_connection = force_connection )
+    def connect_start3(self, end5, type_="intrahelical", force_connection=False):
         if isinstance(end5, SingleStrandedSegment):
             end5 = end5.end5
-        self._connect_ends( self.start3, end5, force_connection = force_connection )
-    def connect_end3(self, end5, force_connection=False):
+        self._connect_ends( self.start3, end5, type_, force_connection = force_connection )
+    def connect_end3(self, end5, type_="intrahelical", force_connection=False):
         if isinstance(end5, SingleStrandedSegment):
             end5 = end5.end5
-        self._connect_ends( self.end3, end5, force_connection = force_connection )
-    def connect_end5(self, end3, force_connection=False):
+        self._connect_ends( self.end3, end5, type_, force_connection = force_connection )
+    def connect_end5(self, end3, type_="intrahelical", force_connection=False):
         if isinstance(end3, SingleStrandedSegment):
             end3 = end3.end3
-        self._connect_ends( self.end5, end3, force_connection = force_connection )
+        self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
 
         
     ## Real work
-    def _connect_ends(self, end1, end2, force_connection):
+    def _connect_ends(self, end1, end2, type_, force_connection):
 
         ## validate the input
         for end in (end1, end2):
@@ -255,7 +270,7 @@ class DoubleStrandedSegment(Segment):
             assert( end.type_ in ("end3","end5") )
         assert( end1.type_ != end2.type_ )
 
-        end1.container._connect( end2.container, Connection( end1, end2, type_="intrahelical" ) )
+        end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ) )
 
     def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
         return (self.num_nts // max_basepairs_per_bead) + 1
@@ -340,6 +355,8 @@ class SingleStrandedSegment(Segment):
         self._connect( other.container, Connection( my_end, other, type_="intrahelical" ) )
 
     def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
+        # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
+        #     pdb.set_trace()
         return (self.num_nts // max_nucleotides_per_bead) + 1
 
     def _generate_one_bead(self, pos, nts, orientation=None):
@@ -361,9 +378,9 @@ class SegmentModel(ArbdModel):
                  dimensions=(1000,1000,1000), temperature=291,
                  timestep=50e-6, cutoff=50, 
                  decompPeriod=10000, pairlistDistance=None, 
-                 nonbondedResolution=0):
-
-
+                 nonbondedResolution=0,DEBUG=0):
+        self.DEBUG = DEBUG
+        if DEBUG > 0: print("Building ARBD Model")
         ArbdModel.__init__(self,segments,
                            dimensions, temperature, timestep, cutoff, 
                            decompPeriod, pairlistDistance=None,
@@ -376,6 +393,9 @@ class SegmentModel(ArbdModel):
         self._bonded_potential = dict() # cache bonded potentials
 
         self._generate_bead_model(segments, max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
+
+        self.useNonbondedScheme( nbDnaScheme )
+
     
     def _get_intrahelical_beads(self):
         ret = []
@@ -462,7 +482,7 @@ class SegmentModel(ArbdModel):
         key = (type_,kSpring,d)
         if key not in self._bonded_potential:
             if type_ == "bond":
-                self._bonded_potential[key] = HarmonicBond(kSpring,d, max_potential=max_potential)
+                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,500), max_potential=max_potential)
             elif type_ == "angle":
                 self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                 # , resolution = 1, maxForce=0.1)
@@ -516,6 +536,7 @@ class SegmentModel(ArbdModel):
 
 
         """ Generate beads """
+        if self.DEBUG: print("Generating beads")
         for s in segments:
             if local_twist:
                 s.local_twist = True
@@ -528,9 +549,11 @@ class SegmentModel(ArbdModel):
                     ...
 
         """ Reassign bead types """
+        if self.DEBUG: print("Assigning bead types")
         beadtype_s = dict()
         for segment in segments:
             for b in segment:
+                b.num_nts = np.around( b.num_nts, decimals=1 )
                 key = (b.type_.name[0].upper(), b.num_nts)
                 if key in beadtype_s:
                     b.type_ = beadtype_s[key]
@@ -545,13 +568,16 @@ class SegmentModel(ArbdModel):
                     else:
                         raise Exception("TODO")
                     # print(t.nts)
-                    t.name = t.name + "%03d" % (100*t.nts)
+                    t.name = t.name + "%03d" % (10*t.nts)
                     beadtype_s[key] = b.type_ = t
 
 
         """ Add intrahelical bond potentials """
+        if self.DEBUG: print("Adding intrahelical bond potentials")
         dists = dict()          # built for later use
-        for b1,b2 in self._get_intrahelical_beads():
+        intra_beads = self._get_intrahelical_beads() 
+        if self.DEBUG: print("  Adding %d bonds" % len(intra_beads))
+        for b1,b2 in intra_beads:
             parent = self._getParent(b1,b2)
             if b1.parent == b2.parent:
                 sep = 0.5*(b1.num_nts+b2.num_nts)
@@ -559,7 +585,7 @@ class SegmentModel(ArbdModel):
                 sep = 1
 
             conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
-            if b1.type_.name[0] == "D" and b1.type_.name[0] == "D":
+            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
                 elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
                 d = 3.4*sep
                 k = conversion*elastic_modulus/d
@@ -568,6 +594,7 @@ class SegmentModel(ArbdModel):
                 elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
                 d = 5*sep
                 k = conversion*elastic_modulus/d
+                # print(sep,d,k)
               
             if b1 not in dists:
                 dists[b1] = []
@@ -576,11 +603,16 @@ class SegmentModel(ArbdModel):
             dists[b1].append([b2,sep])
             dists[b2].append([b1,sep])
 
+
+            # if not (b1.type_.name[0] == "D" and b2.type_.name[0] == "D"):
+            #     continue 
+
             # dists[[b1,b2]] = dists[[b2,b1]] = sep
             bond = self.get_bond_potential(k,d)
             parent.add_bond( b1, b2, bond, exclude=True )
 
         """ Add intrahelical angle potentials """
+        if self.DEBUG: print("Adding intrahelical angle potentials")
         for b1,b2,b3 in self._get_intrahelical_angle_beads():
             
             sep = 0
@@ -607,10 +639,14 @@ class SegmentModel(ArbdModel):
                 ## TODO: get correct number from ssDNA model
                 k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2                
             # k *= 1e-6
+            # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
+            #     pdb.set_trace()
+
             angle = self.get_angle_potential(k,180)
             parent.add_angle( b1, b2, b3, angle )
 
         """ Add intrahelical exclusions """
+        if self.DEBUG: print("Adding intrahelical exclusions")
         beads = dists.keys()
         def _recursively_get_beads_within(b1,d,done=[]):
             ret = []
@@ -626,12 +662,15 @@ class SegmentModel(ArbdModel):
         exclusions = set()
         for b1 in beads:
             exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
+
+        if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
         for b1,b2 in exclusions:
             parent = self._getParent(b1,b2)
             parent.add_exclusion( b1, b2 )
 
         """ Twist potentials """
         if local_twist:
+            if self.DEBUG: print("Adding twist potentials")
 
             ## TODO: decide whether to add bond here
             # """Add bonds between orientation bead and parent"""
@@ -667,12 +706,14 @@ class SegmentModel(ArbdModel):
                     pot = self.get_dihedral_potential(k,angle,max_potential=1)
                     parent.add_dihedral(o1,b1,b2,o2, pot)
 
-                
-                
-        
         """ Add connection potentials """
-        ## TODO
-
+        for s in segments:
+            for c in s.connections:
+                if c.type_ == "terminal_crossover":
+                    if c.A.container == s:
+                        b1,b2 = [loc.particle for loc in (c.A,c.B)]
+                        pot = self.get_bond_potential(4,18.5)
+                        self.add_bond(b1,b2, pot)
 
     # def get_bead(self, location):
     #     if type(location.container) is not list:

segmentmodel_example.py

@@ -14,7 +14,7 @@ PATH=/home/cmaffeo2/anaconda3/bin:$PATH; source activate cadnano
 
 if __name__ == "__main__":
 
-    seg1 = DoubleStrandedSegment("strand", num_nts = 46, local_twist=True)
+    seg1 = DoubleStrandedSegment("strand", num_nts = 46)
 
     seg2 = SingleStrandedSegment("strand", 
                                  start_position = seg1.end_position + np.array((5,0,5)),
@@ -41,5 +41,4 @@ if __name__ == "__main__":
                           dimensions=(5000,5000,5000),
                       )
 
-    model.useNonbondedScheme( nbDnaScheme )
     model.simulate( outputPrefix = 'strand-test', outputPeriod=1e4, numSteps=1e6, gpu=1 )