Write out ENM bonds in blocks according to residue pair type such as stack, (base)pair, etc

mrdna/segmentmodel.py

@@ -3494,7 +3494,8 @@ class SegmentModel(ArbdModel):
 
         with open("%s.exb" % output_name,'w') as fh:
             # natoms=0
-            bonds = []
+            pairtypes = ('pair','stack','cross','paircross')
+            bonds = {k:[] for k in pairtypes}
 
             for seg in self.segments:
                 ## Continue unless dsDNA
@@ -3523,13 +3524,13 @@ class SegmentModel(ArbdModel):
                             if nt2 is not None and strand_piece.is_fwd:
                                 other.append((nt2,'cross'))
 
-                        for nt2,key in other:
+                        for nt2,pairtype in other:
                             """
                             if np.linalg.norm(nt2.position-nt1.position) > 7:
                                 import pdb
                                 pdb.set_trace()
                             """
-                            key = ','.join((key,nt1.sequence[0],nt2.sequence[0]))
+                            key = ','.join((pairtype,nt1.sequence[0],nt2.sequence[0]))
                             for n1, n2, d in enmTemplate[key]:
                                 d = float(d)
                                 k = 0.1
@@ -3542,15 +3543,17 @@ class SegmentModel(ArbdModel):
 
                                 i = nt1._get_atomic_index(name=n1)
                                 j = nt2._get_atomic_index(name=n2)
-                                bonds.append((i,j,k,d))
+                                bonds[pairtype].append((i,j,k,d))
 
             # print("NO STACKS found for:", noStackPrime)
             # print("NO BASEPAIRS found for:", noBasepair)
 
-            if len(bonds) != len(set(bonds)):
-                devlogger.warning("Duplicate ENM bonds")
-            for b in set(bonds):
-                fh.write("bond %d %d %f %.2f\n" % b)
+            for k,bondlist in bonds.items():
+                if len(bondlist) != len(set(bondlist)):
+                    devlogger.warning("Duplicate ENM bonds for pair type {}".format(k))
+                fh.write("# {}\n".format(k.upper()))
+                for b in set(bondlist):
+                    fh.write("bond %d %d %f %.2f\n" % b)
 
         ## Loop dsDNA regions
         push_bonds = []
@@ -3637,6 +3640,7 @@ class SegmentModel(ArbdModel):
         if interhelical_bonds:
             if not self.useTclForces:
                 with open("%s.exb" % output_name, 'a') as fh:
+                    fh.write("# PUSHBONDS\n")
                     for i,j in push_bonds:
                         fh.write("bond %d %d %f %.2f\n" % (i,j,1.0,31))
             else: