Added routines for writing atomic bp restraints and for changing debye length

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: