diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py index f758df5f6ccc6b81afb6a51de0bf3bf54f2e2e4f..54a0b82d19ceb026abe1bb807cf903ce6156798f 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: