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."""