Add bonds to Topology based on CONECT records in PDB file

wrappers/python/simtk/openmm/app/internal/pdbstructure.py

@@ -137,6 +137,14 @@ class PdbStructure(object):
                 self._current_model._current_chain._add_ter_record()
             elif (pdb_line.find("CRYST1") == 0):
                 self._unit_cell_dimensions = (float(pdb_line[6:15]), float(pdb_line[15:24]), float(pdb_line[24:33]))*unit.angstroms
+            elif (pdb_line.find("CONECT") == 0):
+                atoms = [int(pdb_line[7:12])]
+                for pos in (12,17,22,27):
+                    try:
+                        atoms.append(int(pdb_line[pos:pos+5]))
+                    except:
+                        pass
+                self._current_model.connects.append(atoms)
         self._finalize()
 
     def write(self, output_stream=sys.stdout):
@@ -242,6 +250,7 @@ class Model(object):
         self.chains = []
         self._current_chain = None
         self.chains_by_id = {}
+        self.connects = []
 
     def _add_atom(self, atom):
         """

wrappers/python/simtk/openmm/app/pdbfile.py

@@ -47,6 +47,7 @@ class PDBFile(object):
 
         # Build the topology
 
+        atomByNumber = {}
         for chain in pdb.iter_chains():
             c = top.addChain()
             for residue in chain.iter_residues():
@@ -87,7 +88,8 @@ class PDBFile(object):
                             element = elem.get_by_symbol(atomName[0])
                         except KeyError:
                             pass
-                    top.addAtom(atomName, element, r)
+                    newAtom = top.addAtom(atomName, element, r)
+                    atomByNumber[atom.serial_number] = newAtom
                     pos = atom.get_position().value_in_unit(nanometers)
                     coords.append(Vec3(pos[0], pos[1], pos[2]))
         self.positions = coords*nanometers
@@ -96,6 +98,21 @@ class PDBFile(object):
         self.topology.createDisulfideBonds(self.positions)
         self._numpyPositions = None
         
+        # Add bonds based on CONECT records.
+        
+        connectBonds = []
+        for connect in pdb.get_model(0).connects:
+            i = connect[0]
+            for j in connect[1:]:
+                connectBonds.append((atomByNumber[i], atomByNumber[j]))
+        if len(connectBonds) > 0:
+            # Only add bonds that don't already exist.
+            existingBonds = set(top.bonds())
+            for bond in connectBonds:
+                if bond not in existingBonds and (bond[1], bond[0]) not in existingBonds:
+                    top.addBond(bond[0], bond[1])
+                    existingBonds.add(bond)
+
     def getTopology(self):
         """Get the Topology of the model."""
         return self.topology