First pass at adding my CHARMM compatibility layer.

Currently gives exactly the same energy as the same files passed through chamber
(which, in turn, was validated to machine precision against CHARMM itself).

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

+"""
+This package contains code for reading CHARMM structure files for setting up a
+simulation with CHARMM; specifically PSF, PAR, RTF, and STR files
+
+    - PAR : Parameter file (PRM) -- this contains all of the force field
+            parameters (e.g., bond lengths and strengths) for all of the atom
+            types
+
+    - RTF : Residue Topology File -- this file contains the residue
+            connectivity tables as well as a definition of all of the atom
+            types. Also contains an internal coordinate representation of the
+            residues
+
+    - PSF : Protein Structure File -- this is the main file type in CHARMM
+            simulations that defines all of the residues in a system as well as
+            the atom types and connectivity between the atoms
+
+    - STR : Stream file -- Source of additional information and CHARMM commands
+            that can contain RTF and PAR information. Allows users to define
+            additional parameters without 'contaminating' the original force
+            field parameter files
+"""
+
+__all__ = ['psf', 'parameters', 'openmmloader']
+__authors__ = 'Jason Swails'
+__contributors__ = ''
+__license__ = 'GPL v.3'
+__date__ = 'Mar. 31, 2014'
+
+
+__private__ = ['topologyobjects', '_charmmfile']

wrappers/python/simtk/openmm/app/charmm/_charmmfile.py

+"""
+Provides a class for reading CHARMM-style files. The key component to these
+files is that the ! character is a comment character and everything after ! is
+ignored.
+
+Author: Jason M. Swails
+Contributors:
+Date: April 9, 2014
+"""
+from simtk.openmm.app.charmm.exceptions import CharmmFileError
+
+class CharmmFile(object):
+    """
+    A CHARMM file that recognizes the "!" character as a 'comment' token. It
+    can be iterated over and generally treated like a file object, but only
+    spits out strings that have been truncated at its first comment character.
+    
+    There is currently no way to recognize a ! as a _non_ comment character,
+    since allowing an escape character does not seem to be common practice and
+    would likely introduce negative performance implications.
+    """
+
+    def __init__(self, fname, mode='r'):
+        if mode not in ('r', 'w'):
+            raise ValueError('Cannot open CharmmFile with mode "%s"' % mode)
+        if mode == 'r':
+            self.status = 'OLD'
+        else:
+            self.status = 'NEW'
+        try:
+            self._handle = open(fname, mode)
+        except IOError, e:
+            raise CharmmFileError(str(e))
+        self.closed = False
+        self.line_number = 0
+
+    def write(self, *args, **kwargs):
+        return self._handle.write(*args, **kwargs)
+
+    def __iter__(self):
+        # Iterate over the file
+        for line in self._handle:
+            try:
+                idx = line.index('!')
+                end = '\n'
+            except ValueError:
+                # There is no comment...
+                idx = None
+                end = ''
+            yield line[:idx] + end
+
+    def readline(self):
+        self.line_number += 1
+        line = self._handle.readline()
+        try:
+            idx = line.index('!')
+            end = '\n'
+        except ValueError:
+            idx = None
+            end = ''
+        return line[:idx] + end
+
+    def readlines(self):
+        return [line for line in self]
+
+    def read(self):
+        return ''.join(self.readlines())
+
+    def close(self):
+        self._handle.close()
+        self.closed = True
+
+    def rewind(self):
+        """ Return to the beginning of the file """
+        self._handle.seek(0)
+
+    def __del__(self):
+        self.closed or self._handle.close()
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class CharmmStreamFile(object):
+    """
+    The stream file is broken down into sections of commands delimited by the
+    strings:
+        read <section> <options>
+        ....
+        ....
+        end
+    This object provides iterators over those sections and a file-like API for
+    dealing with the text.
+    """
+    def __init__(self, fname):
+        self.lines = CharmmFile(fname, 'r').readlines()
+        self.line_number = 0
+
+    def __iter__(self):
+        return iter(self.lines)
+
+    def rewind(self):
+        """ Return to the beginning of the file """
+        self.line_number = 0
+
+    def next_section(self):
+        """
+        Fast-forwards the file to the next CHARMM command section
+
+        Returns: (str, list)
+            - The first string is the line defining the section that's being
+              returned
+            - The list is a list of all lines contained in the section
+              excluding the "read <blah>" and "end" lines.
+
+        Notes:
+            The line pointer will be set to the line defining the 
+        """
+        lines = []
+        while self.line_number < len(self.lines):
+            line = self.lines[self.line_number].strip()
+            if line[:4].lower() == 'read':
+                title = line.strip()
+                self.line_number += 1
+                line = self.lines[self.line_number]
+                while line and not line.strip().lower().startswith('end'):
+                    lines.append(line)
+                    self.line_number += 1
+                    line = self.lines[self.line_number]
+                return title, lines
+            self.line_number += 1
+        # No sections left
+        return None, None
+
+    def __del__(self):
+        pass

wrappers/python/simtk/openmm/app/charmm/exceptions.py

+"""
+This package contains exceptions that may be raised by the CHARMM components of
+the OpenMM Application layer
+
+Author: Jason M. Swails
+Contributors:
+Date: April 9, 2014
+"""
+
+class CharmmPSFError(CharmmError):
+    """ If there is a problem parsing CHARMM PSF files """
+
+class SplitResidueWarning(CharmmWarning):
+    """ For if a residue with the same number but different names is split """
+
+class ResidueError(CharmmError):
+    """ For when there are problems defining a residue """
+
+class CharmmPSFWarning(CharmmWarning):
+    """ For non-fatal PSF parsing issues """
+
+class CharmmFileError(CharmmError):
+    """ If there is a problem parsing CHARMM files """
+
+class MissingParameter(CharmmError):
+    """ If a parameter is missing from a database """
+
+class CmapError(CharmmError):
+    """ For an error arising from CMAP grid incompatibilities """
+
+class BondError(CharmmError):
+    """ Prevent an atom from bonding to itself """

wrappers/python/simtk/openmm/app/charmm/parameters.py

+"""
+This module contains classes for parsing and processing CHARMM parameter,
+topology, and stream files. It only extracts atom properties from the
+topology files and extracts all parameters from the parameter files
+
+Author: Jason M. Swails
+Contributors:
+Date: April 5, 2014
+"""
+from simtk.openmm.app.charmm._charmmfile import CharmmFile, CharmmStreamFile
+from simtk.openmm.app.charmm.topologyobjects import (AtomType, BondType,
+                    AngleType, DihedralType, ImproperType, CmapType,
+                    UreyBradleyType, NoUreyBradley)
+from simtk.openmm.app.charmm.exceptions import CharmmFileError
+from simtk.openmm.element import Element, get_by_symbol
+from simtk.units import dalton
+import warnings
+
+class ParameterSet(object):
+    """
+    Stores a parameter set defined by CHARMM files. It stores the equivalent of
+    the information found in the MASS section of the CHARMM topology file
+    (TOP/RTF) and all of the information in the parameter files (PAR)
+
+    Attributes:
+        All type lists are dictionaries whose keys are tuples (with however
+        many elements are needed to define that type of parameter). The types
+        that can be in any order are SORTED.
+        
+        - atom_types_str
+        - atom_types_int
+        - atom_types_tuple
+        - bond_types
+        - angle_types
+        - urey_bradley_types
+        - dihedral_types
+        - improper_types
+        - cmap_types
+        - nbfix_types
+
+        The dihedral types can be multiterm, so the values for each dict key is
+        actually a list of DihedralType instances. The atom_types are dicts that
+        match the name (str), number (int), or (name, number) tuple (tuple) to
+        the atom type. The tuple is guaranteed to be the most robust, although
+        when only the integer or string is available the other dictionaries are
+        helpful
+    """
+
+    @staticmethod
+    def _convert(data, type, msg=''):
+        """
+        Converts a data type to a desired type, raising CharmmFileError if it
+        fails
+        """
+        try:
+            return type(data)
+        except ValueError:
+            raise CharmmFileError('Could not convert %s to %s' % (msg, type))
+
+    def __init__(self):
+        # Instantiate the list types
+        self.atom_types_str = dict()
+        self.atom_types_int = dict()
+        self.atom_types_tuple = dict()
+        self.bond_types = dict()
+        self.angle_types = dict()
+        self.urey_bradley_types = dict()
+        self.dihedral_types = dict()
+        self.improper_types = dict()
+        self.cmap_types = dict()
+        self.nbfix_types = dict()
+        self.parametersets = []
+
+    @classmethod
+    def load_set(cls, tfile=None, pfile=None, sfiles=[]):
+        """
+        Instantiates a ParameterSet from a Topology file and a Parameter file
+        (or just a Parameter file if it has all information)
+
+        Parameters:
+            - tfile (str) : Name of the Topology (RTF/TOP) file
+            - pfile (str) : Name of the Parameter (PAR) file
+            - sfiles (list of str) : List or tuple of stream (STR) file names.
+
+        Returns:
+            New ParameterSet populated with the parameters found in the
+            provided files.
+            
+        Notes:
+            The RTF file is read first (if provided), followed by the PAR file,
+            followed by the list of stream files (in the order they are
+            provided). Parameters in each stream file will overwrite those that
+            came before (or simply append to the existing set if they are
+            different)
+        """
+        inst = cls()
+        if tfile is not None:
+            inst.read_topology_file(tfile)
+        if pfile is not None:
+            inst.read_parameter_file(pfile)
+        if isinstance(sfiles, str):
+            # The API docstring requests a list, but allow for users to pass a
+            # string with a single filename instead
+            inst.read_stream_file(sfiles)
+        elif sfiles is not None:
+            for sfile in sfiles:
+                inst.read_stream_file(sfile)
+        return inst
+
+    def read_parameter_file(self, pfile):
+        """
+        Reads all of the parameters from a parameter file. Versions 36 and
+        later of the CHARMM force field files have an ATOMS section defining
+        all of the atom types.  Older versions need to load this information
+        from the RTF/TOP files.
+
+        Parameters:
+            - pfile (str) : Name of the CHARMM PARameter file to read
+
+        Notes: The atom types must all be loaded by the end of this routine.
+        Either supply a PAR file with atom definitions in them or read in a
+        RTF/TOP file first. Failure to do so will result in a raised
+        RuntimeError.
+        """
+        conv = ParameterSet._convert
+        if isinstance(pfile, str):
+            own_handle = True
+            f = CharmmFile(pfile)
+        else:
+            own_handle = False
+            f = pfile
+        # What section are we parsing?
+        section = None
+        # The current cmap we are building (these span multiple lines)
+        current_cmap = None
+        current_cmap_data = []
+        current_cmap_res = 0
+        nonbonded_types = dict() # Holder
+        needs_blank = False
+        parameterset = None
+        for line in f:
+            line = line.strip()
+            if not line:
+                # This is a blank line
+                needs_blank = False
+                continue
+            if parameterset is None and line.strip().startswith('*>>'):
+                parameterset = line.strip()[1:78]
+                continue
+            # Set section if this is a section header
+            if line.startswith('ATOMS'):
+                section = 'ATOMS'
+                continue
+            if line.startswith('BONDS'):
+                section = 'BONDS'
+                continue
+            if line.startswith('ANGLES'):
+                section = 'ANGLES'
+                continue
+            if line.startswith('DIHEDRALS'):
+                section = 'DIHEDRALS'
+                continue
+            if line.startswith('IMPROPER'):
+                section = 'IMPROPER'
+                continue
+            if line.startswith('CMAP'):
+                section = 'CMAP'
+                continue
+            if line.startswith('NONBONDED'):
+                section = 'NONBONDED'
+                needs_blank = True
+                continue
+            if line.startswith('NBFIX'):
+                section = 'NBFIX'
+                continue
+            if line.startswith('HBOND'):
+                section = None
+                continue
+            # If we have no section, skip
+            if section is None: continue
+            # Now handle each section specifically
+            if section == 'ATOMS':
+                if not line.startswith('MASS'): continue # Should this happen?
+                words = line.split()
+                try:
+                    idx = conv(words[1], int, 'atom type')
+                    name = words[2]
+                    mass = conv(words[3], float, 'atom mass')
+                except IndexError:
+                    raise CharmmFileError('Could not parse MASS section.')
+                # The parameter file might or might not have an element name
+                try:
+                    elem = words[4]
+                    atomic_number = get_by_symbol(elem).atomic_number
+                except (IndexError, KeyError):
+                    # Figure it out from the mass
+                    atomic_number = element_by_mass(mass).atomic_number
+                atype = AtomType(name=name, number=idx, mass=mass,
+                                 atomic_number=atomic_number)
+                self.atom_types_str[atype.name] = atype
+                self.atom_types_int[atype.number] = atype
+                self.atom_types_tuple[(atype.name, atype.number)] = atype
+                continue
+            if section == 'BONDS':
+                words = line.split()
+                try:
+                    type1 = words[0]
+                    type2 = words[1]
+                    k = conv(words[2], float, 'bond force constant')
+                    req = conv(words[3], float, 'bond equilibrium dist')
+                except IndexError:
+                    raise CharmmFileError('Could not parse bonds.')
+                key = (min(type1, type2), max(type1, type2))
+                self.bond_types[key] = BondType(k, req)
+                continue
+            if section == 'ANGLES':
+                words = line.split()
+                try:
+                    type1 = words[0]
+                    type2 = words[1]
+                    type3 = words[2]
+                    k = conv(words[3], float, 'angle force constant')
+                    theteq = conv(words[4], float, 'angle equilibrium value')
+                except IndexError:
+                    raise CharmmFileError('Could not parse angles.')
+                key = (min(type1, type3), type2, max(type1, type3))
+                self.angle_types[key] = AngleType(k, theteq)
+                # See if we have a urey-bradley
+                try:
+                    ubk = conv(words[5], float, 'Urey-Bradley force constant')
+                    ubeq = conv(words[6], float, 'Urey-Bradley equil. value')
+                    ubtype = UreyBradleyType(ubk, ubeq)
+                except IndexError:
+                    ubtype = NoUreyBradley
+                self.urey_bradley_types[key] = ubtype
+                continue
+            if section == 'DIHEDRALS':
+                words = line.split()
+                try:
+                    type1 = words[0]
+                    type2 = words[1]
+                    type3 = words[2]
+                    type4 = words[3]
+                    k = conv(words[4], float, 'dihedral force constant')
+                    n = conv(words[5], float, 'dihedral periodicity')
+                    phase = conv(words[6], float, 'dihedral phase')
+                except IndexError:
+                    raise CharmmFileError('Could not parse dihedrals.')
+                # Torsion can be in either direction. Sort by end groups first,
+                # then sort by middle 2
+                if type1 < type4:
+                    key = (type1, type2, type3, type4)
+                elif type1 > type4:
+                    key = (type4, type3, type2, type1)
+                else:
+                    # OK, we need to sort by the middle atoms now
+                    if type2 < type3:
+                        key = (type1, type2, type3, type4)
+                    else:
+                        key = (type4, type3, type2, type1)
+                # See if this is a second (or more) term of the dihedral group
+                # that's already present.
+                dihedral = DihedralType(k, n, phase)
+                if key in self.dihedral_types:
+                    # See if the existing dihedral type list has a term with
+                    # the same periodicity -- If so, replace it
+                    replaced = False
+                    for i, dtype in enumerate(self.dihedral_types[key]):
+                        if dtype.per == dihedral.per:
+                            # Replace. Should we warn here?
+                            warnings.warn('Replacing dihedral %r with %r' % 
+                                          (dtype, dihedral))
+                            self.dihedral_types[key]
+                            replaced = True
+                            break
+                    if not replaced:
+                        self.dihedral_types[key].append(dihedral)
+                else: # key not present
+                    self.dihedral_types[key] = [dihedral]
+                continue
+            if section == 'IMPROPER':
+                words = line.split()
+                try:
+                    type1 = words[0]
+                    type2 = words[1]
+                    type3 = words[2]
+                    type4 = words[3]
+                    k = conv(words[4], float, 'improper force constant')
+                    theteq = conv(words[5], float, 'improper equil. value')
+                except IndexError:
+                    raise CharmmFileError('Could not parse dihedrals.')
+                # If we have a 7th column, that is the real psi0 (and the 6th
+                # is just a dummy 0)
+                try:
+                    tmp = conv(words[6], float, 'improper equil. value')
+                    theteq = tmp
+                except IndexError:
+                    pass # Do nothing
+                # Improper types seem not to have the central atom defined in
+                # the first place, so just have the key a fully sorted list. We
+                # still depend on the PSF having properly ordered improper atoms
+                key = tuple(sorted([type1, type2, type3, type4]))
+                self.improper_types[key] = ImproperType(k, theteq)
+                continue
+            if section == 'CMAP':
+                # This is the most complicated part, since cmap parameters span
+                # many lines. We won't do much error catching here.
+                words = line.split()
+                try:
+                    holder = [float(w) for w in words]
+                    current_cmap_data.extend(holder)
+                except ValueError:
+                    # We assume this is a definition of a new CMAP, so
+                    # terminate the last CMAP if applicable
+                    if current_cmap is not None:
+                        # We have a map to terminate
+                        ty = CmapType(current_cmap_res, current_cmap_data)
+                        self.cmap_types[current_cmap] = ty
+                    try:
+                        type1 = words[0]
+                        type2 = words[1]
+                        type3 = words[2]
+                        type4 = words[3]
+                        type5 = words[4]
+                        type6 = words[5]
+                        type7 = words[6]
+                        type8 = words[7]
+                        res = conv(words[8], int, 'CMAP resolution')
+                    except IndexError:
+                        raise CharmmFileError('Could not parse CMAP data.')
+                    # order the torsions independently
+                    k1 = [type1, type2, type3, type4]
+                    k2 = [type4, type3, type2, type1]
+                    key1 = min(k1, k2)
+                    k1 = [type5, type6, type7, type8]
+                    k2 = [type8, type7, type6, type5]
+                    key2 = min(k1, k2)
+                    current_cmap = tuple(key1 + key2)
+                    current_cmap_res = res
+                    current_cmap_data = []
+                continue
+            if section == 'NONBONDED' and not needs_blank:
+                # Now get the nonbonded values
+                words = line.split()
+                try:
+                    atype = words[0]
+                    # 1st column is ignored
+                    epsilon = conv(words[2], float, 'vdW epsilon term')
+                    rmin = conv(words[3], float, 'vdW Rmin/2 term')
+                except IndexError:
+                    raise CharmmFileError('Could not parse nonbonded terms.')
+                # See if we have 1-4 parameters
+                try:
+                    # 4th column is ignored
+                    eps14 = conv(words[5], float, '1-4 vdW epsilon term')
+                    rmin14 = conv(words[6], float, '1-4 vdW Rmin/2 term')
+                except IndexError:
+                    eps14 = rmin14 = None
+                nonbonded_types[atype] = [epsilon, rmin, eps14, rmin14]
+                continue
+            if section == 'NBFIX':
+                words = line.split()
+                try:
+                    at1 = words[0]
+                    at2 = words[1]
+                    emin = conv(words[2], float, 'NBFIX Emin')
+                    rmin = conv(words[3], float, 'NBFIX Rmin')
+                except IndexError:
+                    raise CharmmFileError('Could not parse NBFIX terms.')
+                self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)
+        # Now we're done. Load the nonbonded types into the relevant AtomType
+        # instances. In order for this to work, all keys in nonbonded_types
+        # must be in the self.atom_types_str dict. Raise a RuntimeError if this
+        # is not satisfied
+        try:
+            for key in nonbonded_types:
+                self.atom_types_str[key].set_lj_params(*nonbonded_types[key])
+        except KeyError:
+            raise RuntimeError('Atom type %s not present in AtomType list' %
+                               key)
+        if parameterset is not None: self.parametersets.append(parameterset)
+        if own_handle: f.close()
+
+    def read_topology_file(self, tfile):
+        """
+        Reads _only_ the atom type definitions from a topology file. This is
+        unnecessary for versions 36 and later of the CHARMM force field.
+
+        Parameters:
+            - tfile (str) : Name of the CHARMM TOPology file to read
+
+        Note: The CHARMM TOPology file is also called a Residue Topology File
+        """
+        conv = ParameterSet._convert
+        if isinstance(tfile, str):
+            own_handle = True
+            f = CharmmFile(tfile)
+        else:
+            own_handle = False
+            f = tfile
+        for line in f:
+            line = line.strip()
+            if line[:4] != 'MASS': continue
+            words = line.split()
+            try:
+                idx = conv(words[1], int, 'atom type')
+                name = words[2]
+                mass = conv(words[3], float, 'atom mass')
+            except IndexError:
+                raise CharmmFileError('Could not parse MASS section of %s' %
+                                      tfile)
+            # The parameter file might or might not have an element name
+            try:
+                elem = words[4]
+                atomic_number = get_by_symbol(elem).atomic_number
+            except (IndexError, KeyError):
+                # Figure it out from the mass
+                atomic_number = element_by_mass(mass).atomic_number
+            atype = AtomType(name=name, number=idx, mass=mass,
+                             atomic_number=atomic_number)
+            self.atom_types_str[atype.name] = atype
+            self.atom_types_int[atype.number] = atype
+            self.atom_types_tuple[(atype.name, atype.number)] = atype
+        if own_handle: f.close()
+
+    def read_stream_file(self, sfile):
+        """
+        Reads RTF and PAR sections from a stream file and dispatches the
+        sections to read_topology_file or read_parameter_file
+
+        Parameters:
+            - sfile (str or CharmmStreamFile) : Stream file to parse
+        """
+        if isinstance(sfile, CharmmStreamFile):
+            f = sfile
+        else:
+            f = CharmmStreamFile(sfile)
+
+        title, section = f.next_section()
+        while title is not None and section is not None:
+            words = title.lower().split()
+            if words[1] == 'rtf':
+                # This is a Residue Topology File section.
+                self.read_topology_file(section)
+            elif words[1].startswith('para'):
+                # This is a Parameter file section
+                self.read_parameter_file(section)
+            title, section = f.next_section()
+
+    def condense(self):
+        """
+        This function goes through each of the parameter type dicts and
+        eliminates duplicate types. After calling this function, every unique
+        bond, angle, dihedral, improper, or cmap type will pair with EVERY key
+        in the type mapping dictionaries that points to the equivalent type
+
+        Returns:
+            - Returns the instance that is being condensed.
+
+        Notes:
+            The return value allows you to condense the types at construction
+            time.
+
+        Example:
+        >>> params = ParameterSet.load_set(pfile='charmm.prm').condense()
+        """
+        # First scan through all of the bond types
+        self._condense_types(self.bond_types)
+        self._condense_types(self.angle_types)
+        self._condense_types(self.urey_bradley_types)
+        self._condense_types(self.improper_types)
+        self._condense_types(self.cmap_types)
+        # Dihedrals have to be handled separately, since each key corresponds to
+        # a list of (potentially multiterm) dihedral terms. Since all terms in a
+        # multiterm dihedral have to have a DIFFERENT periodicity, we don't have
+        # to condense _within_ a single list of torsions assigned to the same
+        # key (they're guaranteed to be different)
+        keylist = self.dihedral_types.keys()
+        for i in range(len(keylist) - 1):
+            key1 = keylist[i]
+            for dihedral in self.dihedral_types[key1]:
+                for j in range(i+1, len(keylist)):
+                    key2 = keylist[j]
+                    for jj, dihedral2 in enumerate(self.dihedral_types[key2]):
+                        if dihedral2 == dihedral:
+                            self.dihedral_types[key2][jj] = dihedral
+        return self
+
+    @staticmethod
+    def _condense_types(typedict):
+        """
+        Loops through the given dict and condenses all types.
+
+        Parameter:
+            - typedict : Type dictionary to condense
+        """
+        keylist = typedict.keys()
+        for i in range(len(keylist) - 1):
+            key1 = keylist[i]
+            for j in range(i+1, len(keylist)):
+                key2 = keylist[j]
+                if typedict[key1] == typedict[key2]:
+                    typedict[key2] = typedict[key1]
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+def element_by_mass(mass):
+    """ Determines what element the given atom is based on its mass """
+
+    diff = mass
+    best_guess = 'EP'
+
+    for key in Element._elements_by_atomic_number:
+        element = Element._elements_by_atomic_number[key]
+        massdiff = abs(element.mass.value_in_unit(u.daltons) - mass)
+        if massdiff < diff:
+            best_guess = element
+            diff = massdiff
+
+    return best_guess

