Merge pull request #406 from swails/master

Add support for native CHARMM files

wrappers/python/CMakeLists.txt

@@ -57,6 +57,8 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.dms"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.top"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*psf"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/charmm22.*"
     )
     foreach(file ${STAGING_INPUT_FILES1})
         set(STAGING_INPUT_FILES ${STAGING_INPUT_FILES} "${file}")

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

@@ -26,6 +26,9 @@ 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/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):
@@ -93,6 +93,27 @@ 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
@@ -117,6 +138,7 @@ class Element(object):
 
 # 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/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 """