Merge pull request #2068 from peastman/insertion

Preserve PDB insertion codes

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

@@ -109,7 +109,7 @@ class Modeller(object):
         for chain in self.topology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 for atom in residue.atoms():
                     newAtom = newTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                     newAtoms[atom] = newAtom
@@ -123,7 +123,7 @@ class Modeller(object):
         for chain in addTopology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 for atom in residue.atoms():
                     newAtom = newTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                     newAtoms[atom] = newAtom
@@ -168,7 +168,7 @@ class Modeller(object):
                                     newChain = newTopology.addChain(chain.id)
                                     needNewChain = False;
                                 if needNewResidue:
-                                    newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                                    newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                                     needNewResidue = False;
                                 newAtom = newTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                                 newAtoms[atom] = newAtom
@@ -210,7 +210,7 @@ class Modeller(object):
         for chain in self.topology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 if residue.name == "HOH":
                     # Copy the oxygen and hydrogens
                     oatom = [atom for atom in residue.atoms() if atom.element == elem.oxygen]
@@ -395,7 +395,7 @@ class Modeller(object):
         for chain in self.topology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 for atom in residue.atoms():
                     newAtom = newTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                     newAtoms[atom] = newAtom
@@ -691,7 +691,7 @@ class Modeller(object):
         for chain in self.topology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 isNTerminal = (residue == chain._residues[0])
                 isCTerminal = (residue == chain._residues[-1])
                 if residue.name in Modeller._residueHydrogens:
@@ -977,7 +977,7 @@ class Modeller(object):
         for chain in self.topology.chains():
             newChain = newTopology.addChain(chain.id)
             for residue in chain.residues():
-                newResidue = newTopology.addResidue(residue.name, newChain, residue.id)
+                newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
 
                 # Look for a matching template.
 
@@ -999,7 +999,7 @@ class Modeller(object):
                     # extra points.
 
                     template = None
-                    residueNoEP = Residue(residue.name, residue.index, residue.chain, residue.id)
+                    residueNoEP = Residue(residue.name, residue.index, residue.chain, residue.id, residue.insertionCode)
                     residueNoEP._atoms = [atom for atom in residue.atoms() if atom.element is not None]
                     if signature in forcefield._templateSignatures:
                         for t in forcefield._templateSignatures[signature]:
@@ -1310,7 +1310,7 @@ class Modeller(object):
                 # Remove the same number of residues from each leaf.
                 skipFromLeaf[lipidLeaf[residue]] -= 1
             else:
-                newResidue = membraneTopology.addResidue(residue.name, lipidChain, residue.id)
+                newResidue = membraneTopology.addResidue(residue.name, lipidChain, residue.id, residue.insertionCode)
                 for atom in residue.atoms():
                     newAtom = membraneTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                     newAtoms[atom] = newAtom
@@ -1322,7 +1322,7 @@ class Modeller(object):
         
         solventChain = membraneTopology.addChain()
         for (residue, pos) in addedWater:
-            newResidue = membraneTopology.addResidue(residue.name, solventChain, residue.id)
+            newResidue = membraneTopology.addResidue(residue.name, solventChain, residue.id, residue.insertionCode)
             for atom in residue.atoms():
                 newAtom = membraneTopology.addAtom(atom.name, atom.element, newResidue, atom.id)
                 newAtoms[atom] = newAtom

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

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2016 Stanford University and the Authors.
+Portions copyright (c) 2012-2018 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -106,7 +106,7 @@ class PDBFile(object):
                 resName = residue.get_name()
                 if resName in PDBFile._residueNameReplacements:
                     resName = PDBFile._residueNameReplacements[resName]
-                r = top.addResidue(resName, c, str(residue.number))
+                r = top.addResidue(resName, c, str(residue.number), residue.insertion_code)
                 if resName in PDBFile._atomNameReplacements:
                     atomReplacements = PDBFile._atomNameReplacements[resName]
                 else:
@@ -352,8 +352,10 @@ class PDBFile(object):
                     resName = res.name
                 if keepIds and len(res.id) < 5:
                     resId = res.id
+                    resIC = res.insertionCode
                 else:
                     resId = "%4d" % ((resIndex+1)%10000)
+                    resIC = " "
                 if res.name in nonHeterogens:
                     recordName = "ATOM  "
                 else:
@@ -370,10 +372,11 @@ class PDBFile(object):
                     else:
                         atomName = atom.name
                     coords = positions[posIndex]
-                    line = "%s%5d %-4s %3s %s%4s    %s%s%s  1.00  0.00          %2s  " % (
-                        recordName, atomIndex%100000, atomName, resName, chainName, resId, _format_83(coords[0]),
+                    line = "%s%5d %-4s %3s %s%4s%1s   %s%s%s  1.00  0.00          %2s  " % (
+                        recordName, atomIndex%100000, atomName, resName, chainName, resId, resIC, _format_83(coords[0]),
                         _format_83(coords[1]), _format_83(coords[2]), symbol)
-                    assert len(line) == 80, 'Fixed width overflow detected'
+                    if len(line) != 80:
+                        raise ValueError('Fixed width overflow detected')
                     print(line, file=file)
                     posIndex += 1
                     atomIndex += 1

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

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2015 Stanford University and the Authors.
+Portions copyright (c) 2015-2018 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors: Jason Swails
 
@@ -84,6 +84,7 @@ class PDBxFile(object):
         atomIdCol = atomData.getAttributeIndex('id')
         resNameCol = atomData.getAttributeIndex('auth_comp_id')
         resNumCol = atomData.getAttributeIndex('auth_seq_id')
+        resInsertionCol = atomData.getAttributeIndex('pdbx_PDB_ins_code')
         chainIdCol = atomData.getAttributeIndex('auth_asym_id')
         elementCol = atomData.getAttributeIndex('type_symbol')
         altIdCol = atomData.getAttributeIndex('label_alt_id')
@@ -117,7 +118,9 @@ class PDBxFile(object):
                     lastResId = None
                 if lastResId != row[resNumCol] or lastChainId != row[chainIdCol] or (lastResId == '.' and row[atomNameCol] in atomsInResidue):
                     # The start of a new residue.
-                    res = top.addResidue(row[resNameCol], chain, None if resNumCol == -1 else row[resNumCol])
+                    resId = (None if resNumCol == -1 else row[resNumCol])
+                    resIC = ('' if resInsertionCol == -1 else row[resInsertionCol])
+                    res = top.addResidue(row[resNameCol], chain, resId, resIC)
                     lastResId = row[resNumCol]
                     atomsInResidue.clear()
                 element = None
@@ -382,16 +385,18 @@ class PDBxFile(object):
             for (resIndex, res) in enumerate(residues):
                 if keepIds:
                     resId = res.id
+                    resIC = (res.insertionCode if len(res.insertionCode) > 0 else '.')
                 else:
                     resId = resIndex + 1
+                    resIC = '.'
                 for atom in res.atoms():
                     coords = positions[posIndex]
                     if atom.element is not None:
                         symbol = atom.element.symbol
                     else:
                         symbol = '?'
-                    line = "ATOM  %5d %-3s %-4s . %-4s %s ? %5s . %10.4f %10.4f %10.4f  0.0  0.0  ?  ?  ?  ?  ?  .  %5s %4s %s %4s %5d"
-                    print(line % (atomIndex, symbol, atom.name, res.name, chainName, resId, coords[0], coords[1], coords[2],
+                    line = "ATOM  %5d %-3s %-4s . %-4s %s ? %5s %s %10.4f %10.4f %10.4f  0.0  0.0  ?  ?  ?  ?  ?  .  %5s %4s %s %4s %5d"
+                    print(line % (atomIndex, symbol, atom.name, res.name, chainName, resId, resIC, coords[0], coords[1], coords[2],
                                   resId, res.name, chainName, atom.name, modelIndex), file=file)
                     posIndex += 1
                     atomIndex += 1

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

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2016 Stanford University and the Authors.
+Portions copyright (c) 2012-2018 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -137,7 +137,7 @@ class Topology(object):
         self._chains.append(chain)
         return chain
 
-    def addResidue(self, name, chain, id=None):
+    def addResidue(self, name, chain, id=None, insertionCode=''):
         """Create a new Residue and add it to the Topology.
 
         Parameters
@@ -149,6 +149,8 @@ class Topology(object):
         id : string=None
             An optional identifier for the residue.  If this is omitted, an id
             is generated based on the residue index.
+        insertionCode: string=''
+            An optional insertion code for the residue.
 
         Returns
         -------
@@ -159,7 +161,7 @@ class Topology(object):
             raise ValueError('All residues within a chain must be contiguous')
         if id is None:
             id = str(self._numResidues+1)
-        residue = Residue(name, self._numResidues, chain, id)
+        residue = Residue(name, self._numResidues, chain, id, insertionCode)
         self._numResidues += 1
         chain._residues.append(residue)
         return residue
@@ -395,7 +397,7 @@ class Chain(object):
 
 class Residue(object):
     """A Residue object represents a residue within a Topology."""
-    def __init__(self, name, index, chain, id):
+    def __init__(self, name, index, chain, id, insertionCode):
         """Construct a new Residue.  You should call addResidue() on the Topology instead of calling this directly."""
         ## The name of the Residue
         self.name = name
@@ -405,6 +407,8 @@ class Residue(object):
         self.chain = chain
         ## A user defined identifier for this Residue
         self.id = id
+        ## A user defined insertion code for this Residue
+        self.insertionCode = insertionCode
         self._atoms = []
 
     def atoms(self):