wrappers/python/simtk/openmm/app/charmm/psf.py

+"""
+Provides a Python class for parsing a PSF file and setting up a system
+structure for it within the OpenMM framework.
+
+Author: Jason M. Swails
+Contributors:
+Date: April 9, 2014
+"""
+from __future__ import division
+
+from functools import wraps
+from simtk.openmm.app.charmm._charmmfile import CharmmFile
+from simtk.openmm.app.charmm.topologyobjects import (ResidueList, AtomList,
+                TrackedList, Bond, Angle, Dihedral, Improper, AcceptorDonor,
+                Group, Cmap, UreyBradley, NoUreyBradley)
+from simtk.openmm.app.charmm.exceptions import (CharmmPSFError, MoleculeError,
+                CharmmPSFWarning, MissingParameter)
+import os
+import warnings
+import simtk.openmm as mm
+from simtk.openmm.vec3 import Vec3
+import simtk.unit as u
+from simtk.openmm.app import (forcefield as ff, Topology, element)
+from simtk.openmm.app.amberprmtopfile import HCT, OBC1, OBC2, GBn, GBn2
+from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce,
+                GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force)
+
+TINY = 1e-8
+WATNAMES = ('WAT', 'HOH', 'TIP3', 'TIP4', 'TIP5', 'SPCE', 'SPC')
+
+
+def _catchindexerror(func):
+    """
+    Protects a function from raising an index error, and replace that exception
+    with a CharmmPSFError instead
+    """
+    @wraps(func)
+    def newfunc(*args, **kwargs):
+        """ Catch the index error """
+        try:
+            return func(*args, **kwargs)
+        except IndexError, e:
+            raise CharmmPSFError('Array is too short: %s' % e)
+
+    return newfunc
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class ProteinStructure(object):
+    """
+    A chemical structure instantiated from CHARMM files. You can instantiate a
+    ProteinStructure from a PSF file using the load_from_psf constructor
+
+    Example:
+    >>> cs = ProteinStructure.load_from_psf("testfiles/test.psf")
+    
+    This structure has numerous attributes that are lists of the elements of
+    this structure, including atoms, bonds, torsions, etc. The attributes are
+        - residue_list
+        - atom_list
+        - bond_list
+        - angle_list
+        - dihedral_list
+        - improper_list
+        - cmap_list
+        - donor_list    # hbonds donors?
+        - acceptor_list # hbond acceptors?
+        - group_list    # list of nonbonded interaction groups
+
+    Additional attribute is available if a ParameterSet is loaded into this
+    structure.
+        
+        - urey_bradley_list
+
+    The lengths of each of these lists gives the pointers (e.g., natom, nres,
+    etc.)
+
+    Example:
+    >>> cs = ProteinStructure.load_from_psf("testfiles/test.psf")
+    >>> len(cs.atom_list)
+    33
+    >>> len(cs.bond_list)
+    32
+    """
+    # Define default force groups for all of the bonded terms. This allows them
+    # to be turned on and off selectively. This is a way to implement per-term
+    # energy decomposition to compare individual components
+
+    BOND_FORCE_GROUP = 0
+    ANGLE_FORCE_GROUP = 1
+    DIHEDRAL_FORCE_GROUP = 2
+    UREY_BRADLEY_FORCE_GROUP = 3
+    IMPROPER_FORCE_GROUP = 4
+    CMAP_FORCE_GROUP = 5
+    NONBONDED_FORCE_GROUP = 6
+    GB_FORCE_GROUP = 6
+    
+    def __init__(self, residues, atoms, bonds=None, angles=None,
+                 dihedrals=None, impropers=None, cmaps=None, donors=None,
+                 acceptors=None, groups=None, title=None, flags=None):
+        self.residue_list = residues
+        self.atom_list = atoms
+        self.bond_list = bonds
+        self.angle_list = angles
+        self.dihedral_list = dihedrals
+        self.improper_list = impropers
+        self.cmap_list = cmaps
+        self.donor_list = donors
+        self.acceptor_list = acceptors
+        self.group_list = groups
+        self.title = title
+        self.flags = flags
+
+    @staticmethod
+    def _convert(string, type, message):
+        """
+        Converts a string to a specific type, making sure to raise
+        CharmmPSFError with the given message in the event of a failure.
+
+        Parameters:
+            - string (str) : Input string to process
+            - type (type) : Type of data to convert to
+            - message (str) : Error message to put in exception if failed
+        """
+        try:
+            return type(string)
+        except ValueError, e:
+            print e
+            raise CharmmPSFError('Could not convert %s' % message)
+
+    @staticmethod
+    def _parse_psf_section(psf, dtype):
+        """
+        This method parses a section of the PSF file
+
+        Parameters:
+            - psf (CharmmFile) : Open file that is pointing to the first line
+                                 of the section that is to be parsed
+            - dtype (type) : The data type to convert all of the data into
+        
+        Returns:
+            (pointers, data)
+            
+            - pointers (int/tuple of ints) : If one pointer is set, pointers is
+                    simply the integer that is value of that pointer. Otherwise
+                    it is a tuple with every pointer value defined in the first
+                    line
+            - data (list) : A list of all data in the parsed section converted
+                    to `dtype'
+        """
+        conv = ProteinStructure._convert
+        words = psf.readline().split()
+        if len(words) == 1:
+            pointers = conv(words[0], int, 'pointer')
+        else:
+            pointers = tuple([conv(w, int, 'pointer') for w in words])
+        line = psf.readline().strip()
+        if not line:
+            # This will correctly handle the NNB section (which has a spurious
+            # blank line) as well as any sections that have 0 members.
+            line = psf.readline().strip()
+        data = []
+        while line:
+            words = line.split()
+            data.extend([conv(w, dtype, 'PSF data') for w in words])
+            line = psf.readline().strip()
+        return pointers, data
+
+    @classmethod
+    @_catchindexerror
+    def load_from_psf(cls, psf_name):
+        """
+        Opens and parses a PSF file, then instantiates a ProteinStructure
+        instance from the data.
+            
+        Parameters:
+            psf_name (str) : Name of the PSF file (it must exist)
+        
+        Exceptions Raised:
+            IOError : If file "psf_name" does not exist
+            CharmmPSFError: If any parsing errors are encountered
+        """
+        conv = ProteinStructure._convert
+        # Make sure the file exists
+        if not os.path.exists(psf_name):
+            raise IOError('Could not find PSF file %s' % psf_name)
+        # Open the PSF and read the first line. It must start with "PSF"
+        psf = CharmmFile(psf_name, 'r')
+        line = psf.readline()
+        if not line.startswith('PSF'):
+            raise CharmmPSFError('Unrecognized PSF file. First line is %s' %
+                                 line.strip())
+        # Store the flags
+        psf_flags = line.split()[1:]
+        # Next line is blank
+        psf.readline()
+        # The next line has one number -- the number of title lines
+        ntitle = conv(psf.readline().strip(), int, 'title count')
+        # store the title
+        title = list()
+        for i in range(ntitle):
+            title.append(psf.readline().rstrip())
+        # Skip the blank line
+        psf.readline()
+        # Next is the number of atoms
+        natom = conv(psf.readline().strip(), int, 'natom')
+        # Parse all of the atoms
+        residue_list = ResidueList()
+        atom_list = AtomList()
+        for i in xrange(natom):
+            words = psf.readline().split()
+            atid = int(words[0])
+            if atid != i + 1:
+                raise CharmmPSFError('Nonsequential atoms detected!')
+            system = words[1]
+            resid = conv(words[2], int, 'residue number')
+            resname = words[3]
+            name = words[4]
+            attype = words[5]
+            # Try to convert the atom type to an integer a la CHARMM
+            try:
+                attype = int(attype)
+            except ValueError:
+                pass
+            charge = conv(words[6], float, 'partial charge')
+            mass = conv(words[7], float, 'atomic mass')
+            props = words[8:]
+            atom = residue_list.add_atom(system, resid, resname, name,
+                                         attype, charge, mass, props)
+            atom_list.append(atom)
+        atom_list.assign_indexes()
+        # Eat the next line
+        psf.readline()
+        # Now get the number of bonds
+        nbond, holder = ProteinStructure._parse_psf_section(psf, int)
+        bond_list = TrackedList()
+        if len(holder) != nbond * 2:
+            raise CharmmPSFError('Got %d indexes for %d bonds' %
+                                 (len(holder), nbond))
+        for i in range(nbond):
+            id1 = holder[2*i  ] - 1
+            id2 = holder[2*i+1] - 1
+            bond_list.append(Bond(atom_list[id1], atom_list[id2]))
+        bond_list.changed = False
+        # Now get the number of angles and the angle list
+        ntheta, holder = ProteinStructure._parse_psf_section(psf, int)
+        angle_list = TrackedList()
+        if len(holder) != ntheta * 3:
+            raise CharmmPSFError('Got %d indexes for %d angles' %
+                                 (len(holder), ntheta))
+        for i in range(ntheta):
+            id1 = holder[3*i  ] - 1
+            id2 = holder[3*i+1] - 1
+            id3 = holder[3*i+2] - 1
+            angle_list.append(
+                    Angle(atom_list[id1], atom_list[id2], atom_list[id3])
+            )
+        angle_list.changed = False
+        # Now get the number of torsions and the torsion list
+        nphi, holder = ProteinStructure._parse_psf_section(psf, int)
+        dihedral_list = TrackedList()
+        if len(holder) != nphi * 4:
+            raise CharmmPSFError('Got %d indexes for %d torsions' %
+                                 (len(holder), nphi))
+        for i in range(nphi):
+            id1 = holder[4*i  ] - 1
+            id2 = holder[4*i+1] - 1
+            id3 = holder[4*i+2] - 1
+            id4 = holder[4*i+3] - 1
+            dihedral_list.append(
+                    Dihedral(atom_list[id1], atom_list[id2], atom_list[id3],
+                             atom_list[id4])
+            )
+        dihedral_list.changed = False
+        # Now get the number of improper torsions
+        nimphi, holder = ProteinStructure._parse_psf_section(psf, int)
+        improper_list = TrackedList()
+        if len(holder) != nimphi * 4:
+            raise CharmmPSFError('Got %d indexes for %d impropers' %
+                                 (len(holder), nimphi))
+        for i in range(nimphi):
+            id1 = holder[4*i  ] - 1
+            id2 = holder[4*i+1] - 1
+            id3 = holder[4*i+2] - 1
+            id4 = holder[4*i+3] - 1
+            improper_list.append(
+                    Improper(atom_list[id1], atom_list[id2], atom_list[id3],
+                             atom_list[id4])
+            )
+        improper_list.changed = False
+        # Now handle the donors (what is this used for??)
+        ndon, holder = ProteinStructure._parse_psf_section(psf, int)
+        donor_list = TrackedList()
+        if len(holder) != ndon * 2:
+            raise CharmmPSFError('Got %d indexes for %d donors' %
+                                 (len(holder), ndon))
+        for i in range(ndon):
+            id1 = holder[2*i  ] - 1
+            id2 = holder[2*i+1] - 1
+            donor_list.append(AcceptorDonor(atom_list[id1], atom_list[id2]))
+        donor_list.changed = False
+        # Now handle the acceptors (what is this used for??)
+        nacc, holder = ProteinStructure._parse_psf_section(psf, int)
+        acceptor_list = TrackedList()
+        if len(holder) != nacc * 2:
+            raise CharmmPSFError('Got %d indexes for %d acceptors' %
+                                 (len(holder), ndon))
+        for i in range(nacc):
+            id1 = holder[2*i  ] - 1
+            id2 = holder[2*i+1] - 1
+            acceptor_list.append(AcceptorDonor(atom_list[id1], atom_list[id2]))
+        acceptor_list.changed = False
+        # Now get the NNB section. Not sure what this section is for or what it
+        # does...
+        nnb, holder = ProteinStructure._parse_psf_section(psf, int)
+        # Now get the group sections
+        pointers, holder = ProteinStructure._parse_psf_section(psf, int)
+        group_list = TrackedList()
+        try:
+            ngrp, nst2 = pointers
+        except ValueError:
+            raise CharmmPSFError('Could not unpack GROUP pointers')
+        group_list.nst2 = nst2
+        # Now handle the groups
+        if len(holder) != ngrp * 3:
+            raise CharmmPSFError('Got %d indexes for %d groups' %
+                                 (len(holder), ngrp))
+        for i in range(ngrp):
+            i1 = holder[3*i  ]
+            i2 = holder[3*i+1]
+            i3 = holder[3*i+2]
+            group_list.append(Group(i1, i2, i3))
+        group_list.changed = False
+        # The next section might be the number of molecules or it might be the
+        # cross-term (cmap) section. The first thing we'll do is determine
+        # molecularity based on the atom connectivity. If every PSF file was
+        # guaranteed to be "correct", we could just compare the MOLNT
+        # section with the one we compute here. However, CHARMM GUI appears
+        # to assign MOLNT as a dummy section (with all 1's), so this
+        # approach will not work. Instead, look at the value of the pointer
+        # and the number of entries in the group. If the # of entries is
+        # NATOM, assume we have MOLNT section. Warn if the MOLNT section is
+        # 'wrong'...
+        pointer, holder = ProteinStructure._parse_psf_section(psf, int)
+
+        # Assign all of the atoms to molecules recursively
+        set_molecules(atom_list)
+        molecule_list = [atom.marked for atom in atom_list]
+        if len(holder) == len(atom_list):
+            if molecule_list != holder:
+                warnings.warn('Detected PSF molecule section that is WRONG. '
+                              'Resetting molecularity.', CharmmPSFWarning)
+            # We have a CHARMM PSF file; now do NUMLP/NUMLPH sections
+            words = psf.readline().split()
+            numlp = conv(words[0], int, 'numlp')
+            numlph = conv(words[1], int, 'numlph')
+            if numlp != 0 or numlph != 0:
+                raise NotImplemented('Cannot currently handle PSFs with lone '
+                                     'pairs defined in the NUMLP/NUMLPH '
+                                     'section.')
+            psf.readline() # blank
+            # Now we get to the cross-term section
+            ncrterm, holder = ProteinStructure._parse_psf_section(psf, int)
+        else:
+            ncrterm = pointer
+        # At this point, ncrterm and holder are both set to the CMAP list for
+        # VMD and non-VMD PSFs.
+        # Now get the cmaps
+        cmap_list = TrackedList()
+        if len(holder) != ncrterm * 8:
+            raise CharmmPSFError('Got %d CMAP indexes for %d cmap terms' %
+                                 (len(holder), ncrterm))
+        for i in range(ncrterm):
+            id1 = holder[8*i  ] - 1
+            id2 = holder[8*i+1] - 1
+            id3 = holder[8*i+2] - 1
+            id4 = holder[8*i+3] - 1
+            id5 = holder[8*i+4] - 1
+            id6 = holder[8*i+5] - 1
+            id7 = holder[8*i+6] - 1
+            id8 = holder[8*i+7] - 1
+            cmap_list.append(
+                    Cmap(atom_list[id1], atom_list[id2], atom_list[id3],
+                         atom_list[id4], atom_list[id5], atom_list[id6],
+                         atom_list[id7], atom_list[id8])
+            )
+        cmap_list.changed = False
+        # Now instantiate the object and return it
+        inst = cls(residues=residue_list, atoms=atom_list, bonds=bond_list,
+                   angles=angle_list, dihedrals=dihedral_list,
+                   impropers=improper_list, donors=donor_list,
+                   acceptors=acceptor_list, groups=group_list, cmaps=cmap_list,
+                   title=title, flags=psf_flags)
+        # Assign some potentially useful attributes
+        return inst
+
+    def write_psf(self, dest, vmd=False):
+        """
+        Writes a PSF file from the stored molecule
+
+        Parameters:
+            - dest (str or file-like) : The place to write the output PSF file.
+                    If it has a "write" attribute, it will be used to print the
+                    PSF file. Otherwise, it will be treated like a string and a
+                    file will be opened, printed, then closed
+            - vmd (bool) : If True, it will write out a PSF in the format that
+                    VMD prints it in (i.e., no NUMLP/NUMLPH or MOLNT sections)
+        Example:
+            >>> cs = ProteinStructure.load_from_psf('testfiles/test.psf')
+            >>> cs.write_psf('testfiles/test2.psf')
+        """
+        # See if this is an extended format
+        ext = 'EXT' in self.flags
+        own_handle = False
+        # Index the atoms and residues
+        self.residue_list.assign_indexes()
+        self.atom_list.assign_indexes()
+        if not hasattr(dest, 'write'):
+            own_handle = True
+            dest = open(dest, 'w')
+
+        # Assign the formats we need to write with
+        if ext:
+            atmfmt1 = ('%10d %-8s %-8i %-8s %-8s %4d %10.6f %13.4f' + 11*' ')
+            atmfmt2 = ('%10d %-8s %-8i %-8s %-8s %-4s %10.6f %13.4f' + 11*' ')
+            intfmt = '%10d' # For pointers
+        else:
+            atmfmt1 = ('%8d %-4s %-4i %-4s %-4s %4d %10.6f %13.4f' + 11*' ')
+            atmfmt2 = ('%8d %-4s %-4i %-4s %-4s %-4s %10.6f %13.4f' + 11*' ')
+            intfmt = '%8d' # For pointers
+
+        # Now print the header then the title
+        dest.write('PSF ' + ' '.join(self.flags) + '\n')
+        dest.write('\n')
+        dest.write(intfmt % len(self.title) + ' !NTITLE\n')
+        dest.write('\n'.join(self.title) + '\n\n')
+        # Now time for the atoms
+        dest.write(intfmt % len(self.atom_list) + ' !NATOM\n')
+        # atmfmt1 is for CHARMM format (i.e., atom types are integers)
+        # atmfmt is for XPLOR format (i.e., atom types are strings)
+        for i, atom in enumerate(self.atom_list):
+            if isinstance(atom.attype, str):
+                fmt = atmfmt2
+            else:
+                fmt = atmfmt1
+            atmstr = fmt % (i+1, atom.system, atom.residue.resnum,
+                            atom.residue.resname, atom.name, atom.attype,
+                            atom.charge, atom.mass)
+            dest.write(atmstr + '   '.join(atom.props) + '\n')
+        dest.write('\n')
+        # Bonds
+        dest.write(intfmt % len(self.bond_list) + ' !NBOND: bonds\n')
+        for i, bond in enumerate(self.bond_list):
+            dest.write((intfmt*2) % (bond.atom1.idx+1, bond.atom2.idx+1))
+            if i % 4 == 3: # Write 4 bonds per line
+                dest.write('\n')
+        # See if we need to terminate
+        if len(self.bond_list) % 4 != 0 or len(self.bond_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # Angles
+        dest.write(intfmt % len(self.angle_list) + ' !NTHETA: angles\n')
+        for i, angle in enumerate(self.angle_list):
+            dest.write((intfmt*3) % (angle.atom1.idx+1, angle.atom2.idx+1,
+                                     angle.atom3.idx+1)
+            )
+            if i % 3 == 2: # Write 3 angles per line
+                dest.write('\n')
+        # See if we need to terminate
+        if len(self.angle_list) % 3 != 0 or len(self.angle_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # Dihedrals
+        dest.write(intfmt % len(self.dihedral_list) + ' !NPHI: dihedrals\n')
+        for i, dih in enumerate(self.dihedral_list):
+            dest.write((intfmt*4) % (dih.atom1.idx+1, dih.atom2.idx+1,
+                                     dih.atom3.idx+1, dih.atom4.idx+1)
+            )
+            if i % 2 == 1: # Write 2 dihedrals per line
+                dest.write('\n')
+        # See if we need to terminate
+        if len(self.dihedral_list) % 2 != 0 or len(self.dihedral_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # Impropers
+        dest.write(intfmt % len(self.improper_list) + ' !NIMPHI: impropers\n')
+        for i, imp in enumerate(self.improper_list):
+            dest.write((intfmt*4) % (imp.atom1.idx+1, imp.atom2.idx+1,
+                                     imp.atom3.idx+1, imp.atom4.idx+1)
+            )
+            if i % 2 == 1: # Write 2 dihedrals per line
+                dest.write('\n')
+        # See if we need to terminate
+        if len(self.improper_list) % 2 != 0 or len(self.improper_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # Donor section
+        dest.write(intfmt % len(self.donor_list) + ' !NDON: donors\n')
+        for i, don in enumerate(self.donor_list):
+            dest.write((intfmt*2) % (don.atom1.idx+1, don.atom2.idx+1))
+            if i % 4 == 3: # 4 donors per line
+                dest.write('\n')
+        if len(self.donor_list) % 4 != 0 or len(self.donor_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # Acceptor section
+        dest.write(intfmt % len(self.acceptor_list) + ' !NACC: acceptors\n')
+        for i, acc in enumerate(self.acceptor_list):
+            dest.write((intfmt*2) % (acc.atom1.idx+1, acc.atom2.idx+1))
+            if i % 4 == 3: # 4 donors per line
+                dest.write('\n')
+        if len(self.acceptor_list) % 4 != 0 or len(self.acceptor_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # NNB section ??
+        dest.write(intfmt % 0 + ' !NNB\n\n')
+        for i in range(len(self.atom_list)):
+            dest.write(intfmt % 0)
+            if i % 8 == 7: # Write 8 0's per line
+                dest.write('\n')
+        if len(self.atom_list) % 8 != 0: dest.write('\n')
+        dest.write('\n')
+        # Group section
+        dest.write((intfmt*2) % (len(self.group_list), self.group_list.nst2))
+        dest.write(' !NGRP NST2\n')
+        for i, gp in enumerate(self.group_list):
+            dest.write((intfmt*3) % (gp.bs, gp.type, gp.move))
+            if i % 3 == 2: dest.write('\n')
+        if len(self.group_list) % 3 != 0 or len(self.group_list) == 0:
+            dest.write('\n')
+        dest.write('\n')
+        # The next two sections are never found in VMD prmtops...
+        if not vmd:
+            # Molecule section; first set molecularity
+            set_molecules(self.atom_list)
+            mollist = [a.marked for a in self.atom_list]
+            dest.write(intfmt % max(mollist) + ' !MOLNT\n')
+            for i, atom in enumerate(self.atom_list):
+                dest.write(intfmt % atom.marked)
+                if i % 8 == 7: dest.write('\n')
+            if len(self.atom_list) % 8 != 0: dest.write('\n')
+            dest.write('\n')
+            # NUMLP/NUMLPH section
+            dest.write((intfmt*2) % (0, 0) + ' !NUMLP NUMLPH\n')
+            dest.write('\n')
+        # CMAP section
+        dest.write(intfmt % len(self.cmap_list) + ' !NCRTERM: cross-terms\n')
+        for i, cmap in enumerate(self.cmap_list):
+            dest.write((intfmt*4) % (cmap.atom1.idx+1, cmap.atom2.idx+1,
+                                     cmap.atom3.idx+1, cmap.atom4.idx+1)
+            )
+            if cmap.consecutive:
+                dest.write((intfmt*4) % (cmap.atom2.idx+1, cmap.atom3.idx+1,
+                                         cmap.atom4.idx+1, cmap.atom5.idx+1)
+                )
+            else:
+                dest.write((intfmt*4) % (cmap.atom5.idx+1, cmap.atom6.idx+1,
+                                         cmap.atom7.idx+1, cmap.atom8.idx+1)
+                )
+            dest.write('\n')
+        # Done!
+        # If we opened our own handle, close it
+        if own_handle:
+            dest.close()
+        
+    def load_parameters(self, parmset):
+        """
+        Loads parameters from a parameter set that was loaded via CHARMM RTF,
+        PAR, and STR files.
+
+        Parameters:
+            - parmset (ParameterSet) : List of all parameters
+
+        Notes:
+            - If any parameters that are necessary cannot be found, a
+              MissingParameter exception is raised.
+
+            - If any dihedral or improper parameters cannot be found, I will try
+              inserting wildcards (at either end for dihedrals and as the two
+              central atoms in impropers) and see if that matches.  Wild-cards
+              will apply ONLY if specific parameters cannot be found.
+
+            - This method will expand the dihedral_list attribute by adding a
+              separate Dihedral object for each term for types that have a
+              multi-term expansion
+        """
+        # First load the atom types
+        types_are_int = False
+        for atom in self.atom_list:
+            try:
+                if isinstance(atom.attype, int):
+                    atype = parmset.atom_types_int[atom.attype]
+                    types_are_int = True # if we have to change back
+                else:
+                    atype = parmset.atom_types_str[atom.attype]
+            except KeyError:
+                raise MissingParameter('Could not find atom type for %s' %
+                                       atom.attype)
+            atom.type = atype
+            # Change to string attype to look up the rest of the parameters
+            atom.type_to_str()
+
+        # Next load all of the bonds
+        for bond in self.bond_list:
+            # Construct the key
+            key = (min(bond.atom1.attype, bond.atom2.attype),
+                   max(bond.atom1.attype, bond.atom2.attype))
+            try:
+                bond.bond_type = parmset.bond_types[key]
+            except KeyError:
+                raise MissingParameter('Missing bond type for %r' % bond)
+        # Next load all of the angles. If a Urey-Bradley term is defined for
+        # this angle, also build the urey_bradley and urey_bradley_type lists
+        self.urey_bradley_list = TrackedList()
+        for ang in self.angle_list:
+            # Construct the key
+            key = (min(ang.atom1.attype, ang.atom3.attype), ang.atom2.attype,
+                   max(ang.atom1.attype, ang.atom3.attype))
+            try:
+                ang.angle_type = parmset.angle_types[key]
+                ubt = parmset.urey_bradley_types[key]
+                if ubt is not NoUreyBradley:
+                    ub = UreyBradley(ang.atom1, ang.atom3, ubt)
+                    self.urey_bradley_list.append(ub)
+            except KeyError:
+                raise MissingParameter('Missing angle type for %r' % ang)
+        # Next load all of the dihedrals. This is a little trickier since we
+        # need to back up the existing dihedral list and replace it with a
+        # longer one that has only one Fourier term per Dihedral instance.
+        dihedral_list = self.dihedral_list
+        self.dihedral_list = TrackedList()
+        for dih in dihedral_list:
+            # Store the atoms
+            a1, a2, a3, a4 = dih.atom1, dih.atom2, dih.atom3, dih.atom4
+            at1, at2, at3, at4 = a1.attype, a2.attype, a3.attype, a4.attype
+            # First see if the exact dihedral is specified
+            key = min((at1,at2,at3,at4), (at4,at3,at2,at1))
+            if not key in parmset.dihedral_types:
+                # Check for wild-cards
+                key = min(('X',at2,at3,'X'), ('X',at3,at2,'X'))
+                if not key in parmset.dihedral_types:
+                    raise MissingParameter('No dihedral parameters found for '
+                                           '%r' % dih)
+            dtlist = parmset.dihedral_types[key]
+            for i, dt in enumerate(dtlist):
+                self.dihedral_list.append(Dihedral(a1, a2, a3, a4, dt))
+                # See if we include the end-group interactions for this
+                # dihedral. We do IFF it is the last or only dihedral term and
+                # it is NOT in the angle/bond partners
+                if i != len(dtlist) - 1:
+                    self.dihedral_list[-1].end_groups_active = False
+                elif a1 in a4.bond_partners or a1 in a4.angle_partners:
+                    self.dihedral_list[-1].end_groups_active = False
+        # Now do the impropers
+        for imp in self.improper_list:
+            # Store the atoms
+            a1, a2, a3, a4 = imp.atom1, imp.atom2, imp.atom3, imp.atom4
+            at1, at2, at3, at4 = a1.attype, a2.attype, a3.attype, a4.attype
+            key = tuple(sorted([at1, at2, at3, at4]))
+            if not key in parmset.improper_types:
+                # Check for wild-cards
+                for anchor in (at2, at3, at4):
+                    key = tuple(sorted([at1, anchor, 'X', 'X']))
+                    if key in parmset.improper_types:
+                        break # This is the right key
+            try:
+                imp.improper_type = parmset.improper_types[key]
+            except KeyError:
+                raise MissingParameter('No improper parameters found for %r' %
+                                       imp)
+        # Now do the cmaps. These will not have wild-cards
+        for cmap in self.cmap_list:
+            # Store the atoms for easy reference
+            if cmap.consecutive:
+                a1, a2, a3, a4 = cmap.atom1, cmap.atom2, cmap.atom3, cmap.atom4
+                a5, a6, a7, a8 = cmap.atom2, cmap.atom3, cmap.atom4, cmap.atom5
+            else:
+                a1, a2, a3, a4 = cmap.atom1, cmap.atom2, cmap.atom3, cmap.atom4
+                a5, a6, a7, a8 = cmap.atom5, cmap.atom6, cmap.atom7, cmap.atom8
+            at1, at2, at3, at4 = a1.attype, a2.attype, a3.attype, a4.attype
+            at5, at6, at7, at8 = a5.attype, a6.attype, a7.attype, a8.attype
+            # Construct the keys
+            k1 = list(min((at1,at2,at3,at4), (at4,at3,at2,at1)))
+            k2 = list(min((at5,at6,at7,at8), (at8,at7,at6,at5)))
+            key = tuple(k1 + k2)
+            try:
+                cmap.cmap_type = parmset.cmap_types[key]
+            except KeyError:
+                raise MissingParameter('No CMAP parameters found for %r' % cmap)
+        # If the types started out as integers, change them back
+        if types_are_int:
+            for atom in self.atom_list: atom.type_to_int()
+
+    def set_coordinates(self, positions, velocities=None):
+        """
+        This method loads the coordinates and velocity information from an
+        external object or passed data.
+
+        Parameters:
+            - positions (list of floats) : A 3-N length iterable with all of the
+                coordinates in the order [x1, y1, z1, x2, y2, z2, ...].
+            - velocities (list of floats) : If not None, is the velocity
+                equivalent of the positions
+        """
+        if len(positions) / 3 != len(self.atom_list):
+            raise ValueError('Coordinates given for %s atoms, but %d atoms '
+                             'exist in this structure.' %
+                             (len(positions)/3, len(self.atom_list)))
+        # Now assign all of the atoms positions
+        for i, atom in enumerate(self.atom_list):
+            atom.xx = positions[3*i  ]
+            atom.xy = positions[3*i+1]
+            atom.xz = positions[3*i+2]
+
+        # Do velocities if given
+        if velocities is not None:
+            if len(velocities) / 3 != len(self.atom_list):
+                raise ValueError('Velocities given for %s atoms, but %d atoms '
+                                 'exist in this structure.' %
+                                 (len(velocities)/3, len(self.atom_list)))
+            for i, atom in enumerate(self.atom_list):
+                atom.vx = velocities[3*i  ]
+                atom.vy = velocities[3*i+1]
+                atom.vz = velocities[3*i+2]
+            self.velocities = velocities
+
+        self.positions = positions
+
+    def set_box(self, box):
+        """
+        Sets the periodic box boundary conditions.
+
+        Parameters:
+            - box (list of 6 floats) : A list of 6 numbers representing a, b, c,
+                alpha, beta, and gamma, respectively (box lengths and angles).
+                If None, the system is assumed to be aperiodic
+        Notes:
+            The box can alternatively be a list of 3 numbers representing the
+            box lengths. In this case, the angles are assumed to be 90 degrees
+
+            The box here is copied via slicing, so changing the box that was
+            passed in will have no effect after set_box is called.
+        """
+        if box is None:
+            self.box = None
+        elif len(box) == 6:
+            self.box = list(box[:])
+        elif len(box) == 3:
+            self.box = list(box[:]) + [90.0, 90.0, 90.0]
+        else:
+            raise ValueError('set_box requires 3 box lengths, 3 box lengths '
+                             'and 3 angles, or None for no box')
+        # If we already have a _topology instance, then we have possibly changed
+        # the existence of box information (whether or not this is a periodic
+        # system), so delete any cached reference to a topology so it's
+        # regenerated with updated information
+        if hasattr(self, '_topology'):
+            del self._topology
+
+    @property
+    def topology(self):
+        """ Create an OpenMM Topology object from the stored bonded network """
+        try:
+            return self._topology
+        except AttributeError:
+            # If none exists, we need to create it
+            pass
+        # Cache the topology for easy returning later
+        self._topology = topology = Topology()
+        
+        last_chain = None
+        last_residue = None
+        # Add each chain (separate 'system's) and residue
+        for atom in self.atom_list:
+            if atom.system != last_chain:
+                chain = topology.addChain()
+                last_chain = atom.system
+            if atom.residue.idx != last_residue:
+                last_residue = atom.residue.idx
+                residue = topology.addResidue(atom.residue.resname, chain)
+            element_name = element_by_mass(atom.mass)
+            elem = element.get_by_symbol(element_name)
+            topology.addAtom(atom.name, elem, residue)
+
+        # Add all of the bonds
+        atoms = list(topology.atoms())
+        # Assign atom indexes to make sure they're current
+        self.atom_list.assign_indexes()
+        for bond in self.bond_list:
+            topology.addBond(atoms[bond.atom1.idx], atoms[bond.atom2.idx])
+
+        # Add the periodic box if there is one
+        if hasattr(self, 'box') and self.box is not None:
+            topology.setUnitCellDimensions(self.box[:3] * u.angstroms)
+
+        return topology
+
+    def _get_gb_params(self, gb_model=HCT):
+        """ Gets the GB parameters. Need this method to special-case GB neck """
+        if gb_model is GBn:
+            screen = [0.5 for atom in self.atom_list]
+            for i, atom in enumerate(self.atom_list):
+                if atom.type.atomic_number == 6:
+                    screen[i] = 0.48435382330
+                elif atom.type.atomic_number == 1:
+                    screen[i] = 1.09085413633
+                elif atom.type.atomic_number == 7:
+                    screen[i] = 0.700147318409
+                elif atom.type.atomic_number == 8:
+                    screen[i] = 1.06557401132
+                elif atom.type.atomic_number == 16:
+                    screen[i] = 0.602256336067
+        elif gb_model is GBn2:
+            # Add non-optimized values as defaults
+            alpha = [1.0 for i in self.atom_list]
+            beta = [0.8 for i in self.atom_list]
+            gamma = [4.85 for i in self.atom_list]
+            screen = [0.5 for i in self.atom_list]
+            for i, atom in enumerate(self.atom_list):
+                if atom.type.atomic_number == 6:
+                    screen[i] = 1.058554
+                    alpha[i] = 0.733756
+                    beta[i] = 0.506378
+                    gamma[i] = 0.205844
+                elif atom.type.atomic_number == 1:
+                    screen[i] = 1.425952
+                    alpha[i] = 0.788440
+                    beta[i] = 0.798699
+                    gamma[i] = 0.437334
+                elif atom.type.atomic_number == 7:
+                    screen[i] = 0.733599
+                    alpha[i] = 0.503364
+                    beta[i] = 0.316828
+                    gamma[i] = 0.192915
+                elif atom.type.atomic_number == 8:
+                    screen[i] = 1.061039
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+                elif atom.type.atomic_number == 16:
+                    screen[i] = -0.703469
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+        else:
+            # Set the default screening parameters
+            for i, atom in enumerate(self.atom_list):
+                if atom.type.atomic_number == 1:
+                    screen[i] = 0.85
+                elif atom.type.atomic_number == 6:
+                    screen[i] = 0.72
+                elif atom.type.atomic_number == 7:
+                    screen[i] = 0.79
+                elif atom.type.atomic_number == 8:
+                    screen[i] = 0.85
+                elif atom.type.atomic_number == 9:
+                    screen[i] = 0.88
+                elif atom.type.atomic_number == 15:
+                    screen[i] = 0.86
+                elif atom.type.atomic_number == 16:
+                    screen[i] = 0.96
+                else:
+                    screen[i] = 0.8
+
+        length_conv = u.angstrom.conversion_factor_to(u.nanometer)
+        radii = [rad * length_conv for rad in self.parm_data['RADII']]
+
+        if gb_model is GBn2:
+            return zip(radii, screen, alpha, beta, gamma)
+        return zip(radii, screen)
+
+    def createSystem(self, nonbondedMethod=ff.NoCutoff,
+                     nonbondedCutoff=1.0*u.nanometer,
+                     switchDistance=0.0*u.nanometer,
+                     constraints=None,
+                     rigidWater=True,
+                     implicitSolvent=None,
+                     implicitSolventKappa=None,
+                     implicitSolventSaltConc=0.0*u.moles/u.liter,
+                     temperature=298.15*u.kelvin,
+                     soluteDielectric=1.0,
+                     solventDielectric=78.5,
+                     removeCMMotion=True,
+                     hydrogenMass=None,
+                     ewaldErrorTolerance=0.0005,
+                     flexibleConstraints=True,
+                     verbose=False):
+        """
+        Construct an OpenMM System representing the topology described by the
+        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.
+         -  switchDistance (distance=0*nanometer) The distance at which the
+               switching function is active for nonbonded interactions. If the
+               switchDistance evaluates to boolean False (if it is 0), no
+               switching function will be used. Illegal values will raise a
+               ValueError
+         -  constraints (object=None) Specifies which bonds or angles should be
+               implemented with constraints. Allowed values are None, HBonds,
+               AllBonds, or HAngles.
+         -  rigidWater (boolean=True) If true, water molecules will be fully
+               rigid regardless of the value passed for the constraints argument
+         -  implicitSolvent (object=None) If not None, the implicit solvent
+               model to use. Allowed values are HCT, OBC1, OBC2, or GBn
+         -  implicitSolventKappa (float=None): Debye screening parameter to
+               model salt concentrations in GB solvent.
+         -  implicitSolventSaltConc (float=0.0*u.moles/u.liter): Salt
+               concentration for GB simulations. Converted to Debye length
+               `kappa'
+         -  temperature (float=298.15*u.kelvin): Temperature used in the salt
+               concentration-to-kappa conversion for GB salt concentration term
+         -  soluteDielectric (float=1.0) The solute dielectric constant to use
+               in the implicit solvent model.
+         -  solventDielectric (float=78.5) The solvent dielectric constant to
+               use in the implicit solvent model.
+         -  removeCMMotion (boolean=True) If true, a CMMotionRemover will be
+               added to the System.
+         -  hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to
+               heavy atoms. Any mass added to a hydrogen is subtracted from the
+               heavy atom to keep their total mass the same.
+         -  ewaldErrorTolerance (float=0.0005) The error tolerance to use if the
+               nonbonded method is Ewald or PME.
+         -  flexibleConstraints (bool=True) Are our constraints flexible or not?
+         -  verbose (bool=False) Optionally prints out a running progress report
+        """
+        hasbox = self.topology.getUnitCellDimensions() is not None
+        # Set the cutoff distance in nanometers
+        cutoff = None
+        if nonbondedMethod is not ff.NoCutoff:
+            cutoff = nonbondedCutoff
+            # Remove units from cutoff
+            if u.is_quantity(cutoff):
+                cutoff = cutoff.value_in_unit(u.nanometers)
+
+        if nonbondedMethod not in (ff.NoCutoff, ff.CutoffNonPeriodic,
+                                   ff.CutoffPeriodic, ff.Ewald, ff.PME):
+            raise ValueError('Illegal value for nonbonded method')
+        if not hasbox and nonbondedMethod in (ff.CutoffPeriodic,
+                                              ff.Ewald, ff.PME):
+            raise ValueError('Illegal nonbonded method for a '
+                             'non-periodic system')
+        if implicitSolvent not in (HCT, OBC1, OBC2, GBn, GBn2, None):
+            raise ValueError('Illegal implicit solvent model choice.')
+        if not constraints in (None, ff.HAngles, ff.HBonds, ff.AllBonds):
+            raise ValueError('Illegal constraints choice')
+      
+        # Define conversion factors
+        length_conv = u.angstrom.conversion_factor_to(u.nanometer)
+        _chmfrc = u.kilocalorie_per_mole/(u.angstrom*u.angstrom)
+        _openmmfrc = u.kilojoule_per_mole/(u.nanometer*u.nanometer)
+        bond_frc_conv = _chmfrc.conversion_factor_to(_openmmfrc)
+        _chmfrc = u.kilocalorie_per_mole/(u.radians*u.radians)
+        _openmmfrc = u.kilojoule_per_mole/(u.radians*u.radians)
+        angle_frc_conv = _chmfrc.conversion_factor_to(_openmmfrc)
+        dihe_frc_conv = u.kilocalorie_per_mole.conversion_factor_to(
+                            u.kilojoule_per_mole)
+        ene_conv = dihe_frc_conv
+      
+        # Create the system
+        system = mm.System()
+        if verbose: print('Adding particles...')
+        for atom in self.atom_list:
+            system.addParticle(atom.mass)
+        # Set up the constraints
+        if verbose and (constraints is not None and not rigidWater):
+            print('Adding constraints...')
+        if constraints in (ff.HBonds, ff.AllBonds, ff.HAngles):
+            for bond in self.bond_list:
+                if (bond.atom1.type.atomic_number != 1 and
+                    bond.atom2.type.atomic_number != 1):
+                    continue
+                system.addConstraint(bond.atom1.idx, bond.atom2.idx,
+                                     bond.bond_type.req*length_conv)
+        if constraints in (ff.AllBonds, ff.HAngles):
+            for bond in self.bond_list:
+                if (bond.atom1.type.atomic_number == 1 or
+                    bond.atom2.type.atomic_number == 1):
+                    continue
+                system.addConstraint(bond.atom1.idx, bond.atom2.idx,
+                                     bond.bond_type.req*length_conv)
+        if rigidWater and constraints is None:
+            for bond in self.bond_list:
+                if (bond.atom1.type.atomic_number != 1 and
+                    bond.atom2.type.atomic_number != 1):
+                    continue
+                if (bond.atom1.residue.resname in WATNAMES and
+                    bond.atom2.residue.resname in WATNAMES):
+                    system.addConstraint(bond.atom1.idx, bond.atom2.idx,
+                                         bond.bond_type.req*length_conv)
+        # Add Bond forces
+        if verbose: print('Adding bonds...')
+        force = mm.HarmonicBondForce()
+        force.setForceGroup(self.BOND_FORCE_GROUP)
+        # See which, if any, energy terms we omit
+        omitall = not flexibleConstraints and constraints is ff.AllBonds
+        omith = omitall or (flexibleConstraints and constraints in
+                            (ff.HBonds, ff.AllBonds, ff.HAngles))
+        for bond in self.bond_list:
+            if omitall: continue
+            if omith and (bond.atom1.type.atomic_number == 1 or
+                          bond.atom2.type.atomic_number == 1):
+                continue
+            force.addBond(bond.atom1.idx, bond.atom2.idx,
+                          bond.bond_type.req*length_conv,
+                          2*bond.bond_type.k*bond_frc_conv)
+        system.addForce(force)
+        # Add Angle forces
+        if verbose: print('Adding angles...')
+        force = mm.HarmonicAngleForce()
+        force.setForceGroup(self.ANGLE_FORCE_GROUP)
+        if constraints is ff.HAngles:
+            num_constrained_bonds = system.getNumConstraints()
+            atom_constraints = [[]] * system.getNumParticles()
+            for i in range(num_constrained_bonds):
+                c = system.getConstraintParameters(i)
+                dist = c[2].value_in_unit(u.nanometer)
+                atom_constraints[c[0]].append((c[1], dist))
+                atom_constraints[c[1]].append((c[0], dist))
+        for angle in self.angle_list:
+            # Only constrain angles including hydrogen here
+            if (angle.atom1.type.atomic_number != 1 and
+                angle.atom2.type.atomic_number != 1 and
+                angle.atom3.type.atomic_number != 1):
+                continue
+            if constraints is ff.HAngles:
+                a1 = angle.atom1.atomic_number
+                a2 = angle.atom2.atomic_number
+                a3 = angle.atom3.atomic_number
+                nh = int(a1==1) + int(a2==1) + int(a3==1)
+                constrained = (nh >= 2 or (nh == 1 and a2 == 8))
+            else:
+                constrained = False # no constraints
+            if constrained:
+                l1 = l2 = None
+                for bond in angle.atom2.bonds:
+                    if bond.atom1 is angle.atom1 or bond.atom2 is angle.atom1:
+                        l1 = bond.bond_type.req * length_conv
+                    elif bond.atom1 is angle.atom3 or bond.atom2 is angle.atom3:
+                        l2 = bond.bond_type.req * length_conv
+                # Compute the distance between the atoms and add a constraint
+                length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*
+                              cos(angle.angle_type.theteq))
+                system.addConstraint(bond.atom1.idx, bond.atom2.idx, length)
+            if flexibleConstraints or not constrained:
+                force.addAngle(angle.atom1.idx, angle.atom2.idx,
+                               angle.atom3.idx, angle.angle_type.theteq*pi/180,
+                               2*angle.angle_type.k*angle_frc_conv)
+        for angle in self.angle_list:
+            # Already did the angles with hydrogen above. So skip those here
+            if (angle.atom1.type.atomic_number == 1 or
+                angle.atom2.type.atomic_number == 1 or
+                angle.atom3.type.atomic_number == 1):
+                continue
+            force.addAngle(angle.atom1.idx, angle.atom2.idx,
+                           angle.atom3.idx, angle.angle_type.theteq*pi/180,
+                           2*angle.angle_type.k*angle_frc_conv)
+        system.addForce(force)
+
+        # Add the urey-bradley terms
+        if verbose: print('Adding Urey-Bradley terms')
+        force = mm.HarmonicBondForce()
+        force.setForceGroup(self.UREY_BRADLEY_FORCE_GROUP)
+        for ub in self.urey_bradley_list:
+            force.addBond(ub.atom1.idx, ub.atom2.idx,
+                          ub.ub_type.req*length_conv,
+                          2*ub.ub_type.k*bond_frc_conv)
+        system.addForce(force)
+
+        # Add dihedral forces
+        if verbose: print('Adding torsions...')
+        force = mm.PeriodicTorsionForce()
+        force.setForceGroup(self.DIHEDRAL_FORCE_GROUP)
+        for tor in self.dihedral_list:
+            force.addTorsion(tor.atom1.idx, tor.atom2.idx, tor.atom3.idx,
+                             tor.atom4.idx, tor.dihedral_type.per,
+                             tor.dihedral_type.phase*pi/180,
+                             tor.dihedral_type.phi_k*dihe_frc_conv)
+        system.addForce(force)
+
+        if verbose: print('Adding impropers...')
+        # Ick. OpenMM does not have an improper torsion class. Need to
+        # construct one from CustomTorsionForce
+        force = mm.CustomTorsionForce('k*(theta-theta0)^2')
+        force.addPerTorsionParameter('k')
+        force.addPerTorsionParameter('theta0')
+        force.setForceGroup(self.IMPROPER_FORCE_GROUP)
+        for imp in self.improper_list:
+            force.addTorsion(imp.atom1.idx, imp.atom2.idx,
+                             imp.atom3.idx, imp.atom4.idx,
+                             (imp.improper_type.k*dihe_frc_conv,
+                              imp.improper_type.phieq*pi/180)
+            )
+        system.addForce(force)
+
+        if hasattr(self, 'cmap_list'):
+            if verbose: print('Adding CMAP coupled torsions...')
+            force = mm.CMAPTorsionForce()
+            force.setForceGroup(self.CMAP_FORCE_GROUP)
+            # First get the list of cmap maps we're going to use. Just store the
+            # IDs so we have simple integer comparisons to do later
+            cmap_type_list = []
+            cmap_map = dict()
+            for cmap in self.cmap_list:
+                if not id(cmap.cmap_type) in cmap_type_list:
+                    ct = cmap.cmap_type
+                    cmap_type_list.append(id(ct))
+                    # Our torsion correction maps need to go from 0 to 360
+                    # degrees
+                    grid = ct.grid.switch_range().T
+                    m = force.addMap(ct.resolution, [x*ene_conv for x in grid])
+                    cmap_map[id(ct)] = m
+            # Now add in all of the cmaps
+            for cmap in self.cmap_list:
+                if cmap.consecutive:
+                    id1, id2 = cmap.atom1.idx, cmap.atom2.idx
+                    id3, id4 = cmap.atom3.idx, cmap.atom4.idx
+                    id5, id6 = cmap.atom2.idx, cmap.atom3.idx
+                    id7, id8 = cmap.atom4.idx, cmap.atom5.idx
+                else:
+                    id1, id2 = cmap.atom1.idx, cmap.atom2.idx
+                    id3, id4 = cmap.atom3.idx, cmap.atom4.idx
+                    id5, id6 = cmap.atom5.idx, cmap.atom6.idx
+                    id7, id8 = cmap.atom7.idx, cmap.atom8.idx
+                force.addTorsion(cmap_map[id(cmap.cmap_type)],
+                                 id1, id2, id3, id4, id5, id6, id7, id8)
+            system.addForce(force)
+        # Add nonbonded terms now
+        if verbose: print('Adding nonbonded interactions...')
+        force = mm.NonbondedForce()
+        force.setForceGroup(self.NONBONDED_FORCE_GROUP)
+        if not hasbox: # non-periodic
+            if nonbondedMethod is ff.NoCutoff:
+                force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+            elif nonbondedMethod is ff.CutoffNonPeriodic:
+                if cutoff is None:
+                    raise ValueError('No cutoff value specified')
+                force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
+                force.setCutoffDistance(cutoff)
+            else:
+                raise ValueError('Illegal nonbonded method for non-periodic '
+                                 'system')
+
+            # See if we need to use a switching function
+            if switchDistance and nonbondedMethod is not ff.NoCutoff:
+                # make sure it's legal
+                if switchDistance >= nonbondedCutoff:
+                    raise ValueError('switchDistance is too large compared '
+                                     'to the cutoff!')
+                    force.setUseSwitchingFunction(True)
+                    force.setSwitchingDistance(switchDistance)
+
+        else: # periodic
+            # Set up box vectors (from inpcrd if available, or fall back to
+            # prmtop definitions
+            system.setDefaultPeriodicBoxVectors(*self.box_vectors)
+
+            # Set cutoff
+            if cutoff is None:
+                # Compute cutoff automatically
+                box = self.box_lengths
+                min_box_width = min((box[0]/u.nanometers,
+                                     box[1]/u.nanometers,
+                                     box[2]/u.nanometers))
+                CLEARANCE_FACTOR = 0.97
+                cutoff = u.Quantity((min_box_width*CLEARANCE_FACTOR)/2.0,
+                                    u.nanometers)
+            if nonbondedMethod is not ff.NoCutoff:
+                force.setCutoffDistance(cutoff)
+
+            # Set nonbonded method.
+            if nonbondedMethod is ff.NoCutoff:
+                force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+            elif nonbondedMethod is ff.CutoffNonPeriodic:
+                force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
+            elif nonbondedMethod is ff.CutoffPeriodic:
+                force.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
+            elif nonbondedMethod is ff.Ewald:
+                force.setNonbondedMethod(mm.NonbondedForce.Ewald)
+            elif nonbondedMethod is ff.PME:
+                force.setNonbondedMethod(mm.NonbondedForce.PME)
+            else:
+                raise ValueError('Cutoff method is not understood')
+
+            # See if we need to use a switching function
+            if switchDistance and nonbondedMethod is not ff.NoCutoff:
+                # make sure it's legal
+                if switchDistance >= nonbondedCutoff:
+                    raise ValueError('switchDistance is too large compared '
+                                     'to the cutoff!')
+                    force.setUseSwitchingFunction(True)
+                    force.setSwitchingDistance(switchDistance)
+
+            if ewaldErrorTolerance is not None:
+                force.setEwaldErrorTolerance(ewaldErrorTolerance)
+
+        # Add per-particle nonbonded parameters (LJ params)
+        sigma_scale = 2**(-1/6) * 2
+        for i, atm in enumerate(self.atom_list):
+            force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv,
+                              abs(atm.type.epsilon*ene_conv))
+
+        # Add 1-4 interactions
+        excluded_atom_pairs = set() # save these pairs so we don't zero them out
+        sigma_scale = 2**(-1/6)
+        for tor in self.dihedral_list:
+            # First check to see if atoms 1 and 4 are already excluded because
+            # they are 1-2 or 1-3 pairs (would happen in 6-member rings or
+            # fewer). Then check that they're not already added as exclusions
+            if tor.atom1 in tor.atom4.bond_partners: continue
+            if tor.atom1 in tor.atom4.angle_partners: continue
+            key = min((tor.atom1.idx, tor.atom4.idx),
+                      (tor.atom4.idx, tor.atom1.idx))
+            if key in excluded_atom_pairs: continue # multiterm...
+            charge_prod = (tor.atom1.charge * tor.atom4.charge)
+            epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv *
+                            abs(tor.atom4.type.epsilon_14) * ene_conv))
+            sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * (
+                     length_conv * sigma_scale)
+            force.addException(tor.atom1.idx, tor.atom4.idx,
+                               charge_prod, sigma, epsilon)
+            excluded_atom_pairs.add(
+                    min((tor.atom1.idx, tor.atom4.idx),
+                        (tor.atom4.idx, tor.atom1.idx))
+            )
+
+        # Add excluded atoms
+        for atom in self.atom_list:
+            # Exclude all bonds and angles
+            for atom2 in atom.bond_partners:
+                if atom2.idx > atom.idx:
+                    force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+            for atom2 in atom.angle_partners:
+                if atom2.idx > atom.idx:
+                    force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+            for atom2 in atom.dihedral_partners:
+                if atom2.idx <= atom.idx: continue
+                if ((atom.idx, atom2.idx) in excluded_atom_pairs):
+                    continue
+                force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+        system.addForce(force)
+
+        # Add GB model if we're doing one
+        if implicitSolvent is not None:
+            if verbose: print('Adding GB parameters...')
+            gb_parms = self._get_gb_params(implicitSolvent)
+
+            # If implicitSolventKappa is None, compute it from salt
+            # concentration
+            if implicitSolventKappa is None:
+                if u.is_quantity(implicitSolventSaltConc):
+                    sc = implicitSolventSaltConc.value_in_unit(u.moles/u.liter)
+                    implicitSolventSaltConc = sc
+                if u.is_quantity(temperature):
+                    temperature = temperature.value_in_unit(u.kelvin)
+                # The constant is 1 / sqrt( epsilon_0 * kB / (2 * NA * q^2 *
+                # 1000) ) where NA is avogadro's number, epsilon_0 is the
+                # permittivity of free space, q is the elementary charge (this
+                # number matches sander/pmemd's kappa conversion factor)
+                implicitSolventKappa = 50.33355 * sqrt(implicitSolventSaltConc /
+                                                solventDielectric / temperature)
+                # Multiply by 0.73 to account for ion exclusions, and multiply
+                # by 10 to convert to 1/nm from 1/angstroms
+                implicitSolventKappa *= 7.3
+            elif implicitSolvent is None:
+                implicitSolventKappa = 0.0
+
+            if u.is_quantity(implicitSolventKappa):
+                implicitSolventKappa = implicitSolventKappa.value_in_unit(
+                                            (1.0/u.nanometer).unit)
+            if implicitSolvent is HCT:
+                gb = GBSAHCTForce(solventDielectric, soluteDielectric, None,
+                                  cutoff, kappa=implicitSolventKappa)
+            elif implicitSolvent is OBC1:
+                gb = GBSAOBC1Force(solventDielectric, soluteDielectric, None,
+                                   cutoff, kappa=implicitSolventKappa)
+            elif implicitSolvent is OBC2:
+                gb = GBSAOBC2Force(solventDielectric, soluteDielectric, None,
+                                   cutoff, kappa=implicitSolventKappa)
+            elif implicitSolvent is GBn:
+                gb = GBSAGBnForce(solventDielectric, soluteDielectric, None,
+                                  cutoff, kappa=implicitSolventKappa)
+            elif implicitSolvent is GBn2:
+                gb = GBSAGBn2Force(solventDielectric, soluteDielectric, None,
+                                   cutoff, kappa=implicitSolventKappa)
+            for i, atom in enumerate(self.atom_list):
+                gb.addParticle([atom.charge] + list(gb_parms[i]))
+            # Set cutoff method
+            if nonbondedMethod is ff.NoCutoff:
+                gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
+            elif nonbondedMethod is ff.CutoffNonPeriodic:
+                gb.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
+                gb.setCutoffDistance(cutoff)
+            elif nonbondedMethod is ff.CutoffPeriodic:
+                gb.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
+                gb.setCutoffDistance(cutoff)
+            else:
+                raise ValueError('Illegal nonbonded method for use with GBSA')
+            gb.setForceGroup(self.GB_FORCE_GROUP)
+            system.addForce(gb)
+            force.setReactionFieldDielectric(1.0) # applies to NonbondedForce
+
+        # See if we repartition the hydrogen masses
+        if hydrogenMass is not None:
+            for bond in self.bond_list:
+                # Only take the ones with at least one hydrogen
+                if (bond.atom1.type.atomic_number != 1 and
+                    bond.atom2.type.atomic_number != 1):
+                    continue
+                atom1, atom2 = bond.atom1, bond.atom2
+                if atom1.type.atomic_number == 1:
+                    atom1, atom2 = atom2, atom1 # now atom2 is hydrogen for sure
+                if atom1.type.atomic_number != 1:
+                    transfer_mass = hydrogenMass - atom2.mass
+                    new_mass1 = (system.getParticleMass(atom1.idx) -
+                                 transfer_mass)
+                    system.setParticleMass(atom2.idx, hydrogenMass)
+                    system.setParticleMass(atom1.idx, new_mass1)
+        # See if we want to remove COM motion
+        if removeCMMotion:
+            system.addForce(mm.CMMotionRemover())
+
+        # Cache our system for easy access
+        self._system = system
+
+        return system
+
+    @property
+    def system(self):
+        """
+        Return the cached system class -- it needs to be initialized via
+        "createSystem" first!
+        """
+        try:
+            return self._system
+        except AttributeError:
+            raise APIError('You must initialize the system with createSystem '
+                           'before accessing the cached object.')
+
+    @property
+    def positions(self):
+        """
+        Return the cached positions or create new ones from the atoms
+        """
+        try:
+            if len(self._positions) == len(self.atom_list):
+                return self._positions
+        except AttributeError:
+            pass
+
+        self._positions = tuple([Vec3(a.xx, a.xy, a.xz)
+                               for a in self.atom_list]) * u.angstroms
+        return self._positions
+
+    @property
+    def velocities(self):
+        """ Same as for positions, but for velocities """
+        try:
+            if len(self._velocities) == len(self.atom_list):
+                return self._velocities
+        except AttributeError:
+            pass
+
+        self._velocities = tuple([Vec3(a.vx, a.vy, a.vz)
+                    for a in self.atom_list]) * (u.angstroms/u.picosecond) 
+        return self._velocities
+
+    @property
+    def box_vectors(self):
+        """ Return tuple of box vectors """
+        if hasattr(self, 'rst7'):
+            box = [x*u.angstrom for x in self.rst7.box[:3]]
+            ang = [self.rst7.box[3], self.rst7.box[4], self.rst7.box[5]]
+            return _box_vectors_from_lengths_angles(box[0], box[1], box[2],
+                                                    ang[0], ang[1], ang[2])
+        else:
+            box = [x*u.angstrom for x in self.parm_data['BOX_DIMENSIONS'][1:]]
+            ang = [self.parm_data['BOX_DIMENSIONS'][0]] * 3
+            return _box_vectors_from_lengths_angles(box[0], box[1], box[2],
+                                                    ang[0], ang[1], ang[2])
+
+    @property
+    def box_lengths(self):
+        """ Return tuple of 3 units """
+        if hasattr(self, 'rst7'):
+            box = [x*u.angstrom for x in self.rst7.box[:3]]
+        else:
+            box = [x*u.angstrom for x in self.parm_data['BOX_DIMENSIONS'][1:]]
+        return tuple(box)
+
+def _box_vectors_from_lengths_angles(a, b, c, alpha, beta, gamma):
+    """
+    This method takes the lengths and angles from a unit cell and creates unit
+    cell vectors.
+
+    Parameters:
+        - a (unit, dimension length): Length of the first vector
+        - b (unit, dimension length): Length of the second vector
+        - c (unit, dimension length): Length of the third vector
+        - alpha (float): Angle between b and c in degrees
+        - beta (float): Angle between a and c in degrees
+        - gamma (float): Angle between a and b in degrees
+
+    Returns:
+        Tuple of box vectors (as Vec3 instances)
+   """
+    if not (u.is_quantity(a) and u.is_quantity(b) and u.is_quantity(c)):
+        raise TypeError('a, b, and c must be units of dimension length')
+    if u.is_quantity(alpha): alpha = alpha.value_in_unit(u.degree)
+    if u.is_quantity(beta): beta = beta.value_in_unit(u.degree)
+    if u.is_quantity(gamma): gamma = gamma.value_in_unit(u.degree)
+    a = a.value_in_unit(u.angstrom)
+    b = b.value_in_unit(u.angstrom)
+    c = c.value_in_unit(u.angstrom)
+
+    if alpha <= 2 * pi and beta <= 2 * pi and gamma <= 2 * pi:
+        raise ValueError('box angles must be given in degrees')
+
+    alpha *= pi / 180
+    beta *= pi / 180
+    gamma *= pi / 180
+
+    av = Vec3(a, 0.0, 0.0) * u.angstrom
+    bx = b * cos(gamma)
+    by = b * sin(gamma)
+    bz = 0.0
+    cx = c * cos(beta)
+    cy = c * (cos(alpha) - cos(beta) * cos(gamma))
+    cz = sqrt(c * c - cx * cx - cy * cy)
+   
+    # Make sure any components that are close to zero are set to zero exactly
+    if abs(bx) < TINY: bx = 0.0
+    if abs(by) < TINY: by = 0.0
+    if abs(cx) < TINY: cx = 0.0
+    if abs(cy) < TINY: cy = 0.0
+    if abs(cz) < TINY: cz = 0.0
+
+    bv = Vec3(bx, by, bz) * u.angstrom
+    cv = Vec3(cx, cy, cz) * u.angstrom
+
+    return (av, bv, cv)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+def set_molecules(atom_list):
+    """
+    Correctly sets the ATOMS_PER_MOLECULE and SOLVENT_POINTERS sections of the
+    topology file.
+    """
+    from sys import setrecursionlimit, getrecursionlimit
+    # Since we use a recursive function here, we make sure that the recursion
+    # limit is large enough to handle the maximum possible recursion depth we'll
+    # need (NATOM). We don't want to shrink it, though, since we use list
+    # comprehensions in list constructors in some places that have an implicit
+    # (shallow) recursion, therefore, reducing the recursion limit too much here
+    # could raise a recursion depth exceeded exception during a _Type/Atom/XList
+    # creation. Therefore, set the recursion limit to the greater of the current
+    # limit or the number of atoms
+    setrecursionlimit(max(len(atom_list), getrecursionlimit()))
+
+    # Unmark all atoms so we can track which molecule each goes into
+    atom_list.unmark()
+
+    # The molecule "ownership" list
+    owner = []
+    # The way I do this is via a recursive algorithm, in which
+    # the "set_owner" method is called for each bonded partner an atom
+    # has, which in turn calls set_owner for each of its partners and 
+    # so on until everything has been assigned.
+    molecule_number = 1 # which molecule number we are on
+    for i in range(len(atom_list)):
+        # If this atom has not yet been "owned", make it the next molecule
+        # However, we only increment which molecule number we're on if 
+        # we actually assigned a new molecule (obviously)
+        if not atom_list[i].marked:
+            tmp = [i]
+            _set_owner(atom_list, tmp, i, molecule_number)
+            # Make sure the atom indexes are sorted
+            tmp.sort()
+            owner.append(tmp)
+            molecule_number += 1
+    return owner
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+def _set_owner(atom_list, owner_array, atm, mol_id):
+    """ Recursively sets ownership of given atom and all bonded partners """
+    atom_list[atm].marked = mol_id
+    for partner in atom_list[atm].bond_partners:
+        if not partner.marked:
+            owner_array.append(partner.idx)
+            _set_owner(atom_list, owner_array, partner.idx, mol_id)
+        elif partner.marked != mol_id:
+            raise MoleculeError('Atom %d in multiple molecules' % 
+                                partner.idx)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()

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

+"""
+This module contains data structures useful for parsing in CHARMM files and
+constructing chemical structures from those files
+
+Author: Jason M. Swails
+Contributors:
+Date: April 9, 2014
+"""
+from simtk.openmm.app.charmm.exceptions import (SplitResidueWarning, BondError,
+                ResidueError, CmapError, MissingParameter)
+import warnings
+
+TINY = 1e-8
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class AtomType(object):
+    """
+    Atom types can either be compared by indexes or names. Can be assigned with
+    a string, integer, (string is not automatically cast) or with both. Create
+    new atom types with the "add" constructor to make sure the registry is
+    filled with only unique types
+
+    Parameters and Attributes:
+        - name (str) : The name of the atom type
+        - number (int) : The integer index of the atom type
+        - mass (float) : The mass of the atom type
+        - atomic_number (int) : The atomic number of the element of the atom
+                                type
+    Attributes:
+        - name (str) : The name of the atom type
+        - number (int) : The integer index of the atom type
+        - _member_number (int, private) : The order in which this atom type
+                was 'added' this is used to make sure that atom types added
+                last have priority in assignment in the generated hash tables
+    Example:
+    >>> at = AtomType('HA', 1, 1.008, 1)
+    >>> at.name, at.number
+    ('HA', 1)
+    >>> at2 = AtomType('CA', 2, 12.01, 6)
+    >>> at2.name, at2.number
+    ('CA', 2)
+    >>> print "%s: %d" % (str(at), int(at))
+    HA: 1
+    >>> print at == WildCard
+    True
+    >>> print at2 == WildCard
+    True
+    """
+
+    def __init__(self, name, number, mass, atomic_number):
+        if number is None and name is not None:
+            # If we were given an integer, assign it to number. Otherwise,
+            # assign it to the name
+            if isinstance(name, int):
+                self.number = name
+                self.name = None
+            else:
+                self.name = name
+                self.number = None
+        else:
+            self.name = name
+            self.number = int(number)
+        self.mass = mass
+        self.atomic_number = atomic_number
+        # We have no LJ parameters as of yet
+        self.epsilon = self.rmin = self.epsilon_14 = self.rmin_14 = None
+
+    def __eq__(self, other):
+        """
+        Compares based on available properties (name and number, just name,
+        or just number)
+        """
+        if other is WildCard: return True # all atom types match wild cards
+        if isinstance(other, AtomType):
+            return self.name == other.name and self.number == other.number
+        if isinstance(other, basestring):
+            return self.name == other
+        if isinstance(other, int):
+            return self.number == other
+        return other == (self.number, self.name)
+
+    def set_lj_params(self, eps, rmin, eps14=None, rmin14=None):
+        """ Sets Lennard-Jones parameters on this atom type """
+        if eps14 is None:
+            eps14 = eps
+        if rmin14 is None:
+            rmin14 = rmin
+        self.epsilon = eps
+        self.rmin = rmin
+        self.epsilon_14 = eps14
+        self.rmin_14 = rmin14
+
+    def __int__(self):
+        """ The integer representation of an AtomType is its index """
+        return self.number
+
+    def __str__(self):
+        return self.name
+
+    # Comparisons are all based on number
+    def __gt__(self, other):
+        return self._member_number > other._member_number
+    def __lt__(self, other):
+        return self._member_number < other._member_number
+    def __ge__(self, other):
+        return self._member_number > other._member_number or self == other
+    def __le__(self, other):
+        return self._member_number < other._member_number or self == other
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class WildCard(AtomType):
+    """
+    This is a wild-card atom type that matches ANY atom type. It is often used
+    in torsion parameters for the end-groups. Some properties:
+
+        - It is a singleton, so seeing if an AtomType is a WildCard should be
+          done using the "is" binary operator
+
+        - It is defined as "equal" to every other object via the "==" operator
+
+        - It is less than everything
+
+        - It is greater than nothing
+    """
+
+    def __init__(self):
+        self._member_number = -1
+        self.name = 'X'
+        self.number = None
+        self.mass = None
+
+    # Define comparison operators
+    def __eq__(self, other): return True
+    def __lt__(self, other): return True
+    def __gt__(self, other): return False
+    def __le__(self, other): return True
+    def __ge__(self, other): return True
+
+WildCard = WildCard() # Turn it into a singleton
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Atom(object):
+    """ An atom in a structure.
+
+    Parameters:
+        system (str) : Name of the system this atom belongs to
+        name (str): name of the atom
+        type (str or int) : Type of the atom
+        charge (float) : Partial atomic charge (elementary charge units)
+        mass (float) : Atomic mass (amu)
+        props (list) : Other properties from the PSF
+
+    Attributes:
+        - attype (str) : Name of the atom type
+        - system (str) : The system name associated with this atom
+        - name (str) : Name of the atom (str)
+        - charge (float) : Partial atomic charge
+        - mass (float) : Mass of the atom (amu)
+        - idx (int) : index of the atom in the system, starting from 0
+        - props (list) : List of extraneous properties parsed from a PSF
+        - type (AtomType) : If assigned, has additional properties like the
+             non-bonded LJ parameters. If None, it has not yet been assigned
+
+    Possible Attributes (SOA == Set of Atom instances)
+        - bond_partners (SOA) : List of all atoms I am bonded to
+        - angle_partners (SOA) : List of all atoms I am angled to
+        - dihedral_partners (SOA) : List of all atoms I am dihedraled to
+        - bonds (list of Bond's) : All bonds to which I belong
+        - angles (list of Angle's) : All angles to which I belong
+        - dihedrals (list of Dihedral's) : All dihedrals to which I belong
+        - impropers (list of Improper's) : All impropers to which I belong
+        - cmaps (list of Cmap's) : All correction maps to which I belong
+    """
+    def __init__(self, system, name, attype, charge, mass, props=None):
+        self.name = name
+        self.attype = attype
+        self.type = None
+        self.charge = charge
+        self.mass = mass
+        self.idx = -1
+        self.props = props
+        self.system = system
+        self.marked = 0 # For recursive molecule determination
+        self._bond_partners = set()
+        self._angle_partners = set()
+        self._dihedral_partners = set()
+        self.bonds = []
+        self.angles = []
+        self.urey_bradleys = []
+        self.dihedrals = []
+        self.impropers = []
+        self.cmaps = []
+
+    def bond_to(self, other):
+        """
+        Register this atom as bonded partners. Cannot bond to itself. If that
+        is attempted, a BondError is raised
+        """
+        if self is other:
+            raise BondError('Cannot bond atom to itself')
+        self._bond_partners.add(other)
+        other._bond_partners.add(self)
+
+    def angle_to(self, other):
+        """
+        Register this atom as angle partners. Cannot angle to itself. If that
+        is attempted, a BondError is raised
+        """
+        if self is other:
+            raise BondError('Cannot angle atom to itself')
+        self._angle_partners.add(other)
+        other._angle_partners.add(self)
+
+    def dihedral_to(self, other):
+        """
+        Register this atom as dihedral partners. Cannot dihedral to itself. If
+        that is attempted, a BondError is raised
+        """
+        if self is other:
+            raise BondError('Cannot dihedral atom to itself')
+        self._dihedral_partners.add(other)
+        other._dihedral_partners.add(self)
+
+    @property
+    def bond_partners(self):
+        return sorted(list(self._bond_partners))
+
+    @property
+    def angle_partners(self):
+        return sorted(list(self._angle_partners - self._bond_partners))
+
+    @property
+    def dihedral_partners(self):
+        return sorted(list(self._dihedral_partners - self._angle_partners -
+                           self._bond_partners))
+
+    def type_to_int(self):
+        """
+        Changes the type to an integer, matching CHARMM conventions. This can
+        only be done if a type mapping has been loaded (i.e., if the type
+        attribute is not None).
+
+        If the type identifier is not already an integer and no mapping is
+        available, MissingParameter is raised.
+        """
+        if isinstance(self.attype, int):
+            return
+        if self.type is None:
+            raise MissingParameter('No type mapping loaded. Cannot change '
+                                   'type identifier to integer for %s' %
+                                   self.attype)
+        self.attype = self.type.number
+
+    def type_to_str(self):
+        """
+        Changes the type to a string, matching XPLOR conventions. This can only
+        be done if a type mapping has been loaded (i.e., if the type attribute
+        is not None).
+
+        If the type identifier is not already a string and no mapping is
+        available, MissingParameter is raised.
+        """
+        if isinstance(self.attype, str):
+            return
+        if self.type is None:
+            raise MissingParameter('No type mapping loaded. Cannot change '
+                                   'type identifier to string for %s' %
+                                   self.attype)
+        self.attype = self.type.name
+
+    def __repr__(self):
+        retstr = '<Atom %d' % self.idx
+        if hasattr(self, 'residue'):
+            retstr += '; %d %s [%s: %s]>' % (
+                        self.residue.idx, self.residue.resname, self.name,
+                        self.attype)
+        else:
+            retstr += '; %s> ' % (self.name)
+        return retstr
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class AtomList(list):
+    """ A list of Atom instances.  """
+
+    def unmark(self):
+        for atom in self: atom.marked = 0
+
+    def assign_indexes(self):
+        self._index_us()
+
+    def _index_us(self):
+        for i, atom in enumerate(self): atom.idx = i
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Residue(object):
+    """ Residue class """
+
+    def __init__(self, resname, idx):
+        self.resname = resname
+        self.idx = idx
+        self.atoms = []
+        self.resnum = None # Numbered based on SYSTEM name
+        self.system = None
+
+    def add_atom(self, atom):
+        if self.system is None:
+            self.system = atom.system
+        else:
+            if self.system != atom.system:
+                raise ResidueError('Added atom has a different system than '
+                                   'the other atoms in the residue!')
+        atom.residue = self
+        self.atoms.append(atom)
+
+    def delete_atom(self, atom):
+        """
+        If an atom is present in this residue, delete it from the list of
+        atoms
+        """
+        self.atoms = [a for a in self.atoms if a is not atom]
+
+    # Implement some container methods over the list of atoms
+    def __contains__(self, thing):
+        """ True if an atom is present in this residue """
+        return thing in self.atoms
+
+    def __len__(self):
+        return len(self.atoms)
+
+    def __getitem__(self, idx):
+        return self.atoms[idx]
+
+    # Equality
+    def __eq__(self, thing):
+        if isinstance(thing, Residue):
+            return self.resname == thing.resname and self.idx == thing.idx
+        if isinstance(thing, tuple) or isinstance(thing, list):
+            # Must be resnum, resname
+            return thing == (self.resname, self.idx)
+        return False # No other type can be equal.
+
+    def __ne__(self, thing):
+        return not self.__eq__(thing)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class ResidueList(list):
+    """
+    A list of residues in a chemical system. Atoms are added to this list and
+    spawn new residues when the residue number or name changes.
+    """
+    def __init__(self):
+        list.__init__(self)
+        self._last_residue = None
+
+    def assign_indexes(self):
+        """ Indexes all residues starting from 1 """
+        last_system = None
+        resnum = 1
+        for i, res in enumerate(self):
+            res.idx = i + 1
+            if res.system != last_system:
+                res.resnum = resnum = 1
+                last_system = res.system
+            else:
+                res.resnum = resnum
+            resnum += 1
+
+    def add_atom(self, system, resnum, resname, name,
+                 attype, charge, mass, props=None):
+        """
+        Adds an atom to the list of residues. If the residue is not the same as
+        the last residue that was created, a new Residue is created and added
+        to this list
+
+        Parameters:
+            - system (str) : The system this atom belongs to
+            - resnum (int) : Residue number
+            - resname (str) : Name of the residue
+            - name (str) : Name of the atom
+            - attype (int or str) : Type of the atom
+            - charge (float) : Partial atomic charge of the atom
+            - mass (float) : Mass (amu) of the atom
+
+        Returns:
+            The Atom instance created and added to the list of residues
+        """
+        if self._last_residue is None:
+            res = self._last_residue = Residue(resname, resnum)
+            list.append(self, res)
+        elif self._last_residue != (resname, resnum):
+            if self._last_residue.idx == resnum:
+                lresname = self._last_residue.resname
+                warnings.warn('Residue %d split into separate residues %s '
+                              'and %s' % (resnum, lresname, resname),
+                              SplitResidueWarning)
+            res = self._last_residue = Residue(resname, resnum)
+            list.append(self, res)
+        else:
+            res = self._last_residue
+        atom = Atom(system, name, attype, float(charge), float(mass), props)
+        res.add_atom(atom)
+        return atom
+    
+    def append(self, thing):
+        raise NotImplemented('Use "add_atom" to build a residue list')
+
+    extend = append
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Bond(object):
+    """
+    A bond object that links 2 atoms
+
+    Parameters:
+        - atom1 (Atom) : First atom included in the bond
+        - atom2 (Atom) : Second atom included in the bond
+        - bond_type (BondType) : Type for the bond (None if unknown)
+    """
+    def __init__(self, atom1, atom2, bond_type=None):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.bond_type = bond_type
+        self.atom1.bond_to(atom2)
+        # Add this bond to the atoms
+        self.atom1.bonds.append(self)
+        self.atom2.bonds.append(self)
+
+    def __repr__(self):
+        return '<Bond; %r -- %r; type=%r>' % (self.atom1, self.atom2,
+                                              self.bond_type)
+
+    def __contains__(self, thing):
+        """ See if an atom is part of this bond """
+        return thing is self.atom1 or thing is self.atom2
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Angle(object):
+    """
+    An angle object that links 3 atoms
+
+    Parameters:
+        - atom1 (Atom) : First atom included in the angle
+        - atom2 (Atom) : Central atom in the valence angle
+        - atom3 (Atom) : Third atom in the valence angle
+        - angle_type (AngleType) : Type for the angle (None if unknown)
+    """
+    def __init__(self, atom1, atom2, atom3, angle_type=None):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.atom3 = atom3
+        self.angle_type = angle_type
+        self.atom1.angle_to(atom3)
+        # Add this angle to the atoms
+        self.atom1.angles.append(self)
+        self.atom2.angles.append(self)
+        self.atom3.angles.append(self)
+
+    def __repr__(self):
+        return '<Angle; %r-%r-%r; type=%r>' % (self.atom1, self.atom2,
+                                               self.atom3, self.angle_type)
+
+    def __contains__(self, thing):
+        """ See if a bond or an atom is in this angle """
+        if isinstance(thing, Bond):
+            return self.atom2 in thing and (self.atom1 in thing or
+                                            self.atom2 in thing)
+        # Otherwise assume it's an atom
+        return self.atom1 is thing or self.atom2 is thing or self.atom3 is thing
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class UreyBradley(object):
+    """
+    A harmonic restraint between two atoms separated by 2 valence bonds (i.e.,
+    involved in a valence angle with each other
+
+    Parameters:
+        - atom1 (Atom) : The first atom included in the Urey-Bradley term
+        - atom2 (Atom) : The second atom included in the Urey-Bradley term
+        - ub_type (UreyBradleyType) : The type for the Urey-Bradley term (None
+                                      if unknown)
+    """
+    def __init__(self, atom1, atom2, ub_type=None):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.ub_type = ub_type
+        # Add this Urey-Bradley to the atoms
+        atom1.urey_bradleys.append(self)
+        atom2.urey_bradleys.append(self)
+
+    def __repr__(self):
+        return '<UreyBradley; %r-%r; type=%r>' % (self.atom1, self.atom2,
+                                                  self.ub_type)
+
+    def __contains__(self, thing):
+        # See the comments in chemistry/amber/topologyobjects.py for what's
+        # being done here (under the same method of the UreyBradley type there)
+        if isinstance(thing, Bond):
+            if not thing.atom1 in self:
+                if not thing.atom2 in self:
+                    return False
+                end1 = thing.atom2
+                cent = thing.atom1
+            else:
+                end1 = thing.atom1
+                cent = thing.atom2
+            for atm in cent.bonds():
+                if atm is end1: continue
+                if atm in self:
+                    return True
+            return False
+        # thing must be an atom
+        return thing is self.atom1 or thing is self.atom2
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Dihedral(object):
+    """
+    A torsion angle object that links 4 atoms
+
+    Parameters:
+        - atom1 (Atom) : First atom included in the torsion
+        - atom2 (Atom) : Second atom included in the torsion
+        - atom3 (Atom) : Third atom included in the torsion
+        - atom4 (Atom) : Fourth atom included in the torsion
+        - dihedral_type (DihedralType) : Type for the torsion (None if unknown)
+    """
+    def __init__(self, atom1, atom2, atom3, atom4, dihedral_type=None):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.atom3 = atom3
+        self.atom4 = atom4
+        self.dihedral_type = dihedral_type
+        self.atom1.dihedral_to(atom4)
+        # Add this torsion to the atoms
+        self.atom1.dihedrals.append(self)
+        self.atom2.dihedrals.append(self)
+        self.atom3.dihedrals.append(self)
+        self.atom4.dihedrals.append(self)
+        # Add a marker for indicating if this dihedral is not the final term in
+        # a multi-term expansion or if the atoms at the end are also 1-2 or 1-3
+        # pairs (this can happen for 5- and 6-member rings, respectively)
+        self.end_groups_active = True
+
+    def __repr__(self):
+        return '<Dihedral; %r-%r-%r-%r; type=%r>' % (
+                    self.atom1, self.atom2, self.atom3, self.atom4,
+                    self.dihedral_type)
+
+    def __contains__(self, thing):
+        """ See if a bond or an atom is in this torsion """
+        if isinstance(thing, Bond):
+            if self.atom1 in thing and self.atom2 in thing:
+                return True
+            if self.atom2 in thing and self.atom3 in thing:
+                return True
+            if self.atom3 in thing and self.atom4 in thing:
+                return True
+            return False
+        # Otherwise assume it's an atom
+        return (thing is self.atom1 or thing is self.atom2 or
+                thing is self.atom3 or thing is self.atom4)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Improper(object):
+    """
+    An improper torsion object. The third atom is bonded to each other atom
+
+    Parameters:
+        - atom1 (Atom) : First atom included in the torsion
+        - atom2 (Atom) : Second atom included in the torsion
+        - atom3 (Atom) : Third atom included in the torsion
+        - atom4 (Atom) : Fourth atom included in the torsion
+        - improper_type (ImproperType) : Type for the improper (None if unknown)
+    """
+    def __init__(self, atom1, atom2, atom3, atom4, improper_type=None):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.atom3 = atom3
+        self.atom4 = atom4
+        self.improper_type = improper_type
+        # Nothing to register -- all bonds should already be formed
+        # Add this improper to the atoms
+        self.atom1.impropers.append(self)
+        self.atom2.impropers.append(self)
+        self.atom3.impropers.append(self)
+        self.atom4.impropers.append(self)
+
+    def __repr__(self):
+        return '<Improper; %r-%r-%r-%r; type=%r>' % (
+                    self.atom1, self.atom2, self.atom3, self.atom4,
+                    self.improper_type)
+
+    def __contains__(self, thing):
+        """
+        See if a bond or an atom is in this improper
+        An improper is defined as shown below
+
+                                A3
+                                |
+                                |
+                       A4 ----- A1 ----- A2
+         
+        So the bonds will either be between atom1 and any other atom
+        """
+        if isinstance(thing, Bond):
+            if self.atom2 in thing and self.atom1 in thing:
+                return True
+            if self.atom3 in thing and self.atom1 in thing:
+                return True
+            if self.atom4 in thing and self.atom1 in thing:
+                return True
+            return False
+        # Otherwise assume it's an atom
+        return (thing is self.atom1 or thing is self.atom2 or
+                thing is self.atom3 or thing is self.atom4)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class AcceptorDonor(object):
+    """
+    Just a holder for donors and acceptors in CHARMM speak
+    
+    Parameters:
+        - atom1 (Atom) : First atom in the donor/acceptor group
+        - atom2 (Atom) : Second atom in the donor/acceptor group
+    """
+    def __init__(self, atom1, atom2):
+        self.atom1 = atom1
+        self.atom2 = atom2
+
+    def __repr__(self):
+        return '<AcceptorDonor; %r %r>' % (self.atom1, self.atom2)
+
+    def __contains__(self, thing):
+        """ See if the atom is in this donor/acceptor """
+        return thing is self.atom1 or thing is self.atom2
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Group(object):
+    """
+    An 'interacting' group defined by the PSF.
+
+    Parameters:
+        - bs (int) : ??
+        - type (int) : The group type
+        - move (int) : If the group moves ??
+
+    Disclaimer: I really don't know what these numbers mean. I'm speculating
+    based on the source code of 'chamber', and this section is simply ignored
+    there.
+    """
+    def __init__(self, bs, type, move):
+        self.bs = bs
+        self.type = type
+        self.move = move
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class Cmap(object):
+    """
+    A correction-map (two coupled torsions). Defined between 8 atoms (2
+    consecutive correction maps). "Consecutive torsions" (i.e., those definable
+    by 5 atoms) will only be recognized if the two torsions have the same order
+
+    Parameters:
+        - atom1 (Atom) : 1st atom of first dihedral
+        - atom2 (Atom) : 2nd atom of first dihedral
+        - atom3 (Atom) : 3rd atom of first dihedral
+        - atom4 (Atom) : 4th atom of first dihedral
+        - atom5 (Atom) : 1st atom of second dihedral
+        - atom6 (Atom) : 2nd atom of second dihedral
+        - atom7 (Atom) : 3rd atom of second dihedral
+        - atom8 (Atom) : 4th atom of second dihedral
+        - cmap_type (CmapType) : Cmap type for this cmap (None if unknown)
+
+    Attributes:
+        - consecutive (bool) : Are the dihedrals consecutive?
+
+     if consecutive:
+        - atom1 (Atom) : 1st atom of 1st dihedral
+        - atom2 (Atom) : 2nd atom of 1st dihedral && 1st atom of 2nd dihedral
+        - atom3 (Atom) : 3rd atom of 1st dihedral && 2nd atom of 2nd dihedral
+        - atom4 (Atom) : 4th atom of 1st dihedral && 3rd atom of 2nd dihedral
+        - atom5 (Atom) :                             4th atom of 2nd dihedral
+     
+     if not consecutive:
+        - atom1 (Atom) : 1st atom of first dihedral
+        - atom2 (Atom) : 2nd atom of first dihedral
+        - atom3 (Atom) : 3rd atom of first dihedral
+        - atom4 (Atom) : 4th atom of first dihedral
+        - atom5 (Atom) : 1st atom of second dihedral
+        - atom6 (Atom) : 2nd atom of second dihedral
+        - atom7 (Atom) : 3rd atom of second dihedral
+        - atom8 (Atom) : 4th atom of second dihedral
+    """
+    def __init__(self, atom1, atom2, atom3, atom4, atom5, atom6, atom7,
+                 atom8, cmap_type=None):
+        self.consecutive = False
+        if atom2 is atom5 and atom3 is atom6 and atom4 is atom7:
+            self.consecutive = True
+        if self.consecutive:
+            self.atom1 = atom1
+            self.atom2 = atom2
+            self.atom3 = atom3
+            self.atom4 = atom4
+            self.atom5 = atom8
+            # Add this cmap to the atoms
+            self.atom1.cmaps.append(self)
+            self.atom2.cmaps.append(self)
+            self.atom3.cmaps.append(self)
+            self.atom4.cmaps.append(self)
+            self.atom5.cmaps.append(self)
+        else:
+            self.atom1 = atom1
+            self.atom2 = atom2
+            self.atom3 = atom3
+            self.atom4 = atom4
+            self.atom5 = atom5
+            self.atom6 = atom6
+            self.atom7 = atom7
+            self.atom8 = atom8
+            # Add this cmap to the atoms
+            self.atom1.cmaps.append(self)
+            self.atom2.cmaps.append(self)
+            self.atom3.cmaps.append(self)
+            self.atom4.cmaps.append(self)
+            self.atom5.cmaps.append(self)
+            self.atom6.cmaps.append(self)
+            self.atom7.cmaps.append(self)
+            self.atom8.cmaps.append(self)
+        self.cmap_type = cmap_type
+
+    def __repr__(self):
+        if self.consecutive:
+            return '<Cmap; %r-%r-%r-%r-%r; cmap_type=%r>' % (
+                        self.atom1, self.atom2, self.atom3, self.atom4,
+                        self.atom5, self.cmap_type)
+        else:
+            return '<Cmap; %r-%r-%r-%r && %r-%r-%r-%r; cmap_type=%r>' % (
+                        self.atom1, self.atom2, self.atom3, self.atom4,
+                        self.atom5, self.atom6, self.atom7, self.atom8,
+                        self.cmap_type)
+
+    def __contains__(self, thing):
+        """ See if a Bond or Atom is inside this torsion """
+        if isinstance(thing, Bond):
+            if self.consecutive:
+                return self._consecutive_has_bond(thing)
+            else:
+                return self._nonconsecutive_has_bond(thing)
+        # Must be an atom
+        if self.consecutive:
+            return (thing is self.atom1 or thing is self.atom2 or
+                    thing is self.atom3 or thing is self.atom4 or
+                    thing is self.atom5)
+        return (thing is self.atom1 or thing is self.atom2 or
+                thing is self.atom3 or thing is self.atom4 or
+                thing is self.atom5 or thing is self.atom6 or
+                thing is self.atom7 or thing is self.atom8)
+
+    def _consecutive_has_bond(self, thing):
+        return ((self.atom1 in thing and self.atom2 in thing) or
+                (self.atom2 in thing and self.atom3 in thing) or
+                (self.atom3 in thing and self.atom4 in thing) or
+                (self.atom4 in thing and self.atom5 in thing))
+
+    def _nonconsecutive_has_bond(self, thing):
+        return ((self.atom1 in thing and self.atom2 in thing) or
+                (self.atom2 in thing and self.atom3 in thing) or
+                (self.atom3 in thing and self.atom4 in thing) or
+                (self.atom5 in thing and self.atom6 in thing) or
+                (self.atom6 in thing and self.atom7 in thing) or
+                (self.atom7 in thing and self.atom8 in thing))
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class BondType(object):
+    """
+    A bond type with an equilibrium length (Angstroms) and force constant
+    (kcal/mol/Angstrom^2)
+
+    Parameters:
+        - k (float) : Force constant (kcal/mol/A^2)
+        - req (float) : Equilibrium distance
+    """
+    def __init__(self, k, req):
+        self.k = k
+        self.req = req
+
+    def __eq__(self, other):
+        return self.k == other.k and self.req == other.req
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class AngleType(object):
+    """
+    An angle type with an equilibrium angle (degrees) and force constant
+    (kcal/mol/radians^2)
+
+    Parameters:
+        - k (float) : Force constant (kcal/mol/radians^2)
+        - theteq (float) : Equilibrium angle value (degrees)
+    """
+    def __init__(self, k, theteq):
+        self.k = k
+        self.theteq = theteq
+
+    def __eq__(self, other):
+        return self.k == other.k and self.theteq == other.theteq
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class UreyBradleyType(BondType):
+    """
+    A harmonic spring connecting two atoms separated by two valence bonds (a
+    valence angle). It is functionally equivalent to a Bond and is actually
+    implemented as a (unaltered) Bond subclass. See BondType documentation.
+    """
+
+# Not all angles have Urey-Bradley terms attached to them. This is a singleton
+# that indicates that there is NO U-B term for this particular type
+NoUreyBradley = UreyBradleyType(None, None)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class DihedralType(object):
+    """
+    A torsion angle type with a force constant (kcal/mol), periodicity (int),
+    and phase (degrees)
+
+    Parameters:
+        - phi_k (float) : Force constant (kcal/mol)
+        - per (int) : Periodicity
+        - phase (float): Phase of the torsion
+    """
+    def __init__(self, phi_k, per, phase):
+        self.phi_k = float(phi_k)
+        self.per = int(per)
+        self.phase = float(phase)
+
+    def __eq__(self, other):
+        return (self.phi_k == other.phi_k and self.per == other.per and
+                self.phase == other.phase)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class ImproperType(object):
+    """
+    An improper torsion angle type with a force constant (kcal/mol) and
+    equilibrium angle (degrees)
+
+    Parameters:
+        - k (float) : Force constant (kcal/mol)
+        - phieq (int) : Equilibrium angle (degrees)
+    """
+    def __init__(self, k, phieq):
+        self.k = k
+        self.phieq = phieq
+
+    def __eq__(self, other):
+        return self.k == other.k and self.phieq == other.phieq
+    
+    def __repr__(self):
+        return '<ImproperType; k=%s; phieq=%s>' % (self.k, self.phieq)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+class CmapType(object):
+    """
+    Contains a correction map interpolation grid
+    
+    Parameters:
+        - resolution (int) : Number of interpolation points for each dihedral
+        - grid (list of floats) : resolution x resolution list of energy values
+                (kcal/mol) for the angles with the 2nd angle changing fastest.
+
+    The grid object is converted to a _CmapGrid instance which can be treated
+    like a normal list, but also has the ability to quickly return a transpose
+    (so the 1st angle changes fastest) and to switch the range from -180 -- 180
+    to 0 -- 360 (and back again). This is particularly helpful because CHARMM
+    defines CMAP tables from -180 -- 180 whereas OpenMM expects them from
+    0 -- 360 (with the 1st angle changing fastest!!)
+    """
+
+    def __init__(self, resolution, grid):
+        self.resolution = resolution
+        self.grid = _CmapGrid(resolution, grid)
+        if len(grid) != self.resolution * self.resolution:
+            raise CmapError('CMAP grid does not match expected resolution')
+
+    def __eq__(self, other):
+        return (self.resolution == other.resolution and
+                all([abs(i - j) < TINY for i, j in zip(self.grid, other.grid)]))
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+# Take the CmapGrid class from the Chamber prmtop topology objects
+class _CmapGrid(object):
+    """
+    A grid object for storing Correction map data. Data can be accessed in one
+    of two ways; either with 1 or 2 indexes. If 2 indexes are given, the index
+    into the flattened array is i*resolution+j. Indexing starts from 0.
+
+    The _CmapGrid usually has ranges for the two angles from -180 to 180. Some
+    places will expect the range to be 0-360 degrees (e.g., OpenMM). The
+    switch_range method returns a _CmapGrid with this range. This method will
+    also work to go backwards.
+
+    Example:
+    >>> g = _CmapGrid(2, [0, 1, 2, 3])
+    >>> print g[0], g[0,0]
+    0 0
+    >>> print g[1], g[0,1]
+    1 1
+    >>> print g[1,0]
+    2
+    >>> g[1,1] = 10
+    >>> print g[3]
+    10
+    >>> print g.switch_range()
+    [10.0000, 2.0000
+     1.0000, 0.0000]
+    """
+
+    def __init__(self, resolution, data=None):
+        self.resolution = resolution
+        if data is None:
+            self._data = [0 for i in range(self.resolution*self.resolution)]
+        else:
+            self._data = data
+
+    @property
+    def transpose(self):
+        """ Returns the transpose of the grid """
+        try:
+            return self._transpose
+        except AttributeError:
+            pass
+        _transpose = []
+        size = len(self._data)
+        for i in range(self.resolution):
+            piece = [self[j] for j in range(i, size, self.resolution)]
+            _transpose += piece
+        self._transpose = _CmapGrid(self.resolution, _transpose)
+        return self._transpose
+
+    T = transpose
+
+    def __getitem__(self, idx):
+        if isinstance(idx, tuple):
+            return self._data[self.resolution*idx[0]+idx[1]]
+        return self._data[idx]
+
+    def __setitem__(self, idx, val):
+        if isinstance(idx, tuple):
+            self._data[self.resolution*idx[0]+idx[1]] = val
+        else:
+            self._data[idx] = val
+
+    def __len__(self):
+        return len(self._data)
+
+    def __iter__(self):
+        return iter(self._data)
+
+    def __repr__(self):
+        return '<_CmapGrid: %dx%d>' % (self.resolution, self.resolution)
+
+    def __str__(self):
+        retstr = '[%.4f,' % self._data[0]
+        fmt = ' %.4f'
+        for i, val in enumerate(self):
+            if i == 0: continue
+            retstr += fmt % val
+            if (i+1) % self.resolution == 0 and i != len(self._data) - 1:
+                retstr += '\n'
+            elif i != len(self) - 1:
+                retstr += ','
+        return retstr + ']'
+
+    def __eq__(self, other):
+        try:
+            if self.resolution != other.resolution:
+                return False
+            for x, y in zip(self, other):
+                if abs(x - y) > TINY:
+                    return False
+            return True
+        except AttributeError:
+            return TypeError('Bad type comparison with _CmapGrid')
+
+    def switch_range(self):
+        """
+        Returns a grid object whose range is 0 to 360 degrees in both dimensions
+        instead of -180 to 180 degrees (or -180 to 180 degrees if the range is
+        already 0 to 360 degrees)
+        """
+        res = self.resolution
+        mid = res // 2
+        newgrid = _CmapGrid(res)
+        for i in range(res):
+            for j in range(res):
+                # Start from the middle
+                newgrid[i, j] = self[(i+mid)%res, (j+mid)%res]
+        return newgrid
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+def _tracking(fcn):
+    """ Decorator to indicate the list has changed """
+    def new_fcn(self, *args):
+        self.changed = True
+        return fcn(self, *args)
+    return new_fcn
+
+class TrackedList(list):
+    """
+    This creates a list type that allows you to see if anything has changed
+    """
+    def __init__(self, arg=[]):
+        self.changed = False
+        list.__init__(self, arg)
+
+    __delitem__ = _tracking(list.__delitem__)
+    append = _tracking(list.append)
+    extend = _tracking(list.extend)
+    __delslice__ = _tracking(list.__delslice__)
+    __setitem__ = _tracking(list.__setitem__)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()