Improve CHARMM parsing when resnums have inscodes

Apparently CHARMM PSF files (particularly those printed by VMD's psfgen) can
have an insertion code tacked on to the end of the residue number, meaning that
casting to an integer will invariably fail.  We can't just ignore the insertion
code, though, since it will differentiate from the residues that came before and
after (which will likely have the same number -- same general idea as insertion
codes from PDB files). So instead we pull the insertion code off of the residue
number and extend the Residue object to accommodate that.

This is the port of the fix from ParmEd: ParmEd/ParmEd#98

I added a test for this case.

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

@@ -37,6 +37,7 @@ from __future__ import division
 from functools import wraps
 from math import pi, cos, sin, sqrt
 import os
+import re
 import simtk.openmm as mm
 from simtk.openmm.vec3 import Vec3
 import simtk.unit as u
@@ -99,6 +100,8 @@ class _ZeroDict(dict):
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+_resre = re.compile(r'(\d+)([a-zA-Z]*)')
+
 class CharmmPsfFile(object):
     """
     A chemical structure instantiated from CHARMM files.
@@ -196,7 +199,12 @@ class CharmmPsfFile(object):
             if atid != i + 1:
                 raise CharmmPSFError('Nonsequential atoms detected!')
             system = words[1]
-            resid = conv(words[2], int, 'residue number')
+            rematch = _resre.match(words[2])
+            if not rematch:
+                raise RuntimeError('Could not parse residue number %s' %
+                                   words[2])
+            resid, inscode = rematch.groups()
+            resid = conv(resid, int, 'residue number')
             resname = words[3]
             name = words[4]
             attype = words[5]
@@ -209,7 +217,7 @@ class CharmmPsfFile(object):
             mass = conv(words[7], float, 'atomic mass')
             props = words[8:]
             atom = residue_list.add_atom(system, resid, resname, name,
-                                         attype, charge, mass, props)
+                            attype, charge, mass, inscode, props)
             atom_list.append(atom)
         atom_list.assign_indexes()
         atom_list.changed = False

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -375,12 +375,13 @@ class AtomList(TrackedList):
 class Residue(object):
     """ Residue class """
 
-    def __init__(self, resname, idx):
+    def __init__(self, resname, idx, inscode=''):
         self.resname = resname
         self.idx = idx
         self.atoms = []
         self.resnum = None # Numbered based on SYSTEM name
         self.system = None
+        self.inscode = inscode
 
     def add_atom(self, atom):
         if self.system is None:
@@ -416,7 +417,7 @@ class Residue(object):
             return self.resname == thing.resname and self.idx == thing.idx
         if isinstance(thing, tuple) or isinstance(thing, list):
             # Must be resnum, resname
-            return thing == (self.resname, self.idx)
+            return thing == (self.resname, self.idx, self.inscode)
         return False # No other type can be equal.
 
     def __ne__(self, thing):
@@ -447,7 +448,7 @@ class ResidueList(list):
             resnum += 1
 
     def add_atom(self, system, resnum, resname, name,
-                 attype, charge, mass, props=None):
+                 attype, charge, mass, inscode, props=None):
         """
         Adds an atom to the list of residues. If the residue is not the same as
         the last residue that was created, a new Residue is created and added
@@ -461,17 +462,19 @@ class ResidueList(list):
             - attype (int or str) : Type of the atom
             - charge (float) : Partial atomic charge of the atom
             - mass (float) : Mass (amu) of the atom
+            - inscode (str) : Insertion code, if it is specified
 
         Returns:
             The Atom instance created and added to the list of residues
         """
         if self._last_residue is None:
-            res = self._last_residue = Residue(resname, resnum)
+            res = self._last_residue = Residue(resname, resnum, inscode)
             list.append(self, res)
-        elif (self._last_residue != (resname, resnum) or
+        elif (self._last_residue != (resname, resnum, inscode) or
               system != self._last_residue.system):
             if (self._last_residue.idx == resnum and
-                self._last_residue.system == system):
+                self._last_residue.system == system and
+                self._last_residue.inscode == inscode):
                 lresname = self._last_residue.resname
                 warnings.warn('Residue %d split into separate residues %s '
                               'and %s' % (resnum, lresname, resname),

wrappers/python/tests/TestCharmmFiles.py

@@ -115,6 +115,13 @@ class TestCharmmFiles(unittest.TestCase):
         ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         self.assertAlmostEqual(ene, 15490.0033559, delta=0.05)
 
+    def test_InsCode(self):
+        """ Test the parsing of PSF files that contain insertion codes in their residue numbers """
+        psf = CharmmPsfFile('systems/4TVP-dmj_wat-ion.psf')
+        self.assertEqual(len(list(psf.topology.atoms())), 66264)
+        self.assertEqual(len(list(psf.topology.residues())), 20169)
+        self.assertEqual(len(list(psf.bonds())), 46634)
+
     def test_ImplicitSolventForces(self):
         """Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
         solventType = [HCT, OBC1, OBC2, GBn, GBn2]