Option and commands to convert crossovers into intrahelical

bin/mrdna

@@ -56,6 +56,9 @@ parser.add_argument('--oxdna-output-period', type=float, default=1e4,
 parser.add_argument('--coarse-bond-cutoff', type=float, default = 0,
                     help='Ignore bonds beyond this cutoff during first step of simulation; a value of 0 implies bonds are not ignored')
 
+parser.add_argument('--crossover-to-intrahelical-cutoff', type=float, default=-1,
+                    help='Set the distance (in Angstroms) beyond which crossovers between helix ends are converted to intrahelical connections; a negative value means no crossovers will be converted')
+
 parser.add_argument('--backbone-scale', type=float, default=1.0,
                     help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
 
@@ -119,6 +122,15 @@ def main():
     else:
         model.dimensions = args.dimensions
 
+    if args.crossover_to_intrahelical_cutoff >= 0:
+        model.convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(
+            args.crossover_to_intrahelical_cutoff)
+
+    if args.sequence_file is not None:
+        from mrdna.model.dna_sequence import read_sequence_file
+        seq = read_sequence_file(args.sequence_file)
+        model.set_sequence( seq )
+
     if args.output_prefix is not None:
         prefix = args.output_prefix
     

mrdna/segmentmodel.py

@@ -8,6 +8,8 @@ from .arbdmodel.nonbonded import *
 from copy import copy, deepcopy
 from .model.nbPot import nbDnaScheme
 
+import re
+
 from scipy.special import erf
 import scipy.optimize as opt
 from scipy import interpolate
@@ -2027,6 +2029,64 @@ class SegmentModel(ArbdModel):
                 tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 )
                 s.quaternion_spline_params = (tck,u)
 
+    def get_segments_matching(self, pattern):
+        """ Returns a list of all segments that match the regular expression 'pattern' """
+
+        re_pattern = re.compile(pattern)
+        return [s for s in self.segments if re_pattern.search(s.name) is not None]
+
+    def get_crossovers_at_ends(self):
+        """
+        Return a list of all "regular" crossovers where both sides of
+        the crossover exist at the start or end of a segment
+        """
+
+        ret_list = []
+        for s in self.segments:
+            for c in filter(lambda c: c.type_ == "crossover", s.connections):
+                seg1 = c.A.container
+                seg2 = c.B.container
+                nt1 = c.A.get_nt_pos()
+                nt2 = c.B.get_nt_pos()
+                if (nt1 < 1 and nt2 > seg2.num_nt-2) or (nt2 < 1 and nt1 > seg1.num_nt-2):
+                    if c not in ret_list:
+                        ret_list.append(c)
+        return ret_list
+
+    def convert_crossover_to_intrahelical(self, crossover):
+        c = crossover
+        seg1 = c.A.container
+        seg2 = c.B.container
+        nt1 = c.A.get_nt_pos()
+        nt2 = c.B.get_nt_pos()
+
+        assert( c.A.on_fwd_strand == c.B.on_fwd_strand )
+        on_fwd_strand = c.A.on_fwd_strand
+
+        c.delete()
+
+        if nt1 < 1 and nt2 > seg2.num_nt-2:
+            if on_fwd_strand:
+                seg2.connect_end3(seg1.start5)
+            else:
+                seg1.connect_start3(seg2.end5)
+
+        elif nt2 < 1 and nt1 > seg1.num_nt-2:
+            if on_fwd_strand:
+                seg1.connect_end3(seg2.start5)
+            else:
+                seg2.connect_start3(seg1.end5)
+
+    def convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(self,cutoff):
+        if cutoff < 0:
+            raise ValueError("Cutoff should be positive!")
+
+        for c in self.get_crossovers_at_ends():
+            r1,r2 = [l.container.contour_to_position(l.address)
+                     for l in (c.A,c.B)]
+            if np.linalg.norm(r2-r1) > cutoff:
+                self.convert_crossover_to_intrahelical(c)
+
     def _generate_bead_model(self,
                              max_basepairs_per_bead = 7,
                              max_nucleotides_per_bead = 4,