Moved OpenMM app into the main OpenMM project

wrappers/python/CMakeLists.txt

@@ -36,6 +36,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.py"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.i"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
     )
     foreach(file ${STAGING_INPUT_FILES1})
         set(STAGING_INPUT_FILES ${STAGING_INPUT_FILES} "${file}")

wrappers/python/setup.py

@@ -81,11 +81,15 @@ def buildKeywordDictionary(major_version_num=MAJOR_VERSION_NUM,
     setupKeywords["download_url"]      = "https://simtk.org/home/openmm"
     setupKeywords["packages"]          = ["simtk",
                                           "simtk.unit",
-                                          "simtk.openmm"]
+                                          "simtk.openmm",
+                                          "simtk.openmm.app",
+                                          "simtk.openmm.app.internal"]
     setupKeywords["data_files"]        = []
     setupKeywords["package_data"]      = {"simtk" : [],
                                           "simtk.unit" : [],
-                                          "simtk.openmm" : []}
+                                          "simtk.openmm" : [],
+                                          "simtk.openmm.app" : ['data/*.xml'],
+                                          "simtk.openmm.app.internal" : []}
     setupKeywords["platforms"]         = ["Linux", "Mac OS X", "Windows"]
     setupKeywords["description"]       = \
     "Python wrapper for OpenMM (a C++ MD package)"

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

+"""
+OpenMM Application
+"""
+__docformat__ = "epytext en"
+
+__author__ = "Peter Eastman"
+__copyright__ = "Copyright 2011, Stanford University and Peter Eastman"
+__credits__ = []
+__license__ = "MIT"
+__maintainer__ = "Peter Eastman"
+__email__ = "peastman@stanford.edu"
+
+from topology import Topology, Chain, Residue, Atom
+from pdbfile import PDBFile
+from forcefield import ForceField
+from simulation import Simulation
+from pdbreporter import PDBReporter
+from amberprmtopfile import AmberPrmtopFile
+from amberinpcrdfile import AmberInpcrdFile
+from dcdfile import DCDFile
+from dcdreporter import DCDReporter
+
+# Enumerated values
+
+NoCutoff = forcefield.NoCutoff
+CutoffNonPeriodic = forcefield.CutoffNonPeriodic
+CutoffPeriodic = forcefield.CutoffPeriodic
+Ewald = forcefield.Ewald
+PME = forcefield.PME
+
+HBonds = forcefield.HBonds
+AllBonds = forcefield.AllBonds
+HAngles = forcefield.HAngles
+
+OBC = amberprmtopfile.OBC

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

+"""
+armberinpcrdfile.py: Used for loading AMBER inpcrd files.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from simtk.openmm.app.internal import amber_file_parser
+
+class AmberInpcrdFile(object):
+    """AmberInpcrdFile parses an AMBER inpcrd file and loads the data stored in it."""
+    
+    def __init__(self, file, loadVelocities=False, loadBoxVectors=False):
+        """Load an inpcrd file.
+        
+        An inpcrd file contains atom positions and, optionally, velocities and periodic box dimensions.
+        Unfortunately, it is sometimes impossible to determine from the file itself exactly what data
+        it contains.  You therefore must specify in advance what data to load.  It is stored into this
+        object's "positions", "velocities", and "boxVectors" fields.
+        
+        Parameters:
+         - file (string) the name of the file to load
+         - loadVelocities (boolean=False) whether to load velocities from the file
+         - loadBoxVectors (boolean=False) whether to load the periodic box vectors
+        """
+        results = amber_file_parser.readAmberCoordinates(file, read_velocities=loadVelocities, read_box=loadBoxVectors)
+        if loadVelocities:
+            self.positions = results[0]
+            if loadBoxVectors:
+                self.boxVectors = results[1]
+                self.velocities = results[2]
+            else:
+                self.velocities = results[1]
+        elif loadBoxVectors:
+            self.positions = results[0]
+            self.boxVectors = results[1]
+        else:
+            self.positions = results

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

+"""
+armberprmtopfile.py: Used for loading AMBER prmtop files.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from simtk.openmm.app import Topology
+from simtk.openmm.app.internal import amber_file_parser
+import forcefield as ff
+import element as elem
+import simtk.unit as unit
+
+# Enumerated values for implicit solvent model
+
+OBC = object()
+
+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()
+        self.topology = top
+        
+        # Load the prmtop file
+        
+        prmtop = amber_file_parser.PrmtopLoader(file)
+        self.prmtop = prmtop
+
+        # Add atoms to the topology
+
+        lastResidue = None
+        c = top.addChain()
+        for index in range(prmtop.getNumAtoms()):
+            resNumber = prmtop.getResidueNumber(index)
+            if resNumber != lastResidue:
+                lastResidue = resNumber
+                resName = prmtop.getResidueLabel(iAtom=index)
+                r = top.addResidue(resName, c)
+            atomName = prmtop.getAtomName(index)
+
+            # Try to guess the element.
+            
+            upper = atomName.upper()
+            if upper.startswith('CL'):
+                element = elem.chlorine
+            elif upper.startswith('NA'):
+                element = elem.sodium
+            elif upper.startswith('MG'):
+                element = elem.magnesium
+            else:
+                try:
+                    element = elem.get_by_symbol(atomName[0])
+                except KeyError:
+                    pass
+            top.addAtom(atomName, element, r)
+        
+        # Add bonds to the topology
+        
+        for bond in prmtop.getBondsWithH():
+            top.addBond(bond[0], bond[1])
+        for bond in prmtop.getBondsNoH():
+            top.addBond(bond[0], bond[1])
+
+    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
+                     constraints=None, rigidWater=True, implicitSolvent=None):
+        """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.
+         - nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions
+         - 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
+        Returns: the newly created System
+        """
+        methodMap = {ff.NoCutoff:'NoCutoff',
+                     ff.CutoffNonPeriodic:'CutoffNonPeriodic',
+                     ff.CutoffPeriodic:'CutoffPeriodic',
+                     ff.Ewald:'Ewald',
+                     ff.PME:'PME'}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal value for nonbonded method')
+        if not self.prmtop.getIfBox() and nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME):
+            raise ValueError('Illegal nonbonded method for a non-periodic system')
+        if nonbondedMethod == ff.NoCutoff:
+            nonbondedCutoff = None
+        constraintMap = {None:None,
+                         ff.HBonds:'h-bonds',
+                         ff.AllBonds:'all-bonds',
+                         ff.HAngles:'h-angles'}
+        if constraints == ff.NoConstraints:
+            constraintString = None
+        elif constraints in constraintMap:
+            constraintString = constraintMap[constraints]
+        else:
+            raise ValueError('Illegal value for constraints')
+        if implicitSolvent is None:
+            implicitString = None
+        elif implicitSolvent == OBC:
+            implicitString = 'OBC'
+        else:
+            raise ValueError('Illegal value for implicit solvent model')
+        return amber_file_parser.readAmberSystem(prmtop_loader=self.prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
+                                                 nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString)
\ No newline at end of file

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

+"""
+dcdfile.py: Used for writing DCD files.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import array
+import os
+import time
+import struct
+from simtk.unit import picoseconds, nanometers, angstroms, is_quantity
+    
+class DCDFile(object):
+    """DCDFile provides methods for creating DCD files.
+    
+    DCD is a file format for storing simulation trajectories.  It is supported by many programs, such
+    as CHARMM, NAMD, and X-PLOR.  Note, however, that different programs produce subtly different
+    versions of the format.  This class generates the CHARMM version.  Also note that there is no
+    standard byte ordering (big-endian or little-endian) for this format.  This class always generates
+    files with little-endian ordering.
+    
+    To use this class, create a DCDFile object, then call writeModel() once for each model in the file."""
+    
+    def __init__(self, file, topology, dt, firstStep=0, interval=1):
+        """Create a DCD file and write out the header.
+        
+        Parameters:
+         - file (file) A file to write to
+         - topology (Topology) The Topology defining the molecular system being written
+         - dt (time) The time step used in the trajectory
+         - firstStep (int=0) The index of the first step in the trajectory
+         - interval (int=1) The frequency (measured in time steps) at which states are written to the trajectory
+        """
+        self._file = file
+        self._topology = topology
+        self._firstStep = firstStep
+        self._interval = interval
+        self._modelCount = 0
+        if is_quantity(dt):
+            dt = dt.value_in_unit(picoseconds)
+        dt /= 0.04888821
+        boxFlag = 0
+        if topology.getUnitCellDimensions() is not None:
+            boxFlag = 1
+        header = struct.pack('<i4c9if', 84, 'C', 'O', 'R', 'D', 0, firstStep, interval, 0, 0, 0, 0, 0, 0, dt)
+        header += struct.pack('<13i', boxFlag, 0, 0, 0, 0, 0, 0, 0, 0, 24, 84, 164, 2)
+        header += struct.pack('<80s', 'Created by OpenMM')
+        header += struct.pack('<80s', 'Created '+time.asctime(time.localtime(time.time())))
+        header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4)
+        file.write(header)
+    
+    def writeModel(self, positions):
+        """Write out a model to the DCD file.
+                
+        Parameters:
+         - positions (list) The list of atomic positions to write
+        """
+        if len(list(self._topology.atoms())) != len(positions):
+            raise ValueError('The number of positions must match the number of atoms') 
+        if is_quantity(positions):
+            positions = positions.value_in_unit(nanometers)
+        file = self._file
+        
+        # Update the header.
+        
+        self._modelCount += 1
+        file.seek(8, os.SEEK_SET)
+        file.write(struct.pack('<i', self._modelCount))
+        file.seek(20, os.SEEK_SET)
+        file.write(struct.pack('<i', self._firstStep+self._modelCount*self._interval))
+        
+        # Write the data.
+        
+        file.seek(0, os.SEEK_END)
+        boxSize = self._topology.getUnitCellDimensions()
+        if boxSize is not None:
+            size = boxSize.value_in_unit(angstroms)
+            file.write(struct.pack('<i6di', 48, size[0], 0, size[1], 0, 0, size[2], 48))
+        length = struct.pack('<i', 4*len(positions))
+        for i in range(3):
+            file.write(length)
+            data = array.array('f', (10*x[i] for x in positions))
+            data.tofile(file)
+            file.write(length)

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

+"""
+dcdreporter.py: Outputs simulation trajectories in DCD format
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import simtk.openmm as mm
+from simtk.openmm.app import DCDFile
+    
+class DCDReporter(object):
+    """DCDReporter outputs a series of frames from a Simulation to a DCD file.
+    
+    To use it, create a DCDReporter, than add it to the Simulation's list of reporters.
+    """
+    
+    def __init__(self, file, reportInterval):
+        """Create a DCDReporter.
+    
+        Parameters:
+         - file (string) The file to write to
+         - reportInterval (int) The interval (in time steps) at which to write frames
+        """
+        self._reportInterval = reportInterval
+        self._out = open(file, 'wb')
+        self._dcd = None
+    
+    def describeNextReport(self, simulation):
+        """Get information about the next report this object will generate.
+        
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+        Returns: A five element tuple.  The first element is the number of steps until the
+        next report.  The remaining elements specify whether that report will require
+        positions, velocities, forces, and energies respectively.
+        """
+        steps = self._reportInterval - simulation.currentStep%self._reportInterval
+        return (steps, True, False, False, False)
+    
+    def report(self, simulation, state):
+        """Generate a report.
+        
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+         - state (State) The current state of the simulation
+        """
+        if self._dcd is None:
+            self._dcd = DCDFile(self._out, simulation.topology, simulation.integrator.getStepSize(), 0, self._reportInterval)
+        self._dcd.writeModel(state.getPositions())
+    
+    def __del__(self):
+        self._out.close()

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

+#!/bin/env python
+
+"""
+element.py: Used for managing elements.
+"""
+__author__ = "Christopher M. Bruns"
+__version__ = "1.0"
+
+
+from simtk.unit import daltons
+
+class Element:
+    elements_by_symbol = {}
+    
+    def __init__(self, number, name, symbol, mass):
+        self.atomic_number = number
+        self.name = name
+        self.symbol = symbol
+        self.mass = mass
+        # Index this element in a global table
+        s = symbol.strip().upper()
+        assert s not in Element.elements_by_symbol
+        Element.elements_by_symbol[s] = self
+
+
+def get_by_symbol(symbol):
+    s = symbol.strip().upper()
+    return Element.elements_by_symbol[s]
+
+
+hydrogen =       Element(  1, "hydrogen", "H", 1.007947*daltons)
+deuterium =      Element(  1, "deuterium", "D", 2.01355321270*daltons)
+helium =         Element(  2, "helium", "He", 4.003*daltons)
+lithium =        Element(  3, "lithium", "Li", 6.9412*daltons)
+beryllium =      Element(  4, "beryllium", "Be", 9.0121823*daltons)
+boron =          Element(  5, "boron", "B", 10.8117*daltons)
+carbon =         Element(  6, "carbon", "C", 12.01078*daltons)
+nitrogen =       Element(  7, "nitrogen", "N", 14.00672*daltons)
+oxygen =         Element(  8, "oxygen", "O", 15.99943*daltons)
+fluorine =       Element(  9, "fluorine", "F", 18.99840325*daltons)
+neon =           Element( 10, "neon", "Ne", 20.17976*daltons)
+sodium =         Element( 11, "sodium", "Na", 22.989769282*daltons)
+magnesium =      Element( 12, "magnesium", "Mg", 24.30506*daltons)
+aluminum =       Element( 13, "aluminum", "Al", 26.98153868*daltons)
+silicon =        Element( 14, "silicon", "Si", 28.08553*daltons)
+phosphorus =     Element( 15, "phosphorus", "P", 30.9737622*daltons)
+sulfur =         Element( 16, "sulfur", "S", 32.0655*daltons)
+chlorine =       Element( 17, "chlorine", "Cl", 35.4532*daltons)
+argon =          Element( 18, "argon", "Ar", 39.9481*daltons)
+potassium =      Element( 19, "potassium", "K", 39.09831*daltons)
+calcium =        Element( 20, "calcium", "Ca", 40.0784*daltons)
+scandium =       Element( 21, "scandium", "Sc", 44.9559126*daltons)
+titanium =       Element( 22, "titanium", "Ti", 47.8671*daltons)
+vanadium =       Element( 23, "vanadium", "V", 50.94151*daltons)
+chromium =       Element( 24, "chromium", "Cr", 51.99616*daltons)
+manganese =      Element( 25, "manganese", "Mn", 54.9380455*daltons)
+iron =           Element( 26, "iron", "Fe", 55.8452*daltons)
+cobalt =         Element( 27, "cobalt", "Co", 58.9331955*daltons)
+nickel =         Element( 28, "nickel", "Ni", 58.69342*daltons)
+copper =         Element( 29, "copper", "Cu", 63.5463*daltons)
+zinc =           Element( 30, "zinc", "Zn", 65.4094*daltons)
+gallium =        Element( 31, "gallium", "Ga", 69.7231*daltons)
+germanium =      Element( 32, "germanium", "Ge", 72.641*daltons)
+arsenic =        Element( 33, "arsenic", "As", 74.921602*daltons)
+selenium =       Element( 34, "selenium", "Se", 78.963*daltons)
+bromine =        Element( 35, "bromine", "Br", 79.9041*daltons)
+krypton =        Element( 36, "krypton", "Kr", 83.7982*daltons)
+rubidium =       Element( 37, "rubidium", "Rb", 85.46783*daltons)
+strontium =      Element( 38, "strontium", "Sr", 87.621*daltons)
+yttrium =        Element( 39, "yttrium", "Y", 88.905852*daltons)
+zirconium =      Element( 40, "zirconium", "Zr", 91.2242*daltons)
+niobium =        Element( 41, "niobium", "Nb", 92.906382*daltons)
+molybdenum =     Element( 42, "molybdenum", "Mo", 95.942*daltons)
+technetium =     Element( 43, "technetium", "Tc", 98*daltons)
+ruthenium =      Element( 44, "ruthenium", "Ru", 101.072*daltons)
+rhodium =        Element( 45, "rhodium", "Rh", 102.905502*daltons)
+palladium =      Element( 46, "palladium", "Pd", 106.421*daltons)
+silver =         Element( 47, "silver", "Ag", 107.86822*daltons)
+cadmium =        Element( 48, "cadmium", "Cd", 112.4118*daltons)
+indium =         Element( 49, "indium", "In", 114.8183*daltons)
+tin =            Element( 50, "tin", "Sn", 118.7107*daltons)
+antimony =       Element( 51, "antimony", "Sb", 121.7601*daltons)
+tellurium =      Element( 52, "tellurium", "Te", 127.603*daltons)
+iodine =         Element( 53, "iodine", "I", 126.904473*daltons)
+xenon =          Element( 54, "xenon", "Xe", 131.2936*daltons)
+cesium =         Element( 55, "cesium", "Cs", 132.90545192*daltons)
+barium =         Element( 56, "barium", "Ba", 137.3277*daltons)
+lanthanum =      Element( 57, "lanthanum", "La", 138.905477*daltons)
+cerium =         Element( 58, "cerium", "Ce", 140.1161*daltons)
+praseodymium =   Element( 59, "praseodymium", "Pr", 140.907652*daltons)
+neodymium =      Element( 60, "neodymium", "Nd", 144.2423*daltons)
+promethium =     Element( 61, "promethium", "Pm", 145*daltons)
+samarium =       Element( 62, "samarium", "Sm", 150.362*daltons)
+europium =       Element( 63, "europium", "Eu", 151.9641*daltons)
+gadolinium =     Element( 64, "gadolinium", "Gd", 157.253*daltons)
+terbium =        Element( 65, "terbium", "Tb", 158.925352*daltons)
+dysprosium =     Element( 66, "dysprosium", "Dy", 162.5001*daltons)
+holmium =        Element( 67, "holmium", "Ho", 164.930322*daltons)
+erbium =         Element( 68, "erbium", "Er", 167.2593*daltons)
+thulium =        Element( 69, "thulium", "Tm", 168.934212*daltons)
+ytterbium =      Element( 70, "ytterbium", "Yb", 173.043*daltons)
+lutetium =       Element( 71, "lutetium", "Lu", 174.9671*daltons)
+hafnium =        Element( 72, "hafnium", "Hf", 178.492*daltons)
+tantalum =       Element( 73, "tantalum", "Ta", 180.947882*daltons)
+tungsten =       Element( 74, "tungsten", "W", 183.841*daltons)
+rhenium =        Element( 75, "rhenium", "Re", 186.2071*daltons)
+osmium =         Element( 76, "osmium", "Os", 190.233*daltons)
+iridium =        Element( 77, "iridium", "Ir", 192.2173*daltons)
+platinum =       Element( 78, "platinum", "Pt", 195.0849*daltons)
+gold =           Element( 79, "gold", "Au", 196.9665694*daltons)
+mercury =        Element( 80, "mercury", "Hg", 200.592*daltons)
+thallium =       Element( 81, "thallium", "Tl", 204.38332*daltons)
+lead =           Element( 82, "lead", "Pb", 207.21*daltons)
+bismuth =        Element( 83, "bismuth", "Bi", 208.980401*daltons)
+polonium =       Element( 84, "polonium", "Po", 209*daltons)
+astatine =       Element( 85, "astatine", "At", 210*daltons)
+radon =          Element( 86, "radon", "Rn", 222.018*daltons)
+francium =       Element( 87, "francium", "Fr", 223*daltons)
+radium =         Element( 88, "radium", "Ra", 226*daltons)
+actinium =       Element( 89, "actinium", "Ac", 227*daltons)
+thorium =        Element( 90, "thorium", "Th", 232.038062*daltons)
+protactinium =   Element( 91, "protactinium", "Pa", 231.035882*daltons)
+uranium =        Element( 92, "uranium", "U", 238.028913*daltons)
+neptunium =      Element( 93, "neptunium", "Np", 237*daltons)
+plutonium =      Element( 94, "plutonium", "Pu", 244*daltons)
+americium =      Element( 95, "americium", "Am", 243*daltons)
+curium =         Element( 96, "curium", "Cm", 247*daltons)
+berkelium =      Element( 97, "berkelium", "Bk", 247*daltons)
+californium =    Element( 98, "californium", "Cf", 251*daltons)
+einsteinium =    Element( 99, "einsteinium", "Es", 252*daltons)
+fermium =        Element(100, "fermium", "Fm", 257*daltons)
+mendelevium =    Element(101, "mendelevium", "Md", 258*daltons)
+nobelium =       Element(102, "nobelium", "No", 259*daltons)
+lawrencium =     Element(103, "lawrencium",     "Lr", 262*daltons)
+rutherfordium =  Element(104, "rutherfordium",  "Rf", 261*daltons)
+dubnium =        Element(105, "dubnium",        "Db", 262*daltons)
+seaborgium =     Element(106, "seaborgium",     "Sg", 266*daltons)
+bohrium =        Element(107, "bohrium",        "Bh", 264*daltons)
+hassium =        Element(108, "hassium",        "Hs", 269*daltons)
+meitnerium =     Element(109, "meitnerium",     "Mt", 268*daltons)
+darmstadtium =   Element(110, "darmstadtium",   "Ds", 281*daltons)
+roentgenium =    Element(111, "roentgenium",    "Rg", 272*daltons)
+ununbium =       Element(112, "ununbium",       "Uub", 285*daltons)
+ununtrium =      Element(113, "ununtrium",      "Uut", 284*daltons)
+ununquadium =    Element(114, "ununquadium",    "Uuq", 289*daltons)
+ununpentium =    Element(115, "ununpentium",    "Uup", 288*daltons)
+ununhexium =     Element(116, "ununhexium",     "Uuh", 292*daltons)

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

+"""
+forcefield.py: Constructs OpenMM System objects based on a Topology and an XML force field description
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import os
+import itertools
+import xml.etree.ElementTree as etree
+from math import sqrt, cos
+import simtk.openmm as mm
+import simtk.unit as unit
+import element as elem
+from simtk.openmm.app import Topology
+
+# Enumerated values for nonbonded method
+
+NoCutoff = object()
+CutoffNonPeriodic = object()
+CutoffPeriodic = object()
+Ewald = object()
+PME = object()
+
+# Enumerated values for constraint type
+
+HBonds = object()
+AllBonds = object()
+HAngles = object()
+
+# A map of functions to parse elements of the XML file.
+
+parsers = {}
+
+class ForceField(object):
+    """A ForceField constructs OpenMM System objects based on a Topology."""
+
+    def __init__(self, *files):
+        """Load one or more XML files and create a ForceField object based on them.
+        
+        Parameters:
+         - A list of XML files defining the force field.  Each entry may be an absolute file path, a path relative to the
+           current working directory, or a path relative to this module's data subdirectory (for built in force fields).
+        """
+        self._atomTypes = {}
+        self._templates = {}
+        self._templateSignatures = {}
+        self._atomClasses = {}
+        self._forces = []
+        for file in files:
+            try:
+                tree = etree.parse(file)
+            except IOError:
+                tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
+            root = tree.getroot()
+            
+            # Load the atom types.
+            
+            if tree.getroot().find('AtomTypes') is not None:
+                for type in tree.getroot().find('AtomTypes').findall('Type'):
+                    self._atomTypes[type.attrib['name']] = (type.attrib['class'], float(type.attrib['mass']), elem.get_by_symbol(type.attrib['element']))
+            
+            # Load the residue templates.
+            
+            if tree.getroot().find('Residues') is not None:
+                for residue in root.find('Residues').findall('Residue'):
+                    resName = residue.attrib['name']
+                    template = ForceField._TemplateData(resName)
+                    self._templates[resName] = template
+                    for atom in residue.findall('Atom'):
+                        template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
+                    for bond in residue.findall('Bond'):
+                        b = (int(bond.attrib['from']), int(bond.attrib['to']))
+                        template.bonds.append(b)
+                        template.atoms[b[0]].bondedTo.append(b[1])
+                        template.atoms[b[1]].bondedTo.append(b[0])
+                    for bond in residue.findall('ExternalBond'):
+                        b = int(bond.attrib['from'])
+                        template.externalBonds.append(b)
+                        template.atoms[b].externalBonds += 1
+            for template in self._templates.values():
+                template.signature = _createResidueSignature([atom.element for atom in template.atoms])
+                sigString = _signatureToString(template.signature)
+                if sigString in self._templateSignatures:
+                    self._templateSignatures[sigString].append(template)
+                else:
+                    self._templateSignatures[sigString] = [template]
+            
+            # Build sets of every atom type belonging to each class
+            
+            for type in self._atomTypes:
+                atomClass = self._atomTypes[type][0]
+                if atomClass in self._atomClasses:
+                    typeSet = self._atomClasses[atomClass]
+                else:
+                    typeSet = set()
+                    self._atomClasses[atomClass] = typeSet
+                typeSet.add(type)
+            self._atomClasses[''] = self._atomTypes.keys()
+            
+            # Load force definitions
+            
+            for child in root:
+                if child.tag in parsers:
+                    parsers[child.tag](child, self)
+
+    def _findAtomTypes(self, node, num):
+        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
+        types = []
+        attrib = node.attrib
+        for i in range(num):
+            if num == 1:
+                suffix = ''
+            else:
+                suffix = str(i+1)
+            classAttrib = 'class'+suffix
+            if classAttrib in attrib:
+                if attrib[classAttrib] not in self._atomClasses:
+                    return None # Unknown atom class
+                types.append(self._atomClasses[attrib[classAttrib]])
+            else:
+                typeAttrib = 'type'+suffix
+                if typeAttrib not in attrib or attrib[typeAttrib] not in self._atomTypes:
+                    return None # Unknown atom type
+                types.append([attrib[typeAttrib]])
+        return types
+
+    def _parseTorsion(self, node):
+        """Parse the node defining a torsion."""
+        types = self._findAtomTypes(node, 4)
+        if types is None:
+            return None
+        torsion = PeriodicTorsion(types)
+        attrib = node.attrib
+        index = 1
+        while 'phase%d'%index in attrib:
+            torsion.periodicity.append(int(attrib['periodicity%d'%index]))
+            torsion.phase.append(float(attrib['phase%d'%index]))
+            torsion.k.append(float(attrib['k%d'%index]))
+            index += 1
+        return torsion
+        
+    class _SystemData:
+        """Inner class used to encapsulate data about the system being created."""
+        def __init__(self):
+            self.atomType = {}
+            self.atoms = []
+            self.bonds = []
+            self.angles = []
+            self.propers = []
+            self.impropers = []
+            self.atomBonds = []
+            self.isAngleConstrained = []
+
+    class _TemplateData:
+        """Inner class used to encapsulate data about a residue template definition."""
+        def __init__(self, name):
+            self.name = name
+            self.atoms = []
+            self.bonds = []
+            self.externalBonds = []
+
+    class _TemplateAtomData:
+        """Inner class used to encapsulate data about an atom in a residue template definition."""
+        def __init__(self, name, type, element):
+            self.name = name
+            self.type = type
+            self.element = element
+            self.bondedTo = []
+            self.externalBonds = 0
+
+    class _BondData:
+        """Inner class used to encapsulate data about a bond."""
+        def __init__(self, atom1, atom2):
+            self.atom1 = atom1
+            self.atom2 = atom2
+            self.isConstrained = False
+            self.length = 0.0
+
+    def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
+                     constraints=None, rigidWater=True, **args):
+        """Construct an OpenMM System representing a Topology with this force field.
+        
+        Parameters:
+         - topology (Topology) The Topology for which to create a System
+         - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions.  Allowed values are
+           NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
+         - nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions
+         - 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
+         - Arbitrary additional keyword arguments may also be specified.  This allows extra parameters to be specified that are specific to
+           particular force fields.
+        Returns: the newly created System
+        """
+        
+        # Record atom indices
+        
+        data = ForceField._SystemData()
+        atomIndices = {}
+        for index, atom in enumerate(topology.atoms()):
+            data.atoms.append(atom)
+            atomIndices[atom] = index
+
+        # Make a list of all bonds
+        
+        for bond in topology.bonds():
+            if bond[0] in atomIndices and bond[1] in atomIndices:
+                data.bonds.append(ForceField._BondData(atomIndices[bond[0]], atomIndices[bond[1]]))
+
+        # Record which atoms are bonded to each other atom
+        
+        bondedToAtom = []
+        for i in range(len(data.atoms)):
+            bondedToAtom.append(set())
+            data.atomBonds.append([])
+        for i in range(len(data.bonds)):
+            bond = data.bonds[i]
+            bondedToAtom[bond.atom1].add(bond.atom2)
+            bondedToAtom[bond.atom2].add(bond.atom1)
+            data.atomBonds[bond.atom1].append(i)
+            data.atomBonds[bond.atom2].append(i)
+
+        # Find the template matching each residue and assign atom types.
+        
+        for chain in topology.chains():
+            for res in chain.residues():
+                signature = _signatureToString(_createResidueSignature([atom.element for atom in res.atoms()]))
+                template = None
+                matches = None
+                if signature in self._templateSignatures:
+                    for t in self._templateSignatures[signature]:
+                        matches = _matchResidue(res, t, bondedToAtom, atomIndices)
+                        if matches is not None:
+                            template = t
+                            break
+                if matches is None:
+                    raise ValueError('No template found for residue %d (%s)' % (res.index+1, res.name))
+                for atom, match in zip(res.atoms(), matches):
+                    data.atomType[atom] = template.atoms[match].type
+
+        # Create the System and add atoms
+        
+        sys = mm.System()
+        for atom in topology.atoms():
+            sys.addParticle(self._atomTypes[data.atomType[atom]][1])
+        
+        # Set periodic boundary conditions.
+        
+        boxSize = topology.getUnitCellDimensions()
+        if boxSize is not None:
+            sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2]))
+        elif nonbondedMethod is not NoCutoff and nonbondedMethod is not CutoffNonPeriodic:
+            raise ValueError('Requested periodic boundary conditions for a Topology that does not specify periodic box dimensions')
+
+        # Make a list of all unique angles
+        
+        uniqueAngles = set()
+        for bond in data.bonds:
+            for atom in bondedToAtom[bond.atom1]:
+                if atom != bond.atom2:
+                    if atom < bond.atom2:
+                        uniqueAngles.add((atom, bond.atom1, bond.atom2))
+                    else:
+                        uniqueAngles.add((bond.atom2, bond.atom1, atom))
+            for atom in bondedToAtom[bond.atom2]:
+                if atom != bond.atom1:
+                    if atom > bond.atom1:
+                        uniqueAngles.add((bond.atom1, bond.atom2, atom))
+                    else:
+                        uniqueAngles.add((atom, bond.atom2, bond.atom1))
+        data.angles = sorted(list(uniqueAngles))
+        
+        # Make a list of all unique proper torsions
+        
+        uniquePropers = set()
+        for angle in data.angles:
+            for atom in bondedToAtom[angle[0]]:
+                if atom != angle[1]:
+                    if atom < angle[2]:
+                        uniquePropers.add((atom, angle[0], angle[1], angle[2]))
+                    else:
+                        uniquePropers.add((angle[2], angle[1], angle[0], atom))
+            for atom in bondedToAtom[angle[2]]:
+                if atom != angle[1]:
+                    if atom > angle[0]:
+                        uniquePropers.add((angle[0], angle[1], angle[2], atom))
+                    else:
+                        uniquePropers.add((atom, angle[2], angle[1], angle[0]))
+        data.propers = sorted(list(uniquePropers))
+        
+        # Make a list of all unique improper torsions
+        
+        for atom in range(len(bondedToAtom)):
+            bondedTo = bondedToAtom[atom]
+            if len(bondedTo) > 2:
+                for subset in itertools.combinations(bondedTo, 3):
+                    data.impropers.append((atom, subset[0], subset[1], subset[2]))
+        
+        # Identify bonds that should be implemented with constraints
+        
+        if constraints == AllBonds or constraints == HAngles:
+            for bond in data.bonds:
+                bond.isConstrained = True
+        elif constraints == HBonds:
+            for bond in data.bonds:
+                atom1 = data.atoms[bond.atom1]
+                atom2 = data.atoms[bond.atom2]
+                bond.isConstrained = atom1.name.startswith('H') or atom2.name.startswith('H')
+        if rigidWater:
+            for bond in data.bonds:
+                atom1 = data.atoms[bond.atom1]
+                atom2 = data.atoms[bond.atom2]
+                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH':
+                    bond.isConstrained = True
+        
+        # Identify angles that should be implemented with constraints
+        
+        if constraints == HAngles:
+            for angle in data.angles:
+                atom1 = data.atoms[angle[0]]
+                atom2 = data.atoms[angle[1]]
+                atom3 = data.atoms[angle[2]]
+                numH = 0
+                if atom1.name.startswith('H'):
+                    numH += 1
+                if atom3.name.startswith('H'):
+                    numH += 1
+                data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.name.startswith('O')))
+        else:
+            data.isAngleConstrained = len(data.angles)*[False]
+        if rigidWater:
+            for i in range(len(data.angles)):
+                angle = data.angles[i]
+                atom1 = data.atoms[angle[0]]
+                atom2 = data.atoms[angle[1]]
+                atom3 = data.atoms[angle[2]]
+                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH' and atom3.residue.name == 'HOH':
+                    data.isAngleConstrained[i] = True
+
+        # Add forces to the System
+        
+        for force in self._forces:
+            force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
+        return sys
+
+
+def _createResidueSignature(elements):
+    """Create a signature for a residue based on the elements of the atoms it contains."""
+    counts = {}
+    for element in elements:
+        if element in counts:
+            counts[element] += 1
+        else:
+            counts[element] = 1
+    sig = []
+    for c in counts:
+        sig.append((c, counts[c]))
+    sig.sort(key=lambda x: -x[0].mass)
+    return sig
+
+
+def _signatureToString(signature):
+    """Convert the signature returned by _createResidueSignature() to a string."""
+    s = ''
+    for element, count in signature:
+        s += element.symbol+str(count)
+    return s
+
+
+def _matchResidue(res, template, bondedToAtom, atomIndices):
+    """Determine whether a residue matches a template and return a list of corresponding atoms.
+    
+    Parameters:
+     - res (Residue) The residue to check
+     - template (_TemplateData) The template to compare it to
+     - bondedToAtom (list) Enumerates which other atoms each atom is bonded to
+     - atomIndices (map) Maps from atoms to their indices in the System
+    Returns: a list specifying which atom of the template each atom of the residue corresponds to,
+    or None if it does not match the template
+    """
+    atoms = list(res.atoms())
+    residueAtomBonds = []
+    templateAtomBonds = []
+    matches = len(atoms)*[0]
+    hasMatch = len(atoms)*[False]
+    
+    # Translate from global to local atom indices, and record the bonds for each atom.
+    
+    renumberAtoms = {}
+    for i in range(len(atoms)):
+        renumberAtoms[atomIndices[atoms[i]]] = i
+    bondedTo = []
+    externalBonds = []
+    for atom in atoms:
+        bonds = [renumberAtoms[x] for x in bondedToAtom[atomIndices[atom]] if x in renumberAtoms]
+        bondedTo.append(bonds)
+        externalBonds.append(len([x for x in bondedToAtom[atomIndices[atom]] if x not in renumberAtoms]))
+    if _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, 0):
+        return matches
+    return None
+
+
+def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, position):
+    """This is called recursively from inside _matchResidue() to identify matching atoms."""
+    if position == len(atoms):
+        return True
+    elem = atoms[position].element
+    for i in range(len(atoms)):
+        atom = template.atoms[i]
+        if atom.element == elem and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
+            # See if the bonds for this identification are consistent
+            
+            allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
+            if allBondsMatch:
+                # This is a possible match, so trying matching the rest of the residue.
+                
+                matches[position] = i
+                hasMatch[i] = True
+                if _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch, position+1):
+                    return True
+                hasMatch[i] = False
+    return False
+
+
+# The following classes are generators that know how to create Force subclasses and add them to a System that is being
+# created.  Each generator class must define two methods: 1) a static method that takes an etree Element and a ForceField,
+# and returns the corresponding generator object; 2) a createForce() method that constructs the Force object and adds it
+# to the System.  The static method should be added to the parsers map.
+
+
+class HarmonicBondGenerator:
+    """A HarmonicBondGenerator constructs a HarmonicBondForce."""
+    
+    def __init__(self):
+        self.types1 = []
+        self.types2 = []
+        self.length = []
+        self.k = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = HarmonicBondGenerator()
+        ff._forces.append(generator)
+        for bond in element.findall('Bond'):
+            types = ff._findAtomTypes(bond, 2)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.length.append(float(bond.attrib['length']))
+                generator.k.append(float(bond.attrib['k']))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.HarmonicBondForce]
+        if len(existing) == 0:
+            force = mm.HarmonicBondForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        for bond in data.bonds:
+            type1 = data.atomType[data.atoms[bond.atom1]]
+            type2 = data.atomType[data.atoms[bond.atom2]]
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                    bond.length = self.length[i]
+                    if bond.isConstrained:
+                        sys.addConstraint(bond.atom1, bond.atom2, self.length[i])
+                    elif self.k[i] != 0:
+                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
+                    break
+
+parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement
+
+
+class HarmonicAngleGenerator:
+    """A HarmonicAngleGenerator constructs a HarmonicAngleForce."""
+    
+    def __init__(self):
+        self.types1 = []
+        self.types2 = []
+        self.types3 = []
+        self.angle = []
+        self.k = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = HarmonicAngleGenerator()
+        ff._forces.append(generator)
+        for angle in element.findall('Angle'):
+            types = ff._findAtomTypes(angle, 3)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+                generator.angle.append(float(angle.attrib['angle']))
+                generator.k.append(float(angle.attrib['k']))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.HarmonicAngleForce]
+        if len(existing) == 0:
+            force = mm.HarmonicAngleForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
+            type1 = data.atomType[data.atoms[angle[0]]]
+            type2 = data.atomType[data.atoms[angle[1]]]
+            type3 = data.atomType[data.atoms[angle[2]]]
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                    if isConstrained:
+                        # Find the two bonds that make this angle.
+                        
+                        bond1 = None
+                        bond2 = None
+                        for bond in data.atomBonds[angle[1]]:
+                            atom1 = data.bonds[bond].atom1
+                            atom2 = data.bonds[bond].atom2
+                            if atom1 == angle[0] or atom2 == angle[0]:
+                                bond1 = bond
+                            elif atom1 == angle[2] or atom2 == angle[2]:
+                                bond2 = bond
+                        
+                        # Compute the distance between atoms and add a constraint
+                        
+                        if bond1 is not None and bond2 is not None:
+                            l1 = data.bonds[bond1].length
+                            l2 = data.bonds[bond2].length
+                            if l1 is not None and l2 is not None:
+                                length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(self.angle[i]))
+                                sys.addConstraint(angle[0], angle[2], length)
+                    elif self.k[i] != 0:
+                        force.addAngle(angle[0], angle[1], angle[2], self.angle[i], self.k[i])
+                    break
+
+parsers["HarmonicAngleForce"] = HarmonicAngleGenerator.parseElement
+
+
+class PeriodicTorsion:
+    """A PeriodicTorsion records the information for a periodic torsion definition."""
+
+    def __init__(self, types):
+        self.types1 = types[0]
+        self.types2 = types[1]
+        self.types3 = types[2]
+        self.types4 = types[3]
+        self.periodicity = []
+        self.phase = []
+        self.k = []
+
+class PeriodicTorsionGenerator:
+    """A PeriodicTorsionGenerator constructs a PeriodicTorsionForce."""
+    
+    def __init__(self):
+        self.proper = []
+        self.improper = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = PeriodicTorsionGenerator()
+        generator.ff = ff
+        ff._forces.append(generator)
+        for torsion in element.findall('Proper'):
+            torsion = ff._parseTorsion(torsion)
+            if torsion is not None:
+                generator.proper.append(torsion)
+        for torsion in element.findall('Improper'):
+            torsion = ff._parseTorsion(torsion)
+            if torsion is not None:
+                generator.improper.append(torsion)
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.PeriodicTorsionForce]
+        if len(existing) == 0:
+            force = mm.PeriodicTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        wildcard = self.ff._atomClasses['']
+        for torsion in data.propers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            match = None
+            for tordef in self.proper:
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
+                    hasWildcard = (wildcard in (types1, types2, types3, types4))
+                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
+                        match = tordef
+                    if not hasWildcard:
+                        break
+            if match is not None:
+                for i in range(len(match.phase)):
+                    if match.k[i] != 0:
+                        force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.periodicity[i], match.phase[i], match.k[i])
+        for torsion in data.impropers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            done = False
+            for tordef in self.improper:
+                if done:
+                    break
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if type1 in types1:
+                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
+                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
+                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
+                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
+                            # to pick the order.
+                            a1 = torsion[t2[1]]
+                            a2 = torsion[t3[1]]
+                            e1 = data.atoms[a1].element
+                            e2 = data.atoms[a2].element
+                            if e1 == e2 and a1 > a2:
+                                (a1, a2) = (a2, a1)
+                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
+                                (a1, a2) = (a2, a1)
+                            for i in range(len(tordef.phase)):
+                                if tordef.k[i] != 0:
+                                    force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                            done = True
+                            break
+
+parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement
+
+
+class RBTorsion:
+    """An RBTorsion records the information for a Ryckaert-Bellemans torsion definition."""
+
+    def __init__(self, types, c):
+        self.types1 = types[0]
+        self.types2 = types[1]
+        self.types3 = types[2]
+        self.types4 = types[3]
+        self.c = c
+
+class RBTorsionGenerator:
+    """An RBTorsionGenerator constructs an RBTorsionForce."""
+    
+    def __init__(self):
+        self.proper = []
+        self.improper = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = RBTorsionGenerator()
+        generator.ff = ff
+        ff._forces.append(generator)
+        for torsion in element.findall('Proper'):
+            types = ff._findAtomTypes(torsion, 4)
+            if types is not None:
+                generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
+        for torsion in element.findall('Improper'):
+            types = ff._findAtomTypes(torsion, 4)
+            if types is not None:
+                generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.RBTorsionForce]
+        if len(existing) == 0:
+            force = mm.RBTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        wildcard = self.ff._atomClasses['']
+        for torsion in data.propers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            match = None
+            for tordef in self.proper:
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
+                    hasWildcard = (wildcard in (types1, types2, types3, types4))
+                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
+                        match = tordef
+                    if not hasWildcard:
+                        break
+            if match is not None:
+                force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.c[0], match.c[1], match.c[2], match.c[3], match.c[4], match.c[5])
+        for torsion in data.impropers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            done = False
+            for tordef in self.improper:
+                if done:
+                    break
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if type1 in types1:
+                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
+                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
+                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
+                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
+                            # to pick the order.
+                            a1 = torsion[t2[1]]
+                            a2 = torsion[t3[1]]
+                            e1 = data.atoms[a1].element
+                            e2 = data.atoms[a2].element
+                            if e1 == e2 and a1 > a2:
+                                (a1, a2) = (a2, a1)
+                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
+                                (a1, a2) = (a2, a1)
+                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
+                            done = True
+                            break
+
+parsers["RBTorsionForce"] = RBTorsionGenerator.parseElement
+
+
+class CMAPTorsion:
+    """A CMAPTorsion records the information for a CMAP torsion definition."""
+
+    def __init__(self, types, map):
+        self.types1 = types[0]
+        self.types2 = types[1]
+        self.types3 = types[2]
+        self.types4 = types[3]
+        self.types5 = types[4]
+        self.map = map
+
+class CMAPTorsionGenerator:
+    """A CMAPTorsionGenerator constructs a CMAPTorsionForce."""
+    
+    def __init__(self):
+        self.torsions = []
+        self.maps = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CMAPTorsionGenerator()
+        generator.ff = ff
+        ff._forces.append(generator)
+        for map in element.findall('Map'):
+            values = [float(x) for x in map.text.split()]
+            size = sqrt(len(values))
+            if size*size != len(values):
+                raise ValueError('CMAP must have the same number of elements along each dimension')
+            generator.maps.append(values)
+        for torsion in element.findall('Torsion'):
+            types = ff._findAtomTypes(torsion, 5)
+            if types is not None:
+                generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map'])))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.CMAPTorsionForce]
+        if len(existing) == 0:
+            force = mm.CMAPTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        for map in self.maps:
+            force.addMap(int(sqrt(len(map))), map)
+        
+        # Find all chains of length 5
+        
+        uniqueTorsions = set()
+        for torsion in data.propers:
+            for bond in (data.bonds[x] for x in data.atomBonds[torsion[0]]):
+                if bond.atom1 == torsion[0]:
+                    atom = bond.atom2
+                else:
+                    atom = bond.atom1
+                if atom != torsion[1]:
+                    uniqueTorsions.add((atom, torsion[0], torsion[1], torsion[2], torsion[3]))
+            for bond in (data.bonds[x] for x in data.atomBonds[torsion[3]]):
+                if bond.atom1 == torsion[3]:
+                    atom = bond.atom2
+                else:
+                    atom = bond.atom1
+                if atom != torsion[2]:
+                    uniqueTorsions.add((torsion[0], torsion[1], torsion[2], torsion[3], atom))
+        torsions = sorted(list(uniqueTorsions))
+        wildcard = self.ff._atomClasses['']
+        for torsion in torsions:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            type5 = data.atomType[data.atoms[torsion[4]]]
+            match = None
+            for tordef in self.torsions:
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                types5 = tordef.types5
+                if (type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4 and type5 in types5) or (type1 in types5 and type2 in types4 and type3 in types3 and type4 in types2 and type5 in types1):
+                    hasWildcard = (wildcard in (types1, types2, types3, types4, types5))
+                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
+                        match = tordef
+                    if not hasWildcard:
+                        break
+            if match is not None:
+                force.addTorsion(match.map, torsion[0], torsion[1], torsion[2], torsion[3], torsion[1], torsion[2], torsion[3], torsion[4])
+
+parsers["CMAPTorsionForce"] = CMAPTorsionGenerator.parseElement
+
+
+class NonbondedGenerator:
+    """A NonbondedGenerator constructs a NonbondedForce."""
+    
+    def __init__(self, coulomb14scale, lj14scale):
+        self.coulomb14scale = coulomb14scale
+        self.lj14scale = lj14scale
+        self.typeMap = {}
+
+    @staticmethod
+    def parseElement(element, ff):
+        existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)]
+        if len(existing) == 0:
+            generator = NonbondedGenerator(float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale']))
+            ff._forces.append(generator)
+        else:
+            # Multiple <NonbondedForce> tags were found, probably in different files.  Simply add more types to the existing one.
+            generator = existing[0]
+            if generator.coulomb14scale != float(element.attrib['coulomb14scale']) or generator.lj14scale != float(element.attrib['lj14scale']):
+                raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales') 
+        for atom in element.findall('Atom'):
+            types = ff._findAtomTypes(atom, 1)
+            if types is not None:
+                values = (float(atom.attrib['charge']), float(atom.attrib['sigma']), float(atom.attrib['epsilon']))
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
+                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
+                     Ewald:mm.NonbondedForce.Ewald,
+                     PME:mm.NonbondedForce.PME}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for NonbondedForce') 
+        force = mm.NonbondedForce()
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                values = self.typeMap[t]
+                force.addParticle(values[0], values[1], values[2])
+            else:
+                raise ValueError('No nonbonded parameters defined for atom type '+t) 
+        bondIndices = []
+        for bond in data.bonds:
+            bondIndices.append((bond.atom1, bond.atom2))
+        force.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        if 'ewaldErrorTolerance' in args:
+            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
+        sys.addForce(force)
+
+parsers["NonbondedForce"] = NonbondedGenerator.parseElement
+
+
+class GBSAOBCGenerator:
+    """A GBSAOBCGenerator constructs a GBSAOBCForce."""
+    
+    def __init__(self):
+        self.typeMap = {}
+
+    @staticmethod
+    def parseElement(element, ff):
+        generator = GBSAOBCGenerator()
+        ff._forces.append(generator)
+        for atom in element.findall('Atom'):
+            types = ff._findAtomTypes(atom, 1)
+            if types is not None:
+                values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['scale']))
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
+                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for GBSAOBCForce') 
+        force = mm.GBSAOBCForce()
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                values = self.typeMap[t]
+                force.addParticle(values[0], values[1], values[2])
+            else:
+                raise ValueError('No GBSAOBC parameters defined for atom type '+t) 
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        sys.addForce(force)
+
+parsers["GBSAOBCForce"] = GBSAOBCGenerator.parseElement
+
+
+class GBVIGenerator:
+
+    """A GBVIGenerator constructs a GBVIForce."""
+    
+    def __init__(self,ff):
+
+        self.ff                                             = ff
+        self.fixedParameters                                = {}
+
+        self.fixedParameters['soluteDielectric']            = 1.0
+        self.fixedParameters['solventDielectric']           = 78.3
+        self.fixedParameters['scalingMethod']               = 1
+        self.fixedParameters['quinticUpperBornRadiusLimit'] = 5.0
+        self.fixedParameters['quinticLowerLimitFactor']     = 0.8
+
+        self.typeMap                                        = {}
+
+    @staticmethod
+    def parseElement(element, ff):
+        generator    = GBVIGenerator(ff)
+        for key in generator.fixedParameters.iterkeys():
+            if( key in element.attrib ):
+                generator.fixedParameters[key] = float( element.attrib[key] )
+
+        ff._forces.append(generator)
+        for atom in element.findall('Atom'):
+            types = ff._findAtomTypes(atom, 1)
+            if types is not None:
+                values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+
+        methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
+                     CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}
+
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for GB/VI Force') 
+
+        # add particles
+
+        force = mm.GBVIForce()
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                values = self.typeMap[t]
+                force.addParticle(values[0], values[1], values[2])
+            else:
+                raise ValueError('No GB/VI parameters defined for atom type '+t) 
+
+        # get HarmonicBond generator -- exit if not found
+
+        hbGenerator = 0
+        for generator in self.ff._forces:
+            if( generator.__class__.__name__ == 'HarmonicBondGenerator' ): 
+               hbGenerator = generator
+               break
+
+        if( hbGenerator == 0 ):
+            raise ValueError('HarmonicBondGenerator not found.') 
+
+        # add bonds
+
+        for bond in data.bonds:
+            type1 = data.atomType[data.atoms[bond.atom1]]
+            type2 = data.atomType[data.atoms[bond.atom2]]
+            hit   = 0
+            for i in range(len(hbGenerator.types1)):
+                types1 = hbGenerator.types1[i]
+                types2 = hbGenerator.types2[i]
+                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                    #bond.length = hbGenerator.length[i]
+                    force.addBond(bond.atom1, bond.atom2, hbGenerator.length[i])
+
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        force.setSolventDielectric( self.fixedParameters['solventDielectric'] )
+        force.setSoluteDielectric( self.fixedParameters['soluteDielectric'] )
+        force.setBornRadiusScalingMethod( self.fixedParameters['scalingMethod'] )
+        force.setQuinticLowerLimitFactor( self.fixedParameters['quinticLowerLimitFactor'] )
+        force.setQuinticUpperBornRadiusLimit( self.fixedParameters['quinticUpperBornRadiusLimit'] )
+        
+        sys.addForce(force)
+
+parsers["GBVIForce"] = GBVIGenerator.parseElement
+
+class CustomBondGenerator:
+    """A CustomBondGenerator constructs a CustomBondForce."""
+    
+    def __init__(self):
+        self.types1 = []
+        self.types2 = []
+        self.globalParams = {}
+        self.perBondParams = []
+        self.paramValues = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomBondGenerator()
+        ff._forces.append(generator)
+        generator.energy = element.attrib['energy']
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerBondParameter'):
+            generator.perBondParams.append(param.attrib['name'])
+        for bond in element.findall('Bond'):
+            types = ff._findAtomTypes(bond, 2)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.paramValues.append([float(bond.attrib[param]) for param in generator.perBondParams])
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        force = mm.CustomBondForce(self.energy)
+        sys.addForce(force)
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perBondParams:
+            force.addPerBondParameter(param)
+        for bond in data.bonds:
+            type1 = data.atomType[data.atoms[bond.atom1]]
+            type2 = data.atomType[data.atoms[bond.atom2]]
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                    force.addBond(bond.atom1, bond.atom2, self.paramValues[i])
+                    break
+
+parsers["CustomBondForce"] = CustomBondGenerator.parseElement
+
+
+class CustomAngleGenerator:
+    """A CustomAngleGenerator constructs a CustomAngleForce."""
+    
+    def __init__(self):
+        self.types1 = []
+        self.types2 = []
+        self.types3 = []
+        self.globalParams = {}
+        self.perAngleParams = []
+        self.paramValues = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomAngleGenerator()
+        ff._forces.append(generator)
+        generator.energy = element.attrib['energy']
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerAngleParameter'):
+            generator.perAngleParams.append(param.attrib['name'])
+        for angle in element.findall('Angle'):
+            types = ff._findAtomTypes(angle, 3)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+                generator.paramValues.append([float(angle.attrib[param]) for param in generator.perAngleParams])
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        force = mm.CustomAngleForce(self.energy)
+        sys.addForce(force)
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perAngleParams:
+            force.addPerAngleParameter(param)
+        for angle in data.angles:
+            type1 = data.atomType[data.atoms[angle[0]]]
+            type2 = data.atomType[data.atoms[angle[1]]]
+            type3 = data.atomType[data.atoms[angle[2]]]
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                    force.addAngle(angle[0], angle[1], angle[2], self.paramValues[i])
+                    break
+
+parsers["CustomAngleForce"] = CustomAngleGenerator.parseElement
+
+
+class CustomTorsion:
+    """A CustomTorsion records the information for a custom torsion definition."""
+
+    def __init__(self, types, paramValues):
+        self.types1 = types[0]
+        self.types2 = types[1]
+        self.types3 = types[2]
+        self.types4 = types[3]
+        self.paramValues = paramValues
+
+class CustomTorsionGenerator:
+    """A CustomTorsionGenerator constructs a CustomTorsionForce."""
+    
+    def __init__(self):
+        self.proper = []
+        self.improper = []
+        self.globalParams = {}
+        self.perTorsionParams = []
+    
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomTorsionGenerator()
+        generator.ff = ff
+        ff._forces.append(generator)
+        generator.energy = element.attrib['energy']
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerTorsionParameter'):
+            generator.perTorsionParams.append(param.attrib['name'])
+        for torsion in element.findall('Proper'):
+            types = ff._findAtomTypes(torsion, 4)
+            if types is not None:
+                generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
+        for torsion in element.findall('Improper'):
+            types = ff._findAtomTypes(torsion, 4)
+            if types is not None:
+                generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        force = mm.CustomTorsionForce(self.energy)
+        sys.addForce(force)
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perTorsionParams:
+            force.addPerTorsionParameter(param)
+        wildcard = self.ff._atomClasses['']
+        for torsion in data.propers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            match = None
+            for tordef in self.proper:
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if (type2 in types2 and type3 in types3 and type4 in types4 and type1 in types1) or (type2 in types3 and type3 in types2 and type4 in types1 and type1 in types4):
+                    hasWildcard = (wildcard in (types1, types2, types3, types4))
+                    if match is None or not hasWildcard: # Prefer specific definitions over ones with wildcards
+                        match = tordef
+                    if not hasWildcard:
+                        break
+            if match is not None:
+                force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.paramValues)
+        for torsion in data.impropers:
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+            done = False
+            for tordef in self.improper:
+                if done:
+                    break
+                types1 = tordef.types1
+                types2 = tordef.types2
+                types3 = tordef.types3
+                types4 = tordef.types4
+                if type1 in types1:
+                    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
+                        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
+                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
+                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
+                            # to pick the order.
+                            a1 = torsion[t2[1]]
+                            a2 = torsion[t3[1]]
+                            e1 = data.atoms[a1].element
+                            e2 = data.atoms[a2].element
+                            if e1 == e2 and a1 > a2:
+                                (a1, a2) = (a2, a1)
+                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
+                                (a1, a2) = (a2, a1)
+                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.paramValues)
+                            done = True
+                            break
+
+parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement
+
+
+class CustomGBGenerator:
+    """A CustomGBGenerator constructs a CustomGBForce."""
+    
+    def __init__(self):
+        self.typeMap = {}
+        self.globalParams = {}
+        self.perParticleParams = []
+        self.paramValues = []
+        self.computedValues = []
+        self.energyTerms = []
+        self.functions = []
+
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomGBGenerator()
+        ff._forces.append(generator)
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerParticleParameter'):
+            generator.perParticleParams.append(param.attrib['name'])
+        for atom in element.findall('Atom'):
+            types = ff._findAtomTypes(atom, 1)
+            if types is not None:
+                values = [float(atom.attrib[param]) for param in generator.perParticleParams]
+                for t in types[0]:
+                    generator.typeMap[t] = values
+        computationMap = {"SingleParticle" : mm.CustomGBForce.SingleParticle,
+                          "ParticlePair" : mm.CustomGBForce.ParticlePair,
+                          "ParticlePairNoExclusions" : mm.CustomGBForce.ParticlePairNoExclusions}
+        for value in element.findall('ComputedValue'):
+            generator.computedValues.append((value.attrib['name'], value.text, computationMap[value.attrib['type']]))
+        for term in element.findall('EnergyTerm'):
+            generator.energyTerms.append((term.text, computationMap[term.attrib['type']]))
+        for function in element.findall("Function"):
+            values = [float(x) for x in function.text.split()]
+            generator.functions.append((function.attrib['name'], values, float(function.attrib['min']), float(function.attrib['max'])))
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        methodMap = {NoCutoff:mm.CustomGBForce.NoCutoff,
+                     CutoffNonPeriodic:mm.CustomGBForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.CustomGBForce.CutoffPeriodic}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for CustomGBForce') 
+        force = mm.CustomGBForce()
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perParticleParams:
+            force.addPerParticleParameter(param)
+        for value in self.computedValues:
+            force.addComputedValue(value[0], value[1], value[2])
+        for term in self.energyTerms:
+            force.addEnergyTerm(term[0], term[1])
+        for function in self.functions:
+            force.addFunction(function[0], function[1], function[2], function[3])
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                values = self.typeMap[t]
+                force.addParticle(self.typeMap[t])
+            else:
+                raise ValueError('No CustomGB parameters defined for atom type '+t) 
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        sys.addForce(force)
+
+parsers["CustomGBForce"] = CustomGBGenerator.parseElement
+
+def getAtomPrint( data, atomIndex ):
+
+    if( atomIndex < len( data.atoms ) ):
+        atom          = data.atoms[atomIndex]
+        returnString  = "%4s %4s %5d" % (atom.name, atom.residue.name, atom.residue.index )
+    else:
+        returnString  = "NA"
+
+    return returnString
+
+#=============================================================================================
+
+def countConstraint( data ):
+
+    bondCount  = 0
+    angleCount = 0
+    for bond in data.bonds:
+        if bond.isConstrained:
+            bondCount += 1
+
+    angleCount = 0
+    for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
+        if( isConstrained ):
+            angleCount += 1
+ 
+    print "Constraints bond=%d angle=%d  total=%d" % (bondCount, angleCount, (bondCount+angleCount) )
+
+class AmoebaHarmonicBondGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaHarmonicBondGenerator constructs a AmoebaHarmonicBondForce."""
+
+    #=============================================================================================
+    
+    def __init__(self, cubic, quartic):
+
+        self.cubic         = cubic
+        self.quartic       = quartic
+        self.types1        = []
+        self.types2        = []
+        self.length        = []
+        self.k             = []
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        # <AmoebaHarmonicBondForce bond-cubic="-25.5" bond-quartic="379.3125">
+        # <Bond class1="1" class2="2" length="0.1437" k="156900.0"/>
+    
+        generator = AmoebaHarmonicBondGenerator( float(element.attrib['bond-cubic']), float(element.attrib['bond-quartic']) )
+        forceField._forces.append(generator)
+        for bond in element.findall('Bond'):
+            types = forceField._findAtomTypes(bond, 2)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.length.append(float(bond.attrib['length']))
+                generator.k.append(float(bond.attrib['k']))
+            else:
+                outputString = "AmoebaHarmonicBondGenerator: error getting types: %s %s" % (
+                                    bond.attrib['class1'],
+                                    bond.attrib['class2'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        verbose            = 0
+        if( self.hasBeenCalled ):
+             return
+
+        if( verbose ): 
+            countConstraint( data )
+
+        self.hasBeenCalled = 1
+
+        existing           = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing           = [f for f in existing if type(f) == mm.AmoebaHarmonicBondForce]
+        if len(existing) == 0:
+            force = mm.AmoebaHarmonicBondForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        force.setAmoebaGlobalHarmonicBondCubic( self.cubic )
+        force.setAmoebaGlobalHarmonicBondQuartic( self.quartic )
+
+        if( verbose ):
+            print "In AmoebaHarmonicBondGenerator bonds=%d " % (len(data.bonds))
+
+        count = 0
+        for bond in data.bonds:
+            type1 = data.atomType[data.atoms[bond.atom1]]
+            type2 = data.atomType[data.atoms[bond.atom2]]
+            hit   = 0
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                    bond.length = self.length[i]
+                    hit  = 1
+                    if bond.isConstrained:
+                        sys.addConstraint(bond.atom1, bond.atom2, self.length[i])
+                        if( verbose ):
+                            atomS1 = getAtomPrint( data, bond.atom1 )
+                            atomS2 = getAtomPrint( data, bond.atom2 )
+                            print "AmoebaHarmonicBondGenerator %5d %5d %5d [%s %s] [%5s %5s] %15.6f %15.6f Constraint" % (count, bond.atom1, bond.atom2, atomS1, atomS2, type1, type2, self.length[i], self.k[i])
+                    elif self.k[i] != 0:
+                        if( verbose ):
+                            atomS1 = getAtomPrint( data, bond.atom1 )
+                            atomS2 = getAtomPrint( data, bond.atom2 )
+                            print "AmoebaHarmonicBondGenerator %5d %5d %5d [%s %s] [%5s %5s] %15.6f %15.6f" % (count, bond.atom1, bond.atom2, atomS1, atomS2, type1, type2, self.length[i], self.k[i])
+                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
+                    break
+            if( hit == 0 ): 
+                print "AmoebaHarmonicBondGenerator missing: %5d types=[%5s %5s] atoms=[%6d %6d] " % (count, type1, type2, bond.atom1, bond.atom2)
+
+            count += 1
+
+parsers["AmoebaHarmonicBondForce"] = AmoebaHarmonicBondGenerator.parseElement
+
+#=============================================================================================
+# Add angle constraint
+#=============================================================================================
+    
+def addAngleConstraint( angle, idealAngle, data, sys ):
+
+    # Find the two bonds that make this angle.
+                    
+    bond1      = None
+    bond2      = None
+    for bond in data.atomBonds[angle[1]]:
+        atom1 = data.bonds[bond].atom1
+        atom2 = data.bonds[bond].atom2
+        if atom1 == angle[0] or atom2 == angle[0]:
+            bond1 = bond
+        elif atom1 == angle[2] or atom2 == angle[2]:
+            bond2 = bond
+                    
+        # Compute the distance between atoms and add a constraint
+                    
+        if bond1 is not None and bond2 is not None:
+            l1 = data.bonds[bond1].length
+            l2 = data.bonds[bond2].length
+            if l1 is not None and l2 is not None:
+                length     = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(idealAngle))
+                sys.addConstraint(angle[0], angle[2], length)
+                return
+
+
+#=============================================================================================
+class AmoebaHarmonicAngleGenerator:
+
+    #=============================================================================================
+    """An AmoebaHarmonicAngleGenerator constructs a AmoebaHarmonicAngleForce."""
+    #=============================================================================================
+    
+    def __init__(self, forceField, cubic, quartic, pentic, sextic):
+
+        self.forceField    = forceField
+        self.cubic         = cubic
+        self.quartic       = quartic
+        self.pentic        = pentic
+        self.sextic        = sextic
+
+        self.types1        = []
+        self.types2        = []
+        self.types3        = []
+
+        self.angle         = []
+        self.k             = []
+
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        # <AmoebaHarmonicAngleForce angle-cubic="-0.014" angle-quartic="5.6e-05" angle-pentic="-7e-07" angle-sextic="2.2e-08">
+        #   <Angle class1="2" class2="1" class3="3" k="0.0637259642196" angle1="122.00"  />
+
+        generator = AmoebaHarmonicAngleGenerator( forceField, float(element.attrib['angle-cubic']), float(element.attrib['angle-quartic']),  float(element.attrib['angle-pentic']), float(element.attrib['angle-sextic'] ) )
+        forceField._forces.append(generator)
+        for angle in element.findall('Angle'):
+            types = forceField._findAtomTypes(angle, 3)
+            if types is not None:
+
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+
+                angleList = []
+                angleList.append( float(angle.attrib['angle1']) )
+
+                try:
+                    angleList.append( float(angle.attrib['angle2']) )
+                    try:
+                        angleList.append( float(angle.attrib['angle3']) )
+                    except:
+                        pass
+                except:
+                    pass
+                generator.angle.append(angleList)
+                generator.k.append(float(angle.attrib['k']))
+            else:
+                outputString = "AmoebaHarmonicAngleGenerator: error getting types: %s %s %s" % (
+                                    angle.attrib['class1'],
+                                    angle.attrib['class2'],
+                                    angle.attrib['class3'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+    # createForce is bypassed here since the AmoebaOutOfPlaneBendForce generator must first execute
+    # and partition angles into in-plane and non-in-plane angles
+    #=============================================================================================
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+        pass
+
+    #=============================================================================================
+    # createForcePostOpBendAngle is called by AmoebaOutOfPlaneBendForce with the list of
+    # non-in-plane angles
+    #=============================================================================================
+    
+    def createForcePostOpBendAngle( self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args ):
+
+        verbose             = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled += 1
+
+        # get force
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaHarmonicAngleForce]
+
+        if len(existing) == 0:
+            force = mm.AmoebaHarmonicAngleForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        # scalars
+
+        force.setAmoebaGlobalHarmonicAngleCubic( self.cubic )
+        force.setAmoebaGlobalHarmonicAngleQuartic( self.quartic )
+        force.setAmoebaGlobalHarmonicAnglePentic( self.pentic )
+        force.setAmoebaGlobalHarmonicAngleSextic( self.sextic )
+
+        if( verbose ):
+            print "In AmoebaHarmonicAngleGenerator angles=%d " % (len(data.angles))
+
+        count = 0
+        for angleDict in angleList:
+            angle          = angleDict['angle']
+            isConstrained  = angleDict['isConstrained']
+
+            type1          = data.atomType[data.atoms[angle[0]]]
+            type2          = data.atomType[data.atoms[angle[1]]]
+            type3          = data.atomType[data.atoms[angle[2]]]
+            hit            = 0
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                    hit      = 1
+                    if isConstrained and self.k[i] != 0.0:
+
+                        angleDict['idealAngle'] = self.angle[i][0]
+                        addAngleConstraint( angle, self.angle[i][0], data, sys )
+
+                    elif self.k[i] != 0:
+#                       print "AmoebaHarmonicAngleGenerator %5d %5d %5d %5d [%5s %5s %5s] %15.6f %15.6f" % (count, angle[0], angle[1], angle[2], type1, type2, type3, self.angle[i], self.k[i])
+                        lenAngle = len( self.angle[i] )
+                        if( lenAngle > 1 ):
+                            # get k-index by counting number of non-angle hydrogens on the central atom
+                            # based on kangle.f
+                            numberOfHydrogens = 0
+                            for bond in data.atomBonds[angle[1]]:
+                                atom1 = data.bonds[bond].atom1
+                                atom2 = data.bonds[bond].atom2
+                                if( atom1 == angle[1] and atom2 != angle[0] and atom2 != angle[2] and (sys.getParticleMass( atom2 )/unit.dalton) < 1.90 ):
+                                    numberOfHydrogens += 1
+                                if( atom2 == angle[1] and atom1 != angle[0] and atom1 != angle[2] and (sys.getParticleMass( atom1 )/unit.dalton) < 1.90 ):
+                                    numberOfHydrogens += 1
+                            if( numberOfHydrogens < lenAngle ):
+                                angleValue =  self.angle[i][numberOfHydrogens]
+                            else:
+                                print "Error: AmoebaHarmonicAngleGenerator angle index=%d is out of range: [0, %5d] " % (numberOfHydrogens, lenAngle)
+                                sys.exit(-1)
+                        else:
+                            angleValue =  self.angle[i][0]
+               
+                        angleDict['idealAngle'] = angleValue
+                        force.addAngle(angle[0], angle[1], angle[2], angleValue, self.k[i] )
+                    break
+            if( hit == 0 ): 
+                print "AmoebaHarmonicAngleGenerator missing: %5d %s %s %s " % (count, type1, type2, type3)
+
+            count += 1
+
+    #=============================================================================================
+    # createForcePostOpBendInPlaneAngle is called by AmoebaOutOfPlaneBendForce with the list of
+    # in-plane angles
+    #=============================================================================================
+    
+    def createForcePostOpBendInPlaneAngle( self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args ):
+
+        verbose             = 0
+        self.hasBeenCalled += 1
+
+        # get force
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaHarmonicInPlaneAngleForce]
+
+        if len(existing) == 0:
+            force = mm.AmoebaHarmonicInPlaneAngleForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        # scalars
+
+        force.setAmoebaGlobalHarmonicInPlaneAngleCubic( self.cubic )
+        force.setAmoebaGlobalHarmonicInPlaneAngleQuartic( self.quartic )
+        force.setAmoebaGlobalHarmonicInPlaneAnglePentic( self.pentic )
+        force.setAmoebaGlobalHarmonicInPlaneAngleSextic( self.sextic )
+
+        if( verbose ):
+            print "In AmoebaHarmonicAngleGenerator angles=%d " % (len(data.angles))
+
+        count = 0
+        for angleDict in angleList:
+ 
+            angle          = angleDict['angle']
+            isConstrained  = angleDict['isConstrained']
+
+            type1          = data.atomType[data.atoms[angle[0]]]
+            type2          = data.atomType[data.atoms[angle[1]]]
+            type3          = data.atomType[data.atoms[angle[2]]]
+
+            hit   = 0
+            for i in range(len(self.types1)):
+
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+
+                if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                    if( verbose ):
+                        print "AmoebaHarmonicInPlaneAngleGenerator %5d %5d %5d %5d len=%d [%5s %5s %5s] %15.6f %15.6f" % (count, angle[0], angle[1], angle[2], len(angle), type1, type2, type3, self.angle[i][0], self.k[i])
+                    hit                     = 1
+                    angleDict['idealAngle'] = self.angle[i][0]
+                    if( isConstrained and self.k[i] != 0.0 ):
+                        addAngleConstraint( angle, self.angle[i][0], data, sys )
+                    else:
+                        force.addAngle(angle[0], angle[1], angle[2], angle[3], self.angle[i][0], self.k[i])
+                    break
+
+            if( hit == 0 ): 
+                print "AmoebaHarmonicInPlaneAngleGenerator missing: %5d %s %s " % (count, type1, type2, type3)
+
+            count += 1
+
+parsers["AmoebaHarmonicAngleForce"] = AmoebaHarmonicAngleGenerator.parseElement
+
+#=============================================================================================
+# Generator for the AmoebaOutOfPlaneBend covalent force; also calls methods in the
+# AmoebaHarmonicAngleGenerator to generate the AmoebaHarmonicAngleForce and
+# AmoebaHarmonicInPlaneAngleForce
+#=============================================================================================
+
+class AmoebaOutOfPlaneBendGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaOutOfPlaneBendGenerator constructs a AmoebaOutOfPlaneBendForce."""
+    
+    #=============================================================================================
+
+    def __init__(self, forceField, type, cubic, quartic, pentic, sextic):
+
+        self.forceField        = forceField
+        self.type              = type
+        self.cubic             = cubic
+        self.quartic           = quartic
+        self.pentic            = pentic
+        self.sextic            = sextic
+
+        self.types1            = []
+        self.types2            = []
+        self.types3            = []
+        self.types4            = []
+
+        self.ks                = []
+        self.hasBeenCalled     = 0
+
+    #=============================================================================================
+    # Local version of findAtomTypes needed since class indices are 0 (i.e., not recognized)
+    # for types3 and 4
+    #=============================================================================================
+    
+    def findAtomTypes(self, forceField, node, num):
+        """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
+        types = []
+        attrib = node.attrib
+        for i in range(num):
+            if num == 1:
+                suffix = ''
+            else:
+                suffix = str(i+1)
+            classAttrib = 'class'+suffix
+            if classAttrib in attrib:
+                if attrib[classAttrib] in forceField._atomClasses:
+                    types.append(forceField._atomClasses[attrib[classAttrib]])
+                else:
+                    types.append(set())
+        return types
+
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaOutOfPlaneBendForce type="ALLINGER" opbend-cubic="-0.014" opbend-quartic="5.6e-05" opbend-pentic="-7e-07" opbend-sextic="2.2e-08">
+        #   <Angle class1="2" class2="1" class3="0" class4="0" k="0.0531474541591"/>
+        #   <Angle class1="3" class2="1" class3="0" class4="0" k="0.0898536095496"/>
+         
+        # get global scalar parameters
+
+        generator = AmoebaOutOfPlaneBendGenerator( forceField, element.attrib['type'],
+                                                   float(element.attrib['opbend-cubic']),
+                                                   float(element.attrib['opbend-quartic']),
+                                                   float(element.attrib['opbend-pentic']),
+                                                   float(element.attrib['opbend-sextic'] ) )
+
+        forceField._forces.append(generator)
+
+        for angle in element.findall('Angle'):
+            types = generator.findAtomTypes( forceField, angle, 4)
+            if types is not None:
+
+                generator.types1.append( types[0] )
+                generator.types2.append( types[1] )
+                generator.types3.append( types[2] )
+                generator.types4.append( types[3] )
+
+                generator.ks.append( float(angle.attrib['k'] ) )
+
+            else:
+                outputString = "AmoebaOutOfPlaneBendGenerator: error getting types: %s %s %s %s." % (
+                                    angle.attrib['class1'],
+                                    angle.attrib['class2'],
+                                    angle.attrib['class3'],
+                                    angle.attrib['class4'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+    # Get middle atom in a angle
+    # return index of middle atom or -1 if no middle is found
+    # This method appears not to be needed since the angle[1] entry appears to always
+    # be the middle atom. However, was unsure if this is guaranteed
+    #=============================================================================================
+    
+    def getMiddleAtom( self, angle, data ):
+
+        # find atom shared by both bonds making up the angle
+
+        middleAtom       = -1
+        for atomIndex in angle: 
+            isMiddle     = 0
+            for bond in data.atomBonds[atomIndex]:
+                atom1 = data.bonds[bond].atom1
+                atom2 = data.bonds[bond].atom2
+                if( atom1 != atomIndex ):
+                    partner = atom1
+                else:
+                    partner = atom2
+                if( partner == angle[0] or partner == angle[1] or partner == angle[2] ): 
+                    isMiddle += 1
+
+            if( isMiddle == 2 ):
+                return atomIndex
+        return -1
+
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        verbose            = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled = 1
+
+        # get force
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaOutOfPlaneBendForce]
+        if len(existing) == 0:
+            force = mm.AmoebaOutOfPlaneBendForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        # set scalars
+
+        force.setAmoebaGlobalOutOfPlaneBendCubic(   self.cubic )
+        force.setAmoebaGlobalOutOfPlaneBendQuartic( self.quartic )
+        force.setAmoebaGlobalOutOfPlaneBendPentic(  self.pentic )
+        force.setAmoebaGlobalOutOfPlaneBendSextic(  self.sextic )
+
+        count                       = 0
+        opBondCount                 = 0
+
+        # this hash is used to insure the out-of-plane-bend bonds
+        # are only added once
+
+        skipAtoms                   = dict()
+
+        # these lists are used in the partitioning of the angles into
+        # angle and inPlane angles
+
+        inPlaneAngles               = []
+        nonInPlaneAngles            = []
+        nonInPlaneAnglesConstrained = []
+        idealAngles                 = []*len(data.angles)
+
+        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
+
+            middleAtom          = self.getMiddleAtom( angle, data )
+            if( middleAtom > -1 ):
+                middleType      = data.atomType[data.atoms[middleAtom]]
+                middleCovalency = len(data.atomBonds[middleAtom])
+            else:
+                middleType      = -1
+                middleCovalency = -1
+
+            # if middle atom has covalency of 3 and 
+            # the types of the middle atom and the partner atom (atom bonded to
+            # middle atom, but not in angle) match types2 and types2, then
+            # three out-of-plane bend angles are generated. Three in-plane angle 
+            # are also generated. If the conditions are not satisfied the angle is generic angle (not a in-plane angle)
+
+            if( middleAtom > -1 and middleCovalency == 3 and middleAtom not in skipAtoms ):
+
+                partners        = []
+                partnerSet      = set()
+                partnerTypes    = []
+                partnerK        = []
+
+                for bond in data.atomBonds[middleAtom]:
+                    atom1 = data.bonds[bond].atom1
+                    atom2 = data.bonds[bond].atom2
+                    if( atom1 != middleAtom ):
+                        partner = atom1
+                    else:
+                        partner = atom2
+
+                    partnerType = data.atomType[data.atoms[partner]]
+                    for i in range(len(self.types1)):
+                        types1 = self.types1[i]
+                        types2 = self.types2[i]
+                        if( middleType in types2 and partnerType in types1 ):
+                            partners.append( partner )
+                            partnerSet.add( partner )
+                            partnerTypes.append( partnerType )
+                            partnerK.append( self.ks[i] )
+             
+                if( len(partners) == 3 ):
+
+                    opBondCount += 3
+                    if( verbose ):
+                        print "%5d Opbend: %d type=%s cov=%d " % (opBondCount, middleAtom, middleType, middleCovalency)
+
+                    force.addOutOfPlaneBend(partners[0], middleAtom, partners[1], partners[2], partnerK[2] )
+                    force.addOutOfPlaneBend(partners[0], middleAtom, partners[2], partners[1], partnerK[1] )
+                    force.addOutOfPlaneBend(partners[1], middleAtom, partners[2], partners[0], partnerK[0] )
+
+                    skipAtoms[middleAtom] = set()
+                    skipAtoms[middleAtom].add( partners[0] )
+                    skipAtoms[middleAtom].add( partners[1] )
+                    skipAtoms[middleAtom].add( partners[2] )
+
+                    angleDict                    = {}
+                    angleList                    = []
+                    angleList.append( angle[0] )
+                    angleList.append( angle[1] )
+                    angleList.append( angle[2] )
+                    angleDict['angle'] = angleList
+
+                    angleDict['isConstrained']    = 0
+
+                    angleSet  = set()
+                    angleSet.add( angle[0] )
+                    angleSet.add( angle[1] )
+                    angleSet.add( angle[2] )
+
+                    for atomIndex in partnerSet:
+                        if( atomIndex not in angleSet ):
+                            angleList.append( atomIndex )
+
+                    if( verbose ):
+                        print "%5d Opbend: middle=%d inPlane1 %d %s" % (opBondCount, middleAtom, len(angleList), str(angleList) )
+
+                    inPlaneAngles.append( angleDict )
+
+                else:
+                    if( verbose ):
+                        print "%5d Opbend: %d type=%s cov=%d  xxx" % (opBondCount, middleAtom, middleType, middleCovalency )
+                        print "%5d Opbend: middle=%d inPlane1 %d %s" % (opBondCount, middleAtom, len(angleList), str(angleList) )
+                    angleDict                     = {}
+                    angleDict['angle']            = angle
+                    angleDict['isConstrained']    = isConstrained
+                    nonInPlaneAngles.append( angleDict )
+            else:
+                if( middleAtom > -1 and middleCovalency == 3 and middleAtom in skipAtoms ):
+
+                    partnerSet                    = skipAtoms[middleAtom]
+                  
+                    angleDict                     = {}
+
+                    angleList                     = []
+                    angleList.append( angle[0] )
+                    angleList.append( angle[1] )
+                    angleList.append( angle[2] )
+                    angleDict['angle']            = angleList
+
+                    angleDict['isConstrained']    = isConstrained
+
+                    angleSet  = set()
+                    angleSet.add( angle[0] )
+                    angleSet.add( angle[1] )
+                    angleSet.add( angle[2] )
+
+                    for atomIndex in partnerSet:
+                        if( atomIndex not in angleSet ):
+                            angleList.append( atomIndex )
+
+                    if( verbose ):
+                        print "%5d Opbend: middle=%d inPlane2 %d angleList=%s partnerSet=%s" % (opBondCount, middleAtom, len(angleList), str(angleList), str(partnerSet) )
+
+                    inPlaneAngles.append( angleDict )
+
+                else:
+                    angleDict                  = {}
+                    angleDict['angle']         = angle
+                    angleDict['isConstrained'] = isConstrained
+                    nonInPlaneAngles.append( angleDict )
+
+            count += 1
+
+        # get AmoebaHarmonicAngleGenerator and add AmoebaHarmonicAngle and AmoebaHarmonicInPlaneAngle forces
+
+        for force in self.forceField._forces:
+            if( force.__class__.__name__ == 'AmoebaHarmonicAngleGenerator' ): 
+                force.createForcePostOpBendAngle( sys, data, nonbondedMethod, nonbondedCutoff, nonInPlaneAngles, args )
+                force.createForcePostOpBendInPlaneAngle( sys, data, nonbondedMethod, nonbondedCutoff, inPlaneAngles, args )
+            if( force.__class__.__name__ == 'AmoebaHarmonicBondGenerator' ): 
+                force.createForce( sys, data, nonbondedMethod, nonbondedCutoff, args )
+
+        for force in self.forceField._forces:
+            if( force.__class__.__name__ == 'AmoebaStretchBendGenerator' ): 
+                for angleDict in inPlaneAngles:
+                    nonInPlaneAngles.append( angleDict )
+                force.createForcePostAmoebaHarmonicBondForce( sys, data, nonbondedMethod, nonbondedCutoff, nonInPlaneAngles, args )
+
+parsers["AmoebaOutOfPlaneBendForce"] = AmoebaOutOfPlaneBendGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaTorsionGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaTorsionGenerator constructs a AmoebaTorsionForce."""
+
+    #=============================================================================================
+
+    def __init__(self, torsionUnit):
+
+        self.torsionUnit   = torsionUnit
+
+        self.types1        = []
+        self.types2        = []
+        self.types3        = []
+        self.types4        = []
+
+        self.t1            = []
+        self.t2            = []
+        self.t3            = []
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaTorsionForce torsionUnit="0.5">
+        #   <Torsion class1="3" class2="1" class3="2" class4="3"   amp1="0.0" angle1="0.0"   amp2="0.0" angle2="3.14159265359"   amp3="0.0" angle3="0.0" />
+        #   <Torsion class1="3" class2="1" class3="2" class4="6"   amp1="0.0" angle1="0.0"   amp2="0.0" angle2="3.14159265359"   amp3="-0.263592" angle3="0.0" />
+         
+        generator = AmoebaTorsionGenerator( float(element.attrib['torsionUnit']) )
+        forceField._forces.append(generator)
+
+        # collect particle classes and t1,t2,t3,
+        # where ti=[amplitude_i,angle_i]
+
+        for torsion in element.findall('Torsion'):
+            types = forceField._findAtomTypes(torsion, 4)
+            if types is not None:
+
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+                generator.types4.append(types[3])
+
+                for ii in range(1,4):
+                    tInfo   = []
+                    suffix  = str(ii)
+                    ampName = 'amp' + suffix
+                    tInfo.append(float(torsion.attrib[ampName]))
+
+                    angName = 'angle' + suffix
+                    tInfo.append(float(torsion.attrib[angName]))
+
+                    if( ii == 1 ):
+                        generator.t1.append( tInfo )
+                    elif( ii == 2 ):
+                        generator.t2.append( tInfo )
+                    elif( ii == 3 ):
+                        generator.t3.append( tInfo )
+
+            else:
+                outputString = "AmoebaTorsionGenerator: error getting types: %s %s %s %s" % (
+                                    stretchBend.attrib['class1'],
+                                    stretchBend.attrib['class2'],
+                                    stretchBend.attrib['class3'],
+                                    stretchBend.attrib['class4'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    def createForce(self, sys, data, nontorsionedMethod, nontorsionedCutoff, args):
+
+        verbose            = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled = 1
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaTorsionForce]
+        if len(existing) == 0:
+            force = mm.AmoebaTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+        if( verbose ):
+            print "In AmoebaTorsionGenerator torsions=%d " % (len(data.propers))
+        count = 0
+        for torsion in data.propers:
+
+            type1 = data.atomType[data.atoms[torsion[0]]]
+            type2 = data.atomType[data.atoms[torsion[1]]]
+            type3 = data.atomType[data.atoms[torsion[2]]]
+            type4 = data.atomType[data.atoms[torsion[3]]]
+
+            hit   = 0
+            for i in range(len(self.types1)):
+
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+                types4 = self.types4[i]
+
+                # match types in forward or reverse direction
+
+                if (type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4) or (type4 in types1 and type3 in types2 and type2 in types3 and type1 in types4 ):
+                    if( verbose ):
+                        print "AmoebaTorsionGenerator %5d [%5d %5d %5d %5d] [%5s %5s %5s %5s]" % (
+                              count, torsion[0], torsion[1], torsion[2], torsion[3],
+                              type1, type2, type3, type4)
+                    hit  = 1
+                    force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], self.t1[i],  self.t2[i], self.t3[i] )
+                    break
+            if( hit == 0 ): 
+                print "AmoebaTorsionGenerator missing: %5d %s %s " % (count, type1, type2, type3, type4)
+
+            count += 1
+
+        if( verbose ):
+            print "AmoebaTorsionGenerator number of torsions added=%d " % (force.getNumTorsions())
+
+parsers["AmoebaTorsionForce"] = AmoebaTorsionGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaPiTorsionGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaPiTorsionGenerator constructs a AmoebaPiTorsionForce."""
+
+    #=============================================================================================
+    
+    def __init__(self, piTorsionUnit):
+        self.piTorsionUnit = piTorsionUnit 
+        self.types1        = []
+        self.types2        = []
+        self.k             = []
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaPiTorsionForce piTorsionUnit="1.0">
+        #   <PiTorsion class1="1" class2="3" k="28.6604" />
+
+        generator = AmoebaPiTorsionGenerator( float(element.attrib['piTorsionUnit']) )
+        forceField._forces.append(generator)
+
+        for piTorsion in element.findall('PiTorsion'):
+            types = forceField._findAtomTypes(piTorsion, 2)
+            if types is not None:
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.k.append(float(piTorsion.attrib['k']))
+            else:
+                outputString = "AmoebaPiTorsionGenerator: error getting types: %s %s " % (
+                                    piTorsion.attrib['class1'],
+                                    piTorsion.attrib['class2'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):
+
+        verbose            = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled = 1
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaPiTorsionForce]
+
+        if len(existing) == 0:
+            force = mm.AmoebaPiTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        count = 0
+        for bond in data.bonds:
+
+            # search for bonds with both atoms in bond having covalency == 3
+ 
+            atom1 = bond.atom1
+            atom2 = bond.atom2
+ 
+            if( len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom1]) == 3 ):
+
+                type1 = data.atomType[data.atoms[atom1]]
+                type2 = data.atomType[data.atoms[atom2]]
+
+                for i in range(len(self.types1)):
+
+                   types1 = self.types1[i]
+                   types2 = self.types2[i]
+
+                   if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                       if( verbose ):
+                            print "AmoebaPiTorsionGenerator %5d %5d %5d [%5s %5s] %15.6f" % (count, atom1, atom2, type1, type2, self.k[i])
+
+                       # piTorsionAtom1, piTorsionAtom2 are the atoms bonded to atom1, excluding atom2 
+                       # piTorsionAtom5, piTorsionAtom6 are the atoms bonded to atom2, excluding atom1 
+
+                       piTorsionAtom1 = -1
+                       piTorsionAtom2 = -1
+                       piTorsionAtom3 = atom1
+
+                       piTorsionAtom4 = atom2
+                       piTorsionAtom5 = -1
+                       piTorsionAtom6 = -1
+
+                       for bond in data.atomBonds[atom1]:
+                           bondedAtom1 = data.bonds[bond].atom1
+                           bondedAtom2 = data.bonds[bond].atom2
+                           if( bondedAtom1 != atom1 ):
+                               b1 = bondedAtom1
+                           else:
+                               b1 = bondedAtom2
+                           if( b1 != atom2 ):
+                               if( piTorsionAtom1 == -1 ):
+                                   piTorsionAtom1 = b1 
+                               else:
+                                   piTorsionAtom2 = b1
+
+                       for bond in data.atomBonds[atom2]:
+                           bondedAtom1 = data.bonds[bond].atom1
+                           bondedAtom2 = data.bonds[bond].atom2
+                           if( bondedAtom1 != atom2 ):
+                               b1 = bondedAtom1
+                           else:
+                               b1 = bondedAtom2
+
+                           if( b1 != atom1 ):
+                               if( piTorsionAtom5 == -1 ):
+                                   piTorsionAtom5 = b1 
+                               else:
+                                   piTorsionAtom6 = b1
+    
+                       force.addPiTorsion( piTorsionAtom1, piTorsionAtom2, piTorsionAtom3, piTorsionAtom4, piTorsionAtom5, piTorsionAtom6, self.k[i])
+
+            count += 1
+        if( verbose ):
+            print "AmoebaPiTorsionGenerator number of pi-torsions added=%d " % (force.getNumPiTorsions())
+
+parsers["AmoebaPiTorsionForce"] = AmoebaPiTorsionGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaTorsionTorsionGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaTorsionTorsionGenerator constructs a AmoebaTorsionTorsionForce."""
+    
+    #=============================================================================================
+
+    def __init__(self ):
+
+        self.types1        = []
+        self.types2        = []
+        self.types3        = []
+        self.types4        = []
+        self.types5        = []
+
+        self.gridIndex     = []
+
+        self.grids         = []
+
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        generator    = AmoebaTorsionTorsionGenerator( )
+        forceField._forces.append(generator)
+        maxGridIndex = -1
+
+        # <AmoebaTorsionTorsionForce >
+        # <TorsionTorsion class1="3" class2="1" class3="2" class4="3" class5="1" grid="0" nx="25" ny="25" />
+
+        for torsionTorsion in element.findall('TorsionTorsion'):
+            types = forceField._findAtomTypes( torsionTorsion, 5)
+            if types is not None:
+
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+                generator.types4.append(types[3])
+                generator.types5.append(types[4])
+
+                gridIndex  = int(torsionTorsion.attrib['grid'])
+                if( gridIndex > maxGridIndex ):
+                    maxGridIndex = gridIndex
+
+                generator.gridIndex.append(gridIndex)
+            else:
+                outputString = "AmoebaTorsionTorsionGenerator: error getting types: %s %s %s %s %s" % (
+                                    torsionTorsion.attrib['class1'],
+                                    torsionTorsion.attrib['class2'],
+                                    torsionTorsion.attrib['class3'],
+                                    torsionTorsion.attrib['class4'],
+                                    torsionTorsion.attrib['class5']  )
+                raise ValueError( outputString ) 
+    
+        # load grid
+
+        # xml source
+
+        # <TorsionTorsionGrid grid="0" nx="25" ny="25" >
+        # <Grid angle1="-180.00" angle2="-180.00" f="0.0" fx="2.31064374824e-05" fy="0.0" fxy="-0.0052801799672" />
+        # <Grid angle1="-165.00" angle2="-180.00" f="-0.66600912" fx="-0.06983370052" fy="-0.075058725744" fxy="-0.0044462732032" />
+
+        # output grid:
+
+        #     grid[x][y][0] = x value 
+        #     grid[x][y][1] = y value 
+        #     grid[x][y][2] = function value 
+        #     grid[x][y][3] = dfdx value 
+        #     grid[x][y][4] = dfdy value 
+        #     grid[x][y][5] = dfd(xy) value 
+
+        maxGridIndex    += 1
+        generator.grids  = maxGridIndex*[]
+        for torsionTorsionGrid in element.findall('TorsionTorsionGrid'):
+
+            gridIndex    = int( torsionTorsionGrid.attrib[ "grid"] )
+            nx           = int( torsionTorsionGrid.attrib[ "nx"] )
+            ny           = int( torsionTorsionGrid.attrib[ "ny"] )
+
+            grid         = []
+            gridCol      = []
+
+            gridColIndex = 0
+
+            for gridEntry in torsionTorsionGrid.findall('Grid'):
+
+                gridRow        = []
+                gridRow.append( float( gridEntry.attrib['angle1'] ) )
+                gridRow.append( float( gridEntry.attrib['angle2'] ) )
+                gridRow.append( float( gridEntry.attrib['f'] ) )
+                gridRow.append( float( gridEntry.attrib['fx'] ) )
+                gridRow.append( float( gridEntry.attrib['fy'] ) )
+                gridRow.append( float( gridEntry.attrib['fxy'] ) )
+                gridCol.append( gridRow )
+
+                gridColIndex  += 1
+                if( gridColIndex == nx ):
+                    grid.append( gridCol )
+                    gridCol      = []
+                    gridColIndex = 0
+
+            
+            if( gridIndex == len(generator.grids) ):
+                generator.grids.append( grid )
+            else:
+                while( len( generator.grids ) < gridIndex ):
+                    generator.grids.append( [] )
+                generator.grids[gridIndex] = grid
+
+    #=============================================================================================
+
+    def getChiralAtomIndex(self, data, sys, atomB, atomC, atomD ):
+
+        chiralAtomIndex = -1
+
+        # if atomC has four bonds, find the
+        # two bonds that do not include atomB and atomD
+        # set chiralAtomIndex to one of these, if they are
+        # not the same atom(type/mass)
+
+        if( len(data.atomBonds[atomC]) == 4 ):
+            atomE = -1
+            atomF = -1
+            for bond in data.atomBonds[atomC]:
+                bondedAtom1 = data.bonds[bond].atom1
+                bondedAtom2 = data.bonds[bond].atom2
+                hit         = -1
+                if(   bondedAtom1 == atomC and bondedAtom2 != atomB and bondedAtom2 != atomD ):
+                    hit = bondedAtom2
+                elif( bondedAtom2 == atomC and bondedAtom1 != atomB and bondedAtom1 != atomD ):
+                    hit = bondedAtom1
+
+                if( hit > -1 ):
+                    if( atomE == -1 ):
+                        atomE = hit
+                    else:
+                        atomF = hit
+       
+            # raise error if atoms E or F not found
+
+            if( atomE == -1 or atomF == -1 ):
+                outputString = "getChiralAtomIndex: error getting bonded partners of atomC=%s %d %s" % (atomC.name, atomC.resiude.index, atomC.resiude.name, )
+                raise ValueError( outputString ) 
+
+            # check for different type/mass between atoms E & F
+
+            typeE = int(data.atomType[data.atoms[atomE]])
+            typeF = int(data.atomType[data.atoms[atomF]])
+            if( typeE > typeF ):
+                chiralAtomIndex = atomE 
+            if( typeF > typeE ):
+                chiralAtomIndex = atomF 
+
+            massE = sys.getParticleMass( atomE )/unit.dalton
+            massF = sys.getParticleMass( atomE )/unit.dalton
+            if( massE > massF ):
+                chiralAtomIndex = massE 
+            if( massF > massE ):
+                chiralAtomIndex = massF 
+
+        return chiralAtomIndex
+
+    #=============================================================================================
+ 
+    def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):
+
+        verbose            = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled = 1
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaTorsionTorsionForce]
+
+        if len(existing) == 0:
+            force = mm.AmoebaTorsionTorsionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        count = 0
+        for angle in data.angles:
+
+            # search for bitorsions; based on TINKER subroutine bitors() 
+ 
+            ib = angle[0]
+            ic = angle[1]
+            id = angle[2]
+
+            for bondIndex in data.atomBonds[ib]:
+                bondedAtom1 = data.bonds[bondIndex].atom1
+                bondedAtom2 = data.bonds[bondIndex].atom2
+                if( bondedAtom1 != ib ):
+                    ia = bondedAtom1
+                else:
+                    ia = bondedAtom2
+
+                if( ia != ic and ia != id ):
+                    for bondIndex in data.atomBonds[id]:
+                        bondedAtom1 = data.bonds[bondIndex].atom1
+                        bondedAtom2 = data.bonds[bondIndex].atom2
+                        if( bondedAtom1 != id ):
+                            ie = bondedAtom1
+                        else:
+                            ie = bondedAtom2
+
+                        if( ie != ic and ie != ib and ie != ia ):
+
+                            # found candidate set of atoms
+                            # check if types match in order or reverse order
+
+                            type1 = data.atomType[data.atoms[ia]]
+                            type2 = data.atomType[data.atoms[ib]]
+                            type3 = data.atomType[data.atoms[ic]]
+                            type4 = data.atomType[data.atoms[id]]
+                            type5 = data.atomType[data.atoms[ie]]
+
+                            for i in range(len(self.types1)):
+
+                                types1 = self.types1[i]
+                                types2 = self.types2[i]
+                                types3 = self.types3[i]
+                                types4 = self.types4[i]
+                                types5 = self.types5[i]
+
+                                # match in order
+
+                                if( type1 in types1 and type2 in types2 and type3 in types3 and type4 in types4 and type5 in types5 ):
+                                    chiralAtomIndex = self.getChiralAtomIndex( data, sys, ib, ic, id )
+                                    force.addTorsionTorsion( ia, ib, ic, id, ie, chiralAtomIndex, self.gridIndex[i])
+
+                                # match in reverse order
+
+                                if( type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5 ):
+                                    chiralAtomIndex = self.getChiralAtomIndex( data, sys, ib, ic, id )
+                                    force.addTorsionTorsion( ie, id, ic, ib, ia, chiralAtomIndex, self.gridIndex[i])
+
+        # set grids
+
+        for (index, grid) in enumerate(self.grids):
+            force.setTorsionTorsionGrid( index, grid )
+ 
+        if( verbose ):
+            print "AmoebaTorsionTorsionGenerator number of bitorsions added=%d " % (force.getNumTorsionTorsions())
+            print "AmoebaTorsionTorsionGenerator number of grids      added=%d " % (force.getNumTorsionTorsionGrids())
+
+parsers["AmoebaTorsionTorsionForce"] = AmoebaTorsionTorsionGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaStretchBendGenerator:
+    """An AmoebaStretchBendGenerator constructs a AmoebaStretchBendForce."""
+    
+    #=============================================================================================
+
+    def __init__(self):
+
+        self.types1        = []
+        self.types2        = []
+        self.types3        = []
+
+        self.k1            = []
+        self.k2            = []
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+        generator = AmoebaStretchBendGenerator( )
+        forceField._forces.append(generator)
+
+        # <AmoebaStretchBendForce stretchBendUnit="1.0">
+        # <StretchBend class1="2" class2="1" class3="3" k1="5.25776946506" k2="5.25776946506" />
+        # <StretchBend class1="2" class2="1" class3="4" k1="3.14005676385" k2="3.14005676385" />
+
+        for stretchBend in element.findall('StretchBend'):
+            types = forceField._findAtomTypes(stretchBend, 3)
+            if types is not None:
+
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+
+                generator.k1.append(float(stretchBend.attrib['k1']))
+                generator.k2.append(float(stretchBend.attrib['k2']))
+
+            else:
+                outputString = "AmoebaStretchBendGenerator : error getting types: %s %s %s" % (
+                                    stretchBend.attrib['class1'],
+                                    stretchBend.attrib['class2'],
+                                    stretchBend.attrib['class3'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    # The setup of this force is dependent on AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce 
+    # having been called since the ideal bond lengths and angle are needed here.
+    # As a conseqeunce, createForce() is not implemented since it is not guaranteed that the generator for
+    # AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce have been called prior to AmoebaStretchBendGenerator(). 
+    # Instead, createForcePostAmoebaHarmonicBondForce() is called 
+    # after the generators for AmoebaHarmonicBondForce and AmoebaHarmonicAngleForce have been called
+
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+        pass
+
+    #=============================================================================================
+
+    # Note: request for constrained bonds is ignored.
+
+    #=============================================================================================
+
+    def createForcePostAmoebaHarmonicBondForce(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
+
+        verbose             = 0
+        if( self.hasBeenCalled ):
+             return
+        self.hasBeenCalled  = 1
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaStretchBendForce]
+        if len(existing) == 0:
+            force = mm.AmoebaStretchBendForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        if( verbose ):
+            print "In AmoebaStretchBendGenerator bonds=%d " % (len(data.bonds))
+
+        count  = 0
+        for angleDict in angleList:
+
+            angle = angleDict['angle']
+            if( 'isConstrained' in angleDict ):
+                isConstrained = angleDict['isConstrained']
+            else:
+                isConstrained = 0
+
+            type1 = data.atomType[data.atoms[angle[0]]]
+            type2 = data.atomType[data.atoms[angle[1]]]
+            type3 = data.atomType[data.atoms[angle[2]]]
+
+            hit     = 0
+            radian  = 57.2957795130
+            for i in range(len(self.types1)):
+
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+
+                # match types
+                # get ideal bond lengths, bondAB, bondCB
+                # get ideal angle
+
+                if( type2 in types2 and ( (type1 in types1 and type3 in types3) or (type3 in types1 and type1 in types3) ) ): 
+                    if isConstrained:
+                        hit = 1 
+                        # sys.addConstraint(bond.atom1, bond.atom2, self.length[i])
+                        # print message
+                    else:
+                        hit     = 1
+                        bondAB  = -1.0
+                        bondCB  = -1.0
+                        swap    = 0
+                        for bond in data.atomBonds[angle[1]]:
+                            atom1  = data.bonds[bond].atom1
+                            atom2  = data.bonds[bond].atom2
+                            length = data.bonds[bond].length
+                            if( atom1 == angle[0] ):
+                                bondAB = length
+                            if( atom1 == angle[2] ):
+                                bondCB = length
+                            if( atom2 == angle[2] ):
+                                bondCB = length
+                            if( atom2 == angle[0] ):
+                                bondAB = length
+                        
+                        # check that ideal angle and bonds are set
+
+                        if( 'idealAngle' not in angleDict ):
+                           outputString  = "AmoebaStretchBendGenerator: ideal angle is not set for following entry:\n"
+                           outputString += "   %5d [%5d %5d %5d] [%5s %5s %5s]" % (count, angle[0], angle[1], angle[2], type1, type2, type3 )
+                           raise ValueError( outputString ) 
+
+                        elif( bondAB < 0 or bondCB < 0 ):
+                           outputString  = "AmoebaStretchBendGenerator: bonds not set: %15.7e %15.7e. for following entry:" % ( bondAB, bondCB)
+                           outputString += "     %5d [%5d %5d %5d] [%5s %5s %5s]" % (count, angle[0], angle[1], angle[2], type1, type2, type3 )
+                           raise ValueError( outputString ) 
+
+                        else:
+                            if( verbose ):
+                                 print "AmoebaStretchBendGenerator %5d [%5d %5d %5d] [%5s %5s %5s] %15.6f %15.6f %15.6f %15.6f" % (count, angle[0], angle[1], angle[2], type1, type2, type3, bondAB, bondCB, angleDict['idealAngle'], self.k1[i])
+                            force.addStretchBend( angle[0], angle[1], angle[2], bondAB, bondCB, angleDict['idealAngle']/radian, self.k1[i])
+                    break
+
+            if( hit == 0 and verbose ): 
+                print "AmoebaStretchBendGenerator missing: %5d missing [%5d %5d %5d] [%5s %5s %5s] " % (count, angle[0], angle[1], angle[2], type1, type2, type3)
+
+            count += 1
+
+parsers["AmoebaStretchBendForce"] = AmoebaStretchBendGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaVdwGenerator:
+
+    """A AmoebaVdwGenerator constructs a AmoebaVdwForce."""
+    
+    #=============================================================================================
+
+    def __init__(self, type, radiusrule, radiustype, radiussize, epsilonrule, vdw13Scale, vdw14Scale, vdw15Scale ):
+
+        self.type           = type 
+
+        self.radiusrule     = radiusrule
+        self.radiustype     = radiustype
+        self.radiussize     = radiussize
+
+        self.epsilonrule    = epsilonrule
+
+        self.vdw13Scale     = vdw13Scale
+        self.vdw14Scale     = vdw14Scale
+        self.vdw15Scale     = vdw15Scale
+
+        self.typeMap        = {}
+
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        # <AmoebaVdwForce type="BUFFERED-14-7" radiusrule="CUBIC-MEAN" radiustype="R-MIN" radiussize="DIAMETER" epsilonrule="HHG" vdw-13-scale="0.0" vdw-14-scale="1.0" vdw-15-scale="1.0" >
+        #   <Vdw class="1" sigma="0.371" epsilon="0.46024" reduction="1.0" /> 
+        #   <Vdw class="2" sigma="0.382" epsilon="0.422584" reduction="1.0" /> 
+         
+        generator = AmoebaVdwGenerator( element.attrib['type'], element.attrib['radiusrule'], element.attrib['radiustype'], element.attrib['radiussize'], element.attrib['epsilonrule'], 
+                                        float(element.attrib['vdw-13-scale']), float(element.attrib['vdw-14-scale']), float(element.attrib['vdw-15-scale']) ) 
+        forceField._forces.append(generator)
+        two_six = 1.122462048309372
+
+        # types[] = [ sigma, epsilon, reductionFactor, class ]
+        # sigma is modified based on radiustype and radiussize
+
+        for atom in element.findall('Vdw'):
+            types = forceField._findAtomTypes(atom, 1)
+            if types is not None:
+
+                values    = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
+
+                classType = atom.attrib['class']
+
+                if( generator.radiustype == 'SIGMA' ):
+                    values[0] *= two_six
+      
+                if( generator.radiussize == 'DIAMETER' ):
+                    values[0] *= 0.5
+
+                values.append( classType )
+
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+            else:
+                outputString = "AmoebaVdwGenerator: error getting type: %s" % ( atom.attrib['class'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    # Return a set containing the indices of particles bonded to particle with index=particleIndex
+
+    #=============================================================================================
+
+    @staticmethod
+    def getBondedParticleSet( particleIndex, data ):
+
+        bondedParticleSet = set()
+
+        for bond in data.atomBonds[particleIndex]:
+            atom1 = data.bonds[bond].atom1
+            atom2 = data.bonds[bond].atom2
+            if( atom1 != particleIndex ):
+                bondedParticleSet.add( atom1 )
+            else:
+                bondedParticleSet.add( atom2 )
+
+        return bondedParticleSet
+          
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        sigmaMap   = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
+        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'HHG':1}
+        verbose    = 0
+
+        # get or create force depending on whether it has already been added to the system
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaVdwForce]
+        if len(existing) == 0:
+            force = mm.AmoebaVdwForce()
+            sys.addForce(force)
+
+            # sigma and epsilon combining rules
+
+            if( 'sigmaCombiningRule' in args ):
+                sigmaRule = args['sigmaCombiningRule'].upper()
+                if( sigmaRule.upper() in sigmaMap ):
+                    force.setSigmaCombiningRule( sigmaRule.upper() )
+                else:
+                    stringList = ' ' . join( str(x) for x in sigmaMap.keys())
+                    print "sigma combining rule %s not recognized; valid values are %s; using default." % (sigmaRule, stringList)  
+            else:
+                force.setSigmaCombiningRule( self.radiusrule )
+
+            if( 'epsilonCombiningRule' in args ):
+                epsilonRule = args['epsilonCombiningRule'].upper()
+                if( epsilonRule.upper() in epsilonMap ):
+                    force.setEpsilonCombiningRule( epsilonRule.upper() )
+                else:
+                    stringList = ' ' . join( str(x) for x in epsilonMap.keys())
+                    print "epsilon combining rule %s not recognized; valid values are %s; using default." % (epsilonRule, stringList)  
+            else:
+                force.setEpsilonCombiningRule( self.epsilonrule )
+               
+            # cutoff
+
+            if( 'vdwCutoff' in args ):
+                force.setCutoff( float( args['vdwCutoff'] ) )
+
+            if( nonbondedMethod == PME ):
+                force.setPBC( 1 )
+                force.setUseNeighborList( 1 )
+           
+        else:
+            force = existing[0]
+
+        # add particles to force
+        # throw error if particle type not available 
+
+        for (i, atom) in enumerate(data.atoms):
+            t = data.atomType[atom]
+            if t in self.typeMap:
+
+                values      = self.typeMap[t]
+                classIndex  = int(values[3])
+   
+                # ivIndex = index of bonded partner for hydrogens; otherwise ivIndex = particle index
+
+                ivIndex     = i
+                mass        = sys.getParticleMass( i )/unit.dalton
+                if( mass < 1.9 and len( data.atomBonds[i]) == 1 ):
+                    bondIndex = data.atomBonds[i][0]
+                    if( data.bonds[bondIndex].atom1 == i ):
+                        ivIndex = data.bonds[bondIndex].atom2
+                    else:
+                        ivIndex = data.bonds[bondIndex].atom1
+
+                force.addParticle( ivIndex, classIndex, values[0], values[1], values[2])
+                if( verbose ):
+                    print "Vdw %5d %5d %2d %15.7e %15.7e %15.7e" % (i, ivIndex, classIndex, values[0], values[1], values[2])
+            else:
+                raise ValueError('No vdw type for atom %s' % (atom.name))
+
+        # set combining rules
+
+        # set particle exclusions: self, 1-2 and 1-3 bonds
+        # (1) collect in bondedParticleSets[i], 1-2 indices for all bonded partners of particle i
+        # (2) add 1-2,1-3 and self to exclusion set
+
+        bondedParticleSets = []
+        for i in range(len(data.atoms)):
+            bondedParticleSets.append( AmoebaVdwGenerator.getBondedParticleSet( i, data ) )
+
+        for (i,atom) in enumerate(data.atoms):
+ 
+            # 1-2 partners
+
+            exclusionSet = bondedParticleSets[i].copy()
+
+            # 1-3 partners
+
+            if( self.vdw13Scale == 0.0 ):
+                for bondedParticle in bondedParticleSets[i]:
+                    exclusionSet = exclusionSet.union( bondedParticleSets[bondedParticle] )
+
+            # self
+
+            exclusionSet.add(i)
+
+            if( verbose ):
+                print "VdwExcl %5d %s || %s" % (i, str(exclusionSet), str(bondedParticleSets[i]))
+                for atomIndex in sorted( exclusionSet ):
+                    atom = data.atoms[atomIndex]
+                    print "   %5d [%s %s %d]" % (atomIndex, atom.name, atom.residue.name, atom.residue.index )
+                print "\n" 
+            
+            force.setParticleExclusions( i, exclusionSet )
+
+parsers["AmoebaVdwForce"] = AmoebaVdwGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaMultipoleGenerator:
+
+    #=============================================================================================
+
+    """A AmoebaMultipoleGenerator constructs a AmoebaMultipoleForce."""
+    
+    #=============================================================================================
+
+    def __init__(self, forceField,
+                       direct11Scale, direct12Scale, direct13Scale, direct14Scale,
+                       mpole12Scale,  mpole13Scale,  mpole14Scale,  mpole15Scale,
+                       mutual11Scale, mutual12Scale, mutual13Scale, mutual14Scale,
+                       polar12Scale,  polar13Scale,  polar14Scale,  polar15Scale ):
+
+        self.forceField              = forceField
+
+        self.direct11Scale           = direct11Scale 
+        self.direct12Scale           = direct12Scale 
+        self.direct13Scale           = direct13Scale 
+        self.direct14Scale           = direct14Scale 
+
+        self.mpole12Scale            = mpole12Scale 
+        self.mpole13Scale            = mpole13Scale 
+        self.mpole14Scale            = mpole14Scale 
+        self.mpole15Scale            = mpole15Scale 
+
+        self.mutual11Scale           = mutual11Scale 
+        self.mutual12Scale           = mutual12Scale 
+        self.mutual13Scale           = mutual13Scale 
+        self.mutual14Scale           = mutual14Scale 
+
+        self.polar12Scale            = polar12Scale 
+        self.polar13Scale            = polar13Scale 
+        self.polar14Scale            = polar14Scale 
+        self.polar15Scale            = polar15Scale 
+
+        self.typeMap                 = {}
+        self.hasBeenCalled           = 0
+
+    #=============================================================================================
+    # Set axis type
+    #=============================================================================================
+
+    @staticmethod
+    def setAxisType( kIndices ):
+
+                # set axis type
+
+                kIndicesLen        = len( kIndices )
+                if( kIndicesLen > 3 ):
+                    ky             = kIndices[3]
+                else:
+                    ky             = 0
+   
+                if( kIndicesLen > 2 ):
+                    kx             = kIndices[2]
+                else:
+                    kx             = 0
+   
+                if( kIndicesLen > 1 ):
+                    kz             = kIndices[1]
+                else:
+                    kz             = 0
+
+                while( len( kIndices ) < 4 ):
+                    kIndices.append( 0 )
+
+                axisType = mm.AmoebaMultipoleForce.ZThenX
+                if( kz == 0 ):
+                    axisType = mm.AmoebaMultipoleForce.NoAxisType
+                if( kz != 0 and kx == 0 ):
+                    axisType = mm.AmoebaMultipoleForce.ZOnly
+                if( kz < 0 or kx < 0 ):
+                    axisType = mm.AmoebaMultipoleForce.Bisector
+                if( kx < 0 and ky < 0 ):
+                    axisType = mm.AmoebaMultipoleForce.ZBisect
+                if( kz < 0 and kx < 0 and ky  < 0 ):
+                    axisType = mm.AmoebaMultipoleForce.ThreeFold
+
+                kIndices[1] = abs( kz )   
+                kIndices[2] = abs( kx )   
+                kIndices[3] = abs( ky )   
+
+                return axisType
+
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #   <AmoebaMultipoleForce  direct11Scale="0.0"  direct12Scale="1.0"  direct13Scale="1.0"  direct14Scale="1.0"  mpole12Scale="0.0"  mpole13Scale="0.0"  mpole14Scale="0.4"  mpole15Scale="0.8"  mutual11Scale="1.0"  mutual12Scale="1.0"  mutual13Scale="1.0"  mutual14Scale="1.0"  polar12Scale="0.0"  polar13Scale="0.0"  polar14Intra="0.5"  polar14Scale="1.0"  polar15Scale="1.0"  > 
+        # <Multipole class="1"    kz="2"    kx="4"    c0="-0.22620" d1="0.08214" d2="0.00000" d3="0.34883" q11="0.11775" q21="0.00000" q22="-1.02185" q31="-0.17555" q32="0.00000" q33="0.90410"  />
+        # <Multipole class="2"    kz="1"    kx="3"    c0="-0.15245" d1="0.19517" d2="0.00000" d3="0.19687" q11="-0.20677" q21="0.00000" q22="-0.48084" q31="-0.01672" q32="0.00000" q33="0.68761"  />
+
+        generator = AmoebaMultipoleGenerator( forceField,
+                                              element.attrib['direct11Scale'],
+                                              element.attrib['direct12Scale'],
+                                              element.attrib['direct13Scale'],
+                                              element.attrib['direct14Scale'],
+
+                                              element.attrib['mpole12Scale'],
+                                              element.attrib['mpole13Scale'],
+                                              element.attrib['mpole14Scale'],
+                                              element.attrib['mpole15Scale'],
+
+                                              element.attrib['mutual11Scale'],
+                                              element.attrib['mutual12Scale'],
+                                              element.attrib['mutual13Scale'],
+                                              element.attrib['mutual14Scale'],
+
+                                              element.attrib['polar12Scale'],
+                                              element.attrib['polar13Scale'],
+                                              element.attrib['polar14Scale'],
+                                              element.attrib['polar15Scale'] )
+
+
+
+        forceField._forces.append(generator)
+
+        # set type map: [ kIndices, multipoles, AMOEBA/OpenMM axis type]
+
+        for atom in element.findall('Multipole'):
+            types = forceField._findAtomTypes(atom, 1)
+            if types is not None:
+
+                # print "Multipole Atom %s types=%s." %( atom.attrib['type'], str(types) )
+
+                # k-indices not provided default to 0
+
+                kIndices = [int(atom.attrib['type'])]
+
+                kStrings = [ 'kz', 'kx', 'ky' ]
+                for kString in kStrings:
+                    try:
+                        if( atom.attrib[kString] ):
+                             kIndices.append( int( atom.attrib[kString] ))
+                    except: 
+                        pass
+
+                # set axis type based on k-Indices 
+
+                axisType     = AmoebaMultipoleGenerator.setAxisType( kIndices )
+
+                # set multipole
+
+                charge       = float(atom.attrib['c0'])
+ 
+                conversion   = 1.0
+                dipole       = [ conversion*float(atom.attrib['d1']), conversion*float(atom.attrib['d2']), conversion*float(atom.attrib['d3'])]
+
+                quadrupole   = []
+                quadrupole.append( conversion*float(atom.attrib['q11']))
+                quadrupole.append( conversion*float(atom.attrib['q21']))
+                quadrupole.append( conversion*float(atom.attrib['q31']))
+                quadrupole.append( conversion*float(atom.attrib['q21']))
+                quadrupole.append( conversion*float(atom.attrib['q22']))
+                quadrupole.append( conversion*float(atom.attrib['q32']))
+                quadrupole.append( conversion*float(atom.attrib['q31']))
+                quadrupole.append( conversion*float(atom.attrib['q32']))
+                quadrupole.append( conversion*float(atom.attrib['q33']))
+
+                for t in types[0]:
+                    if( t not in generator.typeMap ):
+                        generator.typeMap[t] = []
+
+                    valueMap                    = dict()
+                    valueMap['classIndex']      = atom.attrib['type']
+                    valueMap['kIndices']        = kIndices
+                    valueMap['charge']          = charge 
+                    valueMap['dipole']          = dipole
+                    valueMap['quadrupole']      = quadrupole
+                    valueMap['axisType']        = axisType
+                    generator.typeMap[t].append( valueMap )
+    
+            else:
+                outputString  = "AmoebaMultipoleGenerator: error getting type for multipole: %s" % ( atom.attrib['class'] )
+                raise ValueError( outputString ) 
+    
+        # polarization parameters
+ 
+        for atom in element.findall('Polarize'):
+            types = forceField._findAtomTypes(atom, 1)
+            if types is not None:
+
+                classIndex     = atom.attrib['type']
+                polarizability = float( atom.attrib['polarizability'] )
+                thole          = float( atom.attrib['thole'] )
+                if( thole == 0 ):
+                    pdamp = 0
+                else:
+                    pdamp = pow( polarizability, 1.0/6.0 )
+
+                pgrpMap            = dict()
+                for index in range( 1, 7):
+                    pgrp = 'pgrp' + str(index)
+                    if( pgrp in atom.attrib ):
+                        pgrpMap[int(atom.attrib[pgrp])] = -1
+
+                for t in types[0]:
+                    if( t not in generator.typeMap ):
+                        outputString = "AmoebaMultipoleGenerator: polarize type not present: %s" % ( atom.attrib['type'] )
+                        raise ValueError( outputString ) 
+                    else:
+                        typeMapList = generator.typeMap[t]
+                        hit         = 0
+                        for (ii, typeMap ) in enumerate(typeMapList):
+
+                            if( typeMap['classIndex'] == classIndex ):
+                                typeMap['polarizability']    = polarizability
+                                typeMap['thole']             = thole
+                                typeMap['pdamp']             = pdamp
+                                typeMap['pgrpMap']           = pgrpMap 
+                                typeMapList[ii]              = typeMap
+                                hit                          = 1
+                                #print "Adding polarize for %s map=%d len=%d keys=%s" % ( classIndex, ii, len( typeMapList), str( typeMap.keys() ) )
+
+                        if( hit == 0 ):
+                            outputString = "AmoebaMultipoleGenerator: error getting type for polarize: class index=%s not in multipole list?" % ( atom.attrib['class'] )
+                            raise ValueError( outputString ) 
+    
+            else:
+                outputString = "AmoebaMultipoleGenerator: error getting type for polarize: %s" % ( atom.attrib['class'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    def setPolarGroups(self, data, bonded12ParticleSets, force ):
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+
+            # assign multipole parameters via only 1-2 connected atoms
+
+            multipoleDict             = atom.multipoleDict
+            pgrpMap                   = multipoleDict['pgrpMap']
+            bondedAtomIndices         = bonded12ParticleSets[atomIndex]
+            atom.stage                = -1
+            atom.polarizationGroupSet = list()
+            atom.polarizationGroups[atomIndex] = 1
+            for bondedAtomIndex in bondedAtomIndices:
+                bondedAtomType = int(data.atomType[data.atoms[bondedAtomIndex]])
+                bondedAtom     = data.atoms[bondedAtomIndex]
+                if( bondedAtomType in pgrpMap ):
+                    atom.polarizationGroups[bondedAtomIndex] = 1
+                    bondedAtom.polarizationGroups[atomIndex] = 1
+                    
+        # pgrp11
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+
+            if( len(  data.atoms[atomIndex].polarizationGroupSet ) > 0 ):
+                continue
+
+            group      = set()
+            visited    = set()
+            notVisited = set()
+            for pgrpAtomIndex in atom.polarizationGroups:
+                group.add( pgrpAtomIndex )
+                notVisited.add( pgrpAtomIndex )
+            visited.add( atomIndex )
+            while( len( notVisited ) > 0 ):
+                nextAtom = notVisited.pop()
+                if( nextAtom not in visited ):
+                   visited.add( nextAtom )
+                   for ii in data.atoms[nextAtom].polarizationGroups:
+                       group.add( ii )
+                       if( ii not in visited ):
+                           notVisited.add( ii )
+
+            pGroup = group
+            for pgrpAtomIndex in group:
+                data.atoms[pgrpAtomIndex].polarizationGroupSet.append( pGroup )
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+            atom.polarizationGroupSet[0] = sorted( atom.polarizationGroupSet[0] )
+            force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent11, atom.polarizationGroupSet[0] )
+
+        # pgrp12
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+
+            if( len(  data.atoms[atomIndex].polarizationGroupSet ) > 1 ):
+                continue
+
+            pgrp11 = set( atom.polarizationGroupSet[0] )
+            pgrp12 = set()
+            for pgrpAtomIndex in pgrp11:
+                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
+                    pgrp12 = pgrp12.union( data.atoms[bonded12].polarizationGroupSet[0] )
+            pgrp12 = pgrp12 - pgrp11
+            for pgrpAtomIndex in pgrp11:
+                data.atoms[pgrpAtomIndex].polarizationGroupSet.append( pgrp12 )
+                
+        for (atomIndex, atom) in enumerate(data.atoms):
+            atom.polarizationGroupSet[1] = sorted( atom.polarizationGroupSet[1] )
+            force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent12, atom.polarizationGroupSet[1] )
+
+        # pgrp13
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+
+            if( len( data.atoms[atomIndex].polarizationGroupSet ) > 2 ):
+                continue
+
+            pgrp11 = set( atom.polarizationGroupSet[0] )
+            pgrp12 = set( atom.polarizationGroupSet[1] )
+            pgrp13 = set()
+            for pgrpAtomIndex in pgrp12:
+                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
+                    pgrp13 = pgrp13.union( data.atoms[bonded12].polarizationGroupSet[0] )
+            pgrp13 = pgrp13 - pgrp12
+            pgrp13 = pgrp13 - set( pgrp11 )
+            for pgrpAtomIndex in pgrp11:
+                data.atoms[pgrpAtomIndex].polarizationGroupSet.append( pgrp13 )
+                
+        for (atomIndex, atom) in enumerate(data.atoms):
+            atom.polarizationGroupSet[2] = sorted( atom.polarizationGroupSet[2] )
+            force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent13, atom.polarizationGroupSet[2] )
+
+        # pgrp14
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+
+            if( len( data.atoms[atomIndex].polarizationGroupSet ) > 3 ):
+                continue
+
+            pgrp11 = set( atom.polarizationGroupSet[0] )
+            pgrp12 = set( atom.polarizationGroupSet[1] )
+            pgrp13 = set( atom.polarizationGroupSet[2] )
+            pgrp14 = set()
+            for pgrpAtomIndex in pgrp13:
+                for bonded12 in bonded12ParticleSets[pgrpAtomIndex]:
+                    pgrp14 = pgrp14.union( data.atoms[bonded12].polarizationGroupSet[0] )
+
+            pgrp14 = pgrp14 - pgrp13
+            pgrp14 = pgrp14 - pgrp12
+            pgrp14 = pgrp14 - set( pgrp11 )
+
+            for pgrpAtomIndex in pgrp11:
+                data.atoms[pgrpAtomIndex].polarizationGroupSet.append( pgrp14 )
+                
+        for (atomIndex, atom) in enumerate(data.atoms):
+            atom.polarizationGroupSet[3] = sorted( atom.polarizationGroupSet[3] )
+            force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent14, atom.polarizationGroupSet[3] )
+
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        if( self.hasBeenCalled  ):
+            return
+        self.hasBeenCalled = 1
+
+        methodMap = {NoCutoff:mm.AmoebaMultipoleForce.NoCutoff,
+                     PME:mm.AmoebaMultipoleForce.PME}
+        verbose = 0
+
+        # get or create force depending on whether it has already been added to the system
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaMultipoleForce]
+        if len(existing) == 0:
+            force = mm.AmoebaMultipoleForce()
+            sys.addForce(force)
+            if( nonbondedMethod not in methodMap ): 
+                print "Warning: AmoebaMultipoleForce: cutoff method not available using NoCutoff."
+                force.setNonbondedMethod( mm.AmoebaMultipoleForce.NoCutoff )
+            else:
+                force.setNonbondedMethod( methodMap[nonbondedMethod] )
+            force.setCutoffDistance(nonbondedCutoff)
+
+            if( 'ewaldErrorTolerance' in args ):
+                force.setEwaldErrorTolerance(float(args['ewaldErrorTolerance']))
+
+            if( 'polarization' in args ):
+                polarizationType = args['polarization']
+                if( polarizationType.lower() == 'direct' ):
+                    force.setPolarizationType(mm.AmoebaMultipoleForce.Direct)
+                else:
+                    force.setPolarizationType(mm.AmoebaMultipoleForce.Mutual)
+
+            if( 'aEwald' in args ):
+                force.setAEwald( float(args['aEwald']))
+
+            if( 'pmeGridDimensions' in args ):
+                force.setPmeGridDimensions(args['pmeGridDimensions'])
+
+            if( 'mutualInducedMaxIterations' in args ):
+                force.setMutualInducedMaxIterations( int( args['mutualInducedMaxIterations'] ))
+
+            if( 'mutualInducedTargetEpsilon' in args ):
+                force.setMutualInducedTargetEpsilon(float(args['mutualInducedTargetEpsilon']))
+
+        else:
+            force = existing[0]
+
+        #     OpenMM interface call
+
+        #     /** 
+        #      * Add multipole-related info for a particle 
+        #      *
+        #      * @param charge               the particle's charge
+        #      * @param molecularDipole      the particle's molecular dipole (vector of size 3)
+        #      * @param molecularQuadrupole  the particle's molecular quadrupole (vector of size 9)
+        #      * @param axisType             the particle's axis type ( ZThenX, Bisector )
+        #      * @param multipoleAtomZ       index of first atom used in constructing lab<->molecular frames
+        #      * @param multipoleAtomX       index of second atom used in constructing lab<->molecular frames
+        #      * @param multipoleAtomY       index of second atom used in constructing lab<->molecular frames
+        #      * @param thole                Thole parameter
+        #      * @param dampingFactor        dampingFactor parameter
+        #      * @param polarity             polarity parameter
+        #      *
+        #      * @return the index of the particle that was added
+        #      */
+        #     int addParticle( double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType,
+        #                      int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity )
+        # 
+
+        # add particles to force
+        # throw error if particle type not available 
+
+        # get 1-2, 1-3, 1-4, 1-5 bonded sets
+
+        # 1-2
+
+        bonded12ParticleSets = []
+        for i in range(len(data.atoms)):
+            bonded12ParticleSet = AmoebaVdwGenerator.getBondedParticleSet( i, data )
+            bonded12ParticleSet = set( sorted( bonded12ParticleSet ) )
+            bonded12ParticleSets.append( bonded12ParticleSet )
+
+        # 1-3
+
+        bonded13ParticleSets = []
+        for i in range(len(data.atoms)):
+            bonded13Set          = set()
+            bonded12ParticleSet = bonded12ParticleSets[i]
+            for j in bonded12ParticleSet: 
+                bonded13Set = bonded13Set.union( bonded12ParticleSets[j] )
+
+            # remove 1-2 and self from set
+
+            bonded13Set = bonded13Set - bonded12ParticleSet
+            selfSet     = set()
+            selfSet.add(i)
+            bonded13Set = bonded13Set - selfSet
+            bonded13Set = set( sorted( bonded13Set ) )
+            bonded13ParticleSets.append( bonded13Set )
+
+        # 1-4
+
+        bonded14ParticleSets = []
+        for i in range(len(data.atoms)):
+            bonded14Set          = set()
+            bonded13ParticleSet  = bonded13ParticleSets[i]
+            for j in bonded13ParticleSet: 
+                bonded14Set         = bonded14Set.union( bonded12ParticleSets[j] )
+           
+            # remove 1-3, 1-2 and self from set
+
+            bonded14Set = bonded14Set - bonded12ParticleSets[i]
+            bonded14Set = bonded14Set - bonded13ParticleSet
+            selfSet     = set()
+            selfSet.add(i)
+            bonded14Set = bonded14Set - selfSet
+            bonded14Set = set( sorted( bonded14Set ) )
+            bonded14ParticleSets.append( bonded14Set )
+
+        # 1-5
+
+        bonded15ParticleSets = []
+        for i in range(len(data.atoms)):
+            bonded15Set          = set()
+            bonded14ParticleSet  = bonded14ParticleSets[i]
+            for j in bonded14ParticleSet: 
+                bonded15Set         = bonded15Set.union( bonded12ParticleSets[j] )
+
+            # remove 1-4, 1-3, 1-2 and self from set
+
+            bonded15Set = bonded15Set - bonded12ParticleSets[i]
+            bonded15Set = bonded15Set - bonded13ParticleSets[i]
+            bonded15Set = bonded15Set - bonded14ParticleSet
+            selfSet     = set()
+            selfSet.add(i)
+            bonded15Set = bonded15Set - selfSet
+            bonded15Set = set( sorted( bonded15Set ) )
+            bonded15ParticleSets.append( bonded15Set )
+
+        for (atomIndex, atom) in enumerate(data.atoms):
+            t = data.atomType[atom]
+            if t in self.typeMap:
+
+                multipoleList      = self.typeMap[t]
+                hit                = 0
+                savedMultipoleDict = 0
+
+                # assign multipole parameters via only 1-2 connected atoms
+
+                for multipoleDict in multipoleList:
+
+                    if( hit != 0 ):
+                        break
+
+                    kIndices           = multipoleDict['kIndices']
+    
+                    kz                 = kIndices[1]   
+                    kx                 = kIndices[2]
+                    ky                 = kIndices[3]
+    
+                    #print "Atom %s of %s %d k: %d %d %d keys=%s." %( atom.name, atom.residue.name, atom.residue.index, kz, kx, ky, sorted( multipoleDict.keys() ) )
+
+                    # assign multipole parameters
+                    #    (1) get bonded partners
+                    #    (2) match parameter types
+    
+                    bondedAtomIndices = bonded12ParticleSets[atomIndex]
+                    zaxis             = -1
+                    xaxis             = -1
+                    yaxis             = -1
+                    for bondedAtomZIndex in bondedAtomIndices:
+
+                       if( hit != 0 ):
+                           break
+
+                       bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
+                       bondedAtomZ     = data.atoms[bondedAtomZIndex]
+                       if( bondedAtomZType == kz ):
+                          for bondedAtomXIndex in bondedAtomIndices:
+                              if( bondedAtomXIndex == bondedAtomZIndex or hit != 0 ):
+                                  continue
+                              bondedAtomXType = int(data.atomType[data.atoms[bondedAtomXIndex]])
+                              if( bondedAtomXType == kx ):
+                                  if( ky == 0 ):
+                                      zaxis              = bondedAtomZIndex
+                                      xaxis              = bondedAtomXIndex
+                                      savedMultipoleDict = multipoleDict
+                                      hit                = 1
+                                  else:
+                                      for bondedAtomYIndex in bondedAtomIndices:
+                                          if( bondedAtomYIndex == bondedAtomZIndex or bondedAtomYIndex == bondedAtomXIndex or hit != 0 ):
+                                              continue
+                                          bondedAtomYType = int(data.atomType[data.atoms[bondedAtomYIndex]])
+                                          if( bondedAtomYType == ky ):
+                                              zaxis              = bondedAtomZIndex
+                                              xaxis              = bondedAtomXIndex
+                                              yaxis              = bondedAtomYIndex
+                                              savedMultipoleDict = multipoleDict
+                                              hit                = 2
+                                         
+                # assign multipole parameters via 1-2 and 1-3 connected atoms
+
+                for multipoleDict in multipoleList:
+
+                    if( hit != 0 ):
+                        break
+
+                    kIndices            = multipoleDict['kIndices']
+    
+                    kz                  = kIndices[1]   
+                    kx                  = kIndices[2]
+                    ky                  = kIndices[3]
+    
+                    #print "Atom %s of %s %d k: %d %d %d keys=%s." %( atom.name, atom.residue.name, atom.residue.index, kz, kx, ky, sorted( multipoleDict.keys() ) )
+
+                    # assign multipole parameters
+                    #    (1) get bonded partners
+                    #    (2) match parameter types
+    
+                    bondedAtom12Indices = bonded12ParticleSets[atomIndex]
+                    bondedAtom13Indices = bonded13ParticleSets[atomIndex]
+
+                    zaxis               = -1
+                    xaxis               = -1
+                    yaxis               = -1
+
+                    for bondedAtomZIndex in bondedAtom12Indices:
+
+                       if( hit != 0 ):
+                           break
+
+                       bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
+                       bondedAtomZ     = data.atoms[bondedAtomZIndex]
+
+                       if( bondedAtomZType == kz ):
+                          for bondedAtomXIndex in bondedAtom13Indices:
+
+                              if( bondedAtomXIndex == bondedAtomZIndex or hit != 0 ):
+                                  continue
+                              bondedAtomXType = int(data.atomType[data.atoms[bondedAtomXIndex]])
+                              if( bondedAtomXType == kx and bondedAtomZIndex in bonded12ParticleSets[bondedAtomXIndex] ):
+                                  if( ky == 0 ):
+                                      zaxis              = bondedAtomZIndex
+                                      xaxis              = bondedAtomXIndex
+                                      savedMultipoleDict = multipoleDict
+                                      hit                = 3
+                                  else:
+                                      for bondedAtomYIndex in bondedAtom13Indices:
+                                          if( bondedAtomYIndex == bondedAtomZIndex or bondedAtomYIndex == bondedAtomXIndex or hit != 0 ):
+                                              continue
+                                          bondedAtomYType = int(data.atomType[data.atoms[bondedAtomYIndex]])
+                                          if( bondedAtomYType == ky and bondedAtomZIndex in bonded12ParticleSets[bondedAtomYIndex]):
+                                              zaxis              = bondedAtomZIndex
+                                              xaxis              = bondedAtomXIndex
+                                              yaxis              = bondedAtomYIndex
+                                              savedMultipoleDict = multipoleDict
+                                              hit                = 4
+                                         
+                # assign multipole parameters via only a z-defining atom
+
+                for multipoleDict in multipoleList:
+
+                    if( hit != 0 ):
+                        break
+
+                    kIndices            = multipoleDict['kIndices']
+    
+                    kz                  = kIndices[1]   
+                    kx                  = kIndices[2]   
+    
+                    zaxis               = -1
+                    xaxis               = -1
+                    yaxis               = -1
+
+                    for bondedAtomZIndex in bondedAtom12Indices:
+
+                        if( hit != 0 ):
+                            break
+
+                        bondedAtomZType = int(data.atomType[data.atoms[bondedAtomZIndex]])
+                        bondedAtomZ     = data.atoms[bondedAtomZIndex]
+
+                        if( kx == 0 and kz == bondedAtomZType ):
+                            kz                 = bondedAtomZIndex
+                            savedMultipoleDict = multipoleDict
+                            hit                = 5
+
+                # assign multipole parameters via no connected atoms
+
+                for multipoleDict in multipoleList:
+
+                    if( hit != 0 ):
+                        break
+
+                    kIndices            = multipoleDict['kIndices']
+    
+                    kz                  = kIndices[1]   
+    
+                    zaxis               = -1
+                    xaxis               = -1
+                    yaxis               = -1
+
+                    if( kz == 0 ):
+                        savedMultipoleDict = multipoleDict
+                        hit                = 6
+                
+                # add particle if there was a hit
+
+                if( hit != 0 ):
+                    if( verbose ):
+                        print "Multipole hit for Atom %5d %4s of %4s %4d type=%5s hitBranch=%d axes: %4d %4d %4d axisType=%1d." % ( atomIndex,
+                               atom.name, atom.residue.name, atom.residue.index, t, hit, zaxis, xaxis, yaxis, savedMultipoleDict['axisType'] )
+
+                    atom.multipoleDict      = savedMultipoleDict
+                    atom.polarizationGroups = dict()
+                    newIndex                = force.addParticle( savedMultipoleDict['charge'], savedMultipoleDict['dipole'], savedMultipoleDict['quadrupole'], savedMultipoleDict['axisType'],
+                                                                 zaxis, xaxis, yaxis, savedMultipoleDict['thole'], savedMultipoleDict['pdamp'], savedMultipoleDict['polarizability'])
+                    if( atomIndex == newIndex ):
+                        force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.Covalent12, bonded12ParticleSets[atomIndex] )
+                        force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.Covalent13, bonded13ParticleSets[atomIndex] )
+                        force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.Covalent14, bonded14ParticleSets[atomIndex] )
+                        force.setCovalentMap( atomIndex, mm.AmoebaMultipoleForce.Covalent15, bonded15ParticleSets[atomIndex] )
+                    else:
+                        raise ValueError("Atom %s of %s %d is out of synch!." %( atom.name, atom.residue.name, atom.residue.index ) )
+                else:
+                    raise ValueError("Atom %s of %s %d was not assigned." %( atom.name, atom.residue.name, atom.residue.index ) )
+            else:
+                raise ValueError('No multipole type for atom %s %s %d' % ( atom.name, atom.residue.name, atom.residue.index ) )
+
+        # set polar groups
+
+        self.setPolarGroups( data, bonded12ParticleSets, force )
+
+parsers["AmoebaMultipoleForce"] = AmoebaMultipoleGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaWcaDispersionGenerator:
+
+    """A AmoebaWcaDispersionGenerator constructs a AmoebaWcaDispersionForce."""
+    
+    #=========================================================================================
+
+    def __init__(self, epso, epsh, rmino, rminh, awater, slevy, dispoff, shctd ):
+
+        self.epso           = epso 
+        self.epsh           = epsh 
+        self.rmino          = rmino
+        self.rminh          = rminh
+        self.awater         = awater
+        self.slevy          = slevy
+        self.dispoff        = dispoff
+        self.shctd          = shctd 
+
+        self.typeMap        = {}
+
+    #=========================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaWcaDispersionForce epso="0.46024" epsh="0.056484" rmino="0.17025" rminh="0.13275" awater="33.428" slevy="1.0"  dispoff="0.026" shctd="0.81" >
+        #   <WcaDispersion class="1" radius="0.1855" epsilon="0.46024" />
+        #   <WcaDispersion class="2" radius="0.191" epsilon="0.422584" />
+      
+        generator = AmoebaWcaDispersionGenerator( element.attrib['epso'],
+                                                  element.attrib['epsh'],
+                                                  element.attrib['rmino'],
+                                                  element.attrib['rminh'],
+                                                  element.attrib['awater'], 
+                                                  element.attrib['slevy'],
+                                                  element.attrib['dispoff'],
+                                                  element.attrib['shctd'] ) 
+        forceField._forces.append(generator)
+
+        # typeMap[] = [ radius, epsilon ]
+
+        for atom in element.findall('WcaDispersion'):
+            types = forceField._findAtomTypes(atom, 1)
+            if types is not None:
+
+                values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
+                for t in types[0]:
+                    generator.typeMap[t] = values
+            else:
+                outputString = "AmoebaWcaDispersionGenerator: error getting type: %s" % ( atom.attrib['class'] )
+                raise ValueError( outputString ) 
+    
+    #=========================================================================================
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        verbose = 0
+
+        # get or create force depending on whether it has already been added to the system
+
+        existing = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing = [f for f in existing if type(f) == mm.AmoebaWcaDispersionForce]
+        if len(existing) == 0:
+            force = mm.AmoebaWcaDispersionForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        # add particles to force
+        # throw error if particle type not available 
+
+        force.setEpso(    float(self.epso    ) )
+        force.setEpsh(    float(self.epsh    ) )
+        force.setRmino(   float(self.rmino   ) )
+        force.setRminh(   float(self.rminh   ) )
+        force.setDispoff( float(self.dispoff ) )
+        force.setSlevy(   float(self.slevy   ) )
+        force.setAwater(  float(self.awater  ) )
+        force.setShctd(   float(self.shctd   ) )
+
+        for (i, atom) in enumerate(data.atoms):
+            t = data.atomType[atom]
+            if t in self.typeMap:
+
+                values      = self.typeMap[t]
+                force.addParticle( values[0], values[1])
+                if( verbose ):
+                    print "WcaDispersion %5d %15.7e %15.7e" % (i, values[0], values[1])
+            else:
+                raise ValueError('No WcaDispersion type for atom %s of %s %d' % (atom.name, atom.residue.name, atom.residue.index))
+
+parsers["AmoebaWcaDispersionForce"] = AmoebaWcaDispersionGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaGeneralizedKirkwoodGenerator:
+
+    """A AmoebaGeneralizedKirkwoodGenerator constructs a AmoebaGeneralizedKirkwoodForce."""
+    
+    #=========================================================================================
+
+    def __init__(self, forceField, solventDielectric, soluteDielectric, includeCavityTerm, probeRadius, surfaceAreaFactor ):
+
+        self.forceField                = forceField
+        self.solventDielectric         = solventDielectric
+        self.soluteDielectric          = soluteDielectric
+        self.includeCavityTerm         = includeCavityTerm
+        self.probeRadius               = probeRadius
+        self.surfaceAreaFactor         = surfaceAreaFactor
+        self.hasBeenCalled             = 0
+
+        self.radiusTypeMap             = {}
+        self.radiusTypeMap['Bondi']    = {}
+        bondiMap                       = self.radiusTypeMap['Bondi'] 
+        rscale                         = 1.03
+
+        bondiMap[0]                    = 0.00
+        bondiMap[1]                    = 0.12*rscale
+        bondiMap[2]                    = 0.14*rscale
+        bondiMap[5]                    = 0.18*rscale
+
+        bondiMap[6]                    = 0.170*rscale
+        bondiMap[7]                    = 0.155*rscale
+        bondiMap[8]                    = 0.152*rscale
+        bondiMap[9]                    = 0.147*rscale
+
+        bondiMap[10]                   = 0.154*rscale
+        bondiMap[14]                   = 0.210*rscale
+        bondiMap[15]                   = 0.180*rscale
+        bondiMap[16]                   = 0.180*rscale
+
+        bondiMap[17]                   = 0.175 *rscale
+        bondiMap[18]                   = 0.188*rscale
+        bondiMap[34]                   = 0.190*rscale
+        bondiMap[35]                   = 0.185*rscale
+
+        bondiMap[36]                   = 0.202*rscale
+        bondiMap[53]                   = 0.198*rscale
+        bondiMap[54]                   = 0.216*rscale
+
+    #=========================================================================================
+
+    def getObcShct( self, data, atomIndex):
+
+        atom         = data.atoms[atomIndex]
+        atomicNumber = atom.element.atomic_number
+        shct         = -1.0
+
+        # shct
+ 
+        if( atomicNumber == 1 ):                 # H(1)
+            shct = 0.85       
+        elif( atomicNumber == 6 ):               # C(6)
+            shct = 0.72         
+        elif( atomicNumber == 7 ):               # N(7)
+            shct = 0.79        
+        elif( atomicNumber == 8 ):               # O(8)
+            shct = 0.85       
+        elif( atomicNumber == 9 ):               # F(9)
+            shct = 0.88    
+        elif( atomicNumber == 15 ):              # P(15)              
+            shct = 0.86 
+        elif( atomicNumber == 16 ):              # S(16)
+            shct = 0.96
+        elif( atomicNumber == 26 ):              # Fe(26)
+            shct = 0.88
+
+        if( shct < 0.0 ): 
+            shct = 0.80
+            print "getObcShct: Warning no GK overlap scale factor for atom %s of %s %d using default value=%f" % ( atom.name, atom.residue.name, atom.residue.index, shct)
+ 
+        return shct 
+
+    #=========================================================================================
+
+    def getAmoebaTypeRadius( self, data, bondedAtomIndices, atomIndex):
+
+        atom         = data.atoms[atomIndex]
+        atomicNumber = atom.element.atomic_number
+        #print "CCC atomIndex %d ele=%d %s %s %d" % ( atomIndex, atom.element.atomic_number, atom.name, atom.residue.name, atom.residue.index)
+        radius       = -1.0
+
+        if( atomicNumber == 1 ):                  # H(1)
+ 
+            radius          = 0.132
+ 
+            if( len( bondedAtomIndices ) < 1 ):
+                 outputString = "AmoebaGeneralizedKirkwoodGenerator: error getting atom bonded to %s of %s %d " % ( atom.name, atom.residue.name, atom.residue.index )
+                 raise ValueError( outputString ) 
+ 
+            for bondedAtomIndex in bondedAtomIndices:
+                bondedAtomAtomicNumber  = data.atoms[bondedAtomIndex].element.atomic_number
+
+            if( bondedAtomAtomicNumber == 7 ):
+                radius = 0.11
+            if( bondedAtomAtomicNumber == 8 ):
+                radius = 0.105
+ 
+        elif( atomicNumber == 3 ):               # Li(3)
+            radius = 0.15
+        elif( atomicNumber == 6 ):               # C(6)
+            
+            radius = 0.20
+            if( len( bondedAtomIndices ) == 3 ):
+                radius = 0.205
+
+            elif( len( bondedAtomIndices ) == 4 ):
+                for bondedAtomIndex in bondedAtomIndices:
+                   bondedAtomAtomicNumber = data.atoms[bondedAtomIndex].element.atomic_number
+                   if( bondedAtomAtomicNumber == 7 or bondedAtomAtomicNumber == 8 ):
+                       radius = 0.175
+
+        elif( atomicNumber == 7 ):               # N(7)
+            radius = 0.16
+        elif( atomicNumber == 8 ):               # O(8)
+            radius = 0.155
+            if( len( bondedAtomIndices ) == 2 ):
+                radius = 0.145
+        elif( atomicNumber == 9 ):               # F(9)
+            radius = 0.154
+        elif( atomicNumber == 10 ):              
+            radius = 0.146
+        elif( atomicNumber == 11 ):              
+            radius = 0.209
+        elif( atomicNumber == 12 ):              
+            radius = 0.179
+        elif( atomicNumber == 14 ):              
+            radius = 0.189
+        elif( atomicNumber == 15 ):              # P(15)              
+            radius = 0.196
+        elif( atomicNumber == 16 ):              # S(16)
+            radius = 0.186
+        elif( atomicNumber == 17 ):              
+            radius = 0.182
+        elif( atomicNumber == 18 ):              
+            radius = 0.179
+        elif( atomicNumber == 19 ):              
+            radius = 0.223
+        elif( atomicNumber == 20 ):              
+            radius = 0.191
+        elif( atomicNumber == 35 ):         
+            radius = 2.00
+        elif( atomicNumber == 36 ):   
+            radius = 0.190
+        elif( atomicNumber == 37 ):              
+            radius = 0.226
+        elif( atomicNumber == 53 ):              
+            radius = 0.237
+        elif( atomicNumber == 54 ):              
+            radius = 0.207
+        elif( atomicNumber == 55 ):              
+            radius = 0.263
+        elif( atomicNumber == 56 ):         
+            radius = 0.230
+
+        if( radius < 0.0 ): 
+            radius = 2.0
+            print "Warning no GK radius for atom %s of %s %d using default value=%f" % ( atom.name, atom.residue.name, atom.residue.index, radius)
+ 
+        return radius
+
+    #=========================================================================================
+
+    def getBondiTypeRadius( self, data, bondedAtomIndices, atomIndex):
+
+        bondiMap     = self.radiusTypeMap['Bondi'] 
+        atom         = data.atoms[atomIndex]
+        atomicNumber = atom.element.atomic_number
+        if( atomicNumber in bondiMap ): 
+            radius = bondiMap[atomicNumber]
+        else:
+            radius = 0.206
+            print "Warning no Bondi radius for atom %s of %s %d using default value=%f" % ( atom.name, atom.residue.name, atom.residue.index, radius)
+ 
+        return radius
+
+    #=========================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
+        #   <GeneralizedKirkwood type="1" charge="-0.22620" shct="0.79"  />
+        #   <GeneralizedKirkwood type="2" charge="-0.15245" shct="0.72"  />
+        
+        generator = AmoebaGeneralizedKirkwoodGenerator( forceField, element.attrib['solventDielectric'], element.attrib['soluteDielectric'],
+                                                        element.attrib['includeCavityTerm'], 
+                                                        element.attrib['probeRadius'], element.attrib['surfaceAreaFactor'] ) 
+        forceField._forces.append(generator)
+
+    #=========================================================================================
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        if( self.hasBeenCalled  ):
+            return
+
+        verbose = 0
+      
+        # check if AmoebaMultipoleForce exists since charges needed
+        # if it has not been created, raise an error
+
+        existing                 = [sys.getForce(i) for i in range(sys.getNumForces())]
+        amoebaMultipoleForceList = [f for f in existing if type(f) == mm.AmoebaMultipoleForce]
+        if( len( amoebaMultipoleForceList ) > 0 ):
+            amoebaMultipoleForce = amoebaMultipoleForceList[0]
+        else:
+            # call AmoebaMultipoleForceGenerator.createForce() to ensure charges have been set
+
+            for force in self.forceField._forces:
+                if( force.__class__.__name__ == 'AmoebaMultipoleGenerator' ): 
+                    force.createForce( sys, data, nonbondedMethod, nonbondedCutoff, args )
+
+        # get or create force depending on whether it has already been added to the system
+
+        existing = [f for f in existing if type(f) == mm.AmoebaGeneralizedKirkwoodForce]
+        if len(existing) == 0:
+
+            force = mm.AmoebaGeneralizedKirkwoodForce()
+            sys.addForce(force)
+ 
+            if( 'solventDielectric' in args ):
+                force.setSolventDielectric( float( args['solventDielectric'] ) )
+            else:
+                force.setSolventDielectric(    float(self.solventDielectric ) )
+
+            if( 'soluteDielectric' in args ):
+                force.setSoluteDielectric( float( args['soluteDielectric'] ) )
+            else:
+                force.setSoluteDielectric(     float(self.soluteDielectric ) )
+
+            if( 'includeCavityTerm' in args ):
+                force.setIncludeCavityTerm( int( args['includeCavityTerm'] ) )
+            else:
+               force.setIncludeCavityTerm(    int(self.includeCavityTerm) )
+
+        else:
+            print "AmoebaGeneralizedKirkwoodForce exists"
+            force = existing[0]
+
+        self.hasBeenCalled = 1
+
+        # add particles to force
+        # throw error if particle type not available 
+
+        force.setProbeRadius(          float(self.probeRadius) )
+        force.setSurfaceAreaFactor(    float(self.surfaceAreaFactor) )
+
+        # 1-2
+
+        bonded12ParticleSets = []
+        for i in range(len(data.atoms)):
+            bonded12ParticleSet = AmoebaVdwGenerator.getBondedParticleSet( i, data )
+            bonded12ParticleSet = set( sorted( bonded12ParticleSet ) )
+            bonded12ParticleSets.append( bonded12ParticleSet )
+
+        radiusType = 'Bondi'
+        for atomIndex in range( 0, amoebaMultipoleForce.getNumMultipoles() ):
+            multipoleParameters = amoebaMultipoleForce.getMultipoleParameters( atomIndex )
+            if( radiusType == 'Amoeba' ):
+                radius              = self.getAmoebaTypeRadius( data, bonded12ParticleSets[atomIndex], atomIndex)
+            else:
+                radius              = self.getBondiTypeRadius( data, bonded12ParticleSets[atomIndex], atomIndex)
+            shct                = self.getObcShct( data, atomIndex)
+            force.addParticle( multipoleParameters[0], radius, shct )
+            if( verbose ):
+               print "GeneralizedKirkwood %5d %15.7e %15.7e" % ( atomIndex, multipoleParameters[0], radius, shct)
+
+parsers["AmoebaGeneralizedKirkwoodForce"] = AmoebaGeneralizedKirkwoodGenerator.parseElement
+
+#=============================================================================================
+
+class AmoebaUreyBradleyGenerator:
+
+    #=============================================================================================
+
+    """An AmoebaUreyBradleyGenerator constructs a AmoebaUreyBradleyForce."""
+
+    #=============================================================================================
+    
+    def __init__(self, cubic, quartic):
+
+        self.cubic         = cubic
+        self.quartic       = quartic
+
+        self.types1        = []
+        self.types2        = []
+        self.types3        = []
+
+        self.length        = []
+        self.k             = []
+
+        self.hasBeenCalled = 0
+    
+    #=============================================================================================
+
+    @staticmethod
+    def parseElement(element, forceField):
+
+        #  <AmoebaUreyBradleyForce cubic="0.0" quartic="0.0"  >
+        #   <UreyBradley class1="74" class2="73" class3="74" k="16003.8" d="0.15537" /> 
+
+        generator = AmoebaUreyBradleyGenerator( float(element.attrib['cubic']), float(element.attrib['quartic']) )
+        forceField._forces.append(generator)
+        for bond in element.findall('UreyBradley'):
+            types = forceField._findAtomTypes(bond, 3)
+            if types is not None:
+
+                generator.types1.append(types[0])
+                generator.types2.append(types[1])
+                generator.types3.append(types[2])
+
+                generator.length.append(float(bond.attrib['d']))
+                generator.k.append(float(bond.attrib['k']))
+
+            else:
+                outputString = "AmoebaUreyBradleyGenerator: error getting types: %s %s %s" % (
+                                    bond.attrib['class1'], bond.attrib['class2'], bond.attrib['class3'] )
+                raise ValueError( outputString ) 
+    
+    #=============================================================================================
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args ):
+
+        verbose            = 0
+
+        if( self.hasBeenCalled ):
+             return
+ 
+        self.hasBeenCalled = 1
+
+        existing           = [sys.getForce(i) for i in range(sys.getNumForces())]
+        existing           = [f for f in existing if type(f) == mm.AmoebaUreyBradleyForce]
+
+        if len(existing) == 0:
+            force          = mm.AmoebaUreyBradleyForce()
+            sys.addForce(force)
+        else:
+            force = existing[0]
+
+        force.setAmoebaGlobalUreyBradleyCubic( self.cubic )
+        force.setAmoebaGlobalUreyBradleyQuartic( self.quartic )
+
+        for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
+            if( isConstrained ):
+                continue
+            type1 = data.atomType[data.atoms[angle[0]]]
+            type2 = data.atomType[data.atoms[angle[1]]]
+            type3 = data.atomType[data.atoms[angle[2]]]
+            for i in range(len(self.types1)):
+                types1 = self.types1[i]
+                types2 = self.types2[i]
+                types3 = self.types3[i]
+                if( (type1 in types1 and type2 in types2 and type3 in types3) or (type3 in types1 and type2 in types2 and type1 in types3) ):
+
+                    if( verbose ):
+                        print "AmoebaUreyBradleyGenerator %5d %5d %5d [%5s %5s %5s] %15.6f %15.6f" % (angle[0], angle[1], angle[2], type1, type2, type3, self.length[i], self.k[i])
+
+                    force.addUreyBradley(angle[0], angle[2], self.length[i], self.k[i])
+                    break
+
+parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
+
+#=============================================================================================

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

+#!/bin/env python
+
+#=============================================================================================
+# MODULE DOCSTRING
+#=============================================================================================
+
+"""
+Tools for constructing systems from AMBER prmtop/crd files.
+
+This module incorporates parts of 'zander', but Randall J. Radmer.
+Licensing to be worked out.
+
+@author Randall J. Radmer 
+@author John D. Chodera <jchodera@gmail.com>
+
+"""
+
+#=============================================================================================
+# GLOBAL IMPORTS
+#=============================================================================================
+
+import os
+import os.path
+import copy
+import re
+import math
+
+try:
+    import numpy
+except:
+    pass
+
+import simtk.unit as units
+import simtk.openmm
+
+#=============================================================================================
+# AMBER parmtop loader (from 'zander', by Randall J. Radmer)
+#=============================================================================================
+
+# A regex for extracting print format info from the FORMAT lines.
+FORMAT_RE_PATTERN=re.compile("([0-9]+)([a-zA-Z]+)([0-9]+)\.?([0-9]*)")
+
+# Pointer labels which map to pointer numbers at top of prmtop files
+POINTER_LABELS  = """
+              NATOM,  NTYPES, NBONH,  MBONA,  NTHETH, MTHETA,
+              NPHIH,  MPHIA,  NHPARM, NPARM,  NEXT,   NRES,
+              NBONA,  NTHETA, NPHIA,  NUMBND, NUMANG, NPTRA,
+              NATYP,  NPHB,   IFPERT, NBPER,  NGPER,  NDPER,
+              MBPER,  MGPER,  MDPER,  IFBOX,  NMXRS,  IFCAP
+"""
+
+# Pointer labels (above) as a list, not string.
+POINTER_LABEL_LIST = POINTER_LABELS.replace(',', '').split()
+
+class PrmtopLoader(object):
+    """Parsed AMBER prmtop file.
+
+    ParmtopLoader reads, parses and manages content from a AMBER prmtop file.
+
+    EXAMPLES
+
+    Parse a prmtop file of alanine dipeptide in implicit solvent.
+
+    >>> import os, os.path
+    >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
+    >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
+    >>> prmtop = PrmtopLoader(prmtop_filename)
+
+    Parse a prmtop file of alanine dipeptide in explicit solvent.
+
+    >>> 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)    
+
+    """
+    def __init__(self, inFilename):
+        """
+        Create a PrmtopLoader object from an AMBER prmtop file.
+
+        ARGUMENTS
+
+        inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap        
+
+        """
+
+        self._prmtopVersion=None
+        self._flags=[]
+        self._raw_format={}
+        self._raw_data={}
+
+        fIn=open(inFilename)
+        for line in fIn:
+            if line.startswith('%VERSION'):
+                tag, self._prmtopVersion = line.rstrip().split(None, 1)
+            elif line.startswith('%FLAG'):
+                tag, flag = line.rstrip().split(None, 1)
+                self._flags.append(flag)
+                self._raw_data[flag] = []
+            elif line.startswith('%FORMAT'):
+                format = line.rstrip()
+                index0=format.index('(')
+                index1=format.index(')')
+                self._raw_format[self._flags[-1]] = format[index0+1:index1]
+            elif self._flags \
+                 and 'TITLE'==self._flags[-1] \
+                 and not self._raw_data['TITLE']:
+                self._raw_data['TITLE'] = line.rstrip()
+            else:
+                flag=self._flags[-1]
+                (format, numItems, itemType,
+                 itemLength, itemPrecision) = self._getFormat(flag)
+                iLength=int(itemLength)
+                for index in range(0, len(line), iLength):
+                    item = line.rstrip()[index:index+iLength]
+                    if item:
+                        self._raw_data[flag].append(item)
+        fIn.close()
+
+    def _getFormat(self, flag=None):
+        if not flag:
+            flag=self._flags[-1]
+        format=self._raw_format[flag]
+        m=FORMAT_RE_PATTERN.search(format)
+        return (format, m.group(1),
+                m.group(2), m.group(3), m.group(4))
+
+    def _getPointerValue(self, pointerLabel):
+        """Return pointer value given pointer label
+
+           Parameter:
+            - pointerLabel: a string matching one of the following:
+
+            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
+            NTHETH : number of angles containing hydrogen
+            MTHETA : number of angles not containing hydrogen
+            NPHIH  : number of dihedrals containing hydrogen
+            MPHIA  : number of dihedrals not containing hydrogen
+            NHPARM : currently not used
+            NPARM  : currently not used
+            NEXT   : number of excluded atoms
+            NRES   : number of residues
+            NBONA  : MBONA + number of constraint bonds
+            NTHETA : MTHETA + number of constraint angles
+            NPHIA  : MPHIA + number of constraint dihedrals
+            NUMBND : number of unique bond types
+            NUMANG : number of unique angle types
+            NPTRA  : number of unique dihedral types
+            NATYP  : number of atom types in parameter file, see SOLTY below
+            NPHB   : number of distinct 10-12 hydrogen bond pair types
+            IFPERT : set to 1 if perturbation info is to be read in
+            NBPER  : number of bonds to be perturbed
+            NGPER  : number of angles to be perturbed
+            NDPER  : number of dihedrals to be perturbed
+            MBPER  : number of bonds with atoms completely in perturbed group
+            MGPER  : number of angles with atoms completely in perturbed group
+            MDPER  : number of dihedrals with atoms completely in perturbed groups
+            IFBOX  : set to 1 if standard periodic box, 2 when truncated octahedral
+            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) 
+        return float(self._raw_data['POINTERS'][index])
+
+    def getNumAtoms(self):
+        """Return the number of atoms in the system"""
+        return int(self._getPointerValue('NATOM'))
+
+    def getNumTypes(self):
+        """Return the number of AMBER atoms types in the system"""
+        return int(self._getPointerValue('NTYPES'))
+
+    def getIfBox(self):
+        """Return True if the system was build with periodic boundary conditions (PBC)"""
+        return int(self._getPointerValue('IFBOX'))
+
+    def getIfCap(self):
+        """Return True if the system was build with the cap option)"""
+        return int(self._getPointerValue('IFCAP'))
+
+    def getIfPert(self):
+        """Return True if the system was build with the perturbation parameters)"""
+        return int(self._getPointerValue('IFPERT'))
+
+    def getMasses(self):
+        """Return a list of atomic masses in the system"""
+        try:
+            return self._massList
+        except AttributeError:
+            pass
+
+        self._massList=[]
+        raw_masses=self._raw_data['MASS']
+        for ii in range(self.getNumAtoms()):
+            self._massList.append(float(raw_masses[ii]))
+        self._massList = units.Quantity(self._massList, units.amu)
+        return self._massList
+
+    def getCharges(self):
+        """Return a list of atomic charges in the system"""
+        try:
+            return self._chargeList
+        except AttributeError:
+            pass
+
+        self._chargeList=[]
+        raw_charges=self._raw_data['CHARGE']
+        for ii in range(self.getNumAtoms()):
+            self._chargeList.append(float(raw_charges[ii])/18.2223)
+        self._chargeList = units.Quantity(self._chargeList, units.elementary_charge)
+        return self._chargeList
+
+    def getAtomName(self, iAtom):
+        """Return the atom name for iAtom"""
+        atomNames = self.getAtomNames()
+        return atomNames[iAtom]
+
+    def getAtomNames(self):
+        """Return the list of the system atom names"""
+        return self._raw_data['ATOM_NAME']
+
+    def _getAtomTypeIndexes(self):
+        try:
+            return self._atomTypeIndexes
+        except AttributeError:
+            pass
+        self._atomTypeIndexes=[]
+        for atomTypeIndex in  self._raw_data['ATOM_TYPE_INDEX']:
+            self._atomTypeIndexes.append(int(atomTypeIndex))
+        return self._atomTypeIndexes
+
+    def getAtomType(self, iAtom):
+        """Return the AMBER atom type for iAtom"""
+        atomTypes=self.getAtomTypes()
+        return atomTypes[iAtom]
+
+    def getAtomTypes(self):
+        """Return the list of the AMBER atom types"""
+        return self._raw_data['AMBER_ATOM_TYPE']
+
+    def getResidueNumber(self, iAtom):
+        """Return iAtom's residue number"""
+        return self._getResiduePointer(iAtom)+1
+
+    def getResidueLabel(self, iAtom=None, iRes=None):
+        """Return residue label for iAtom OR iRes"""
+        if iRes==None and iAtom==None:
+            raise Exception("only specify iRes or iAtom, not both")
+        if iRes!=None and iAtom!=None:
+            raise Exception("iRes or iAtom must be set")
+        if iRes!=None:
+            return self._raw_data['RESIDUE_LABEL'][iRes]
+        else:
+            return self.getResidueLabel(iRes=self._getResiduePointer(iAtom))
+
+    def _getResiduePointer(self, iAtom):
+        try:
+            return self.residuePointerDict[iAtom]
+        except:
+            pass
+        resPointers=self._raw_data['RESIDUE_POINTER']
+        iRes=len(resPointers)
+        for ii in range(1, len(resPointers)):
+            firstAtom=int(resPointers[ii])-1
+            if firstAtom>iAtom:
+                iRes=ii
+                break
+        try:
+            self.residuePointerDict[iAtom]=iRes-1
+        except AttributeError:
+            self.residuePointerDict={iAtom:iRes-1}
+        return self.residuePointerDict[iAtom]
+
+    def getNonbondTerms(self):
+        """Return list of all rVdw, epsilon pairs for each atom"""
+        try:
+            return self._nonbondTerms
+        except AttributeError:
+            pass
+        self._nonbondTerms=[]
+        for iAtom in range(self.getNumAtoms()):
+            numTypes=self.getNumTypes()
+            atomTypeIndexes=self._getAtomTypeIndexes()
+            index=(numTypes+1)*(atomTypeIndexes[iAtom]-1)
+            nbIndex=int(self._raw_data['NONBONDED_PARM_INDEX'][index])-1
+            if nbIndex<0:
+                raise Exception("10-12 interactions are not supported")
+            acoef = float(self._raw_data['LENNARD_JONES_ACOEF'][nbIndex])
+            bcoef = float(self._raw_data['LENNARD_JONES_BCOEF'][nbIndex])
+            try:
+                rMin = (2*acoef/bcoef)**(1/6.0)
+                epsilon = 0.25*bcoef*bcoef/acoef
+            except ZeroDivisionError:
+                rMin = 1.0
+                epsilon = 0.0
+            rVdw =    units.Quantity(rMin/2.0, units.angstrom)
+            epsilon = units.Quantity(epsilon, units.kilocalorie_per_mole)
+            self._nonbondTerms.append( (rVdw, epsilon) )
+        return self._nonbondTerms
+
+    def _getBonds(self, bondPointers):
+        forceConstant=self._raw_data["BOND_FORCE_CONSTANT"]
+        bondEquil=self._raw_data["BOND_EQUIL_VALUE"]
+        returnList=[]
+        for ii in range(0,len(bondPointers),3):
+             if int(bondPointers[ii])<0 or \
+                int(bondPointers[ii+1])<0:
+                 raise Exception("Found negative bonded atom pointers %s"
+                                 % ((bondPointers[ii],
+                                     bondPointers[ii+1]),))
+             iType=int(bondPointers[ii+2])-1
+             forceConstUnit=units.kilocalorie_per_mole/(units.angstrom*units.angstrom)
+             returnList.append((int(bondPointers[ii])/3,
+                                int(bondPointers[ii+1])/3,
+                                units.Quantity(float(forceConstant[iType]), forceConstUnit),
+                                units.Quantity(float(bondEquil[iType]), units.angstrom)))
+        return returnList
+
+    def getBondsWithH(self):
+        """Return list of bonded atom pairs, K, and Rmin for each bond with a hydrogen"""
+        try:
+            return self._bondListWithH
+        except AttributeError:
+            pass
+        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"""
+        try:
+            return self._bondListNoH
+        except AttributeError:
+            pass
+        bondPointers=self._raw_data["BONDS_WITHOUT_HYDROGEN"]
+        self._bondListNoH = self._getBonds(bondPointers)
+        return self._bondListNoH
+
+    def getAngles(self):
+        """Return list of atom triplets, K, and ThetaMin for each bond angle"""
+        try:
+            return self._angleList
+        except AttributeError:
+            pass
+        forceConstant=self._raw_data["ANGLE_FORCE_CONSTANT"]
+        angleEquil=self._raw_data["ANGLE_EQUIL_VALUE"]
+        anglePointers = self._raw_data["ANGLES_INC_HYDROGEN"] \
+                       +self._raw_data["ANGLES_WITHOUT_HYDROGEN"]
+        self._angleList=[]
+        for ii in range(0,len(anglePointers),4):
+             if int(anglePointers[ii])<0 or \
+                int(anglePointers[ii+1])<0 or \
+                int(anglePointers[ii+2])<0:
+                 raise Exception("Found negative angle atom pointers %s"
+                                 % ((anglePointers[ii],
+                                     anglePointers[ii+1],
+                                     anglePointers[ii+2]),))
+             iType=int(anglePointers[ii+3])-1
+             forceConstUnit=units.kilocalorie_per_mole/(units.radian*units.radian)
+             self._angleList.append((int(anglePointers[ii])/3,
+                                int(anglePointers[ii+1])/3,
+                                int(anglePointers[ii+2])/3,
+                                units.Quantity(float(forceConstant[iType]), forceConstUnit),
+                                units.Quantity(float(angleEquil[iType])*180/math.pi, units.degree)))
+        return self._angleList
+
+    def getDihedrals(self):
+        """Return list of atom quads, K, phase and periodicity for each dihedral angle"""
+        try:
+            return self._dihedralList
+        except AttributeError:
+            pass
+        forceConstant=self._raw_data["DIHEDRAL_FORCE_CONSTANT"]
+        phase=self._raw_data["DIHEDRAL_PHASE"]
+        periodicity=self._raw_data["DIHEDRAL_PERIODICITY"]
+        dihedralPointers = self._raw_data["DIHEDRALS_INC_HYDROGEN"] \
+                          +self._raw_data["DIHEDRALS_WITHOUT_HYDROGEN"]
+        self._dihedralList=[]
+        for ii in range(0,len(dihedralPointers),5):
+             if int(dihedralPointers[ii])<0 or int(dihedralPointers[ii+1])<0:
+                 raise Exception("Found negative dihedral atom pointers %s"
+                                 % ((dihedralPointers[ii],
+                                    dihedralPointers[ii+1],
+                                    dihedralPointers[ii+2],
+                                    dihedralPointers[ii+3]),))
+             iType=int(dihedralPointers[ii+4])-1
+             self._dihedralList.append((int(dihedralPointers[ii])/3,
+                                int(dihedralPointers[ii+1])/3,
+                                abs(int(dihedralPointers[ii+2]))/3,
+                                abs(int(dihedralPointers[ii+3]))/3,
+                                units.Quantity(float(forceConstant[iType]), units.kilocalorie_per_mole),
+                                units.Quantity(float(phase[iType])*180/math.pi, units.degree),
+                                int(0.5+float(periodicity[iType]))))
+        return self._dihedralList
+
+    def get14Interactions(self):
+        """Return list of atom pairs, chargeProduct, rMin and epsilon for each 1-4 interaction"""
+        dihedralPointers = self._raw_data["DIHEDRALS_INC_HYDROGEN"] \
+                          +self._raw_data["DIHEDRALS_WITHOUT_HYDROGEN"]
+        returnList=[]
+        charges=self.getCharges()
+        nonbondTerms = self.getNonbondTerms()
+        for ii in range(0,len(dihedralPointers),5):
+             if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
+                 iAtom = int(dihedralPointers[ii])/3
+                 lAtom = int(dihedralPointers[ii+3])/3
+                 chargeProd = charges[iAtom]*charges[lAtom]
+                 (rVdwI, epsilonI) = nonbondTerms[iAtom]
+                 (rVdwL, epsilonL) = nonbondTerms[lAtom]
+                 rMin = (rVdwI+rVdwL)
+                 epsilon = units.sqrt(epsilonI*epsilonL)
+                 returnList.append((iAtom, lAtom, chargeProd, rMin, epsilon))
+        return returnList
+
+    def getExcludedAtoms(self):
+        """Return list of lists, giving all pairs of atoms that should have no non-bond interactions"""
+        try:
+            return self._excludedAtoms
+        except AttributeError:
+            pass
+        self._excludedAtoms=[]
+        numExcludedAtomsList=self._raw_data["NUMBER_EXCLUDED_ATOMS"]
+        excludedAtomsList=self._raw_data["EXCLUDED_ATOMS_LIST"]
+        total=0
+        for iAtom in range(self.getNumAtoms()):
+            index0=total
+            n=int(numExcludedAtomsList[iAtom])
+            total+=n
+            index1=total
+            atomList=[]
+            for jAtom in excludedAtomsList[index0:index1]:
+                j=int(jAtom)
+                if j>0:
+                    atomList.append(j-1)
+            self._excludedAtoms.append(atomList)
+        return self._excludedAtoms
+
+    def getGBSA_OBC(self):
+        """Return list giving GB params, Radius and screening factor"""
+        try:
+            return self._gbsa_obcList
+        except AttributeError:
+            pass
+        self._gbsa_obcList=[]
+        radii=self._raw_data["RADII"]
+        screen=self._raw_data["SCREEN"]
+        for iAtom in range(len(radii)):
+            self._gbsa_obcList.append( (units.Quantity(float(radii[iAtom]),
+                                                 units.angstrom),
+                                        units.Quantity(float(screen[iAtom]),
+                                                 units.dimensionless)) )
+        return self._gbsa_obcList
+
+    def getBoxBetaAndDimensions(self):
+        """Return periodic boundary box beta angle and dimensions"""
+        beta=float(self._raw_data["BOX_DIMENSIONS"][0])
+        x=float(self._raw_data["BOX_DIMENSIONS"][1])
+        y=float(self._raw_data["BOX_DIMENSIONS"][2])
+        z=float(self._raw_data["BOX_DIMENSIONS"][3])
+        return (units.Quantity(beta, units.degree),
+                units.Quantity(x, units.angstrom),
+                units.Quantity(y, units.angstrom),
+                units.Quantity(z, units.angstrom))
+
+#=============================================================================================
+# AMBER System builder (based on, but not identical to, systemManager from 'zander')
+#=============================================================================================
+
+def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=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)
+      nonbondedCutoff (float) - if specified, will set nonbondedCutoff (default: None)
+      scnb (float) - 1-4 Lennard-Jones scaling factor (default: 1.2)
+      scee (float) - 1-4 electrostatics scaling factor (default: 2.0)
+      mm - if specified, this module will be used in place of pyopenmm (default: None)
+      verbose (boolean) - if True, print out information on progress (default: False)
+      flexibleConstraints (boolean) - if True, flexible bonds will be added in addition ot constrained bonds
+
+    NOTES
+
+    Even if bonds are SHAKEn, their harmonic stretch terms are still included in the potential.
+
+    TODO
+
+    Should these option names be changed to reflect their 'sander' counterparts?
+
+    EXAMPLES
+
+    Create a system of alanine dipeptide in implicit solvent.
+
+    >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
+    >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
+    >>> system = readAmberSystem(prmtop_filename)
+
+    Parse a prmtop file of alanine dipeptide in explicit solvent.
+
+    >>> 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)    
+
+    """
+    
+    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:
+        raise Exception("Cannot specify both a filename and a loader")
+    if prmtop_filename is not None:
+        # Load prmtop file.
+        if verbose: print "Reading prmtop file '%s'..." % prmtop_filename
+        prmtop = PrmtopLoader(prmtop_filename)
+    else:
+        prmtop = prmtop_loader
+
+    if prmtop.getIfCap()>0:
+        raise Exception("CAP option not currently supported")
+
+    if prmtop.getIfPert()>0:
+        raise Exception("perturbation not currently supported")
+        
+    if prmtop.getIfBox()>1:
+        raise Exception("only standard periodic boxes are currently supported")
+
+    # Use pyopenmm implementation of OpenMM by default.
+    if mm is None:
+        mm = simtk.openmm
+
+    # Create OpenMM System.
+    if verbose: print "Creating OpenMM system..."
+    system = mm.System()
+
+    # Populate system with atomic masses.
+    if verbose: print "Adding particles..."
+    for mass in prmtop.getMasses():
+        system.addParticle(mass)
+
+    # Add constraints.
+    if shake in ('h-bonds', 'all-bonds', 'h-angles'):
+        for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
+            system.addConstraint(iAtom, jAtom, rMin)
+    if shake in ('all-bonds', 'h-angles'):
+        for (iAtom, jAtom, k, rMin) in prmtop.getBondsNoH():
+            system.addConstraint(iAtom, jAtom, rMin)
+            
+    # Add harmonic 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():
+            force.addBond(iAtom, jAtom, rMin, 2*k)                            
+    if flexibleConstraints or (shake not in ('all-bonds', 'h-angles')):
+        for (iAtom, jAtom, k, rMin) in prmtop.getBondsNoH():
+            force.addBond(iAtom, jAtom, rMin, 2*k)
+    system.addForce(force)
+
+    # Add harmonic angles.
+    if verbose: print "Adding angles..."    
+    force = mm.HarmonicAngleForce()
+    if shake == 'h-angles':
+        numConstrainedBonds = system.getNumConstraints()
+        atomConstraints = [[]]*system.getNumParticles()
+        for i in range(system.getNumConstraints()):
+            c = system.getConstraintParameters(i)
+            atomConstraints[c[0]].append((c[1], c[2]))
+            atomConstraints[c[1]].append((c[0], c[2]))
+    for (iAtom, jAtom, kAtom, k, aMin) in prmtop.getAngles():
+        if shake == 'h-angles':
+            type1 = prmtop.getAtomType(iAtom)
+            type2 = prmtop.getAtomType(jAtom)
+            type3 = prmtop.getAtomType(kAtom)
+            numH = len([type for type in (type1, type3) if type.startswith('H')])
+            constrained = (numH == 2 or (numH == 1 and type2.startswith('O')))
+        else:
+            constrained = False
+        if constrained:
+            # Find the two bonds that make this angle.
+            l1 = None
+            l2 = None
+            for bond in atomConstraints[jAtom]:
+                if bond[0] == iAtom:
+                    l1 = bond[1]
+                elif bond[0] == jAtom:
+                    l2 = bond[1]
+            
+            # Compute the distance between atoms and add a constraint
+            length = units.sqrt(l1*l1 + l2*l2 - 2*l1*l2*units.cos(aMin))
+            system.addConstraint(iAtom, kAtom, length)
+        if flexibleConstraints or not constrained:
+            force.addAngle(iAtom, jAtom, kAtom, aMin, 2*k)
+    system.addForce(force)
+
+    # Add 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..."    
+    force = mm.NonbondedForce()
+    if (prmtop.getIfBox() == 0):
+        # System is non-periodic.
+        if nonbondedMethod == 'NoCutoff':
+            force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            if nonbondedCutoff is None:
+                raise Exception("No cutoff value specified")
+            force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)            
+            force.setCutoffDistance(nonbondedCutoff)
+        else:
+            raise Exception("Illegal nonbonded method for a non-periodic system")
+    else:
+        # System is periodic. 
+        # Set periodic box vectors for periodic system
+        (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
+        d0 = units.Quantity(0.0, units.angstroms)
+        xVec = units.Quantity((boxX, d0,   d0))
+        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            
+            nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers)
+        force.setCutoffDistance(nonbondedCutoff)
+        
+        # Set nonbonded method.
+        if nonbondedMethod == 'NoCutoff':
+            force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
+        elif nonbondedMethod == 'CutoffPeriodic':
+            force.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
+        elif nonbondedMethod == 'Ewald':
+            force.setNonbondedMethod(mm.NonbondedForce.Ewald)
+        elif nonbondedMethod == 'PME':
+            force.setNonbondedMethod(mm.NonbondedForce.PME)
+        else:
+            raise Exception("Cutoff method not understood.")
+
+        if EwaldErrorTolerance is not None:
+            force.setEwaldErrorTolerance(EwaldErrorTolerance)
+
+    # Add per-particle nonbonded parameters.
+    for (charge, (rVdw, epsilon)) in zip(prmtop.getCharges(), prmtop.getNonbondTerms()):
+        sigma = rVdw * 2**(-1./6.) * 2.0
+        force.addParticle(charge, sigma, epsilon)
+
+    # Add 1-4 Interactions
+    excludedAtomPairs = set()
+    for (iAtom, lAtom, chargeProd, rMin, epsilon) in prmtop.get14Interactions():
+        chargeProd /= scee
+        epsilon /= scnb
+        sigma = rMin * 2**(-1./6.)
+        force.addException(iAtom, lAtom, chargeProd, sigma, epsilon)
+        excludedAtomPairs.add(min((iAtom, lAtom), (lAtom, iAtom)))
+
+    # Add Excluded Atoms
+    excludedAtoms=prmtop.getExcludedAtoms()
+    for iAtom in range(prmtop.getNumAtoms()):
+        for jAtom in excludedAtoms[iAtom]:
+            if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue            
+            force.addException(iAtom, jAtom, 0.0 * units.elementary_charge**2, 1.0 * units.angstroms, 0.0 * units.kilocalories_per_mole)
+
+    system.addForce(force)
+
+    # Add GBSA-OBC model.
+    if gbmodel == 'OBC':
+        if verbose: print "Adding GB parameters..."            
+        gb = mm.GBSAOBCForce()
+        charges = prmtop.getCharges()
+        gbsa_obc = prmtop.getGBSA_OBC()
+        #for charge, radius, scalingFactor in prmtop.getGBSA_OBC():
+        for iAtom in range(prmtop.getNumAtoms()):
+            gb.addParticle(charges[iAtom], gbsa_obc[iAtom][0], gbsa_obc[iAtom][1])
+        system.addForce(gb)
+        if nonbondedMethod == 'NoCutoff':
+            gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            gb.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
+            gb.setCutoffDistance(nonbondedCutoff)
+        elif nonbondedMethod == 'CutoffPeriodic':
+            gb.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
+            gb.setCutoffDistance(nonbondedCutoff)
+        else:
+            raise Exception("Illegal nonbonded method for use with GBSA")
+
+    # TODO: Add GBVI terms?
+
+    return system
+
+#=============================================================================================
+# AMBER INPCRD loader
+#=============================================================================================
+
+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.    
+
+    ARGUMENTS
+
+    filename (string) - name of Amber coordinates file to be read in
+    system (simtk.openmm.System) - System object for which coordinates are to be read
+
+    OPTIONAL ARGUMENTS
+
+    verbose (boolean) - if True, will print out verbose information about the file being read
+    asNumpy (boolean) - if True, results will be returned as Numpy arrays instead of lists of Vec3s
+
+    EXAMPLES
+
+    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')    
+    >>> 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')    
+    >>> [coordinates, box_vectors] = readAmberCoordinates(crd_filename, read_box=True)
+
+    """
+
+    # Open coordinate file for reading.
+    infile = open(filename, 'r')
+
+    # Read title
+    title = infile.readline().strip()
+    if verbose: print "title: '%s'" % title
+
+    # Read number of atoms
+    natoms = int(infile.readline().strip())
+    if verbose: print "%d atoms" % natoms
+
+    # Allocate storage for coordinates
+    coordinates = []
+
+    # Read coordinates
+    mm = simtk.openmm
+    natoms_read = 0
+    while (natoms_read < natoms):
+        line = infile.readline().strip()
+        elements = line.split()
+        while (len(elements) > 0):
+            coordinates.append(mm.Vec3(float(elements.pop(0)), float(elements.pop(0)), float(elements.pop(0))))
+            natoms_read += 1
+    if asNumpy:
+        newcoords = numpy.zeros([natoms,3], numpy.float32)
+        for i in range(len(coordinates)):
+            for j in range(3):
+                newcoords[i,j] = coordinates[i][j]
+        coordinates = newcoords
+    # Assign units.
+    coordinates = units.Quantity(coordinates, units.angstroms)
+
+    # Read velocities if requested.
+    velocities = None
+    if (read_velocities):
+        # Read velocities
+        velocities = []
+        natoms_read = 0
+        while (natoms_read < natoms):
+            line = infile.readline().strip()
+            elements = line.split()
+            while (len(elements) > 0):
+                velocities.append(mm.Vec3(float(elements.pop(0)), float(elements.pop(0)), float(elements.pop(0))))
+            natoms_read += 1
+        if asNumpy:
+            newvel = numpy.zeros([natoms,3], numpy.float32)
+            for i in range(len(velocities)):
+                for j in range(3):
+                    newvel[i,j] = velocities[i][j]
+            velocities = newvel
+        # Assign units.
+        velocities = units.Quantity(velocities, units.angstroms)
+            
+    # Read box size if present
+    box_vectors = None
+    if (read_box):
+        line = infile.readline().strip()
+        elements = line.split()
+        nelements = len(elements)
+        box_dimensions = [0.0]*nelements
+        for i in range(nelements):
+            box_dimensions[i] = float(elements[i])
+        # TODO: Deal with non-standard box sizes.
+        if nelements == 6:
+            if asNumpy:
+                a = units.Quantity(numpy.array([box_dimensions[0], 0.0, 0.0]), units.angstroms)
+                b = units.Quantity(numpy.array([0.0, box_dimensions[1], 0.0]), units.angstroms)
+                c = units.Quantity(numpy.array([0.0, 0.0, box_dimensions[2]]), units.angstroms)
+            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)            
+            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
+
+#=============================================================================================
+# MAIN AND TESTS
+#=============================================================================================
+
+if __name__ == "__main__":
+    import doctest
+    doctest.testmod()
+
+    

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

+#!/bin/env python
+
+"""
+pdbstructure.py: Used for managing PDB formated files.
+"""
+__author__ = "Christopher M. Bruns"
+__version__ = "1.0"
+
+
+from simtk.openmm.vec3 import Vec3
+import simtk.unit as unit
+from .. import element
+import warnings
+import sys
+
+class PdbStructure(object):
+    """
+    PdbStructure object holds a parsed Protein Data Bank format file.
+
+    Examples:
+
+    Load a pdb structure from a file:
+    > pdb = PdbStructure(open("1ARJ.pdb"))
+
+    Fetch the first atom of the structure:
+    > print pdb.iter_atoms().next()
+    ATOM      1  O5'   G N  17      13.768  -8.431  11.865  1.00  0.00           O
+
+    Loop over all of the atoms of the structure
+    > for atom in pdb.iter_atoms():
+    >     print atom
+    ATOM      1  O5'   G N  17      13.768  -8.431  11.865  1.00  0.00           O
+    ...
+
+    Get a list of all atoms in the structure:
+    > atoms = list(pdb.iter_atoms())
+
+    also:
+    residues = list(pdb.iter_residues())
+    positions = list(pdb.iter_positions())
+    chains = list(pdb.iter_chains())
+    models = list(pdb.iter_models())
+
+    Fetch atomic coordinates of first atom:
+    > print pdb.iter_positions().next()
+    [13.768, -8.431, 11.865] A
+
+     or
+
+    > print pdb.iter_atoms().next().position
+    [13.768, -8.431, 11.865] A
+
+    Strip the length units from an atomic position:
+    > import simtk.unit
+    > pos = pdb.iter_positions().next()
+    > print pos
+    [13.768, -8.431, 11.865] A
+    > print pos / simtk.unit.angstroms
+    [13.768, -8.431, 11.865]
+    > print pos / simtk.unit.nanometers
+    [1.3768, -0.8431, 1.1865]
+
+
+    The hierarchical structure of the parsed PDB structure is as follows:
+    PdbStructure
+      Model
+        Chain
+          Residue
+            Atom
+              Location
+
+    Model - A PDB structure consists of one or more Models.  Each model corresponds to one version of
+    an NMR structure, or to one frame of a molecular dynamics trajectory.
+
+    Chain - A Model contains one or more Chains.  Each chain corresponds to one molecule, although multiple
+    water molecules are frequently included in the same chain.
+
+    Residue - A Chain contains one or more Residues.  One Residue corresponds to one of the repeating
+    unit that constitutes a polymer such as protein or DNA.  For non-polymeric molecules, one Residue
+    represents one molecule.
+
+    Atom - A Residue contains one or more Atoms.  Atoms are chemical atoms.
+
+    Location - An atom can sometimes have more that one position, due to static disorder in X-ray
+    crystal structures.  To see all of the atom positions, use the atom.iter_positions() method,
+    or pass the parameter "include_alt_loc=True" to one of the other iter_positions() methods.
+
+    > for pos in pdb.iter_positions(include_alt_loc=True):
+    >   ...
+
+    Will loop over all atom positions, including multiple alternate locations for atoms that have
+    multiple positions.  The default value of include_alt_loc is False for the iter_positions()
+    methods.
+    """
+
+
+    def __init__(self, input_stream, load_all_models = False):
+        """Create a PDB model from a PDB file stream.
+
+        Parameters:
+         - self (PdbStructure) The new object that is created.
+         - input_stream (stream) An input file stream, probably created with
+             open().
+         - load_all_models (bool) Whether to load every model of an NMR
+             structure or trajectory, or just load the first model, to save memory.
+        """
+        # initialize models
+        self.load_all_models = load_all_models
+        self.models = []
+        self._current_model = None
+        self.default_model = None
+        self.models_by_number = {}
+        self._unit_cell_dimensions = None
+        # read file
+        self._load(input_stream)
+
+    def _load(self, input_stream):
+        # Read one line at a time
+        for pdb_line in input_stream:
+            # Look for atoms
+            if (pdb_line.find("ATOM  ") == 0) or (pdb_line.find("HETATM") == 0):
+                self._add_atom(Atom(pdb_line))
+            # Notice MODEL punctuation, for the next level of detail
+            # in the structure->model->chain->residue->atom->position hierarchy
+            elif (pdb_line.find("MODEL") == 0):
+                model_number = int(pdb_line[10:14])
+                self._add_model(Model(model_number))
+            elif (pdb_line.find("ENDMDL") == 0):
+                self._current_model._finalize()
+                if not self.load_all_models:
+                    break
+            elif (pdb_line.find("END") == 0):
+                self._current_model._finalize()
+                if not self.load_all_models:
+                    break
+            elif (pdb_line.find("TER   ") == 0):
+                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
+        self._finalize()
+
+    def write(self, output_stream=sys.stdout):
+        """Write out structure in PDB format"""
+        for model in self.models:
+            if len(model.chains) == 0:
+                continue
+            if len(self.models) > 1:
+                print >>output_stream, "MODEL     %4d" % (model.number)
+            model.write(output_stream)
+            if len(self.models) > 1:
+                print >>output_stream, "ENDMDL"
+        print >>output_stream, "END"
+
+    def _add_model(self, model):
+        if self.default_model == None:
+            self.default_model = model
+        self.models.append(model)
+        self._current_model = model
+        if model.number not in self.models_by_number:
+            self.models_by_number[model.number] = model
+
+    def get_model(self, model_number):
+        return self.models_by_number[model_number]
+
+    def model_numbers(self):
+        return self.models_by_number.keys()
+
+    def __contains__(self, model_number):
+        return self.models_by_number.__contains__(model_number)
+
+    def __getitem__(self, model_number):
+        return self.models_by_number[model_number]
+
+    def __iter__(self):
+        for model in self.models:
+                yield model
+
+    def iter_models(self, use_all_models=False):
+        if use_all_models:
+            for model in self:
+                yield model
+        elif len(self.models) > 0:
+            yield self.models[0]
+
+    def iter_chains(self, use_all_models=False):
+        for model in self.iter_models(use_all_models):
+            for chain in model.iter_chains():
+                yield chain
+
+    def iter_residues(self, use_all_models=False):
+        for model in self.iter_models(use_all_models):
+            for res in model.iter_residues():
+                yield res
+
+    def iter_atoms(self, use_all_models=False):
+        for model in self.iter_models(use_all_models):
+            for atom in model.iter_atoms():
+                yield atom
+
+    def iter_positions(self, use_all_models=False, include_alt_loc=False):
+        """
+        Iterate over atomic positions.
+
+        Parameters
+         - use_all_models (bool=False) Get positions from all models or just the first one.
+         - include_alt_loc (bool=False) Get all positions for each atom, or just the first one.
+        """
+        for model in self.iter_models(use_all_models):
+            for loc in model.iter_positions(include_alt_loc):
+                yield loc
+
+    def __len__(self):
+        return len(self.models)
+
+    def _add_atom(self, atom):
+        """
+        """
+        if self._current_model == None:
+            self._add_model(Model(0))
+        atom.model_number = self._current_model.number
+        # Atom might be alternate position for existing atom
+        self._current_model._add_atom(atom)
+
+    def _finalize(self):
+        """Establish first and last residues, atoms, etc."""
+        for model in self.models:
+            model._finalize()
+    
+    def get_unit_cell_dimensions(self):
+        """Get the dimensions of the crystallographic unit cell (may be None)."""
+        return self._unit_cell_dimensions
+
+
+class Model(object):
+    """Model holds one model of a PDB structure.
+
+    NMR structures usually have multiple models.  This represents one
+    of them.
+    """
+    def __init__(self, model_number=1):
+        self.number = model_number
+        self.chains = []
+        self._current_chain = None
+        self.chains_by_id = {}
+
+    def _add_atom(self, atom):
+        """
+        """
+        if len(self.chains) == 0:
+            self._add_chain(Chain(atom.chain_id))
+        # Create a new chain if the chain id has changed
+        if self._current_chain.chain_id != atom.chain_id:
+            self._add_chain(Chain(atom.chain_id))
+        # Create a new chain after TER record, even if ID is the same
+        elif self._current_chain.has_ter_record:
+            self._add_chain(Chain(atom.chain_id))
+        self._current_chain._add_atom(atom)
+
+    def _add_chain(self, chain):
+        self.chains.append(chain)
+        self._current_chain = chain
+        if not chain.chain_id in self.chains_by_id:
+            self.chains_by_id[chain.chain_id] = chain
+
+    def get_chain(self, chain_id):
+        return self.chains_by_id[chain_id]
+
+    def chain_ids(self):
+        return self.chains_by_id.keys()
+
+    def __contains__(self, chain_id):
+        return self.chains_by_id.__contains__(chain_id)
+
+    def __getitem__(self, chain_id):
+        return self.chains_by_id[chain_id]
+
+    def __iter__(self):
+        return iter(self.chains)
+
+    def iter_chains(self):
+        for chain in self:
+            yield chain
+
+    def iter_residues(self):
+        for chain in self:
+            for res in chain.iter_residues():
+                yield res
+
+    def iter_atoms(self):
+        for chain in self:
+            for atom in chain.iter_atoms():
+                yield atom
+
+    def iter_positions(self, include_alt_loc=False):
+        for chain in self:
+            for loc in chain.iter_positions(include_alt_loc):
+                yield loc
+
+    def __len__(self):
+        return len(self.chains)
+
+    def write(self, output_stream=sys.stdout):
+        # Start atom serial numbers at 1
+        sn = Model.AtomSerialNumber(1)
+        for chain in self.chains:
+            chain.write(sn, output_stream)
+
+    def _finalize(self):
+        for chain in self.chains:
+            chain._finalize()
+
+
+    class AtomSerialNumber(object):
+        """pdb.Model inner class for pass-by-reference incrementable serial number"""
+        def __init__(self, val):
+            self.val = val
+
+        def increment(self):
+            self.val += 1
+
+
+class Chain(object):
+    def __init__(self, chain_id=' '):
+        self.chain_id = chain_id
+        self.residues = []
+        self.has_ter_record = False
+        self._current_residue = None
+        self.residues_by_num_icode = {}
+        self.residues_by_number = {}
+
+    def _add_atom(self, atom):
+        """
+        """
+        # Create a residue if none have been created
+        if len(self.residues) == 0:
+            self._add_residue(Residue(atom.residue_name_with_spaces, atom.residue_number, atom.insertion_code, atom.alternate_location_indicator))
+        # Create a residue if the residue information has changed
+        elif self._current_residue.number != atom.residue_number:
+            self._add_residue(Residue(atom.residue_name_with_spaces, atom.residue_number, atom.insertion_code, atom.alternate_location_indicator))
+        elif self._current_residue.insertion_code != atom.insertion_code:
+            self._add_residue(Residue(atom.residue_name_with_spaces, atom.residue_number, atom.insertion_code, atom.alternate_location_indicator))
+        elif self._current_residue.name_with_spaces == atom.residue_name_with_spaces:
+            # This is a normal case: number, name, and iCode have not changed
+            pass
+        elif atom.alternate_location_indicator != ' ':
+            # OK - this is a point mutation, Residue._add_atom will know what to do
+            pass
+        else: # Residue name does not match
+            # Only residue name does not match
+            warnings.warn("WARNING: two consecutive residues with same number (%s, %s)" % (atom, self._current_residue.atoms[-1]))
+            self._add_residue(Residue(atom.residue_name_with_spaces, atom.residue_number, atom.insertion_code, atom.alternate_location_indicator))
+        self._current_residue._add_atom(atom)
+
+    def _add_residue(self, residue):
+        if len(self.residues) == 0:
+            residue.is_first_in_chain = True
+        self.residues.append(residue)
+        self._current_residue = residue
+        key = str(residue.number) + residue.insertion_code
+        # only store the first residue with a particular key
+        if key not in self.residues_by_num_icode:
+            self.residues_by_num_icode[key] = residue
+        if residue.number not in self.residues_by_number:
+            self.residues_by_number[residue.number] = residue
+
+    def write(self, next_serial_number, output_stream=sys.stdout):
+        for residue in self.residues:
+            residue.write(next_serial_number, output_stream)
+        if self.has_ter_record:
+            r = self.residues[-1]
+            print >>output_stream, "TER   %5d      %3s %1s%4d%1s" % (next_serial_number.val, r.name_with_spaces, self.chain_id, r.number, r.insertion_code)
+            next_serial_number.increment()
+
+    def _add_ter_record(self):
+        self.has_ter_record = True
+        self._finalize()
+
+    def get_residue(self, residue_number, insertion_code=' '):
+        return residues_by_num_icode[str(residue_number) + insertion_code]
+
+    def __contains__(self, residue_number):
+        return self.residues_by_number.__contains__(residue_number)
+
+    def __getitem__(self, residue_number):
+        """Returns the FIRST residue in this chain with a particular residue number"""
+        return self.residues_by_number[residue_number]
+
+    def __iter__(self):
+        for res in self.residues:
+            yield res
+
+    def iter_residues(self):
+        for res in self:
+            yield res
+
+    def iter_atoms(self):
+        for res in self:
+            for atom in res:
+                yield atom;
+
+    def iter_positions(self, include_alt_loc=False):
+        for res in self:
+            for loc in res.iter_positions(include_alt_loc):
+                yield loc
+
+    def __len__(self):
+        return len(self.residues)
+
+    def _finalize(self):
+        self.residues[0].is_first_in_chain = True
+        self.residues[-1].is_final_in_chain = True
+        for residue in self.residues:
+            residue._finalize()
+
+
+class Residue(object):
+    def __init__(self, name, number, insertion_code=' ', primary_alternate_location_indicator=' '):
+        alt_loc = primary_alternate_location_indicator
+        self.primary_location_id = alt_loc
+        self.locations = {}
+        self.locations[alt_loc] = Residue.Location(alt_loc, name)
+        self.name_with_spaces = name
+        self.number = number
+        self.insertion_code = insertion_code
+        self.atoms = []
+        self.atoms_by_name = {}
+        self.is_first_in_chain = False
+        self.is_final_in_chain = False
+        self._current_atom = None
+
+    def _add_atom(self, atom):
+        """
+        """
+        alt_loc = atom.alternate_location_indicator
+        if not self.locations.has_key(alt_loc):
+            self.locations[alt_loc] = Residue.Location(alt_loc, atom.residue_name_with_spaces)
+        assert atom.residue_number == self.number
+        assert atom.insertion_code == self.insertion_code
+        # Check whether this is an existing atom with another position
+        if (atom.name_with_spaces in self.atoms_by_name):
+            old_atom = self.atoms_by_name[atom.name_with_spaces]
+            # Unless this is a duplicated atom (warn about file error)
+            if atom.alternate_location_indicator in old_atom.locations:
+                warnings.warn("WARNING: duplicate atom (%s, %s)" % (atom, old_atom._pdb_string(old_atom.serial_number, atom.alternate_location_indicator)))
+            else:
+                for alt_loc, position in atom.locations.items():
+                    old_atom.locations[alt_loc] = position
+                return # no new atom added
+        # actually use new atom
+        self.atoms_by_name[atom.name] = atom
+        self.atoms_by_name[atom.name_with_spaces] = atom
+        self.atoms.append(atom)
+        self._current_atom = atom
+
+    def write(self, next_serial_number, output_stream=sys.stdout, alt_loc = "*"):
+        for atom in self.atoms:
+            atom.write(next_serial_number, output_stream, alt_loc)
+
+    def _finalize(self):
+        if len(self.atoms) > 0:
+            self.atoms[0].is_first_atom_in_chain = self.is_first_in_chain
+            self.atoms[-1].is_final_atom_in_chain = self.is_final_in_chain
+            for atom in self.atoms:
+                atom.is_first_residue_in_chain = self.is_first_in_chain
+                atom.is_final_residue_in_chain = self.is_final_in_chain
+
+    def set_name_with_spaces(self, name, alt_loc=None):
+        # Gromacs ffamber PDB files can have 4-character residue names
+        # assert len(name) == 3
+        if alt_loc == None:
+            alt_loc = self.primary_location_id
+        loc = self.locations[alt_loc]
+        loc.name_with_spaces = name
+        loc.name = name.strip()
+    def get_name_with_spaces(self, alt_loc=None):
+        if alt_loc == None:
+            alt_loc = self.primary_location_id
+        loc = self.locations[alt_loc]
+        return loc.name_with_spaces
+    name_with_spaces = property(get_name_with_spaces, set_name_with_spaces, doc='four-character residue name including spaces')
+
+    def get_name(self, alt_loc=None):
+        if alt_loc == None:
+            alt_loc = self.primary_location_id
+        loc = self.locations[alt_loc]
+        return loc.name
+    name = property(get_name, doc='residue name')
+
+    def get_atom(self, atom_name):
+        return self.atoms_by_name[atom_name]
+
+    def __contains__(self, atom_name):
+        return self.atoms_by_name.__contains__(atom_name)
+
+    def __getitem__(self, atom_name):
+        """Returns the FIRST atom in this residue with a particular atom name"""
+        return self.atoms_by_name[atom_name]
+
+    def __iter__(self):
+        """
+        >>> pdb_lines = [ \
+                "ATOM    188  N   CYS A  42      40.714  -5.292  12.123  1.00 11.29           N",\
+                "ATOM    189  CA  CYS A  42      39.736  -5.883  12.911  1.00 10.01           C",\
+                "ATOM    190  C   CYS A  42      40.339  -6.654  14.087  1.00 22.28           C",\
+                "ATOM    191  O   CYS A  42      41.181  -7.530  13.859  1.00 13.70           O",\
+                "ATOM    192  CB  CYS A  42      38.949  -6.825  12.002  1.00  9.67           C",\
+                "ATOM    193  SG  CYS A  42      37.557  -7.514  12.922  1.00 20.12           S"]
+        >>> res = Residue("CYS", 42)
+        >>> for l in pdb_lines:
+        ...     res._add_atom(Atom(l))
+        ...
+        >>> for atom in res:
+        ...     print atom
+        ATOM    188  N   CYS A  42      40.714  -5.292  12.123  1.00 11.29           N  
+        ATOM    189  CA  CYS A  42      39.736  -5.883  12.911  1.00 10.01           C  
+        ATOM    190  C   CYS A  42      40.339  -6.654  14.087  1.00 22.28           C  
+        ATOM    191  O   CYS A  42      41.181  -7.530  13.859  1.00 13.70           O  
+        ATOM    192  CB  CYS A  42      38.949  -6.825  12.002  1.00  9.67           C  
+        ATOM    193  SG  CYS A  42      37.557  -7.514  12.922  1.00 20.12           S  
+        """
+        for atom in self.iter_atoms():
+            yield atom
+
+    # Three possibilities: primary alt_loc, certain alt_loc, or all alt_locs
+    def iter_atoms(self, alt_loc=None):
+        if alt_loc == None:
+            locs = [self.primary_location_id]
+        elif alt_loc == "":
+            locs = [self.primary_location_id]
+        elif alt_loc == "*":
+            locs = None
+        else:
+            locs = list(alt_loc)
+        # If an atom has any location in alt_loc, emit the atom
+        for atom in self.atoms:
+            use_atom = False # start pessimistic
+            for loc2 in atom.locations.keys():
+                # print "#%s#%s" % (loc2,locs)
+                if locs == None: # means all locations
+                    use_atom = True
+                elif loc2 in locs:
+                    use_atom = True
+            if use_atom:
+                yield atom
+
+    def iter_positions(self, include_alt_loc=False):
+        """
+        Returns one position per atom, even if an individual atom has multiple positions.
+
+        >>> pdb_lines = [ \
+                         "ATOM    188  N   CYS A  42      40.714  -5.292  12.123  1.00 11.29           N",\
+                         "ATOM    189  CA  CYS A  42      39.736  -5.883  12.911  1.00 10.01           C",\
+                         "ATOM    190  C   CYS A  42      40.339  -6.654  14.087  1.00 22.28           C",\
+                         "ATOM    191  O   CYS A  42      41.181  -7.530  13.859  1.00 13.70           O",\
+                         "ATOM    192  CB  CYS A  42      38.949  -6.825  12.002  1.00  9.67           C",\
+                         "ATOM    193  SG  CYS A  42      37.557  -7.514  12.922  1.00 20.12           S"]
+        >>> res = Residue("CYS", 42)
+        >>> for l in pdb_lines: res._add_atom(Atom(l))
+        >>> for c in res.iter_positions:
+        ...     print c
+        Traceback (most recent call last):
+          File "<stdin>", line 1, in <module>
+        TypeError: 'instancemethod' object is not iterable
+        >>> for c in res.iter_positions():
+        ...     print c
+        [40.714, -5.292, 12.123] A
+        [39.736, -5.883, 12.911] A
+        [40.339, -6.654, 14.087] A
+        [41.181, -7.53, 13.859] A
+        [38.949, -6.825, 12.002] A
+        [37.557, -7.514, 12.922] A
+        """
+        for atom in self:
+            if include_alt_loc:
+                for loc in atom.iter_positions():
+                    yield loc
+            else:
+                yield atom.position
+
+    def __len__(self):
+        return len(self.atoms)
+
+    # Residues can have multiple locations, based on alt_loc indicator
+    class Location:
+        """
+        Inner class of residue to allow different residue names for different alternate_locations.
+        """
+        def __init__(self, alternate_location_indicator, residue_name_with_spaces):
+            self.alternate_location_indicator = alternate_location_indicator
+            self.residue_name_with_spaces = residue_name_with_spaces
+
+
+class Atom(object):
+    """Atom represents one atom in a PDB structure.
+    """
+    def __init__(self, pdb_line):
+        """Create a new pdb.Atom from an ATOM or HETATM line.
+
+        Example line:
+        ATOM   2209  CB  TYR A 299       6.167  22.607  20.046  1.00  8.12           C
+        00000000011111111112222222222333333333344444444445555555555666666666677777777778
+        12345678901234567890123456789012345678901234567890123456789012345678901234567890
+
+        ATOM line format description from
+          http://deposit.rcsb.org/adit/docs/pdb_atom_format.html:
+
+        COLUMNS        DATA TYPE       CONTENTS
+        --------------------------------------------------------------------------------
+         1 -  6        Record name     "ATOM  "
+         7 - 11        Integer         Atom serial number.
+        13 - 16        Atom            Atom name.
+        17             Character       Alternate location indicator.
+        18 - 20        Residue name    Residue name.
+        22             Character       Chain identifier.
+        23 - 26        Integer         Residue sequence number.
+        27             AChar           Code for insertion of residues.
+        31 - 38        Real(8.3)       Orthogonal coordinates for X in Angstroms.
+        39 - 46        Real(8.3)       Orthogonal coordinates for Y in Angstroms.
+        47 - 54        Real(8.3)       Orthogonal coordinates for Z in Angstroms.
+        55 - 60        Real(6.2)       Occupancy (Default = 1.0).
+        61 - 66        Real(6.2)       Temperature factor (Default = 0.0).
+        73 - 76        LString(4)      Segment identifier, left-justified.
+        77 - 78        LString(2)      Element symbol, right-justified.
+        79 - 80        LString(2)      Charge on the atom.
+
+        """
+        # We might modify first/final status during _finalize() methods
+        self.is_first_atom_in_chain = False
+        self.is_final_atom_in_chain = False
+        self.is_first_residue_in_chain = False
+        self.is_final_residue_in_chain = False
+        # Start parsing fields from pdb line
+        self.record_name = pdb_line[0:6].strip()
+        self.serial_number = int(pdb_line[6:11])
+        self.name_with_spaces = pdb_line[12:16]
+        alternate_location_indicator = pdb_line[16]
+
+        self.residue_name_with_spaces = pdb_line[17:20]
+        # In some MD codes, notably ffamber in gromacs, residue name has a fourth character in
+        # column 21
+        possible_fourth_character = pdb_line[20:21]
+        if possible_fourth_character != " ":
+            # Fourth character should only be there if official 3 are already full
+            if len(self.residue_name_with_spaces.strip()) != 3:
+                raise ValueError('Misaligned residue name: %s' % pdb_line)
+            self.residue_name_with_spaces += possible_fourth_character
+        self.residue_name = self.residue_name_with_spaces.strip()
+
+        self.chain_id = pdb_line[21]
+        self.residue_number = int(pdb_line[22:26])
+        self.insertion_code = pdb_line[26]
+        # coordinates, occupancy, and temperature factor belong in Atom.Location object
+        x = float(pdb_line[30:38])
+        y = float(pdb_line[38:46])
+        z = float(pdb_line[46:54])
+        try:
+            occupancy = float(pdb_line[54:60])
+        except:
+            occupancy = 1.0
+        try:
+            temperature_factor = float(pdb_line[60:66]) * unit.angstroms * unit.angstroms
+        except:
+            temperature_factor = 0.0
+        self.locations = {}
+        loc = Atom.Location(alternate_location_indicator, Vec3(x,y,z) * unit.angstroms, occupancy, temperature_factor, self.residue_name_with_spaces)
+        self.locations[alternate_location_indicator] = loc
+        self.default_location_id = alternate_location_indicator
+        # segment id, element_symbol, and formal_charge are not always present
+        self.segment_id = pdb_line[72:76].strip()
+        self.element_symbol = pdb_line[76:78].strip()
+        try: self.formal_charge = int(pdb_line[78:80])
+        except ValueError: self.formal_charge = None
+        # figure out atom element
+        try:
+            # First try to find a sensible element symbol from columns 76-77
+            self.element = element.get_by_symbol(self.element_symbol)
+        except KeyError:
+            # otherwise, deduce element from first two characters of atom name
+            # remove digits found in some hydrogen atom names
+            symbol = self.name_with_spaces[0:2].strip().lstrip("0123456789")
+            try:
+                # Some molecular dynamics PDB files, such as gromacs with ffamber force
+                # field, include 4-character hydrogen atom names beginning with "H".
+                # Hopefully elements like holmium (Ho) and mercury (Hg) will have fewer than four
+                # characters in the atom name.  This problem is the fault of molecular
+                # dynamics code authors who feel the need to make up their own atom
+                # nomenclature because it is too tedious to read that provided by the PDB.
+                # These are the same folks who invent their own meanings for biochemical terms
+                # like "dipeptide".  Clowntards.
+                if len(self.name) == 4 and self.name[0:1] == "H":
+                    self.element = element.hydrogen
+                else:
+                    self.element = element.get_by_symbol(symbol)
+            except KeyError:
+                # OK, I give up
+                self.element = None
+
+    def iter_locations(self):
+        """
+        Iterate over Atom.Location objects for this atom, including primary location.
+
+        >>> atom = Atom("ATOM   2209  CB  TYR A 299       6.167  22.607  20.046  1.00  8.12           C")
+        >>> for c in atom.iter_locations():
+        ...     print c
+        ...
+        [6.167, 22.607, 20.046] A
+        """
+        for alt_loc in self.locations:
+            yield self.locations[alt_loc]
+
+    def iter_positions(self):
+        """
+        Iterate over atomic positions.  Returns Quantity(Vec3(), unit) objects, unlike
+        iter_locations, which returns Atom.Location objects.
+        """
+        for loc in self.iter_locations():
+            yield loc.position
+
+    def iter_coordinates(self):
+        """
+        Iterate over x, y, z values of primary atom position.
+
+        >>> atom = Atom("ATOM   2209  CB  TYR A 299       6.167  22.607  20.046  1.00  8.12           C")
+        >>> for c in atom.iter_coordinates():
+        ...     print c
+        ...
+        6.167 A
+        22.607 A
+        20.046 A
+        """
+        for coord in self.position:
+            yield coord
+
+    # Hide existence of multiple alternate locations to avoid scaring casual users
+    def get_location(self, location_id=None):
+        id = location_id
+        if (id == None):
+            id = self.default_location_id
+        return self.locations[id]
+    def set_location(self, new_location, location_id=None):
+        id = location_id
+        if (id == None):
+            id = self.default_location_id
+        self.locations[id] = new_location
+    location = property(get_location, set_location, doc='default Atom.Location object')
+
+    def get_position(self):
+        return self.location.position
+    def set_position(self, coords):
+        self.location.position = coords
+    position = property(get_position, set_position, doc='orthogonal coordinates')
+
+    def get_alternate_location_indicator(self):
+        return self.location.alternate_location_indicator
+    alternate_location_indicator = property(get_alternate_location_indicator)
+
+    def get_occupancy(self):
+        return self.location.occupancy
+    occupancy = property(get_occupancy)
+
+    def get_temperature_factor(self):
+        return self.location.temperature_factor
+    temperature_factor = property(get_temperature_factor)
+
+    def get_x(self): return self.position[0]
+    x = property(get_x)
+
+    def get_y(self): return self.position[1]
+    y = property(get_y)
+
+    def get_z(self): return self.position[2]
+    z = property(get_z)
+
+    def _pdb_string(self, serial_number=None, alternate_location_indicator=None):
+        """
+        Produce a PDB line for this atom using a particular serial number and alternate location
+        """
+        if serial_number == None:
+            serial_number = self.serial_number
+        if alternate_location_indicator == None:
+            alternate_location_indicator = self.alternate_location_indicator
+        # produce PDB line in three parts: names, numbers, and end
+        # Accomodate 4-character residue names that use column 21
+        long_res_name = self.residue_name_with_spaces
+        if len(long_res_name) == 3:
+            long_res_name += " "
+        assert len(long_res_name) == 4
+        names = "%-6s%5d %4s%1s%4s%1s%4d%1s   " % (
+            self.record_name, serial_number, \
+            self.name_with_spaces, alternate_location_indicator, \
+            long_res_name, self.chain_id, \
+            self.residue_number, self.insertion_code)
+        numbers = "%8.3f%8.3f%8.3f%6.2f%6.2f      " % (
+            self.x.value_in_unit(unit.angstroms), \
+            self.y.value_in_unit(unit.angstroms), \
+            self.z.value_in_unit(unit.angstroms), \
+            self.occupancy, \
+            self.temperature_factor.value_in_unit(unit.angstroms * unit.angstroms))
+        end =  "%-4s%2s" % (\
+            self.segment_id, self.element_symbol)
+        formal_charge = "  "
+        if (self.formal_charge != None): formal_charge = "%+2d" % self.formal_charge
+        return names+numbers+end+formal_charge
+
+    def __str__(self):
+        return self._pdb_string(self.serial_number, self.alternate_location_indicator)
+
+    def write(self, next_serial_number, output_stream=sys.stdout, alt_loc = "*"):
+        """
+        alt_loc = "*" means write all alternate locations
+        alt_loc = None means write just the primary location
+        alt_loc = "AB" means write locations "A" and "B"
+        """
+        if alt_loc == None:
+            locs = [self.default_location_id]
+        elif alt_loc == "":
+            locs = [self.default_location_id]
+        elif alt_loc == "*":
+            locs = self.locations.keys()
+            locs.sort()
+        else:
+            locs = list(alt_loc)
+        for loc_id in locs:
+            print >>output_stream, self._pdb_string(next_serial_number.val, loc_id)
+            next_serial_number.increment()
+
+    def set_name_with_spaces(self, name):
+        assert len(name) == 4
+        self._name_with_spaces = name
+        self._name = name.strip()
+    def get_name_with_spaces(self):
+        return self._name_with_spaces
+    name_with_spaces = property(get_name_with_spaces, set_name_with_spaces, doc='four-character residue name including spaces')
+
+    def get_name(self):
+        return self._name
+    name = property(get_name, doc='residue name')
+
+    class Location(object):
+        """
+        Inner class of Atom for holding alternate locations
+        """
+        def __init__(self, alt_loc, position, occupancy, temperature_factor, residue_name):
+            self.alternate_location_indicator = alt_loc
+            self.position = position
+            self.occupancy = occupancy
+            self.temperature_factor = temperature_factor
+            self.residue_name = residue_name
+
+        def __iter__(self):
+            """
+            Examples
+
+            >>> from simtk.openmm.vec3 import Vec3
+            >>> import simtk.unit as unit
+            >>> l = Atom.Location(' ', Vec3(1,2,3)*unit.angstroms, 1.0, 20.0*unit.angstroms**2, "XXX")
+            >>> for c in l:
+            ...     print c
+            ...
+            1 A
+            2 A
+            3 A
+            """
+            for coord in self.position:
+                yield coord
+
+        def __str__(self):
+            return str(self.position)
+
+
+# run module directly for testing
+if __name__=='__main__':
+    # Test the examples in the docstrings
+    import doctest, sys
+    doctest.testmod(sys.modules[__name__])
+
+    import sys
+    import os
+    import gzip
+    import re
+    import time
+
+    # Test atom line parsing
+    pdb_line = "ATOM   2209  CB  TYR A 299       6.167  22.607  20.046  1.00  8.12           C"
+    a = Atom(pdb_line)
+    assert a.record_name == "ATOM"
+    assert a.serial_number == 2209
+    assert a.name == "CB"
+    assert a.name_with_spaces == " CB "
+    assert a.residue_name == "TYR"
+    assert a.residue_name_with_spaces == "TYR"
+    assert a.chain_id == "A"
+    assert a.residue_number == 299
+    assert a.insertion_code == " "
+    assert a.alternate_location_indicator == " "
+    assert a.x ==  6.167 * unit.angstroms
+    assert a.y == 22.607 * unit.angstroms
+    assert a.z == 20.046 * unit.angstroms
+    assert a.occupancy == 1.00
+    assert a.temperature_factor == 8.12 * unit.angstroms * unit.angstroms
+    assert a.segment_id == ""
+    assert a.element_symbol == "C"
+
+    # print pdb_line
+    # print str(a)
+    assert str(a).rstrip() == pdb_line.rstrip()
+
+    a = Atom("ATOM   2209  CB  TYR A 299       6.167  22.607  20.046  1.00  8.12           C")
+
+    # misaligned residue name - bad
+    try:
+        a = Atom("ATOM   2209  CB   TYRA 299       6.167  22.607  20.046  1.00  8.12           C")
+        assert(False)
+    except ValueError: pass
+    # four character residue name -- not so bad
+    a = Atom("ATOM   2209  CB  NTYRA 299       6.167  22.607  20.046  1.00  8.12           C")
+
+    atom_count = 0
+    residue_count = 0
+    chain_count = 0
+    model_count = 0
+    structure_count = 0
+
+    def parse_one_pdb(pdb_file_name):
+        global atom_count, residue_count, chain_count, model_count, structure_count
+        print pdb_file_name
+        if pdb_file_name[-3:] == ".gz":
+            fh = gzip.open(pdb_file_name)
+        else:
+            fh = open(pdb_file_name)
+        pdb = PdbStructure(fh, load_all_models=True)
+        # print "  %d atoms found" % len(pdb.atoms)
+        atom_count += len(list(pdb.iter_atoms()))
+        residue_count += len(list(pdb.iter_residues()))
+        chain_count += len(list(pdb.iter_chains()))
+        model_count += len(list(pdb.iter_models()))
+        structure_count += 1
+        fh.close
+        return pdb
+
+    # Parse one file
+    pdb_file_name = "/home/Christopher Bruns/Desktop/1ARJ.pdb"
+    if os.path.exists(pdb_file_name):
+        parse_one_pdb(pdb_file_name)
+
+    # try parsing the entire PDB
+    pdb_dir = "/cygdrive/j/pdb/data/structures/divided/pdb"
+    if os.path.exists(pdb_dir):
+        parse_entire_pdb = False
+        parse_one_division = False
+        parse_one_file = False
+        start_time = time.time()
+        if parse_one_file:
+            pdb_id = "2aed"
+            middle_two = pdb_id[1:3]
+            full_pdb_file = os.path.join(pdb_dir, middle_two, "pdb%s.ent.gz" % pdb_id)
+            parse_one_pdb(full_pdb_file)
+        if parse_one_division:
+            subdir = "ae"
+            full_subdir = os.path.join(pdb_dir, subdir)
+            for pdb_file in os.listdir(full_subdir):
+                if not re.match("pdb.%2s.\.ent\.gz" % subdir , pdb_file):
+                    continue
+                full_pdb_file = os.path.join(full_subdir, pdb_file)
+                parse_one_pdb(full_pdb_file)
+        if parse_entire_pdb:
+            for subdir in os.listdir(pdb_dir):
+                if not len(subdir) == 2: continue
+                full_subdir = os.path.join(pdb_dir, subdir)
+                if not os.path.isdir(full_subdir):
+                    continue
+                for pdb_file in os.listdir(full_subdir):
+                    if not re.match("pdb.%2s.\.ent\.gz" % subdir , pdb_file):
+                        continue
+                    full_pdb_file = os.path.join(full_subdir, pdb_file)
+                    parse_one_pdb(full_pdb_file)
+
+        end_time = time.time()
+        elapsed = end_time - start_time
+        minutes = elapsed / 60
+        seconds = elapsed % 60
+        hours = minutes / 60
+        minutes = minutes % 60
+        print "%dh:%02dm:%02ds elapsed" % (hours, minutes, seconds)
+
+        print "%d atoms found" % atom_count
+        print "%d residues found" % residue_count
+        print "%d chains found" % chain_count
+        print "%d models found" % model_count
+        print "%d structures found" % structure_count

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

+"""
+pdbfile.py: Used for loading PDB files.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import os
+import sys
+import xml.etree.ElementTree as etree
+from copy import copy
+from simtk.openmm import Vec3
+from simtk.openmm.app.internal.pdbstructure import PdbStructure
+from simtk.openmm.app import Topology
+from simtk.unit import nanometers, angstroms, is_quantity
+import element as elem
+    
+class PDBFile(object):
+    """PDBFile parses a Protein Data Bank (PDB) file and constructs a Topology and a set of atom positions from it.
+    
+    This class also provides methods for creating PDB files.  To write a file containing a single model, call
+    writeFile().  You also can create files that contain multiple models.  To do this, first call writeHeader(),
+    then writeModel() once for each model in the file, and finally writeFooter() to complete the file."""
+    
+    _residueNameReplacements = {}
+    _atomNameReplacements = {}
+    
+    def __init__(self, file):
+        """Load a PDB file.
+        
+        The atom positions and Topology are stored into this object's "positions" and "topology" fields.
+        
+        Parameters:
+         - file (string) the name of the file to load
+        """
+        top = Topology()
+        coords = [];
+        self.topology = top
+        
+        # Load the PDB file
+        
+        pdb = PdbStructure(open(file))
+        PDBFile._loadNameReplacementTables()
+
+        # Build the topology
+
+        for chain in pdb.iter_chains():
+            c = top.addChain()
+            for residue in chain.iter_residues():
+                resName = residue.get_name()
+                if resName in PDBFile._residueNameReplacements:
+                    resName = PDBFile._residueNameReplacements[resName]
+                r = top.addResidue(resName, c)
+                if resName in PDBFile._atomNameReplacements:
+                    atomReplacements = PDBFile._atomNameReplacements[resName]
+                else:
+                    atomReplacements = {}
+                for atom in residue.iter_atoms():
+                    atomName = atom.get_name()
+                    if atomName in atomReplacements:
+                        atomName = atomReplacements[atomName]
+                    atomName = atomName.strip()
+                    element = None
+
+                    # Try to guess the element.
+                    
+                    upper = atomName.upper()
+                    if upper.startswith('CL'):
+                        element = elem.chlorine
+                    elif upper.startswith('NA'):
+                        element = elem.sodium
+                    elif upper.startswith('MG'):
+                        element = elem.magnesium
+                    elif upper.startswith('BE'):
+                        element = elem.beryllium
+                    elif upper.startswith('LI'):
+                        element = elem.lithium
+                    elif upper.startswith('K'):
+                        element = elem.potassium
+                    elif( len( residue ) == 1 and upper.startswith('CA') ):
+                        element = elem.calcium
+                    else:
+                        try:
+                            element = elem.get_by_symbol(atomName[0])
+                        except KeyError:
+                            pass
+                    top.addAtom(atomName, element, r)
+                    pos = atom.get_position().value_in_unit(nanometers)
+                    coords.append(Vec3(pos[0], pos[1], pos[2]))
+        self.positions = coords*nanometers
+        self.topology.setUnitCellDimensions(pdb.get_unit_cell_dimensions())
+        self.topology.createStandardBonds()
+        self.topology.createDisulfideBonds(self.positions)
+
+    @staticmethod
+    def _loadNameReplacementTables():
+        """Load the list of atom and residue name replacements."""
+        if len(PDBFile._residueNameReplacements) == 0:
+            tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', 'pdbNames.xml'))
+            allResidues = {}
+            proteinResidues = {}
+            nucleicAcidResidues = {}
+            for residue in tree.getroot().findall('Residue'):
+                name = residue.attrib['name']
+                if name == 'All':
+                    PDBFile._parseResidueAtoms(residue, allResidues)
+                elif name == 'Protein':
+                    PDBFile._parseResidueAtoms(residue, proteinResidues)
+                elif name == 'Nucleic':
+                    PDBFile._parseResidueAtoms(residue, nucleicAcidResidues)
+            for atom in allResidues:
+                proteinResidues[atom] = allResidues[atom]
+                nucleicAcidResidues[atom] = allResidues[atom]
+            for residue in tree.getroot().findall('Residue'):
+                name = residue.attrib['name']
+                for id in residue.attrib:
+                    if id == 'name' or id.startswith('alt'):
+                        PDBFile._residueNameReplacements[residue.attrib[id]] = name
+                if 'type' not in residue.attrib:
+                    atoms = copy(allResidues)
+                elif residue.attrib['type'] == 'Protein':
+                    atoms = copy(proteinResidues)
+                elif residue.attrib['type'] == 'Nucleic':
+                    atoms = copy(nucleicAcidResidues)
+                else:
+                    atoms = copy(allResidues)
+                PDBFile._parseResidueAtoms(residue, atoms)
+                PDBFile._atomNameReplacements[name] = atoms
+
+    @staticmethod
+    def _parseResidueAtoms(residue, map):
+        for atom in residue.findall('Atom'):
+            name = atom.attrib['name']
+            for id in atom.attrib:
+                map[atom.attrib[id]] = name
+    
+    @staticmethod
+    def writeFile(topology, positions, file=sys.stdout, modelIndex=None):
+        """Write a PDB file containing a single model.
+                
+        Parameters:
+         - topology (Topology) The Topology defining the model to write
+         - positions (list) The list of atomic positions to write
+         - file (file=stdout) A file to write to
+        """
+        PDBFile.writeHeader(topology, file)
+        PDBFile.writeModel(topology, positions, file)
+        PDBFile.writeFooter(topology, file)
+
+    @staticmethod
+    def writeHeader(topology, file=sys.stdout):
+        """Write out the header for a PDB file.
+        
+        Parameters:
+         - topology (Topology) The Topology defining the molecular system being written
+         - file (file=stdout) A file to write the file to
+        """
+        boxSize = topology.getUnitCellDimensions()
+        if boxSize is not None:
+            size = boxSize.value_in_unit(angstroms)
+            print >>file, "CRYST1%9.3f%9.3f%9.3f  90.00  90.00  90.00 P 1           1 " % size
+    
+    @staticmethod
+    def writeModel(topology, positions, file=sys.stdout, modelIndex=None):
+        """Write out a model to a PDB file.
+                
+        Parameters:
+         - topology (Topology) The Topology defining the model to write
+         - positions (list) The list of atomic positions to write
+         - file (file=stdout) A file to write the model to
+         - modelIndex (int=None) If not None, the model will be surrounded by MODEL/ENDMDL records with this index
+        """
+        if len(list(topology.atoms())) != len(positions):
+            raise ValueError('The number of positions must match the number of atoms') 
+        if is_quantity(positions):
+            positions = positions.value_in_unit(angstroms)
+        atomIndex = 1
+        posIndex = 0
+        if modelIndex is not None:
+            print >>file, "MODEL     %4d" % modelIndex
+        for (chainIndex, chain) in enumerate(topology.chains()):
+            chainName = chr(ord('A')+chainIndex)
+            residues = list(chain.residues())
+            for (resIndex, res) in enumerate(residues):
+                if len(res.name) > 3:
+                    resName = res.name[:3]
+                else:
+                    resName = res.name
+                for atom in res.atoms():
+                    if len(atom.name) < 4:
+                        atomName = ' '+atom.name
+                    elif len(atom.name) > 4:
+                        atomName = atom.name[:4]
+                    else:
+                        atomName = atom.name
+                    coords = positions[posIndex]
+                    print >>file, "ATOM  %5d %-4s %3s %s%4d    %8.3f%8.3f%8.3f  1.00  0.00" % (atomIndex, atomName, resName, chainName, resIndex+1, coords[0], coords[1], coords[2])
+                    posIndex += 1
+                    atomIndex += 1
+                if resIndex == len(residues)-1:
+                    print >>file, "TER   %5d      %3s %s%4d" % (atomIndex, resName, chainName, resIndex+1)
+                    atomIndex += 1
+        if modelIndex is not None:
+            print >>file, "ENDMDL"
+
+    @staticmethod
+    def writeFooter(topology, file=sys.stdout):
+        """Write out the footer for a PDB file.
+        
+        Parameters:
+         - topology (Topology) The Topology defining the molecular system being written
+         - file (file=stdout) A file to write the file to
+        """
+        print >>file, "END"
+

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

