Ensure no duplicate bonds in ENM, logging violations; thanks to Elija Feigl...

Ensure no duplicate bonds in ENM, logging violations; thanks to Elija Feigl for bringing duplicate bonds to our attention

mrdna/segmentmodel.py

@@ -3491,8 +3491,10 @@ class SegmentModel(ArbdModel):
         ## TODO: allow ENM to be created without first building atomic model
         noStackPrime = 0
         noBasepair = 0
+
         with open("%s.exb" % output_name,'w') as fh:
             # natoms=0
+            bonds = []
 
             for seg in self.segments:
                 ## Continue unless dsDNA
@@ -3540,11 +3542,16 @@ class SegmentModel(ArbdModel):
 
                                 i = nt1._get_atomic_index(name=n1)
                                 j = nt2._get_atomic_index(name=n2)
-                                fh.write("bond %d %d %f %.2f\n" % (i,j,k,d))
+                                bonds.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)
+
         ## Loop dsDNA regions
         push_bonds = []
         processed_segs = set()
@@ -3622,6 +3629,10 @@ class SegmentModel(ArbdModel):
                                 # print("WARNING: could not find 'P' atom in {}:{} or {}:{}".format( segI, ntI_idx, segJ, ntJ_idx ))
                                 ...
 
+        if len(push_bonds) != len(set(push_bonds)):
+            devlogger.warning("Duplicate enrgmd push bonds")
+        push_bonds = list(set( push_bonds ))
+
         # print("PUSH BONDS:", len(push_bonds))
         if interhelical_bonds:
             if not self.useTclForces: