From dc3feb95bdcbfbcb169aa43153b881ba3f21801f Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Thu, 13 Sep 2018 11:55:32 -0500
Subject: [PATCH] Improved support for short ssDNA segments

Previous commit fixed an issue for from_pdb reader, which used "sscrossover" connections, but broke models from the vHelix reader which used "intrahelical" connections. Now we convert "sscrossover" connections at the ends of segments into "intrahelical" connections.
---
 mrdna/segmentmodel.py | 28 ++++++++++++++++++++++------
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index f642d6e..d7eb84d 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -34,7 +34,7 @@ TODO:
  - rework Location class 
  - remove recursive calls
  - document
- - add unit test of helices connected to themselves
+ - develop unit test suite
 """
 class CircularDnaError(Exception):
     pass
@@ -747,8 +747,9 @@ class Segment(ConnectableElement, Group):
         existing_beads0 = {l.particle for l in self.locations if l.particle is not None}
         existing_beads = sorted( list(existing_beads0), key=lambda b: b.get_contour_position(self) )
 
-        if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
-            return
+        # if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
+        #     pdb.set_trace()
+        #     return
 
         for b in existing_beads:
             assert(b.parent is not None)
@@ -995,9 +996,13 @@ class SingleStrandedSegment(Segment):
     def connect_end3(self, end5, force_connection=False):
         self._connect_end( end5,  _5_to_3 = True, force_connection = force_connection )
 
-    def connect_5end(self, end3, force_connection=False): # TODO: change name or possibly deprecate
+    def connect_start5(self, end3, force_connection=False):
         self._connect_end( end3,  _5_to_3 = False, force_connection = force_connection )
 
+    def connect_5end(self, end3, force_connection=False): # TODO: change name or possibly deprecate
+        print("WARNING: 'connect_5end' will be deprecated")
+        return self.connect_start5( end3, force_connection=False)
+
     def _connect_end(self, other, _5_to_3, force_connection):
         assert( isinstance(other, Location) )
         if _5_to_3 == True:
@@ -1012,7 +1017,7 @@ class SingleStrandedSegment(Segment):
             seg1 = other.container
             seg2 = self
             end1 = other
-            end2 = self.end5
+            end2 = self.start
             # assert( other.type_ == "end3" )
             if (other.type_ is not "end3"):
                 print("WARNING: code does not prevent connecting 3prime to 3prime, etc")
@@ -1035,13 +1040,21 @@ class SingleStrandedSegment(Segment):
         ##   TODO
        
         ## TODO: fix direction
+        if nt in (0,1,self.num_nt) and other_nt in (0,1,other.num_nt):
+            if nt_on_5prime == True:
+                other_end = other.start5 if strands_fwd[1] else other.end5
+                self.connect_end3( other_end )
+            else:
+                other_end = other.end3 if strands_fwd[1] else other.start3
+                self.connect_start5( other_end )
+            return
 
         # c1 = self.nt_pos_to_contour(nt)
         # # TODOTODO
         # ## Ensure connections occur at ends, otherwise the structure doesn't make sense
         # # assert(np.isclose(c1,0) or np.isclose(c1,1))
         # assert(np.isclose(nt,0) or np.isclose(nt,self.num_nt-1))
-        if nt == 0:
+        if nt == 0 and (nt_on_5prime or self.num_nt > 1):
             c1 = 0
         elif nt == self.num_nt-1:
             c1 = 1
@@ -1825,6 +1838,7 @@ class SegmentModel(ArbdModel):
             # assert( not np.isclose( np.linalg.norm(b1.collapsedPosition() - b2.collapsedPosition()), 0 ) )
             if np.linalg.norm(b1.collapsedPosition() - b2.collapsedPosition()) < 1:
                 print("WARNING: some beads are very close")
+
             parent = self._getParent(b1,b2)
 
             ## TODO: could be sligtly smarter about sep
@@ -1851,6 +1865,8 @@ class SegmentModel(ArbdModel):
             dists[b1][b2] = sep
             dists[b2][b1] = sep
 
+            if b1 is b2: continue
+
             # dists[[b1,b2]] = dists[[b2,b1]] = sep
             bond = self.get_bond_potential(k,d)
             parent.add_bond( b1, b2, bond, exclude=True )
-- 
GitLab