+"""
+pdbreporter.py: Outputs simulation trajectories in PDB format
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import simtk.openmm as mm
+from simtk.openmm.app import PDBFile
+    
+class PDBReporter(object):
+    """PDBReporter outputs a series of frames from a Simulation to a PDB file.
+    
+    To use it, create a PDBReporter, than add it to the Simulation's list of reporters.
+    """
+    
+    def __init__(self, file, reportInterval):
+        """Create a PDBReporter.
+    
+        Parameters:
+         - file (string) The file to write to
+         - reportInterval (int) The interval (in time steps) at which to write frames
+        """
+        self._reportInterval = reportInterval
+        self._out = open(file, 'w')
+        self._topology = None
+        self._nextModel = 0
+    
+    def describeNextReport(self, simulation):
+        """Get information about the next report this object will generate.
+        
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+        Returns: A five element tuple.  The first element is the number of steps until the
+        next report.  The remaining elements specify whether that report will require
+        positions, velocities, forces, and energies respectively.
+        """
+        steps = self._reportInterval - simulation.currentStep%self._reportInterval
+        return (steps, True, False, False, False)
+    
+    def report(self, simulation, state):
+        """Generate a report.
+        
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+         - state (State) The current state of the simulation
+        """
+        if self._nextModel == 0:
+            PDBFile.writeHeader(simulation.topology, self._out)
+            self._topology = simulation.topology
+            self._nextModel += 1
+        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
+        self._nextModel += 1
+    
+    def __del__(self):
+        PDBFile.writeFooter(self._topology, self._out)
+        self._out.close()

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

+"""
+simulation.py: Provides a simplified API for running simulations.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import simtk.openmm as mm
+import simtk.unit as unit
+    
+class Simulation(object):
+    """Simulation provides a simplified API for running simulations with OpenMM and reporting results.
+    
+    A Simulation ties together various objects used for running a simulation: a Topology, System,
+    Integrator, and Context.  To use it, you provide the Topology, System, and Integrator, and it
+    creates the Context automatically.
+    
+    Simulation also maintains a list of "reporter" objects that record or analyze data as the simulation
+    runs, such as writing coordinates to files or displaying structures on the screen.  For example,
+    the following line will cause a file called "output.pdb" to be created, and a structure written to
+    it every 1000 time steps:
+    
+    simulation.reporters.append(PDBReporter('output.pdb', 1000))
+    """
+    
+    def __init__(self, topology, system, integrator, platform=None, platformProperties=None):
+        """Create a Simulation.
+        
+        Parameters:
+         - topology (Topology) A Topology describing the the system to simulation
+         - system (System) The OpenMM System object to simulate
+         - integrator (Integrator) The OpenMM Integrator to use for simulating the System
+         - platform (Platform=None) If not None, the OpenMM Platform to use
+         - platformProperties (map=None) If not None, a set of platform-specific properties to pass
+           to the Context's constructor
+        """
+        self.topology = topology
+        self.system = system
+        self.integrator = integrator
+        self.currentStep = 0
+        self.reporters = []
+        if platform is None:
+            self.context = mm.Context(system, integrator)
+        elif platformProperties is None:
+            self.context = mm.Context(system, integrator, platform)
+        else:
+            self.context = mm.Context(system, integrator, platform, platformProperties)
+    
+    def minimizeEnergy(self, tolerance=1*unit.kilojoule/unit.mole, maxIterations=0):
+        """Perform a local energy minimization on the system.
+        
+        Parameters:
+         - tolerance (energy=1*kilojoule/mole) The energy tolerance to which the system should be minimized
+         - maxIterations (int=0) The maximum number of iterations to perform.  If this is 0, minimization is continued
+           until the results converge without regard to how many iterations it takes.
+        """
+        mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations)
+       
+    def step(self, steps):
+        """Advance the simulation by integrating a specified number of time steps."""
+        stepTo = self.currentStep+steps
+        nextReport = [None]*len(self.reporters)
+        while self.currentStep < stepTo:
+            nextSteps = stepTo-self.currentStep
+            anyReport = False
+            for i, reporter in enumerate(self.reporters):
+                nextReport[i] = reporter.describeNextReport(self)
+                if nextReport[i][0] > 0 and nextReport[i][0] <= nextSteps:
+                    nextSteps = nextReport[i][0]
+                    anyReport = True
+            self.integrator.step(nextSteps)
+            self.currentStep += nextSteps
+            if anyReport:
+                getPositions = False
+                getVelocities = False
+                getForces = False
+                getEnergy = False
+                for reporter, next in zip(self.reporters, nextReport):
+                    if next[0] == nextSteps:
+                        if next[1]:
+                            getPositions = True
+                        if next[2]:
+                            getVelocities = True
+                        if next[3]:
+                            getForces = True
+                        if next[4]:
+                            getEnergy = True
+                state = self.context.getState(getPositions=getPositions, getVelocities=getVelocities, getForces=getForces, getEnergy=getEnergy, getParameters=True)
+                for reporter, next in zip(self.reporters, nextReport):
+                    if next[0] == nextSteps:
+                        reporter.report(self, state)
+

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

