Remove tri-sulfite bond formation. (#3203)

* Remove tri-sulfite bond formation.

Pdbfixer attempts to add di-sulfide bonds to any pair of CYS "S" atoms that are less than 3 A. This procedure fails when multiple such bonds are assigned to a given atom. 

Here we rectify this issue by a greedy algorithm: Atom I is chosen to be bonded to `argmin(distance([:I-1, I])` over all candidate di-sulfide bonds. This procedure will not generally result in a set of minimal distance di-sulfide bonds, but will fix the "tri-sulfide bond" problem. Fixing the more general case would require a solution to global constraint satisfaction problem.

**Legal**
This patch is provided under the CC0 license (https://creativecommons.org/share-your-work/public-domain/cc0/).
This is not an officially-supported Google product. It is provided ‘as-is’ without any warranty of any kind, whether expressed or implied.

* Update topology.py

Remove a redundant comment.

wrappers/python/openmm/app/topology.py

@@ -353,19 +353,33 @@ class Topology(object):
         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.
+        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."""