From 4526e407beb33371b6b8343ef9dd8b17f4330e8a Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Wed, 22 Apr 2020 12:51:32 -0500
Subject: [PATCH] Added routines for writing atomic bp restraints and for
changing debye length
---
mrdna/segmentmodel.py | 40 +++++++++++++++++++++++++++++++++++++++-
1 file changed, 39 insertions(+), 1 deletion(-)
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index f758df5..54a0b82 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -1675,10 +1675,21 @@ class SegmentModel(ArbdModel):
if debye_length is not None:
nbDnaScheme.debye_length = debye_length
+ self.debye_length = debye_length
self.useNonbondedScheme( nbDnaScheme )
self.useTclForces = False
+ def set_debye_length(self, value):
+ if value <= 0:
+ raise ValueError("The Debye length must be positive")
+ for s,tA,tB in self.nbSchemes:
+ try:
+ s.debye_length = value
+ except:
+ ...
+ self.debye_length = value
+
def get_connections(self,type_=None,exclude=()):
""" Find all connections in model, without double-counting """
added=set()
@@ -2125,7 +2136,7 @@ class SegmentModel(ArbdModel):
""" 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]
+ return [s for s in self.segments if re_pattern.match(s.name) is not None]
def get_crossovers_at_ends(self):
"""
@@ -3274,6 +3285,33 @@ class SegmentModel(ArbdModel):
nt1.basepair = nt2
nt2.basepair = nt1
+ def write_atomic_bp_restraints(self, output_name, spring_constant=1.0):
+ ## TODO: ensure atomic model was generated already
+ ## TODO: allow ENM to be created without first building atomic model
+
+ with open("%s.exb" % output_name,'w') as fh:
+ for seg in self.segments:
+ ## Continue unless dsDNA
+ if not isinstance(seg,DoubleStrandedSegment): continue
+ for strand_piece in seg.strand_pieces['fwd']:
+ assert( strand_piece.is_fwd )
+ for nt1 in strand_piece.children:
+ nt2 = nt1.basepair
+ if nt1.resname == 'ADE':
+ names = (('N1','N3'),('N6','O4'))
+ elif nt1.resname == 'GUA':
+ names = (('N2','O2'),('N1','N3'),('O6','N4'))
+ elif nt1.resname == 'CYT':
+ names = (('O2','N2'),('N3','N1'),('N4','O6'))
+ elif nt1.resname == 'THY':
+ names = (('N3','N1'),('O4','N6'))
+ else:
+ raise Exception("Unrecognized nucleotide!")
+ for n1, n2 in names:
+ i = nt1._get_atomic_index(name=n1)
+ j = nt2._get_atomic_index(name=n2)
+ fh.write("bond %d %d %f %.2f\n" % (i,j,spring_constant,2.8))
+
def write_atomic_ENM(self, output_name, lattice_type=None):
## TODO: ensure atomic model was generated already
if lattice_type is None:
--
GitLab