Added routines for inserting and removing DNA; made it so atomic simulation...

Added routines for inserting and removing DNA; made it so atomic simulation doesn't resize model dimensions

mrdna/model/arbdmodel.py

@@ -990,7 +990,7 @@ component "data" value 3
                 item = tuple(int(p.idx) for p in ex)
                 fh.write("EXCLUDE %d %d\n" % item)
 
-    def set_dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
+    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
         ## TODO: cache coordinates using numpy arrays for quick min/max
         raise(NotImplementedError)
 
@@ -1015,7 +1015,7 @@ tclForcesScript $prefix.forces.tcl
             format_data['tcl_forces'] = ""
 
         if update_dimensions:
-            self.set_dimensions_from_structure()
+            format_data['dimensions'] = self.dimensions_from_structure()
 
         for k,v in zip('XYZ', self.dimensions):
             format_data['origin'+k] = -v*0.5

mrdna/segmentmodel.py

@@ -409,6 +409,59 @@ class Segment(ConnectableElement, Group):
         tck, u = self.position_spline_params
         return np.mean(self.contour_to_position(u), axis=0)
 
+    def _get_location_positions(self):
+        return [self.contour_to_nt_pos(l.address) for l in self.locations]
+
+    def insert_dna(self, at_nt: int, num_nt: int):
+        assert(np.isclose(np.around(num_nt),num_nt))
+        if at_nt < 0:
+            raise ValueError("Attempted to insert DNA into {} at a negative location".format(self))
+        if at_nt > self.num_nt-1:
+            raise ValueError("Attempted to insert DNA into {} at beyond the end of the Segment".format(self))
+        if num_nt < 0:
+            raise ValueError("Attempted to insert DNA a negative amount of DNA into {}".format(self))
+
+        num_nt = np.around(num_nt)
+        nt_positions = self._get_location_positions()
+        new_nt_positions = [p if p <= at_nt else p+num_nt for p in nt_positions]
+
+        self.num_nt = self.num_nt+num_nt
+
+        for l,p in zip(self.locations, new_nt_positions):
+            l.address = self.nt_pos_to_contour(p)
+
+    def remove_dna(self, first_nt: int, last_nt: int):
+        """ Removes nucleotides between first_nt and last_nt """
+        assert(np.isclose(np.around(first_nt),first_nt))
+        assert(np.isclose(np.around(last_nt),last_nt))
+        tmp = min((first_nt,last_nt))
+        last_nt = max((first_nt,last_nt))
+        fist_nt = tmp
+
+        if first_nt < 0 or first_nt > self.num_nt-2:
+            raise ValueError("Attempted to remove DNA from {} starting at an invalid location {}".format(self, first_nt))
+        if last_nt < 1 or last_nt > self.num_nt-1:
+            raise ValueError("Attempted to remove DNA from {} ending at an invalid location {}".format(self, last_nt))
+        if first_nt == last_nt:
+            return
+
+        first_nt = np.around(first_nt)
+        last_nt = np.around(last_nt)
+
+        nt_positions = self._get_location_positions()
+
+        bad_locations = list(filter(lambda p: p >= first_nt and p <= last_nt, nt_positions))
+        if len(bad_locations) > 0:
+            raise Exception("Attempted to remove DNA containing locations {} from {} between {} and {}".format(bad_locations,self,first_nt,last_nt))
+
+        num_nt = last_nt-first_nt
+        new_nt_positions = [p if p <= last_nt else p-num_nt for p in nt_positions]
+
+        self.num_nt = self.num_nt-num_nt
+
+        for l,p in zip(self.locations, new_nt_positions):
+            l.address = self.nt_pos_to_contour(p)
+
     def __filter_contours(contours, positions, position_filter, contour_filter):
         u = contours
         r = positions
@@ -1527,11 +1580,17 @@ class SegmentModel(ArbdModel):
     def _clear_beads(self):
         ## TODO: deprecate
         for s in self.segments:
-            s.clear_all()
+            try:
+                s.clear_all()
+            except:
+                ...
         self.clear_all(keep_children=True)
 
-        if len(self.strands[0].children[0].children) > 0:
-            self.clear_atomic()
+        try:
+            if len(self.strands[0].children[0].children) > 0:
+                self.clear_atomic()
+        except:
+            ...
 
         ## Check that it worked
         assert( len([b for b in self]) == 0 )
@@ -2658,10 +2717,10 @@ class SegmentModel(ArbdModel):
                     elif dot1 < -0.5 and dot2 < -0.5:
                         ## TODO, reverse
                         ...
-                        print("Warning: {} and {} are on antiparallel helices (not yet implemented)... skipping".format(A1,B1))
+                        # print("Warning: {} and {} are on antiparallel helices (not yet implemented)... skipping".format(A1,B1))
                         continue
                     else:
-                        print("Warning: {} and {} are on helices that do not point in similar direction... skipping".format(A1,B1))
+                        # print("Warning: {} and {} are on helices that do not point in similar direction... skipping".format(A1,B1))
                         continue
 
                     ## Go through each nucleotide between the two
@@ -2732,7 +2791,7 @@ proc calcforces {} {
 }
 """)
 
-    def set_dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
+    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
         positions = []
         for s in self.segments:
             positions.append(s.contour_to_position(0))
@@ -2742,7 +2801,7 @@ proc calcforces {} {
         dx,dy,dz = [(np.max(positions[:,i])-np.min(positions[:,i])+30)*padding_factor for i in range(3)]
         if isotropic:
             dx = dy = dz = max((dx,dy,dz))
-        self.dimensions = [dx,dy,dz]
+        return [dx,dy,dz]
 
     def add_grid_potential(self, grid_file, scale=1, per_nucleotide=True):
         grid_file = Path(grid_file)