From 385f358ce2b4e95d01f5c30b72eca5628a6c6643 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Fri, 21 Feb 2020 16:57:00 -0600
Subject: [PATCH] Option and commands to convert crossovers into intrahelical
---
bin/mrdna | 12 +++++++++
mrdna/segmentmodel.py | 60 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 72 insertions(+)
diff --git a/bin/mrdna b/bin/mrdna
index 0c7d04e..e43df23 100755
--- a/bin/mrdna
+++ b/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
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 8867774..6a8c6a6 100644
--- a/mrdna/segmentmodel.py
+++ b/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,
--
GitLab