Changed gb settings to use strings rather than named constants.

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

@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
-Permission is hereby granted, free of charge, to any person obtaining a 
+Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
 to deal in the Software without restriction, including without limitation
 the rights to use, copy, modify, merge, publish, distribute, sublicense,
@@ -40,23 +40,23 @@ import simtk.unit as unit
 import simtk.openmm as mm
 
 # Enumerated values for implicit solvent model
-
-HCT = object()
-OBC1 = object()
-OBC2 = object()
-GBn = object()
+# Prefer to use strings now, but these are here for backwards compatibility
+HCT = 'HCT'
+OBC1 = 'OBC1'
+OBC2 = 'OBC2'
+GBn = 'GBn'
 
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
-    
+
     def __init__(self, file):
         """Load a prmtop file."""
         top = Topology()
         ## The Topology read from the prmtop file
         self.topology = top
-        
+
         # Load the prmtop file
-        
+
         prmtop = amber_file_parser.PrmtopLoader(file)
         self._prmtop = prmtop
 
@@ -82,7 +82,7 @@ class AmberPrmtopFile(object):
                 atomName = atomReplacements[atomName]
 
             # Try to guess the element.
-            
+
             upper = atomName.upper()
             if upper.startswith('CL'):
                 element = elem.chlorine
@@ -96,17 +96,17 @@ class AmberPrmtopFile(object):
                 except KeyError:
                     element = None
             top.addAtom(atomName, element, r)
-        
+
         # Add bonds to the topology
-        
+
         atoms = list(top.atoms())
         for bond in prmtop.getBondsWithH():
             top.addBond(atoms[bond[0]], atoms[bond[1]])
         for bond in prmtop.getBondsNoH():
             top.addBond(atoms[bond[0]], atoms[bond[1]])
-        
+
         # Set the periodic box size.
-        
+
         if prmtop.getIfBox():
             top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)
 
@@ -114,7 +114,7 @@ class AmberPrmtopFile(object):
                      constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True,
                      ewaldErrorTolerance=0.0005):
         """Construct an OpenMM System representing the topology described by this prmtop file.
-        
+
         Parameters:
          - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions.  Allowed values are
            NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
@@ -122,7 +122,8 @@ class AmberPrmtopFile(object):
          - constraints (object=None) Specifies which bonds angles should be implemented with constraints.
            Allowed values are None, HBonds, AllBonds, or HAngles.
          - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
-         - implicitSolvent (object=None) If not None, the implicit solvent model to use.  Allowed values are HCT, OBC1, OBC2, or GBn.
+         - implicitSolvent (object=None) If not None, the implicit solvent
+           model to use.  Allowed values are 'HCT', 'OBC1', 'OBC2', or 'GBn'.
          - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
@@ -148,24 +149,15 @@ class AmberPrmtopFile(object):
             constraintString = constraintMap[constraints]
         else:
             raise ValueError('Illegal value for constraints')
-        if implicitSolvent is None:
-            implicitString = None
-        elif implicitSolvent == HCT:
-            implicitString = 'HCT'
-        elif implicitSolvent == OBC1:
-            implicitString = 'OBC1'
-        elif implicitSolvent == OBC2:
-            implicitString = 'OBC2'
-        elif implicitSolvent == GBn:
-            implicitString = 'GBn'
-        else:
-            raise ValueError('Illegal value for implicit solvent model')
-        sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
-                                                 nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
-                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater)
+
+        sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString,
+                                                nonbondedCutoff=nonbondedCutoff, nonbondedMethod=methodMap[nonbondedMethod],
+                                                flexibleConstraints=False, gbmodel=implicitSolvent,
+                                                soluteDielectric=soluteDielectric, solventDielectric=solventDielectric,
+                                                rigidWater=rigidWater)
         for force in sys.getForces():
             if isinstance(force, mm.NonbondedForce):
                 force.setEwaldErrorTolerance(ewaldErrorTolerance)
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
-        return sys
+        return sys
\ No newline at end of file

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

@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
-Permission is hereby granted, free of charge, to any person obtaining a 
+Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
 to deal in the Software without restriction, including without limitation
 the rights to use, copy, modify, merge, publish, distribute, sublicense,
@@ -45,11 +45,11 @@ HBonds = ff.HBonds
 AllBonds = ff.AllBonds
 HAngles = ff.HAngles
 
-OBC2 = prmtop.OBC2
+OBC2 = 'OBC2'
 
 class GromacsTopFile(object):
     """GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it."""
-    
+
     class _MoleculeType(object):
         """Inner class to store information about a molecule type."""
         def __init__(self):
@@ -60,7 +60,7 @@ class GromacsTopFile(object):
             self.exclusions = []
             self.pairs = []
             self.cmaps = []
-    
+
     def _processFile(self, file):
         append = ''
         for line in open(file):
@@ -69,7 +69,7 @@ class GromacsTopFile(object):
             else:
                 self._processLine(append+' '+line, file)
                 append = ''
-    
+
     def _processLine(self, line, file):
         """Process one line from a file."""
         if ';' in line:
@@ -79,13 +79,13 @@ class GromacsTopFile(object):
         if stripped.startswith('*') or len(stripped) == 0:
             # A comment or empty line.
             return
-        
+
         elif stripped.startswith('[') and not ignore:
             # The start of a category.
             if not stripped.endswith(']'):
                 raise ValueError('Illegal line in .top file: '+line)
             self._currentCategory = stripped[1:-1].strip()
-                
+
         elif stripped.startswith('#'):
             # A preprocessor command.
             fields = stripped.split()
@@ -127,7 +127,7 @@ class GromacsTopFile(object):
                 if len(self._ifStack) == 0:
                     raise ValueError('Unexpected line in .top file: '+line)
                 del(self._ifStack[-1])
-                
+
         elif not ignore:
             # A line of data for the current category
             if self._currentCategory is None:
@@ -166,7 +166,7 @@ class GromacsTopFile(object):
                 self._processPairType(line)
             elif self._currentCategory == 'cmaptypes':
                 self._processCmapType(line)
-    
+
     def _processDefaults(self, line):
         """Process the [ defaults ] line."""
         fields = line.split()
@@ -179,7 +179,7 @@ class GromacsTopFile(object):
         if fields[2].lower() == 'no':
             raise ValueError('gen_pairs=no is not supported')
         self._defaults = fields
-    
+
     def _processMoleculeType(self, line):
         """Process a line in the [ moleculetypes ] category."""
         fields = line.split()
@@ -188,14 +188,14 @@ class GromacsTopFile(object):
         type = GromacsTopFile._MoleculeType()
         self._moleculeTypes[fields[0]] = type
         self._currentMoleculeType = type
-    
+
     def _processMolecule(self, line):
         """Process a line in the [ molecules ] category."""
         fields = line.split()
         if len(fields) < 2:
             raise ValueError('Too few fields in [ molecules ] line: '+line);
         self._molecules.append((fields[0], int(fields[1])))
-    
+
     def _processAtom(self, line):
         """Process a line in the [ atoms ] category."""
         if self._currentMoleculeType is None:
@@ -204,7 +204,7 @@ class GromacsTopFile(object):
         if len(fields) < 5:
             raise ValueError('Too few fields in [ atoms ] line: '+line);
         self._currentMoleculeType.atoms.append(fields)
-    
+
     def _processBond(self, line):
         """Process a line in the [ bonds ] category."""
         if self._currentMoleculeType is None:
@@ -215,7 +215,7 @@ class GromacsTopFile(object):
         if fields[2] != '1':
             raise ValueError('Unsupported function type in [ bonds ] line: '+line);
         self._currentMoleculeType.bonds.append(fields)
-    
+
     def _processAngle(self, line):
         """Process a line in the [ angles ] category."""
         if self._currentMoleculeType is None:
@@ -226,7 +226,7 @@ class GromacsTopFile(object):
         if fields[3] not in ('1', '5'):
             raise ValueError('Unsupported function type in [ angles ] line: '+line);
         self._currentMoleculeType.angles.append(fields)
-    
+
     def _processDihedral(self, line):
         """Process a line in the [ dihedrals ] category."""
         if self._currentMoleculeType is None:
@@ -237,7 +237,7 @@ class GromacsTopFile(object):
         if fields[4] not in ('1', '2', '3', '4', '9'):
             raise ValueError('Unsupported function type in [ dihedrals ] line: '+line);
         self._currentMoleculeType.dihedrals.append(fields)
-    
+
     def _processExclusion(self, line):
         """Process a line in the [ exclusions ] category."""
         if self._currentMoleculeType is None:
@@ -246,7 +246,7 @@ class GromacsTopFile(object):
         if len(fields) < 2:
             raise ValueError('Too few fields in [ exclusions ] line: '+line);
         self._currentMoleculeType.exclusions.append(fields)
-    
+
     def _processPair(self, line):
         """Process a line in the [ pairs ] category."""
         if self._currentMoleculeType is None:
@@ -257,7 +257,7 @@ class GromacsTopFile(object):
         if fields[2] != '1':
             raise ValueError('Unsupported function type in [ pairs ] line: '+line);
         self._currentMoleculeType.pairs.append(fields)
-    
+
     def _processCmap(self, line):
         """Process a line in the [ cmaps ] category."""
         if self._currentMoleculeType is None:
@@ -266,14 +266,14 @@ class GromacsTopFile(object):
         if len(fields) < 6:
             raise ValueError('Too few fields in [ pairs ] line: '+line);
         self._currentMoleculeType.cmaps.append(fields)
-    
+
     def _processAtomType(self, line):
         """Process a line in the [ atomtypes ] category."""
         fields = line.split()
         if len(fields) < 7:
             raise ValueError('Too few fields in [ atomtypes ] line: '+line);
         self._atomTypes[fields[0]] = fields
-    
+
     def _processBondType(self, line):
         """Process a line in the [ bondtypes ] category."""
         fields = line.split()
@@ -282,7 +282,7 @@ class GromacsTopFile(object):
         if fields[2] != '1':
             raise ValueError('Unsupported function type in [ bondtypes ] line: '+line);
         self._bondTypes[tuple(fields[:2])] = fields
-    
+
     def _processAngleType(self, line):
         """Process a line in the [ angletypes ] category."""
         fields = line.split()
@@ -291,7 +291,7 @@ class GromacsTopFile(object):
         if fields[3] not in ('1', '5'):
             raise ValueError('Unsupported function type in [ angletypes ] line: '+line);
         self._angleTypes[tuple(fields[:3])] = fields
-    
+
     def _processDihedralType(self, line):
         """Process a line in the [ dihedraltypes ] category."""
         fields = line.split()
@@ -305,14 +305,14 @@ class GromacsTopFile(object):
             self._dihedralTypes[key].append(fields)
         else:
             self._dihedralTypes[key] = [fields]
-    
+
     def _processImplicitType(self, line):
         """Process a line in the [ implicit_genborn_params ] category."""
         fields = line.split()
         if len(fields) < 6:
             raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line);
         self._implicitTypes[fields[0]] = fields
-    
+
     def _processPairType(self, line):
         """Process a line in the [ pairtypes ] category."""
         fields = line.split()
@@ -321,7 +321,7 @@ class GromacsTopFile(object):
         if fields[2] != '1':
             raise ValueError('Unsupported function type in [ pairtypes ] line: '+line);
         self._pairTypes[tuple(fields[:2])] = fields
-    
+
     def _processCmapType(self, line):
         """Process a line in the [ cmaptypes ] category."""
         fields = line.split()
@@ -330,10 +330,10 @@ class GromacsTopFile(object):
         if fields[5] != '1':
             raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
         self._cmapTypes[tuple(fields[:5])] = fields
-    
+
     def __init__(self, file, unitCellDimensions=None, includeDir='/usr/local/gromacs/share/gromacs/top', defines={}):
         """Load a top file.
-        
+
         Parameters:
          - file (string) the name of the file to load
          - unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
@@ -343,9 +343,9 @@ class GromacsTopFile(object):
          """
         self._includeDirs = (os.path.dirname(file), includeDir)
         self._defines = defines
-        
+
         # Parse the file.
-        
+
         self._currentCategory = None
         self._ifStack = []
         self._moleculeTypes = {}
@@ -359,9 +359,9 @@ class GromacsTopFile(object):
         self._pairTypes = {}
         self._cmapTypes = {}
         self._processFile(file)
-        
+
         # Create the Topology from it.
-        
+
         top = Topology()
         ## The Topology read from the prmtop file
         self.topology = top
@@ -371,9 +371,9 @@ class GromacsTopFile(object):
             if moleculeName not in self._moleculeTypes:
                 raise ValueError("Unknown molecule type: "+moleculeName)
             moleculeType = self._moleculeTypes[moleculeName]
-            
+
             # Create the specified number of molecules of this type.
-            
+
             for i in range(moleculeCount):
                 atoms = []
                 lastResidue = None
@@ -393,9 +393,9 @@ class GromacsTopFile(object):
                     atomName = fields[4]
                     if atomName in atomReplacements:
                         atomName = atomReplacements[atomName]
-        
+
                     # Try to guess the element.
-                    
+
                     upper = atomName.upper()
                     if upper.startswith('CL'):
                         element = elem.chlorine
@@ -409,16 +409,16 @@ class GromacsTopFile(object):
                         except KeyError:
                             element = None
                     atoms.append(top.addAtom(atomName, element, r))
-                
+
                 # Add bonds to the topology
-                
+
                 for fields in moleculeType.bonds:
                     top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True):
         """Construct an OpenMM System representing the topology described by this prmtop file.
-        
+
         Parameters:
          - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions.  Allowed values are
            NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
@@ -426,7 +426,7 @@ class GromacsTopFile(object):
          - constraints (object=None) Specifies which bonds and angles should be implemented with constraints.
            Allowed values are None, HBonds, AllBonds, or HAngles.
          - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
-         - implicitSolvent (object=None) If not None, the implicit solvent model to use.  The only allowed value is OBC2.
+         - implicitSolvent (string=None) If not None, the implicit solvent model to use.  The only allowed value is 'OBC2'.
          - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
@@ -434,7 +434,7 @@ class GromacsTopFile(object):
         Returns: the newly created System
         """
         # Create the System.
-        
+
         sys = mm.System()
         boxSize = self.topology.getUnitCellDimensions()
         if boxSize is not None:
@@ -443,7 +443,7 @@ class GromacsTopFile(object):
             raise ValueError('Illegal nonbonded method for a non-periodic system')
         nb = mm.NonbondedForce()
         sys.addForce(nb)
-        if implicitSolvent is OBC2:
+        if implicitSolvent == 'OBC2':
             gb = mm.GBSAOBCForce()
             gb.setSoluteDielectric(soluteDielectric)
             gb.setSolventDielectric(solventDielectric)
@@ -461,33 +461,33 @@ class GromacsTopFile(object):
         topologyAtoms = list(self.topology.atoms())
         exceptions = []
         fudgeQQ = float(self._defaults[4])
-        
+
         # Loop over molecules and create the specified number of each type.
-        
+
         for moleculeName, moleculeCount in self._molecules:
             moleculeType = self._moleculeTypes[moleculeName]
             for i in range(moleculeCount):
-                
+
                 # Record the types of all atoms.
-                
+
                 baseAtomIndex = sys.getNumParticles()
                 atomTypes = [atom[1] for atom in moleculeType.atoms]
                 try:
                     [self._atomTypes[t][1] for t in atomTypes]
                 except KeyError as e:
                     raise ValueError('Unknown atom type: '+e.message)
-                
+
                 # Add atoms.
-                
+
                 for fields in moleculeType.atoms:
                     if len(fields) >= 8:
                         mass = float(fields[7])
                     else:
                         mass = float(self._atomTypes[fields[1]][2])
                     sys.addParticle(mass)
-                
+
                 # Add bonds.
-                
+
                 atomBonds = [{} for x in range(len(moleculeType.atoms))]
                 for fields in moleculeType.bonds:
                     atoms = [int(x)-1 for x in fields[:2]]
@@ -522,9 +522,9 @@ class GromacsTopFile(object):
                     # Record information that will be needed for constraining angles.
                     atomBonds[atoms[0]][atoms[1]] = length
                     atomBonds[atoms[1]][atoms[0]] = length
-                
+
                 # Add angles.
-                
+
                 degToRad = math.pi/180
                 for fields in moleculeType.angles:
                     atoms = [int(x)-1 for x in fields[:3]]
@@ -571,7 +571,7 @@ class GromacsTopFile(object):
                                 bonds.addBond(baseAtomIndex+atoms[0], baseAtomIndex+atoms[2], float(params[2]), k)
 
                 # Add torsions.
-                
+
                 for fields in moleculeType.dihedrals:
                     atoms = [int(x)-1 for x in fields[:4]]
                     types = tuple(atomTypes[i] for i in atoms)
@@ -617,7 +617,7 @@ class GromacsTopFile(object):
                                     rb = mm.RBTorsionForce()
                                     sys.addForce(rb)
                                 rb.addTorsion(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c[0], c[1], c[2], c[3], c[4], c[5])
-                
+
                 # Add CMAP terms.
 
                 for fields in moleculeType.cmaps:
@@ -648,7 +648,7 @@ class GromacsTopFile(object):
                                  baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], baseAtomIndex+atoms[4])
 
                 # Set nonbonded parameters for particles.
-                
+
                 for fields in moleculeType.atoms:
                     params = self._atomTypes[fields[1]]
                     if len(fields) > 6:
@@ -656,7 +656,7 @@ class GromacsTopFile(object):
                     else:
                         q = float(params[3])
                     nb.addParticle(q, float(params[5]), float(params[6]))
-                    if implicitSolvent is OBC2:
+                    if implicitSolvent == 'OBC2':
                         if fields[1] not in self._implicitTypes:
                             raise ValueError('No implicit solvent parameters specified for atom type: '+fields[1])
                         gbparams = self._implicitTypes[fields[1]]
@@ -664,9 +664,9 @@ class GromacsTopFile(object):
                 for fields in moleculeType.bonds:
                     atoms = [int(x)-1 for x in fields[:2]]
                     bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
-                
+
                 # Record nonbonded exceptions.
-                
+
                 for fields in moleculeType.pairs:
                     atoms = [int(x)-1 for x in fields[:2]]
                     types = tuple(atomTypes[i] for i in atoms)
@@ -681,15 +681,15 @@ class GromacsTopFile(object):
                     atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
                     atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
                     exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
-        
+
         # Create nonbonded exceptions.
-        
+
         nb.createExceptionsFromBonds(bondIndices, fudgeQQ, float(self._defaults[3]))
         for exception in exceptions:
             nb.addException(exception[0], exception[1], exception[2], float(exception[3]), float(exception[4]), True)
 
         # Finish configuring the NonbondedForce.
-        
+
         methodMap = {ff.NoCutoff:mm.NonbondedForce.NoCutoff,
                      ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
                      ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,

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

@@ -16,7 +16,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
 Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
 Contributors: Christoph Klein, Michael R. Shirts
 
-Permission is hereby granted, free of charge, to any person obtaining a 
+Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
 to deal in the Software without restriction, including without limitation
 the rights to use, copy, modify, merge, publish, distribute, sublicense,
@@ -92,7 +92,7 @@ class PrmtopLoader(object):
     >>> import os, os.path
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
     >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
-    >>> prmtop = PrmtopLoader(prmtop_filename)    
+    >>> prmtop = PrmtopLoader(prmtop_filename)
 
     """
     def __init__(self, inFilename):
@@ -101,7 +101,7 @@ class PrmtopLoader(object):
 
         ARGUMENTS
 
-        inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap        
+        inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap
 
         """
 
@@ -152,7 +152,7 @@ class PrmtopLoader(object):
            Parameter:
             - pointerLabel: a string matching one of the following:
 
-            NATOM  : total number of atoms 
+            NATOM  : total number of atoms
             NTYPES : total number of distinct atom types
             NBONH  : number of bonds containing hydrogen
             MBONA  : number of bonds not containing hydrogen
@@ -183,7 +183,7 @@ class PrmtopLoader(object):
             NMXRS  : number of atoms in the largest residue
             IFCAP  : set to 1 if the CAP option from edit was specified
         """
-        index = POINTER_LABEL_LIST.index(pointerLabel) 
+        index = POINTER_LABEL_LIST.index(pointerLabel)
         return float(self._raw_data['POINTERS'][index])
 
     def getNumAtoms(self):
@@ -350,7 +350,7 @@ class PrmtopLoader(object):
         bondPointers=self._raw_data["BONDS_INC_HYDROGEN"]
         self._bondListWithH = self._getBonds(bondPointers)
         return self._bondListWithH
-        
+
 
     def getBondsNoH(self):
         """Return list of bonded atom pairs, K, and Rmin for each bond with no hydrogen"""
@@ -509,11 +509,11 @@ class PrmtopLoader(object):
 def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
     """
     Create an OpenMM System from an Amber prmtop file.
-    
+
     ARGUMENTS (specify  one or the other, but not both)
       prmtop_filename (String) - name of Amber prmtop file (new-style only)
       prmtop_loader (PrmtopLoader) - the loaded prmtop file
-      
+
     OPTIONAL ARGUMENTS
       shake (String) - if 'h-bonds', will SHAKE all bonds to hydrogen and water; if 'all-bonds', will SHAKE all bonds and water (default: None)
       gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None)
@@ -547,10 +547,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
     >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
-    >>> system = readAmberSystem(prmtop_filename)    
+    >>> system = readAmberSystem(prmtop_filename)
 
     """
-    
+
     if prmtop_filename is None and prmtop_loader is None:
         raise Exception("Must specify a filename or loader")
     if prmtop_filename is not None and prmtop_loader is not None:
@@ -567,7 +567,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     if prmtop.getIfPert()>0:
         raise Exception("perturbation not currently supported")
-        
+
     if prmtop.getIfBox()>1:
         raise Exception("only standard periodic boxes are currently supported")
 
@@ -596,9 +596,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
             if isWater[iAtom] and isWater[jAtom]:
                 system.addConstraint(iAtom, jAtom, rMin)
-            
+
     # Add harmonic bonds.
-    if verbose: print "Adding bonds..."    
+    if verbose: print "Adding bonds..."
     force = mm.HarmonicBondForce()
     if flexibleConstraints or (shake not in ('h-bonds', 'all-bonds', 'h-angles')):
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
@@ -610,7 +610,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     system.addForce(force)
 
     # Add harmonic angles.
-    if verbose: print "Adding angles..."    
+    if verbose: print "Adding angles..."
     force = mm.HarmonicAngleForce()
     if shake == 'h-angles':
         numConstrainedBonds = system.getNumConstraints()
@@ -638,7 +638,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
                     l1 = bond[1]
                 elif bond[0] == kAtom:
                     l2 = bond[1]
-            
+
             # Compute the distance between atoms and add a constraint
             length = math.sqrt(l1*l1 + l2*l2 - 2*l1*l2*math.cos(aMin))
             system.addConstraint(iAtom, kAtom, length)
@@ -647,14 +647,14 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     system.addForce(force)
 
     # Add torsions.
