openmm.patch

Index: simtk/openmm/app/topology.py
===================================================================
--- simtk.orig/openmm/app/topology.py
+++ simtk/openmm/app/topology.py
@@ -356,19 +356,35 @@
         def isCyx(res):
             names = [atom.name for atom in res._atoms]
             return 'SG' in names and 'HG' not in names
+        # This function is used to prevent multiple di-sulfide bonds from being
+        # assigned to a given atom. This is a DeepMind modification.
+        def isDisulfideBonded(atom):
+            for b in self._bonds:
+                if (atom in b and b[0].name == 'SG' and
+                    b[1].name == 'SG'):
+                    return True
+
+            return False
 
         cyx = [res for res in self.residues() if res.name == 'CYS' and isCyx(res)]
         atomNames = [[atom.name for atom in res._atoms] for res in cyx]
         for i in range(len(cyx)):
             sg1 = cyx[i]._atoms[atomNames[i].index('SG')]
             pos1 = positions[sg1.index]
+            candidate_distance, candidate_atom = 0.3*nanometers, None
             for j in range(i):
                 sg2 = cyx[j]._atoms[atomNames[j].index('SG')]
                 pos2 = positions[sg2.index]
                 delta = [x-y for (x,y) in zip(pos1, pos2)]
                 distance = sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2])
-                if distance < 0.3*nanometers:
-                    self.addBond(sg1, sg2)
+                if distance < candidate_distance and not isDisulfideBonded(sg2):
+                    candidate_distance = distance
+                    candidate_atom = sg2
+            # Assign bond to closest pair.
+            if candidate_atom:
+                self.addBond(sg1, candidate_atom)
+
+
 
 class Chain(object):
     """A Chain object represents a chain within a Topology."""