+"""
+topology.py: Used for storing topological information about a system.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import os
+import xml.etree.ElementTree as etree
+from simtk.unit import nanometers, sqrt
+
+class Topology(object):
+    """Topology stores the topological information about a system.
+    
+    The structure of a Topology object is similar to that of a PDB file.  It consists of a set of Chains
+    (often but not always corresponding to polymer chains).  Each Chain contains a set of Residues,
+    and each Residue contains a set of Atoms.  In addition, the Topology stores a list of which atom
+    pairs are bonded to each other, and the dimensions of the crystallographic unit cell.
+    
+    Atom and residue names should follow the PDB 3.0 nomenclature for all molecules for which one exists.
+    """
+    
+    _standardBonds = {}
+    
+    def __init__(self):
+        """Create a new Topology object"""
+        self._chains = []
+        self._numResidues = 0
+        self._numAtoms = 0
+        self._bonds = []
+        self._unitCellDimensions = None
+        
+    def addChain(self):
+        """Create a new Chain and add it to the Topology.
+        
+        Returns: the newly created Chain
+        """
+        chain = Chain(len(self._chains), self)
+        self._chains.append(chain)
+        return chain
+    
+    def addResidue(self, name, chain):
+        """Create a new Residue and add it to the Topology.
+        
+        Parameters:
+         - name (string) The name of the residue to add
+         - chain (Chain) The Chain to add it to
+        Returns: the newly created Residue
+        """
+        residue = Residue(name, self._numResidues, chain)
+        self._numResidues += 1
+        chain._residues.append(residue)
+        return residue
+    
+    def addAtom(self, name, element, residue):
+        """Create a new Atom and add it to the Topology.
+        
+        Parameters:
+         - name (string) The name of the atom to add
+         - element (Element) The element of the atom to add
+         - residue (Residue) The Residue to add it to
+        Returns: the newly created Atom
+        """
+        atom = Atom(name, element, self._numAtoms, residue)
+        self._numAtoms += 1
+        residue._atoms.append(atom)
+        return atom
+    
+    def addBond(self, atom1, atom2):
+        """Create a new bond and add it to the Topology.
+        
+        Parameters:
+         - atom1 (Atom) The first Atom connected by the bond
+         - atom2 (Atom) The second Atom connected by the bond
+        """
+        self._bonds.append((atom1, atom2))
+    
+    def chains(self):
+        """Iterate over all Chains in the Topology."""
+        return iter(self._chains)
+    
+    def residues(self):
+        """Iterate over all Residues in the Topology."""
+        for chain in self._chains:
+            for residue in chain._residues:
+                yield residue
+    
+    def atoms(self):
+        """Iterate over all Atoms in the Topology."""
+        for chain in self._chains:
+            for residue in chain._residues:
+                for atom in residue._atoms:
+                    yield atom
+    
+    def bonds(self):
+        """Iterate over all bonds (each represented as a tuple of two Atoms) in the Topology."""
+        return iter(self._bonds)
+        
+    def getUnitCellDimensions(self):
+        """Get the dimensions of the crystallographic unit cell.
+        
+        The return value may be None if this Topology does not represent a periodic structure.
+        """
+        return self._unitCellDimensions
+        
+    def setUnitCellDimensions(self, dimensions):
+        """Set the dimensions of the crystallographic unit cell."""
+        self._unitCellDimensions = dimensions
+    
+    def createStandardBonds(self):
+        """Create bonds based on the atom and residue names for all standard residue types."""
+        if len(Topology._standardBonds) == 0:
+            # Load the standard bond defitions.
+            
+            tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', 'residues.xml'))
+            for residue in tree.getroot().findall('Residue'):
+                bonds = []
+                Topology._standardBonds[residue.attrib['name']] = bonds
+                for bond in residue.findall('Bond'):
+                    bonds.append((bond.attrib['from'], bond.attrib['to']))
+        for chain in self._chains:
+            # First build a map of atom names to atoms.
+            
+            atomMaps = []
+            for residue in chain._residues:
+                atomMap = {}
+                atomMaps.append(atomMap)
+                for atom in residue._atoms:
+                    atomMap[atom.name] = atom
+            
+            # Loop over residues and construct bonds.
+            
+            for i in range(len(chain._residues)):
+                name = chain._residues[i].name
+                if name in Topology._standardBonds:
+                    for bond in Topology._standardBonds[name]:
+                        if bond[0].startswith('-') and i > 0:
+                            fromResidue = i-1
+                            fromAtom = bond[0][1:]
+                        elif bond[0].startswith('+') and i <len(chain._residues):
+                            fromResidue = i+1
+                            fromAtom = bond[0][1:]
+                        else:
+                            fromResidue = i
+                            fromAtom = bond[0]
+                        if bond[1].startswith('-') and i > 0:
+                            toResidue = i-1
+                            toAtom = bond[1][1:]
+                        elif bond[1].startswith('+') and i <len(chain._residues):
+                            toResidue = i+1
+                            toAtom = bond[1][1:]
+                        else:
+                            toResidue = i
+                            toAtom = bond[1]
+                        if fromAtom in atomMaps[fromResidue] and toAtom in atomMaps[toResidue]:
+                            self.addBond(atomMaps[fromResidue][fromAtom], atomMaps[toResidue][toAtom])
+                            
+    def createDisulfideBonds(self, positions):
+        """Identify disulfide bonds based on proximity and add them to the Topology.
+        
+        Parameters:
+         - positions (list) The list of atomic positions based on which to identify bonded atoms
+        """
+        def isCyx(res):
+            names = [atom.name for atom in res._atoms]
+            return 'SG' in names and 'HG' not in names
+        
+        cyx = [res for res in self.residues() if res.name == 'CYS' and isCyx(res)]
+        atomNames = [[atom.name for atom in res._atoms] for res in cyx]
+        for i in range(len(cyx)):
+            sg1 = cyx[i]._atoms[atomNames[i].index('SG')]
+            pos1 = positions[sg1.index]
+            for j in range(i):
+                sg2 = cyx[j]._atoms[atomNames[j].index('SG')]
+                pos2 = positions[sg2.index]
+                delta = [x-y for (x,y) in zip(pos1, pos2)]
+                distance = sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2])
+                if distance < 0.3*nanometers:
+                    self.addBond(sg1, sg2)
+
+class Chain(object):
+    def __init__(self, index, topology):
+        """Construct a new Chain.  You should call addChain() on the Topology instead of calling this directly."""
+        self.index = index
+        self.topology = topology
+        self._residues = []
+    
+    def residues(self):
+        """Iterate over all Residues in the Chain."""
+        return iter(self._residues)
+
+    def atoms(self):
+        """Iterate over all Atoms in the Chain."""
+        for residue in self._residues:
+            for atom in residue._atoms:
+                yield atom
+
+class Residue(object):
+    def __init__(self, name, index, chain):
+        """Construct a new Residue.  You should call addResidue() on the Topology instead of calling this directly."""
+        self.name = name
+        self.index = index
+        self.chain = chain
+        self._atoms = []
+    
+    def atoms(self):
+        """Iterate over all Atoms in the Residue."""
+        return iter(self._atoms)
+
+class Atom(object):
+    def __init__(self, name, element, index, residue):
+        """Construct a new Atom.  You should call addAtom() on the Topology instead of calling this directly."""
+        self.name = name
+        self.element = element
+        self.index = index
+        self.residue = residue
+