-    if verbose: print "Adding torsions..."    
+    if verbose: print "Adding torsions..."
     force = mm.PeriodicTorsionForce()
     for (iAtom, jAtom, kAtom, lAtom, forceConstant, phase, periodicity) in prmtop.getDihedrals():
         force.addTorsion(iAtom, jAtom, kAtom, lAtom, periodicity, phase, forceConstant)
     system.addForce(force)
 
     # Add nonbonded interactions.
-    if verbose: print "Adding nonbonded interactions..."    
+    if verbose: print "Adding nonbonded interactions..."
     force = mm.NonbondedForce()
     if (prmtop.getIfBox() == 0):
         # System is non-periodic.
@@ -663,12 +663,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         elif nonbondedMethod == 'CutoffNonPeriodic':
             if nonbondedCutoff is None:
                 raise Exception("No cutoff value specified")
-            force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)            
+            force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
             force.setCutoffDistance(nonbondedCutoff)
         else:
             raise Exception("Illegal nonbonded method for a non-periodic system")
     else:
-        # System is periodic. 
+        # System is periodic.
         # Set periodic box vectors for periodic system
         (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
         d0 = units.Quantity(0.0, units.angstroms)
@@ -676,16 +676,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         yVec = units.Quantity((d0,   boxY, d0))
         zVec = units.Quantity((d0,   d0,   boxZ))
         system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
-        
+
         # Set cutoff.
         if nonbondedCutoff is None:
             # Compute cutoff automatically.
             min_box_width = min([boxX / units.nanometers, boxY / units.nanometers, boxZ / units.nanometers])
-            CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length            
+            CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length
             nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers)
         if nonbondedMethod != 'NoCutoff':
             force.setCutoffDistance(nonbondedCutoff)
-        
+
         # Set nonbonded method.
         if nonbondedMethod == 'NoCutoff':
             force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
@@ -724,7 +724,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     excludeParams = (0.0, 0.1, 0.0)
     for iAtom in range(prmtop.getNumAtoms()):
         for jAtom in excludedAtoms[iAtom]:
-            if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue            
+            if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue
             force.addException(iAtom, jAtom, excludeParams[0], excludeParams[1], excludeParams[2])
 
     system.addForce(force)
@@ -812,7 +812,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             if gbmodel == 'OBC2':
                 gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1])
             else:
-                gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])        
+                gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
         system.addForce(gb)
         if nonbondedMethod == 'NoCutoff':
             gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
@@ -835,7 +835,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
 def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbose=False, asNumpy=False):
     """
-    Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.    
+    Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.
 
     ARGUMENTS
 
@@ -852,13 +852,13 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
     Read coordinates in vacuum.
 
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
-    >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')    
+    >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
     >>> coordinates = readAmberCoordinates(crd_filename)
 
     Read coordinates in solvent.
 
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
-    >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')    
+    >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
     >>> [coordinates, box_vectors] = readAmberCoordinates(crd_filename, read_box=True)
 
     """
@@ -921,7 +921,7 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
             velocities = newvel
         # Assign units.
         velocities = units.Quantity(velocities, units.angstroms/units.picoseconds)
-            
+
     # Read box size if present
     box_vectors = None
     if (read_box):
@@ -943,20 +943,20 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
             else:
                 a = units.Quantity(mm.Vec3(box_dimensions[0], 0.0, 0.0), units.angstroms)
                 b = units.Quantity(mm.Vec3(0.0, box_dimensions[1], 0.0), units.angstroms)
-                c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)            
+                c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)
             box_vectors = [a,b,c]
         else:
             raise Exception("Don't know what to do with box vectors: %s" % line)
 
     # Close file
     infile.close()
-    
+
     if box_vectors and velocities:
         return (coordinates, box_vectors, velocities)
     if box_vectors:
         return (coordinates, box_vectors)
     if velocities:
-        return (coordinates, velocities)    
+        return (coordinates, velocities)
     return coordinates
 
 #=============================================================================================