Improvements to Gromacs file processing

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

@@ -122,7 +122,7 @@ class GromacsGroFile(object):
                 na = int(line.strip())
             elif _is_gro_coord(line):
                 if frame == 0: # Create the list of residues, atom names etc. only if it's the first frame.
-                    (thisresnum, thisresname, thisatomname, thisatomnum) = [line[i*5:i*5+5].strip() for i in range(4)]
+                    (thisresnum, thisresname, thisatomname) = [line[i*5:i*5+5].strip() for i in range(3)]
                     resname.append(thisresname)
                     resid.append(int(thisresnum))
                     atomname.append(thisatomname)
@@ -133,7 +133,10 @@ class GromacsGroFile(object):
                             elements.append(elem.get_by_symbol(thiselem))
                         except KeyError:
                             elements.append(None)
-                pos = [float(line[20+i*8:28+i*8]) for i in range(3)]
+                firstDecimalPos = line.index('.', 20)
+                secondDecimalPos = line.index('.', firstDecimalPos+1)
+                digits = secondDecimalPos-firstDecimalPos
+                pos = [float(line[20+i*digits:20+(i+1)*digits]) for i in range(3)]
                 xyz.append(Vec3(pos[0], pos[1], pos[2]))
             elif _is_gro_box(line) and ln == na + 2:
                 sline = line.split()

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

@@ -270,8 +270,19 @@ class GromacsTopFile(object):
     def _processAtomType(self, line):
         """Process a line in the [ atomtypes ] category."""
         fields = line.split()
-        if len(fields) < 7:
+        if len(fields) < 6:
             raise ValueError('Too few fields in [ atomtypes ] line: '+line);
+        if len(fields[3]) == 1:
+            # Bonded type and atomic number are both missing.
+            fields.insert(1, None)
+            fields.insert(1, None)
+        elif len(fields[4]) == 1 and len(fields[5]) > 1:
+            if fields[1][0].isalpha():
+                # Atomic number is missing.
+                fields.insert(2, None)
+            else:
+                # Bonded type is missing.
+                fields.insert(1, None)
         self._atomTypes[fields[0]] = fields
 
     def _processBondType(self, line):
@@ -492,9 +503,10 @@ class GromacsTopFile(object):
                 baseAtomIndex = sys.getNumParticles()
                 atomTypes = [atom[1] for atom in moleculeType.atoms]
                 try:
-                    [self._atomTypes[t][1] for t in atomTypes]
+                    bondedTypes = [self._atomTypes[t][1] for t in atomTypes]
                 except KeyError as e:
                     raise ValueError('Unknown atom type: '+e.message)
+                bondedTypes = [b if b is not None else a for a, b in zip(atomTypes, bondedTypes)]
 
                 # Add atoms.
 
@@ -502,7 +514,7 @@ class GromacsTopFile(object):
                     if len(fields) >= 8:
                         mass = float(fields[7])
                     else:
-                        mass = float(self._atomTypes[fields[1]][2])
+                        mass = float(self._atomTypes[fields[1]][3])
                     sys.addParticle(mass)
 
                 # Add bonds.
@@ -510,7 +522,7 @@ class GromacsTopFile(object):
                 atomBonds = [{} for x in range(len(moleculeType.atoms))]
                 for fields in moleculeType.bonds:
                     atoms = [int(x)-1 for x in fields[:2]]
-                    types = tuple(atomTypes[i] for i in atoms)
+                    types = tuple(bondedTypes[i] for i in atoms)
                     if len(fields) >= 5:
                         params = fields[3:5]
                     elif types in self._bondTypes:
@@ -547,7 +559,7 @@ class GromacsTopFile(object):
                 degToRad = math.pi/180
                 for fields in moleculeType.angles:
                     atoms = [int(x)-1 for x in fields[:3]]
-                    types = tuple(atomTypes[i] for i in atoms)
+                    types = tuple(bondedTypes[i] for i in atoms)
                     if len(fields) >= 6:
                         params = fields[4:]
                     elif types in self._angleTypes:
@@ -593,7 +605,7 @@ class GromacsTopFile(object):
 
                 for fields in moleculeType.dihedrals:
                     atoms = [int(x)-1 for x in fields[:4]]
-                    types = tuple(atomTypes[i] for i in atoms)
+                    types = tuple(bondedTypes[i] for i in atoms)
                     dihedralType = fields[4]
                     reversedTypes = types[::-1]+(dihedralType,)
                     types = types+(dihedralType,)
@@ -645,7 +657,7 @@ class GromacsTopFile(object):
 
                 for fields in moleculeType.cmaps:
                     atoms = [int(x)-1 for x in fields[:5]]
-                    types = tuple(atomTypes[i] for i in atoms)
+                    types = tuple(bondedTypes[i] for i in atoms)
                     if len(fields) >= 8 and len(fields) >= 8+int(fields[6])*int(fields[7]):
                         params = fields
                     elif types in self._cmapTypes:
@@ -677,8 +689,8 @@ class GromacsTopFile(object):
                     if len(fields) > 6:
                         q = float(fields[6])
                     else:
-                        q = float(params[3])
-                    nb.addParticle(q, float(params[5]), float(params[6]))
+                        q = float(params[4])
+                    nb.addParticle(q, float(params[6]), float(params[7]))
                     if implicitSolvent is OBC2:
                         if fields[1] not in self._implicitTypes:
                             raise ValueError('No implicit solvent parameters specified for atom type: '+fields[1])