Various changes to make it easy to run circular structures

mrdna/segmentmodel.py

@@ -199,8 +199,7 @@ class ConnectableElement():
         self.connections.append(connection)
         if other is not self:
             other.connections.append(connection)
-        else:
-            raise NotImplementedError("Segments cannot yet be connected to themselves; if you are attempting to make a circular object, try breaking the object into multiple segments")
+
         l = A.container.locations
         if A not in l: l.append(A)
         l = B.container.locations
@@ -977,7 +976,6 @@ class Segment(ConnectableElement, Group):
 
 
     def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
-
         """ Generate beads (positions, types, etc) and bonds, angles, dihedrals, exclusions """
         ## TODO: decide whether to remove bead_model argument
         ##       (currently unused)
@@ -1003,7 +1001,7 @@ class Segment(ConnectableElement, Group):
 
         ## Add ends if they don't exist yet
         ## TODOTODO: test 1 nt segments?
-        if len(existing_beads) == 0 or existing_beads[0][0].get_nt_position(self,0) >= 0.5:
+        if len(existing_beads) == 0 or (existing_beads[0][0].get_nt_position(self,0) >= 0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__):
             # if len(existing_beads) > 0:            
             #     assert(existing_beads[0].get_nt_position(self) >= 0.5)
             c = self.nt_pos_to_contour(0)
@@ -1011,7 +1009,7 @@ class Segment(ConnectableElement, Group):
             b = self._generate_one_bead(c, 0)
             existing_beads = [(b,0)] + existing_beads
 
-        if existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 or len(existing_beads)==1:
+        if (existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__) or len(existing_beads)==1:
             c = self.nt_pos_to_contour(self.num_nt-1)
             if self.num_nt == 1: c += 0.4
             b = self._generate_one_bead(c, 0)
@@ -2141,6 +2139,7 @@ class SegmentModel(ArbdModel):
         """ Generate beads at intrahelical junctions """
         if self.DEBUG: print( "Adding intrahelical beads at junctions" )
 
+
         ## Loop through all connections, generating beads at appropriate locations
         for c,A,B in self.get_connections("intrahelical"):
             s1,s2 = [l.container for l in (A,B)]
@@ -2155,21 +2154,22 @@ class SegmentModel(ArbdModel):
                 assert( np.isclose(a,0) or np.isclose(a,1) )
             
             ## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
-            """ Search to see whether bead at location is already found """
-            b = None
-            if isinstance(s1,DoubleStrandedSegment):
-                b = s1.get_nearest_bead(a1) 
+
+            def _get_nearest(seg,address):
+                b = seg.get_nearest_bead(address)
                 if b is not None:
-                    assert( b.parent is s1 )
+                    assert( b.parent is seg )
                     """ if above assertion is true, no problem here """
-                    if np.abs(b.get_nt_position(s1) - s1.contour_to_nt_pos(a1)) > 0.5:
-                        b = None
-            if b is None and isinstance(s2,DoubleStrandedSegment):
-                b = s2.get_nearest_bead(a2)
-                if b is not None:
-                    if np.abs(b.get_nt_position(s2) - s2.contour_to_nt_pos(a2)) > 0.5:
+                    if np.abs(b.get_nt_position(seg) - seg.contour_to_nt_pos(address)) > 0.5:
                         b = None
+                return b
 
+            """ Search to see whether bead at location is already found """
+            b = None
+            if isinstance(s1,DoubleStrandedSegment):
+                b = _get_nearest(s1,a1)
+            if b is None and isinstance(s2,DoubleStrandedSegment):
+                b = _get_nearest(s2,a2)
 
             if b is not None and b.parent not in (s1,s2):
                 b = None
@@ -2183,6 +2183,8 @@ class SegmentModel(ArbdModel):
             A.particle = B.particle = b
             b.is_intrahelical = True
             b.locations.extend([A,B])
+            if s1 is s2:
+                b.is_circular_bead = True
 
         # pdb.set_trace()
 
@@ -2242,12 +2244,27 @@ class SegmentModel(ArbdModel):
         # pdb.set_trace()
 
         """ Add intrahelical neighbors at connections """
-
+        special_seps = dict()
         for c,A,B in self.get_connections("intrahelical"):
             b1,b2 = [l.particle for l in (A,B)]
             if b1 is b2:
-                ## already handled by Segment._generate_beads
-                pass
+                ## already handled by Segment._generate_beads, unless A,B are on same segment
+                if A.container == B.container:
+                    sorted_sites = sorted([(abs(L.address-b1.contour_position),L)
+                                       for L in (A,B)])
+                    close_site, far_site = [s[1] for s in sorted_sites]
+
+                    b3 = far_site.container.get_nearest_bead(far_site.address)
+                    b1.intrahelical_neighbors.append(b3)
+                    b3.intrahelical_neighbors.append(b1)
+
+                    if far_site.address > 0.5:
+                        sep = b1.contour_position + (1-b3.contour_position)
+                    else:
+                        sep = b3.contour_position + (1-b1.contour_position)
+                    sep = A.container.num_nt * sep
+                    special_seps[(b1,b3)] = special_seps[(b3,b1)] = sep
+
             else:
                 for b in (b1,b2): assert( b is not None )
                 b1.make_intrahelical_neighbor(b2)
@@ -2358,11 +2375,13 @@ class SegmentModel(ArbdModel):
 
             parent = self._getParent(b1,b2)
 
+            if (b1,b2) in special_seps:
+                sep = special_seps[(b1,b2)]
+            else:
                 seg = b2.parent
                 c0 = b2.contour_position
                 sep = np.abs(b1.get_nt_position(seg,c0)-b2.get_nt_position(seg))
 
-
             is_dsdna = b1.type_.name[0] == "D" and b2.type_.name[0] == "D"
             if is_dsdna:
                 d = 3.4*sep
@@ -2426,7 +2445,7 @@ class SegmentModel(ArbdModel):
             parent = self._getParent(b1,b2,b3)
             seg = b2.parent
             c0 = b2.contour_position
-            sep = np.abs(b1.get_nt_position(seg,c0)-b2.get_nt_position(seg)) + np.abs(b3.get_nt_position(seg,c0)-b2.get_nt_position(seg))
+            sep = dists[b2][b1] + dists[b2][b3]
 
 
             if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":

mrdna/simulate.py

@@ -10,6 +10,8 @@ import numpy as np
 
 def _remove_bonds_beyond(model, cutoff):
     bonds = model.get_bonds()
+    for seg in model.segments:
+        bonds.extend(seg.get_bonds())
     new_bonds = []
     for bond in bonds:
         n1,n2,b,ex = bond