Merge branch 'master' into qc

wrappers/python/setup.py

@@ -144,7 +144,8 @@ def buildKeywordDictionary(major_version_num=MAJOR_VERSION_NUM,
                                           "simtk.unit",
                                           "simtk.openmm",
                                           "simtk.openmm.app",
-                                          "simtk.openmm.app.internal"]
+                                          "simtk.openmm.app.internal",
+                                          "simtk.openmm.app.internal.charmm"]
     setupKeywords["data_files"]        = []
     setupKeywords["package_data"]      = {"simtk" : [],
                                           "simtk.unit" : [],

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

@@ -25,6 +25,10 @@ from modeller import Modeller
 from statedatareporter import StateDataReporter
 from element import Element
 from desmonddmsfile import DesmondDMSFile
+from checkpointreporter import CheckpointReporter
+from charmmcrdfiles import CharmmCrdFile, CharmmRstFile
+from charmmparameterset import CharmmParameterSet
+from charmmpsffile import CharmmPsfFile
 
 # Enumerated values
 

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

+"""
+Provides a class for parsing CHARMM-style coordinate files, namely CHARMM .crd
+(coordinate) files and CHARMM .rst (restart) file. Uses CharmmFile class in
+_charmmfile.py for reading files 
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason Deckman
+Contributors: Jason M. Swails
+Date: April 19, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+
+from simtk.openmm.app.internal.charmm.exceptions import CharmmFileError
+import simtk.unit as u
+from simtk.openmm.vec3 import Vec3
+
+CHARMMLEN = 22
+TIMESCALE = 4.888821E-14 * 1e12 # AKMA time units to picoseconds
+ONE_TIMESCALE = 1 / TIMESCALE
+
+class CharmmCrdFile(object):
+    """
+    Reads and parses a CHARMM coordinate file (.crd) into its components,
+    namely the coordinates, CHARMM atom types, resid, resname, etc.
+
+    Main attributes:
+        - natom (int) : Number of atoms in the system
+        - resname (list) : Names of all residues
+        - positions (list) : All cartesian coordinates [x1, y1, z1, x2, ...]
+
+    Example:
+    >>> chm = CharmmCrdFile('testfiles/1tnm.crd')
+    >>> print '%d atoms; %d coords' % (chm.natom, len(chm.positions))
+    1414 atoms; 1414 coords
+    """
+
+    def __init__(self, fname):
+        self.atomno = []                   # Atom number
+        self.resno = []                    # Residue number
+        self.resname = []                  # Residue name
+        self.resid = []                    # Residue ID
+        self.attype = []                   # Atom type
+        self.positions = []                # 3N atomic coordinates
+        self.title = []                    # .crd file title block
+        self.segid = []                    # Segment ID
+        self.weighting = []                # Atom weighting
+
+        self.natom = 0                     # Number of atoms specified in file
+        self._parse(fname)
+
+    def _parse(self, fname):
+
+        crdfile = open(fname, 'r')
+        line = crdfile.readline()
+
+        while len(line.strip()) == 0:      # Skip whitespace, as a precaution
+            line = crdfile.readline()
+
+        intitle = True
+
+        while intitle:
+            self.title.append(line.strip())
+            line = crdfile.readline()
+            if len(line.strip()) == 0:
+                intitle = False
+            elif line.strip()[0] != '*':
+                intitle = False
+            else: 
+                intitle = True
+
+        while len(line.strip()) == 0:      # Skip whitespace
+            line = crdfile.readline()
+        
+        try:
+            self.natom = int(line.strip().split()[0])
+            
+            for row in range(self.natom):
+                line = crdfile.readline().strip().split()
+                self.atomno.append(int(line[0]))
+                self.resno.append(int(line[1]))
+                self.resname.append(line[2])
+                self.attype.append(line[3])
+                pos = Vec3(float(line[4]), float(line[5]), float(line[6]))
+                self.positions.append(pos * u.angstroms)
+                self.segid.append(line[7])
+                self.resid.append(int(line[8]))
+                self.weighting.append(float(line[9]))
+
+            if self.natom != len(self.positions):
+                raise CharmmFileError("Error parsing CHARMM .crd file: %d "
+                                      "atoms requires %d positions (not %d)" %
+                                      (self.natom, self.natom,
+                                       len(self.positions))
+                )
+
+        except (ValueError, IndexError), e:
+            raise CharmmFileError('Error parsing CHARMM coordinate file')
+
+class CharmmRstFile(object):
+    """
+    Reads and parses data, velocities and coordinates from a CHARMM restart
+    file (.rst) of file name 'fname' into class attributes
+
+    Main attributes:
+        - natom (int) : Number of atoms in the system
+        - resname (list) : Names of all residues
+        - positions (list) : All cartesian coordinates [x1, y1, z1, x2, ...]
+        - positionsold (list) : Old cartesian coordinates
+        - velocities (list) : List of all cartesian velocities
+
+    Example:
+    >>> chm = CharmmRstFile('testfiles/sample-charmm.rst')
+    >>> print chm.header[0]
+    REST    37     1
+    >>> natom, nc, nco = chm.natom, len(chm.positions), len(chm.positionsold)
+    >>> nv = len(chm.velocities)
+    >>> print '%d atoms; %d crds; %d old crds; %d vels' % (natom, nc, nco, nv)
+    256 atoms; 256 crds; 256 old crds; 256 vels
+    """
+
+    def __init__(self, fname):
+        self.header = []
+        self.title = []
+        self.enrgstat = []
+        self.positionsold = []
+        self.positions = []
+        self.velocities = []
+        
+        self.ff_version = 0
+        self.natom = 0
+        self.npriv = 0
+        self.nstep = 0
+        self.nsavc = 0
+        self.nsavv = 0
+        self.jhstrt = 0
+
+        self._parse(fname)
+
+    def _parse(self, fname):
+
+        crdfile = open(fname, 'r')
+        readingHeader = True 
+
+        while readingHeader:
+            line = crdfile.readline()
+            if not len(line):
+                raise CharmmFileError('Premature end of file')
+            line = line.strip()
+            words = line.split()
+            if len(line) != 0:  
+                if words[0] == 'ENERGIES' or words[0] == '!ENERGIES':
+                    readingHeader = False
+                else:
+                    self.header.append(line.strip())
+            else:
+                self.header.append(line.strip())
+
+        for row in range(len(self.header)):
+            if len(self.header[row].strip()) != 0:
+                line = self.header[row].strip().split()
+                if line[0][0:5] == 'NATOM' or line[0][0:6] == '!NATOM':
+                    try:
+                        line = self.header[row+1].strip().split()
+                        self.natom = int(line[0])     
+                        self.npriv = int(line[1])     # num. previous steps
+                        self.nstep = int(line[2])     # num. steps in file  
+                        self.nsavc = int(line[3])     # coord save frequency 
+                        self.nsavv = int(line[4])     # velocities "
+                        self.jhstrt = int(line[5])    # Num total steps?
+                        break
+                   
+                    except (ValueError, IndexError), e:
+                        raise CharmmFileError('Problem parsing CHARMM restart')
+
+        self._scan(crdfile, '!XOLD')
+        self._get_formatted_crds(crdfile, self.positionsold)
+
+        self._scan(crdfile, '!VX')
+        self._get_formatted_crds(crdfile, self.velocities)
+
+        self._scan(crdfile, '!X')
+        self._get_formatted_crds(crdfile, self.positions)
+
+        # Convert velocities to angstroms/ps
+        self.velocities = [v * ONE_TIMESCALE for v in self.velocities]
+
+        # Add units to positions and velocities
+        self.positions *= u.angstroms
+        self.velocities *= u.angstroms / u.picoseconds
+
+    def _scan(self, handle, str, r=0): # read lines in file until str is found
+        scanning = True
+
+        if(r): handle.seek(0)
+
+        while scanning:
+            line = handle.readline()
+            if not line:
+                raise CharmmFileError('Premature end of file')
+
+            if len(line.strip()) != 0:
+                if line.strip().split()[0][0:len(str)] == str:
+                    scanning = False
+
+
+    def _get_formatted_crds(self, crdfile, crds):
+        for row in range(self.natom):
+            line = crdfile.readline()
+
+            if not line:
+                raise CharmmFileError('Premature end of file')
+
+            if len(line) < 3*CHARMMLEN:
+                raise CharmmFileError("Less than 3 coordinates present in "
+                                      "coordinate row or positions may be "
+                                      "truncated.") 
+
+            line = line.replace('D','E')     # CHARMM uses 'D' for exponentials
+
+            # CHARMM uses fixed format (len = CHARMMLEN = 22) for crds in .rst's
+
+            c = Vec3(float(line[0:CHARMMLEN]), float(line[CHARMMLEN:2*CHARMMLEN]),
+                     float(line[2*CHARMMLEN:3*CHARMMLEN]))
+            crds.append(c)
+
+
+    def printcoords(self, crds):
+        for crd in range(len(crds)):
+            print crds[crd],
+            if not (crd+1) % 3:
+                print '\n',
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()

wrappers/python/simtk/openmm/app/charmmparameterset.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
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason M. Swails
+Contributors:
+Date: April 18, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+import os
+from simtk.openmm.app.internal.charmm._charmmfile import (
+            CharmmFile, CharmmStreamFile)
+from simtk.openmm.app.internal.charmm.topologyobjects import (
+            AtomType, BondType, AngleType, DihedralType, ImproperType, CmapType,
+            UreyBradleyType, NoUreyBradley)
+from simtk.openmm.app.internal.charmm.exceptions import CharmmFileError
+from simtk.openmm.app.element import Element, get_by_symbol
+import simtk.unit as u
+import warnings
+
+class CharmmParameterSet(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)
+
+    Parameters:
+        - filenames : List of topology, parameter, and stream files to load into
+          the parameter set. The following file type suffixes are recognized.
+          Unrecognized file types raise a TypeError
+            .rtf, .top -- Residue topology file
+            .par, .prm -- Parameter file
+            .str -- Stream file
+            .inp -- If "par" is in the file name, it is a parameter file. If
+                    "top" is in the file name, it is a topology file. Otherwise,
+                    raise TypeError
+
+    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
+
+    Example:
+    >>> params = CharmmParameterSet('charmm22.top', 'charmm22.par', 'file.str')
+    """
+
+    @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, *args):
+        # 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 = []
+
+        # Load all of the files
+        tops, pars, strs = [], [], []
+        for arg in args:
+            if arg.endswith('.rtf') or arg.endswith('.top'):
+                tops.append(arg)
+            elif arg.endswith('.par') or arg.endswith('.prm'):
+                pars.append(arg)
+            elif arg.endswith('.str'):
+                strs.append(arg)
+            elif arg.endswith('.inp'):
+                fname = os.path.split(arg)[1]
+                if 'par' in fname:
+                    pars.append(arg)
+                elif 'top' in fname:
+                    tops.append(arg)
+                else:
+                    raise TypeError('Unrecognized file type: %s' % arg)
+            else:
+                raise TypeError('Unrecognized file type: %s' % arg)
+        for top in tops: self.readTopologyFile(top)
+        for par in pars: self.readParameterFile(par)
+        for strf in strs: self.readStreamFile(strf)
+
+    @classmethod
+    def loadSet(cls, tfile=None, pfile=None, sfiles=[]):
+        """
+        Instantiates a CharmmParameterSet 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 CharmmParameterSet 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.readTopologyFile(tfile)
+        if pfile is not None:
+            inst.readParameterFile(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.readStreamFile(sfiles)
+        elif sfiles is not None:
+            for sfile in sfiles:
+                inst.readStreamFile(sfile)
+        return inst
+
+    def readParameterFile(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 = CharmmParameterSet._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
+                    masselem = Element.getByMass(mass)
+                    if masselem is None:
+                        atomic_number = 0 # Extra point or something
+                    else:
+                        atomic_number = masselem.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 readTopologyFile(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 = CharmmParameterSet._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
+                masselem = Element.getByMass(mass)
+                if masselem is None:
+                    atomic_number = 0 # Extra point or something
+                else:
+                    atomic_number = masselem.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 readStreamFile(self, sfile):
+        """
+        Reads RTF and PAR sections from a stream file and dispatches the
+        sections to readTopologyFile or readParameterFile
+
+        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.readTopologyFile(section)
+            elif words[1].startswith('para'):
+                # This is a Parameter file section
+                self.readParameterFile(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 = CharmmParameterSet('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]
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

+"""
+Provides a Python class for parsing a PSF file and setting up a system
+structure for it within the OpenMM framework.
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason M. Swails
+Contributors:
+Date: April 18, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+from __future__ import division
+
+from functools import wraps
+from math import pi, cos, sin, sqrt
+import os
+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)
+# CHARMM imports
+from simtk.openmm.app.internal.charmm._charmmfile import CharmmFile
+from simtk.openmm.app.internal.charmm.topologyobjects import (
+                ResidueList, AtomList, TrackedList, Bond, Angle, Dihedral,
+                Improper, AcceptorDonor, Group, Cmap, UreyBradley,
+                NoUreyBradley)
+from simtk.openmm.app.internal.charmm.exceptions import (
+                CharmmPSFError, MoleculeError, CharmmPSFWarning,
+                MissingParameter)
+import warnings
+
+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 CharmmPsfFile(object):
+    """
+    A chemical structure instantiated from CHARMM files.
+
+    Example:
+    >>> cs = CharmmPsfFile("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 CharmmParameterSet 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 = CharmmPsfFile("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
+    
+    @_catchindexerror
+    def __init__(self, psf_name):
+        """
+        Opens and parses a PSF file, then instantiates a CharmmPsfFile 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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._parse_psf_section(psf, int)
+        # Now get the group sections
+        pointers, holder = CharmmPsfFile._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 = CharmmPsfFile._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 = CharmmPsfFile._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 set the instance attributes
+        self.residue_list = residue_list
+        self.atom_list = atom_list
+        self.bond_list = bond_list
+        self.angle_list = angle_list
+        self.dihedral_list = dihedral_list
+        self.improper_list = improper_list
+        self.cmap_list = cmap_list
+        self.donor_list = donor_list
+        self.acceptor_list = acceptor_list
+        self.group_list = group_list
+        self.title = title
+        self.flags = psf_flags
+        self.box_vectors = None
+
+    @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 = CharmmPsfFile._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
+
+    def writePsf(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 = CharmmPsfFile('testfiles/test.psf')
+            >>> cs.writePsf('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 loadParameters(self, parmset):
+        """
+        Loads parameters from a parameter set that was loaded via CHARMM RTF,
+        PAR, and STR files.
+
+        Parameters:
+            - parmset (CharmmParameterSet) : 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 setCoordinates(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 setBox(self, a, b, c, alpha=90.0*u.degrees, beta=90.0*u.degrees,
+               gamma=90.0*u.degrees):
+        """
+        Sets the periodic box boundary conditions.
+
+        Parameters:
+            - a, b, c (floats) : Lengths of the periodic cell
+            - alpha, beta, gamma (floats, optional) : Angles between the
+                periodic cells.
+        """
+        self.box_vectors = _box_vectors_from_lengths_angles(a, b, c,
+                                                            alpha, beta, gamma)
+        # 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)
+            if atom.type is not None:
+                # This is the most reliable way of determining the element
+                atomic_num = atom.type.atomic_number
+                elem = element.Element.getByAtomicNumber(atomic_num)
+            else:
+                # Figure it out from the mass
+                elem = element.Element.getByMass(atom.mass)
+            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 self.box_vectors is not None:
+            topology.setUnitCellDimensions(self.boxLengths)
+
+        return topology
+
+    def _get_gb_params(self, gb_model=HCT):
+        """ Gets the GB parameters. Need this method to special-case GB neck """
+        screen = [0 for atom in self.atom_list]
+        if gb_model is GBn:
+            radii = _bondi_radii(self.atom_list)
+            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:
+            radii = _mbondi3_radii(self.atom_list)
+            # 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
+            # Determine which radii set we need
+            if gb_model is OBC1 or gb_model is OBC2:
+                radii = _mbondi2_radii(self.atom_list)
+            elif gb_model is HCT:
+                radii = _mbondi_radii(self.atom_list)
+
+        length_conv = u.angstrom.conversion_factor_to(u.nanometer)
+        radii = [x * length_conv for x in radii]
+
+        if gb_model is GBn2:
+            return zip(radii, screen, alpha, beta, gamma)
+        return zip(radii, screen)
+
+    def createSystem(self, params, 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. You MUST have loaded a parameter set into this PSF before
+        calling createSystem. If not, AttributeError will be raised. ValueError
+        is raised for illegal input.
+
+        Parameters:
+         -  params (CharmmParameterSet) The parameter set to use to parametrize
+               this molecule
+         -  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
+        """
+        # back up the dihedral list
+        dihedral_list = self.dihedral_list
+        # Load the parameter set
+        self.loadParameters(params.condense())
+        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.boxLengths
+                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
+
+        # Restore the dihedral list to allow reparametrization later
+        self.dihedral_list = dihedral_list
+
+        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 AttributeError('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
+
+    @positions.setter
+    def positions(self, stuff):
+        """
+        Replace the cached positions and the positions of each atom. If no units
+        are applied to "stuff", it is assumed to be Angstroms.
+        """
+        if not u.is_quantity(stuff):
+            # Assume this is Angstroms
+            stuff *= u.angstroms
+
+        # If we got a 1-D array, reshape it into an natom list of Vec3's
+        if len(stuff) == len(self.atom_list) * 3:
+            stuff = [Vec3(stuff[i*3], stuff[i*3+1], stuff[i*3+2])
+                     for i in range(len(self.atom_list))]
+        self._positions = stuff
+        for atom, pos in zip(self.atom_list, stuff):
+            atom.xx, atom.xy, atom.xz = pos.value_in_unit(u.angstrom)
+
+    @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 boxLengths(self):
+        """ Return tuple of 3 units """
+        if self.box_vectors is not None:
+            return (self.box_vectors[0][0], self.box_vectors[0][1],
+                    self.box_vectors[0][2])
+        return None
+
+    @boxLengths.setter
+    def boxLengths(self, stuff):
+        raise RuntimeError('Use setBox to set a box with lengths and angles '
+                           'or set the boxVectors attribute with box vectors')
+    
+    @property
+    def boxVectors(self):
+        """ Return the box vectors """
+        return self.box_vectors
+
+    @boxVectors.setter
+    def boxVectors(self, stuff):
+        """ Sets the box vectors """
+        self.box_vectors = stuff
+
+    def deleteCmap(self):
+        """ Deletes the CMAP terms from the CHARMM PSF """
+        self.cmap_list = TrackedList()
+
+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 molecularity of the system based on connectivity
+    """
+    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)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+# Routines that set the necessary radii lists based on a list of atoms with
+# their connectivities
+
+def _bondi_radii(atom_list):
+    """ Sets the bondi radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        if atom.type.atomic_number == 6:
+            radii[i] = 1.7
+        elif atom.type.atomic_number == 1:
+            radii[i] = 1.2
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi_radii(atom_list):
+    """ Sets the mbondi radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        # Radius of H atom depends on element it is bonded to
+        if atom.type.atomic_number == 1:
+            bondeds = list(atom.bond_partners)
+            if bondeds[0].type.atomic_number in (6, 7): # C or N
+                radii[i] = 1.3
+            elif bondeds[0].type.atomic_number in (8, 16): # O or S
+                radii[i] = 0.8
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.type.atomic_number == 6:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi2_radii(atom_list):
+    """ Sets the mbondi2 radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        # Radius of H atom depends on element it is bonded to
+        if atom.type.atomic_number == 1:
+            if atom.bond_partners[0].type.atomic_number == 7:
+                radii[i] = 1.3
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.type.atomic_number == 6:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # Converted to nanometers above
+
+def _mbondi3_radii(atom_list):
+    """ Sets the mbondi3 radii """
+    radii = _mbondi2_radii(atom_list)
+    for i, atom in enumerate(atom_list):
+        # Adjust OE (GLU), OD (ASP) and HH/HE (ARG)
+        if atom.residue.resname in ('GLU', 'ASP'):
+            if atom.name.startswith('OE') or atom.name.startswith('OD'):
+                radii[i] = 1.4
+        elif atom.residue.resname == 'ARG':
+            if atom.name.startswith('HH') or atom.name.startswith('HE'):
+                radii[i] = 1.17
+        # Adjust carboxylate O radii on C-termini
+        if atom.name == 'OXT':
+            radii[i] = 1.4
+    return radii  # Converted to nanometers above
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()

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

+"""
+checkpointreporter.py: Saves checkpoint files for a simulation
+
+This is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of
+Biological Structures at Stanford, funded under the NIH Roadmap for
+Medical Research, grant U54 GM072970. See https://simtk.org.
+
+Portions copyright (c) 2014 Stanford University and the Authors.
+Authors: Robert McGibbon
+Contributors: 
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+__author__ = "Robert McGibbon"
+__version__ = "1.0"
+
+import simtk.openmm as mm
+__all__ = ['CheckpointReporter']
+
+
+class CheckpointReporter(object):
+    """CheckpointReporter saves periodic checkpoints of a simulation.
+    The checkpoints will overwrite one another -- only the last checkpoint
+    will be saved in the file.
+
+    To use it, create a CheckpointReporter, then add it to the Simulation's
+    list of reporters. To load a checkpoint file and continue a simulation,
+    use the following recipe:
+
+    >>> with open('checkput.chk', 'rb') as f:
+    >>>     simulation.context.loadCheckpoint(f.read())
+
+    Notes:
+    A checkpoint contains not only publicly visible data such as the particle
+    positions and velocities, but also internal data such as the states of
+    random number generators.  Ideally, loading a checkpoint should restore the
+    Context to an identical state to when it was written, such that continuing
+    the simulation will produce an identical trajectory.  This is not strictly
+    guaranteed to be true, however, and should not be relied on.  For most
+    purposes, however, the internal state should be close enough to be
+    reasonably considered equivalent.
+
+    A checkpoint contains data that is highly specific to the Context from
+    which it was created. It depends on the details of the System, the
+    Platform being used, and the hardware and software of the computer it was
+    created on.  If you try to load it on a computer with different hardware,
+    or for a System that is different in any way, loading is likely to fail.
+    Checkpoints created with different versions of OpenMM are also often
+    incompatible.  If a checkpoint cannot be loaded, that is signaled by
+    throwing an exception.
+
+    """
+    def __init__(self, file, reportInterval):
+        """Create a CheckpointReporter.
+
+        Parameters:
+         - file (string or open file object) The file to write to. Any current
+           contents will be overwritten.
+         - reportInterval (int) The interval (in time steps) at which to write checkpoints
+        """
+
+        self._reportInterval = reportInterval
+        if isinstance(file, basestring):
+            self._own_handle = true
+            self._out = open(file, 'w+b', 0)
+        else:
+            self._out = file
+            self._own_handle = False
+
+    def describeNextReport(self, simulation):
+        """Get information about the next report this object will generate.
+
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+
+        Returns: A five element tuple.  The first element is the number of steps until the
+        next report.  The remaining elements specify whether that report will require
+        positions, velocities, forces, and energies respectively.
+        """
+        steps = self._reportInterval - simulation.currentStep%self._reportInterval
+        return (steps, False, False, False, False)
+
+    def report(self, simulation, state):
+        """Generate a report.
+
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+         - state (State) The current state of the simulation
+        """
+        self._out.seek(0)
+        chk = simulation.context.createCheckpoint()
+        self._out.write(chk)
+        self._out.truncate()
+        self._out.flush()
+
+    def __del__(self):
+        if self._own_handle:
+            self._out.close()
+

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

@@ -34,7 +34,7 @@ __author__ = "Christopher M. Bruns"
 __version__ = "1.0"
 
 
-from simtk.unit import daltons
+from simtk.unit import daltons, is_quantity
 import copy_reg
 
 class Element(object):
@@ -42,22 +42,34 @@ class Element(object):
 
     The simtk.openmm.app.element module contains objects for all the standard chemical elements,
     such as element.hydrogen or element.carbon.  You can also call the static method Element.getBySymbol() to
-    look up the Element with a particular chemical symbol."""
+    look up the Element with a particular chemical symbol.
+
+    Element objects should be considered immutable
+    """
 
     _elements_by_symbol = {}
     _elements_by_atomic_number = {}
 
     def __init__(self, number, name, symbol, mass):
+        """Create a new element
+
+        Parameters:
+        number (int) The atomic number of the element
+        name (string) The name of the element
+        symbol (string) The chemical symbol of the element
+        mass (float) The atomic mass of the element
+        """
         ## The atomic number of the element
-        self.atomic_number = number
+        self._atomic_number = number
         ## The name of the element
-        self.name = name
+        self._name = name
         ## The chemical symbol of the element
-        self.symbol = symbol
+        self._symbol = symbol
         ## The atomic mass of the element
-        self.mass = mass
+        self._mass = mass
         # Index this element in a global table
         s = symbol.strip().upper()
+
         assert s not in Element._elements_by_symbol
         Element._elements_by_symbol[s] = self
         if number in Element._elements_by_atomic_number:
@@ -81,8 +93,52 @@ class Element(object):
     def getByAtomicNumber(atomic_number):
         return Element._elements_by_atomic_number[atomic_number]
 
+    @staticmethod
+    def getByMass(mass):
+        """
+        Get the element whose mass is CLOSEST to the requested mass. This method
+        should not be used for repartitioned masses
+        """
+        # Assume masses are in daltons if they are not units
+        if not is_quantity(mass):
+            mass = mass * daltons
+        diff = mass
+        best_guess = None
+
+        for key in Element._elements_by_atomic_number:
+            element = Element._elements_by_atomic_number[key]
+            massdiff = abs(element.mass - mass)
+            if massdiff < diff:
+                best_guess = element
+                diff = massdiff
+
+        return best_guess
+
+    @property
+    def atomic_number(self):
+        return self._atomic_number
+
+    @property
+    def name(self):
+        return self._name
+
+    @property
+    def symbol(self):
+        return self._symbol
+
+    @property
+    def mass(self):
+        return self._mass
+
+    def __str__(self):
+        return '<Element %s>' % self.name
+
+    def __repr__(self):
+        return '<Element %s>' % self.name
+
 # This is for backward compatibility.
 def get_by_symbol(symbol):
+    """ Get the element with a particular chemical symbol. """
     s = symbol.strip().upper()
     return Element._elements_by_symbol[s]
 

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

@@ -267,9 +267,14 @@ class ForceField(object):
             elif self.type == 'outOfPlane':
                 self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
                 self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
+            elif self.type == 'localCoords':
+                self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
+                self.originWeights = [float(attrib['wo1']), float(attrib['wo2']), float(attrib['wo3'])]
+                self.xWeights = [float(attrib['wx1']), float(attrib['wx2']), float(attrib['wx3'])]
+                self.yWeights = [float(attrib['wy1']), float(attrib['wy2']), float(attrib['wy3'])]
+                self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
             else:
                 raise ValueError('Unknown virtual site type: %s' % self.type)
-            self.atoms = [x-self.index for x in self.atoms]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
@@ -326,11 +331,12 @@ class ForceField(object):
                             break
                 if matches is None:
                     raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
+                matchAtoms = dict(zip(matches, res.atoms()))
                 for atom, match in zip(res.atoms(), matches):
                     data.atomType[atom] = template.atoms[match].type
                     for site in template.virtualSites:
                         if match == site.index:
-                            data.virtualSites[atom] = site
+                            data.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms])
 
         # Create the System and add atoms
 
@@ -445,14 +451,20 @@ class ForceField(object):
         # Add virtual sites
 
         for atom in data.virtualSites:
-            site = data.virtualSites[atom]
+            (site, atoms) = data.virtualSites[atom]
             index = atom.index
             if site.type == 'average2':
-                sys.setVirtualSite(index, mm.TwoParticleAverageSite(index+site.atoms[0], index+site.atoms[1], site.weights[0], site.weights[1]))
+                sys.setVirtualSite(index, mm.TwoParticleAverageSite(atoms[0], atoms[1], site.weights[0], site.weights[1]))
             elif site.type == 'average3':
-                sys.setVirtualSite(index, mm.ThreeParticleAverageSite(index+site.atoms[0], index+site.atoms[1], index+site.atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+                sys.setVirtualSite(index, mm.ThreeParticleAverageSite(atoms[0], atoms[1], atoms[2], site.weights[0], site.weights[1], site.weights[2]))
             elif site.type == 'outOfPlane':
-                sys.setVirtualSite(index, mm.OutOfPlaneSite(index+site.atoms[0], index+site.atoms[1], index+site.atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+                sys.setVirtualSite(index, mm.OutOfPlaneSite(atoms[0], atoms[1], atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+            elif site.type == 'localCoords':
+                sys.setVirtualSite(index, mm.LocalCoordinatesSite(atoms[0], atoms[1], atoms[2],
+                                                                  mm.Vec3(site.originWeights[0], site.originWeights[1], site.originWeights[2]),
+                                                                  mm.Vec3(site.xWeights[0], site.xWeights[1], site.xWeights[2]),
+                                                                  mm.Vec3(site.yWeights[0], site.yWeights[1], site.yWeights[2]),
+                                                                  mm.Vec3(site.localPos[0], site.localPos[1], site.localPos[2])))
 
         # Add forces to the System
 
@@ -563,7 +575,7 @@ def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch
     name = atoms[position].name
     for i in range(len(atoms)):
         atom = template.atoms[i]
-        if (atom.element == elem or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
+        if ((atom.element is not None and atom.element == elem) or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
             # See if the bonds for this identification are consistent
 
             allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))

wrappers/python/simtk/openmm/app/internal/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']
+__authors__ = 'Jason Swails'
+__contributors__ = ''
+__license__ = 'MIT'
+__copyright__ = 'Copyright 2014, Jason Swails'
+__date__ = 'Apr. 18, 2014'
+
+
+__private__ = ['topologyobjects', '_charmmfile']

wrappers/python/simtk/openmm/app/internal/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.
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason M. Swails
+Contributors:
+Date: April 18, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+from simtk.openmm.app.internal.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'):
+        self.closed = False
+        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.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):
+        try:
+            self.closed or self._handle.close()
+        except AttributeError:
+            # It didn't make it out of the constructor
+            pass
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+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/internal/charmm/exceptions.py

+"""
+This package contains exceptions that may be raised by the CHARMM components of
+the OpenMM Application layer
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason M. Swails
+Contributors:
+Date: April 18, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+
+class CharmmError(Exception):
+    """ Base class for all exceptions raised in this package """
+
+class CharmmWarning(Warning):
+    """ Base class for all warnings emitted in this package """
+
+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 """
+
+class MoleculeError(CharmmError):
+    """ For (impossibly) messed up connectivity """

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

+"""
+This module contains data structures useful for parsing in CHARMM files and
+constructing chemical structures from those files
+
+This file is part of the OpenMM molecular simulation toolkit originating from
+Simbios, the NIH National Center for Physics-Based Simulation of Biological
+Structures at Stanford, funded under the NIH Roadmap for Medical Research,
+grant U54 GM072970. See https://simtk.org.  This code was originally part of
+the ParmEd program and was ported for use with OpenMM.
+
+Copyright (c) 2014 the Authors
+
+Author: Jason M. Swails
+Contributors:
+Date: April 18, 2014
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+from simtk.openmm.app.internal.charmm.exceptions import (
+                SplitResidueWarning, BondError, ResidueError, CmapError,
+                MissingParameter)
+import simtk.unit as u
+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 * u.daltons
+        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 * u.daltons
+        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()

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

@@ -39,6 +39,7 @@ from simtk.openmm.app.topology import Residue
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LocalEnergyMinimizer
 from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, degree, sqrt, acos, is_quantity, dot, norm
+import simtk.unit as unit
 import element as elem
 import os
 import random
@@ -811,7 +812,7 @@ class Modeller(object):
         else:
             context = Context(system, VerletIntegrator(0.0), platform)
         context.setPositions(newPositions)
-        LocalEnergyMinimizer.minimize(context)
+        LocalEnergyMinimizer.minimize(context, 1.0, 50)
         self.topology = newTopology
         self.positions = context.getState(getPositions=True).getPositions()
         del context
@@ -843,7 +844,9 @@ class Modeller(object):
                     if atom.element is not None:
                         newIndex[i] = index
                         index += 1
-                        newTemplate.atoms.append(ForceField._TemplateAtomData(atom.name, atom.type, atom.element))
+                        newAtom = ForceField._TemplateAtomData(atom.name, atom.type, atom.element)
+                        newAtom.externalBonds = atom.externalBonds
+                        newTemplate.atoms.append(newAtom)
                 for b1, b2 in template.bonds:
                     if b1 in newIndex and b2 in newIndex:
                         newTemplate.bonds.append((newIndex[b1], newIndex[b2]))
@@ -968,7 +971,7 @@ class Modeller(object):
                                 # and hope that energy minimization will fix it.
 
                                 knownPositions = [x for x in templateAtomPositions if x is not None]
-                                position = sum(knownPositions)/len(knownPositions)
+                                position = unit.sum(knownPositions)/len(knownPositions)
                             newPositions.append(position*nanometer)
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:

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

@@ -377,5 +377,5 @@ def _format_83(f):
         return '%8.3f' % f
     if -9999999 < f < 99999999:
         return ('%8.3f' % f)[:8]
-    raise ValueError('coordinate "%s" could not be represnted '
+    raise ValueError('coordinate "%s" could not be represented '
                      'in a width-8 field' % f)

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

@@ -31,8 +31,16 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
-import bz2
-import gzip
+try:
+    import bz2
+    have_bz2 = True
+except: have_bz2 = False
+
+try:
+    import gzip
+    have_gzip = True
+except: have_gzip = False
+
 import simtk.openmm as mm
 import simtk.unit as unit
 import math
@@ -83,8 +91,12 @@ class StateDataReporter(object):
             # Detect the desired compression scheme from the filename extension
             # and open all files unbuffered
             if file.endswith('.gz'):
+                if not have_gzip:
+                    raise RuntimeError("Cannot write .gz file because Python could not import gzip library")
                 self._out = gzip.GzipFile(fileobj=open(file, 'wb', 0))
             elif file.endswith('.bz2'):
+                if not have_bz2:
+                    raise RuntimeError("Cannot write .bz2 file because Python could not import bz2 library")
                 self._out = bz2.BZ2File(file, 'w', 0)
             else:
                 self._out = open(file, 'w', 0)
@@ -189,9 +201,12 @@ class StateDataReporter(object):
         if self._density:
             values.append((self._totalMass/volume).value_in_unit(unit.gram/unit.item/unit.milliliter))
         if self._speed:
-            elapsedDays = (clockTime-self._initialClockTime)/86400
+            elapsedDays = (clockTime-self._initialClockTime)/86400.0
             elapsedNs = (state.getTime()-self._initialSimulationTime).value_in_unit(unit.nanosecond)
+            if elapsedDays > 0.0:
                 values.append('%.3g' % (elapsedNs/elapsedDays))
+            else:
+                values.append('--')
         if self._remainingTime:
             elapsedSeconds = clockTime-self._initialClockTime
             elapsedSteps = simulation.currentStep-self._initialSteps

wrappers/python/src/swig_doxygen/OpenMM.i

@@ -13,10 +13,10 @@ See https://simtk.org/home/pyopenmm for details"
 
 %module (docstring=DOCSTRING) openmm
 
-%include "typemaps.i"
 %include "factory.i"
 %include "std_string.i"
 %include "std_iostream.i"
+%include "typemaps.i"
 
 %include "std_map.i"
 %include "std_pair.i"
@@ -64,7 +64,7 @@ using namespace OpenMM;
 
 %include OpenMM_docstring.i
 
-%include OpenMM_headers.i
+%include OpenMMSwigHeaders.i
 
 %pythoncode %{
   # when we import * from the python module, we only want to import the

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -191,6 +191,10 @@ UNITS = {
 ("*", "getWeight12") : (None, ()),
 ("*", "getWeight13") : (None, ()),
 ("*", "getWeightCross") : (None, ()),
+("LocalCoordinatesSite", "getOriginWeights") : (None, ()),
+("LocalCoordinatesSite", "getXWeights") : (None, ()),
+("LocalCoordinatesSite", "getYWeights") : (None, ()),
+("LocalCoordinatesSite", "getLocalPosition") : ("unit.nanometer", ()),
 ("SerializationNode", "getChildren") : (None, ()),
 ("SerializationNode", "getChildNode") : (None, ()),
 ("SerializationNode", "getProperties") : (None, ()),
@@ -297,6 +301,7 @@ UNITS = {
 ("AmoebaWcaDispersionForce",              "getShctd")                                      :  ( None, ()),
 
 ("Context", "getParameter") : (None, ()),
+("Context", "getMolecules") : (None, ()),
 ("CMAPTorsionForce", "getMapParameters") : (None, ()),
 ("CMAPTorsionForce", "getTorsionParameters") : (None, ()),
 ("CMMotionRemover", "getFrequency") : (None, ()),
@@ -398,6 +403,7 @@ UNITS = {
            'unit.kilojoule_per_mole/(unit.nanometer*unit.nanometer)')),
 ("MonteCarloBarostat", "getFrequency") : (None, ()),
 ("MonteCarloAnisotropicBarostat", "getFrequency") : (None, ()),
+("NonbondedForce", "getPMEParameters") : (None, ('1/unit.nanometer', None, None, None)),
 ("NonbondedForce", "getExceptionParameters")
  : (None, (None, None,
            'unit.elementary_charge*unit.elementary_charge',
@@ -424,5 +430,6 @@ UNITS = {
 ("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()),
 ("DrudeSCFIntegrator", "getMinimizationErrorTolerance") : ("unit.kilojoules_per_mole/unit.nanometer", ()),
 ("RPMDIntegrator", "getContractions") : (None, ()),
+("RPMDIntegrator", "getTotalEnergy") : ("unit.kilojoules_per_mole", ()),
 }
 

wrappers/python/src/swig_doxygen/swig_lib/python/typemaps.i

@@ -117,6 +117,14 @@
   Py_DECREF(args);
 }
 
+%typemap(out) const Vec3& {
+  PyObject* mm = PyImport_AddModule("simtk.openmm");
+  PyObject* vec3 = PyObject_GetAttrString(mm, "Vec3");
+  PyObject* args = Py_BuildValue("(d,d,d)", (*$1)[0], (*$1)[1], (*$1)[2]);
+  $result = PyObject_CallObject(vec3, args);
+  Py_DECREF(args);
+}
+
 /* Convert C++ (Vec3&, Vec3&, Vec3&) object to python tuple or tuples */
 %typemap(argout) (Vec3& a, Vec3& b, Vec3& c) {
   // %typemap(argout) (Vec3& a, Vec3& b, Vec3& c)
@@ -156,3 +164,36 @@
   $3 = &tempC;
 }
 
+%typemap(out) std::string OpenMM::Context::createCheckpoint{
+    // createCheckpoint returns a bytes object
+    $result = PyBytes_FromStringAndSize($1.c_str(), $1.length());
+}
+
+
+%typemap(in) std::string {
+    // if we have a C++ method that takes in a std::string, we're most happy to accept
+    // a python bytes object. But if the user passes in a unicode object we'll try
+    // to recover by encoding it to UTF-8 bytes
+    PyObject* temp = NULL;
+    char* c_str = NULL;
+    Py_ssize_t len = 0;
+
+    if (PyUnicode_Check($input)) {
+        temp = PyUnicode_AsUTF8String($input);
+        if (temp == NULL) {
+            SWIG_exception_fail(SWIG_TypeError, "'utf-8' codec can't decode byte");
+        }
+        PyBytes_AsStringAndSize(temp, &c_str, &len);
+        Py_XDECREF(temp);
+    } else if (PyBytes_Check($input)) {
+        PyBytes_AsStringAndSize($input, &c_str, &len);
+    } else {
+         SWIG_exception_fail(SWIG_TypeError, "argument must be str or bytes");
+    }
+
+    if (c_str == NULL) {
+        SWIG_exception_fail(SWIG_TypeError, "argument must be str or bytes");
+    }
+
+    $1 = std::string(c_str, len);
+}
\ No newline at end of file

wrappers/python/tests/TestBytes.py

+import unittest
+import simtk.openmm as mm
+
+
+class TestBytes(unittest.TestCase):
+    def test_createCheckpoint(self):
+        system = mm.System()
+        system.addParticle(1.0)
+        refPositions = [(0,0,0)]
+
+        context = mm.Context(system, mm.VerletIntegrator(0))
+        context.setPositions(refPositions)
+        chk = context.createCheckpoint()
+        # check that the return value of createCheckpoint is of type bytes (non-unicode)
+        assert isinstance(chk, bytes)
+
+        # set the positions to something random then reload the checkpoint, and
+        # make sure that the positions get restored correctly
+        context.setPositions([(12345, 12345, 123451)])
+        context.loadCheckpoint(chk)
+        newPositions = context.getState(getPositions=True).getPositions()._value
+
+        assert newPositions == refPositions
+
+        # try encoding the checkpoint in utf-8. OpenMM should be able to handle this too
+        context.setPositions([(12345, 12345, 123451)])
+        context.loadCheckpoint(chk.decode('utf-8'))
+        newPositions = context.getState(getPositions=True).getPositions()._value
+
+        assert newPositions == refPositions
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+    

wrappers/python/tests/TestCharmmFiles.py

+import unittest
+from validateConstraints import *
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+import simtk.openmm.app.element as elem
+
+class TestCharmmFiles(unittest.TestCase):
+
+    """Test the GromacsTopFile.createSystem() method."""
+ 
+    def setUp(self):
+        """Set up the tests by loading the input files."""
+
+        # alanine tripeptide; no waters
+        self.psf_c = CharmmPsfFile('systems/ala_ala_ala.psf')
+        self.psf_x = CharmmPsfFile('systems/ala_ala_ala.xpsf')
+        self.psf_v = CharmmPsfFile('systems/ala_ala_ala.vpsf')
+        self.params = CharmmParameterSet(
+                            'systems/charmm22.rtf', 'systems/charmm22.par')
+        self.pdb = PDBFile('systems/ala_ala_ala.pdb')
+
+    def test_NonbondedMethod(self):
+        """Test both non-periodic methods for the systems"""
+
+        methodMap = {NoCutoff:NonbondedForce.NoCutoff, 
+                     CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic}
+        for top in (self.psf_c, self.psf_x, self.psf_v):
+            for method in methodMap:
+                system = top.createSystem(self.params, nonbondedMethod=method)
+                forces = system.getForces()
+                self.assertTrue(any(isinstance(f, NonbondedForce) and 
+                                    f.getNonbondedMethod()==methodMap[method] 
+                                    for f in forces))
+
+    def test_Cutoff(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for top in (self.psf_c, self.psf_x, self.psf_v):
+            for method in [CutoffNonPeriodic]:
+                system = top.createSystem(self.params, nonbondedMethod=method,
+                                          nonbondedCutoff=2*nanometer, 
+                                          constraints=HBonds)
+                cutoff_distance = 0.0*nanometer
+                cutoff_check = 2.0*nanometer
+                for force in system.getForces():
+                    if isinstance(force, NonbondedForce):
+                        cutoff_distance = force.getCutoffDistance()
+                self.assertEqual(cutoff_distance, cutoff_check)
+
+    def test_RemoveCMMotion(self):
+        """Test both options (True and False) for the removeCMMotion parameter."""
+
+        for b in [True, False]:
+            system = self.psf_c.createSystem(self.params, removeCMMotion=b)
+            self.assertEqual(any(isinstance(f, CMMotionRemover) for f in system.getForces()), b)
+
+    def test_ImplicitSolvent(self):
+        """Test implicit solvent using the implicitSolvent parameter.
+
+        """
+        system = self.psf_v.createSystem(self.params, implicitSolvent=OBC2)
+        self.assertTrue(any(isinstance(f, CustomGBForce) for f in system.getForces()))
+
+    def test_ImplicitSolventParameters(self):
+        """Test that solventDielectric and soluteDielectric are passed correctly.
+
+        """
+        system = self.psf_x.createSystem(self.params, implicitSolvent=GBn,
+                                         solventDielectric=50.0, 
+                                         soluteDielectric = 0.9)
+        found_matching_solvent_dielectric=False
+        found_matching_solute_dielectric=False
+        for force in system.getForces():
+            if isinstance(force, CustomGBForce):
+                for i in range(force.getNumGlobalParameters()):
+                    if force.getGlobalParameterName(i) == 'solventDielectric':
+                        if force.getGlobalParameterDefaultValue(i) == 50.0:
+                            found_matching_solvent_dielectric = True
+                    elif force.getGlobalParameterName(i) == 'soluteDielectric':
+                        if force.getGlobalParameterDefaultValue(i) == 0.9:
+                            found_matching_solute_dielectric = True
+            if isinstance(force, NonbondedForce):
+                self.assertEqual(force.getReactionFieldDielectric(), 1.0)
+        self.assertTrue(found_matching_solvent_dielectric and 
+                        found_matching_solute_dielectric)
+
+    def test_HydrogenMass(self):
+        """Test that altering the mass of hydrogens works correctly."""
+        
+        topology = self.psf_v.topology
+        hydrogenMass = 4*amu
+        system1 = self.psf_v.createSystem(self.params)
+        system2 = self.psf_v.createSystem(self.params, hydrogenMass=hydrogenMass)
+        for atom in topology.atoms():
+            if atom.element == elem.hydrogen:
+                self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                self.assertEqual(hydrogenMass, system2.getParticleMass(atom.index))
+        totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
+        totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
+        self.assertAlmostEqual(totalMass1, totalMass2)
+
+if __name__ == '__main__':
+    unittest.main()
+