Rewrite Amber inpcrd/restart file handling. Among the improvements:

    - existence of velocities and box vectors is auto-detected, so
      loadVelocities and loadBoxVectors are unnecessary
    - recognizes NetCDF-format Amber restart files (which can be written by
      tleap, sander, pmemd, and cpptraj).
    - Restart files coordinate/velocity/box fields have a hard width of 12
      characters and does not ensure that fields are whitespace delimited. This
      commit will correctly parse (valid) Amber restart files in which a
      particular coordinate is <= -100 A or >= 1000 A (not uncommon for restart
      files from sander/pmemd).
    - File format (ASCII vs. NetCDF) is auto-detected.

I deprecated the loadVelocities and loadBoxVectors arguments, since the behavior
does not change (and is correct) whether these variables are set or not. A
DeprecationWarning is emitted if users specify these variables.

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

@@ -6,9 +6,9 @@ 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) 2012 Stanford University and the Authors.
+Portions copyright (c) 2012-2014 Stanford University and the Authors.
 Authors: Peter Eastman
-Contributors:
+Contributors: Jason Swails
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -31,51 +31,55 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
+from functools import wraps
 from simtk.openmm.app.internal import amber_file_parser
 from simtk.unit import Quantity, nanometers, picoseconds
+import warnings
 try:
-    import numpy
+    import numpy as np
 except:
-    pass
+    np = None
+
+def numpy_protector(func):
+    """
+    Decorator to emit useful error messages if users try to request numpy
+    processing if numpy is not available. Raises ImportError if numpy could not
+    be found
+    """
+    @wraps(func)
+    def wrapper(self, asNumpy=False):
+        if asNumpy and np is None:
+            raise ImportError('Could not import numpy. Cannot set asNumpy=True')
+        return func(self, asNumpy=asNumpy)
+    return wrapper
 
 class AmberInpcrdFile(object):
     """AmberInpcrdFile parses an AMBER inpcrd file and loads the data stored in it."""
 
-    def __init__(self, file, loadVelocities=False, loadBoxVectors=False):
+    def __init__(self, file, loadVelocities=None, loadBoxVectors=None):
         """Load an inpcrd file.
 
         An inpcrd file contains atom positions and, optionally, velocities and periodic box dimensions.
-        Unfortunately, it is sometimes impossible to determine from the file itself exactly what data
-        it contains.  You therefore must specify in advance what data to load.  It is stored into this
-        object's "positions", "velocities", and "boxVectors" fields.
 
         Parameters:
          - file (string) the name of the file to load
-         - loadVelocities (boolean=False) whether to load velocities from the file
-         - loadBoxVectors (boolean=False) whether to load the periodic box vectors
+         - loadVelocities (boolean=None) deprecated. Velocities are loaded automatically if present
+         - loadBoxVectors (boolean=None) deprecated. Box vectors are loaded automatically if present
         """
-        results = amber_file_parser.readAmberCoordinates(file, read_velocities=loadVelocities, read_box=loadBoxVectors)
-        if loadVelocities:
-            ## The atom positions read from the inpcrd file
-            self.positions = results[0]
-            if loadBoxVectors:
-                ## The periodic box vectors read from the inpcrd file
-                self.boxVectors = results[1]
-                ## The atom velocities read from the inpcrd file
-                self.velocities = results[2]
-            else:
-                self.velocities = results[1]
-        elif loadBoxVectors:
-            self.positions = results[0]
-            self.boxVectors = results[1]
-        else:
-            self.positions = results
+        self.file = file
+        if loadVelocities is not None or loadBoxVectors is not None:
+            warnings.warn('loadVelocities and loadBoxVectors have been '
+                          'deprecated. velocities and box information '
+                          'is loaded automatically if the inpcrd file contains '
+                          'them.', DeprecationWarning)
+        results = amber_file_parser.readAmberCoordinates(file)
+        self.positions, self.velocities, self.boxVectors = results
+        # Cached numpy arrays
         self._numpyPositions = None
-        if loadVelocities:
         self._numpyVelocities = None
-        if loadBoxVectors:
         self._numpyBoxVectors = None
 
+    @numpy_protector
     def getPositions(self, asNumpy=False):
         """Get the atomic positions.
 
@@ -84,34 +88,40 @@ class AmberInpcrdFile(object):
         """
         if asNumpy:
             if self._numpyPositions is None:
-                self._numpyPositions = Quantity(numpy.array(self.positions.value_in_unit(nanometers)), nanometers)
+                self._numpyPositions = Quantity(np.array(self.positions.value_in_unit(nanometers)), nanometers)
             return self._numpyPositions
         return self.positions
 
+    @numpy_protector
     def getVelocities(self, asNumpy=False):
         """Get the atomic velocities.
 
         Parameters:
          - asNumpy (boolean=False) if true, the vectors are returned as numpy arrays instead of Vec3s
         """
+        if self.velocities is None:
+            raise AttributeError('velocities not found in %s' % self.file)
         if asNumpy:
             if self._numpyVelocities is None:
-                self._numpyVelocities = Quantity(numpy.array(self.velocities.value_in_unit(nanometers/picoseconds)), nanometers/picoseconds)
+                self._numpyVelocities = Quantity(np.array(self.velocities.value_in_unit(nanometers/picoseconds)), nanometers/picoseconds)
             return self._numpyVelocities
         return self.velocities
 
+    @numpy_protector
     def getBoxVectors(self, asNumpy=False):
         """Get the periodic box vectors.
 
         Parameters:
          - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
         """
+        if self.boxVectors is None:
+            raise AttributeError('Box information not found in %s' % self.file)
         if asNumpy:
             if self._numpyBoxVectors is None:
                 self._numpyBoxVectors = []
-                self._numpyBoxVectors.append(Quantity(numpy.array(self.boxVectors[0].value_in_unit(nanometers)), nanometers))
-                self._numpyBoxVectors.append(Quantity(numpy.array(self.boxVectors[1].value_in_unit(nanometers)), nanometers))
-                self._numpyBoxVectors.append(Quantity(numpy.array(self.boxVectors[2].value_in_unit(nanometers)), nanometers))
+                self._numpyBoxVectors.append(Quantity(np.array(self.boxVectors[0].value_in_unit(nanometers)), nanometers))
+                self._numpyBoxVectors.append(Quantity(np.array(self.boxVectors[1].value_in_unit(nanometers)), nanometers))
+                self._numpyBoxVectors.append(Quantity(np.array(self.boxVectors[2].value_in_unit(nanometers)), nanometers))
             return self._numpyBoxVectors
         return self.boxVectors
 

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

@@ -40,15 +40,14 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
 #=============================================================================================
 
 import os
-import os.path
 import re
-import math
+from math import ceil, cos, sin, asin, sqrt, pi
 import warnings
 
 try:
-    import numpy
+    import numpy as np
 except:
-    pass
+    np = None
 
 import simtk.unit as units
 import simtk.openmm
@@ -74,6 +73,9 @@ POINTER_LABELS  = """
 # Pointer labels (above) as a list, not string.
 POINTER_LABEL_LIST = POINTER_LABELS.replace(',', '').split()
 
+VELSCALE = 20.455 # velocity conversion factor to angstroms/picosecond
+TINY = 1.0e-8
+
 class PrmtopLoader(object):
     """Parsed AMBER prmtop file.
 
@@ -83,14 +85,14 @@ class PrmtopLoader(object):
 
     Parse a prmtop file of alanine dipeptide in implicit solvent.
 
-    >>> import os, os.path
+    >>> import os
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
     >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
     >>> prmtop = PrmtopLoader(prmtop_filename)
 
     Parse a prmtop file of alanine dipeptide in explicit solvent.
 
-    >>> import os, os.path
+    >>> import os
     >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
     >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
     >>> prmtop = PrmtopLoader(prmtop_filename)
@@ -437,7 +439,7 @@ class PrmtopLoader(object):
                 (rVdwI, epsilonI) = nonbondTerms[iAtom]
                 (rVdwL, epsilonL) = nonbondTerms[lAtom]
                 rMin = (rVdwI+rVdwL)
-                epsilon = math.sqrt(epsilonI*epsilonL)
+                epsilon = sqrt(epsilonI*epsilonL)
                 try:
                     iScee = float(self._raw_data["SCEE_SCALE_FACTOR"][iidx])
                 except KeyError:
@@ -708,7 +710,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
                     l2 = bond[1]
 
             # Compute the distance between atoms and add a constraint
-            length = math.sqrt(l1*l1 + l2*l2 - 2*l1*l2*math.cos(aMin))
+            length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(aMin))
             system.addConstraint(iAtom, kAtom, length)
         if flexibleConstraints or not constrained:
             force.addAngle(iAtom, jAtom, kAtom, aMin, 2*k)
@@ -848,13 +850,13 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             if len(waterO[res]) == 1 and len(waterH[res]) == 2:
                 if len(waterEP[res]) == 1:
                     # Four point water
-                    weightH = distOE[res]/math.sqrt(distOH[res]**2-(0.5*distHH[res])**2)
+                    weightH = distOE[res]/sqrt(distOH[res]**2-(0.5*distHH[res])**2)
                     system.setVirtualSite(waterEP[res][0], mm.ThreeParticleAverageSite(waterO[res][0], waterH[res][0], waterH[res][1], 1-weightH, weightH/2, weightH/2))
                 elif len(waterEP[res]) == 2:
                     # Five point water
-                    weightH = cosOOP*distOE[res]/math.sqrt(distOH[res]**2-(0.5*distHH[res])**2)
-                    angleHOH = 2*math.asin(0.5*distHH[res]/distOH[res])
-                    lenCross = (distOH[res]**2)*math.sin(angleHOH)
+                    weightH = cosOOP*distOE[res]/sqrt(distOH[res]**2-(0.5*distHH[res])**2)
+                    angleHOH = 2*asin(0.5*distHH[res]/distOH[res])
+                    lenCross = (distOH[res]**2)*sin(angleHOH)
                     weightCross = sinOOP*distOE[res]/lenCross
                     system.setVirtualSite(waterEP[res][0], mm.OutOfPlaneSite(waterO[res][0], waterH[res][0], waterH[res][1], weightH/2, weightH/2, weightCross))
                     system.setVirtualSite(waterEP[res][1], mm.OutOfPlaneSite(waterO[res][0], waterH[res][0], waterH[res][1], weightH/2, weightH/2, -weightCross))
@@ -917,23 +919,338 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     return system
 
 #=============================================================================================
-# AMBER INPCRD loader
+# AMBER INPCRD loader classes
 #=============================================================================================
 
-def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbose=False, asNumpy=False):
+class AmberAsciiRestart(object):
+    """
+    Class responsible for parsing Amber coordinates in the ASCII format.
+    Automatically detects the presence of velocities or box parameters in the
+    file.
+
+    Parameters
+    ----------
+    filename : str
+        Name of the restart file
+    asNumpy : bool (False)
+        Load the coordinates, velocities, and box as numpy ndarray objects
+
+    Attributes
+    ----------
+    coordinates : natom x 3 array, Quantity
+        Particle positions with units of length
+    velocities : natom x 3 array, Quantity
+        Particle velocities with units of length per time (None if velocities
+        are not present in the inpcrd file)
+    boxVectors : 3 x 3 array, Quantity
+        Box vectors with units of length (None if no box is present in the
+        inpcrd file)
+    time : float, Quantity
+        Simulation time (None if not present) with units of time
+    title : str
+        Title of the inpcrd file
+    filename : str
+        Name of the file we are parsing
+    natom : int
+        Number of atoms in the inpcrd file
+
+    Raises
+    ------
+        `IOError' if the file does not exist
+        `TypeError' if the format of the file is not recognized
+        `ValueError' if not all fields are numbers (for example, if a field is
+                     filled with ****'s)
+        `IndexError' if the file is empty
+        `ImportError' if numpy is requested but could not be imported
+    Example
+    -------
+    >>> f = AmberAsciiRestart('alanine-dipeptide.inpcrd')
+    >>> coordinates = f.coordinates
+    """
+
+    def __init__(self, filename, asNumpy=False):
+        # Make sure numpy is available if requested
+        if asNumpy and np is None:
+            raise ImportError('asNumpy=True: numpy is not available')
+        self._asNumpy = asNumpy
+        self.filename = filename
+        with open(filename, 'r') as f:
+            lines = f.readlines()
+            # Get rid of trailing blank lines
+            while lines and not lines[-1].strip():
+                lines.pop()
+            self._parse(lines)
+
+    def __str__(self):
+        return self.filename
+
+    def _parse(self, lines):
+        """ Parses through the inpcrd file """
+        global VELSCALE
+        self.title = lines[0].strip()
+        self.time = None
+
+        try:
+            words = lines[1].split()
+            self.natom = int(words[0])
+        except (IndexError, ValueError):
+            raise TypeError('Unrecognized file type [%s]' % self.filename)
+        
+        if len(words) >= 2:
+            self.time = float(words[1]) * units.picoseconds
+
+        if len(lines) == int(ceil(self.natom / 2.0) + 2):
+            hasbox = hasvels = False
+            self.boxVectors = self.velocities = None
+        elif self.natom in (1, 2) and len(lines) == 4:
+            # This is the _only_ case where line counting does not work -- there
+            # is either 1 or 2 atoms and there are 4 lines. The 1st 3 lines are
+            # the title, natom/time, and coordinates. The 4th are almost always
+            # velocities since Amber does not make it easy to make a periodic
+            # system with only 2 atoms. If natom is 1, the 4th line is either a
+            # velocity (3 #'s) or a box (6 #'s). If natom is 2, it is a bit
+            # ambiguous. However, velocities (which are scaled by 20.445) have a
+            # ~0% chance of being 60+, so we can pretty easily tell if the last
+            # line has box dimensions and angles or velocities. I cannot
+            # envision a _plausible_ scenario where the detection here will fail
+            # in real life.
+            line = lines[3]
+            if self.natom == 1:
+                tmp = [line[i:i+12] for i in range(0, 72, 12) if line[i:i+12]]
+                if len(tmp) == 3:
+                    hasvels = True
+                    hasbox = False
+                    self.boxVectors = False
+                elif len(tmp) == 6:
+                    hasbox = True
+                    hasvels = False
+                    self.velocities = None
+                else:
+                    raise TypeError('Unrecognized line in restart file %s' %
+                                    self.filename)
+            else:
+                # Ambiguous case
+                tmp = [float(line[i:i+12]) >= 60.0 for i in range(0, 72, 12)]
+                if any(tmp):
+                    hasbox = True
+                    hasvels = False
+                    self.velocities = False
+                else:
+                    hasvels = True
+                    hasbox = False
+                    self.boxVectors = False
+        elif len(lines) == int(ceil(self.natom / 2.0) + 3):
+            hasbox = True
+            hasvels = False
+            self.velocities = None
+        elif len(lines) == int(2 * ceil(self.natom / 2.0) + 2):
+            hasbox = False
+            self.boxVectors = None
+            hasvels = True
+        elif len(lines) == int(2 * ceil(self.natom / 2.0) + 3):
+            hasbox = hasvels = True
+        else:
+            raise TypeError('Badly formatted restart file. Has %d lines '
+                            'for %d atoms.' % (len(self.lines), self.natom))
+
+        if self._asNumpy:
+            coordinates = np.zeros((self.natom, 3), np.float32)
+            if hasvels:
+                velocities = np.zeros((self.natom, 3), np.float32)
+            if hasbox:
+                boxVectors = np.zeros((3, 3), np.float32)
+        else:
+            coordinates = [[0.0, 0.0, 0.0] for i in range(self.natom)]
+            if hasvels:
+                velocities = [[0.0, 0.0, 0.0] for i in range(self.natom)]
+            if hasbox:
+                boxVectors = [[0.0, 0.0, 0.0] for i in range(3)]
+
+        # Now it's time to parse.  Coordinates first
+        startline = 2
+        endline = startline + int(ceil(self.natom / 2.0))
+        idx = 0
+        for i in range(startline, endline):
+            line = lines[i]
+            coordinates[idx][0] = float(line[ 0:12])
+            coordinates[idx][1] = float(line[12:24])
+            coordinates[idx][2] = float(line[24:36])
+            idx += 1
+            if idx < self.natom:
+                coordinates[idx][0] = float(line[36:48])
+                coordinates[idx][1] = float(line[48:60])
+                coordinates[idx][2] = float(line[60:72])
+                idx += 1
+        self.coordinates = units.Quantity(coordinates, units.angstroms)
+        startline = endline
+        # Now it's time to parse velocities if we have them
+        if hasvels:
+            endline = startline + int(ceil(self.natom / 2.0))
+            idx = 0
+            for i in range(startline, endline):
+                line = lines[i]
+                velocities[idx][0] = float(line[ 0:12]) * VELSCALE
+                velocities[idx][1] = float(line[12:24]) * VELSCALE
+                velocities[idx][2] = float(line[24:36]) * VELSCALE
+                idx += 1
+                if idx < self.natom:
+                    velocities[idx][0] = float(line[36:48]) * VELSCALE
+                    velocities[idx][1] = float(line[48:60]) * VELSCALE
+                    velocities[idx][2] = float(line[60:72]) * VELSCALE
+                    idx += 1
+            startline = endline
+            self.velocities = units.Quantity(velocities,
+                                             units.angstroms/units.picoseconds)
+        if hasbox:
+            line = lines[startline]
+            try:
+                tmp = [float(line[i:i+12]) for i in range(0, 72, 12)]
+            except (IndexError, ValueError):
+                raise ValueError('Could not parse box line in %s' %
+                                 self.filename)
+            lengths = tmp[:3]
+            angles = tmp[3:]
+            _box_vectors_from_lengths_angles(lengths, angles, boxVectors)
+            self.boxVectors = units.Quantity(boxVectors, units.angstroms)
+
+class AmberNetcdfRestart(object):
+    """
+    Amber restart/inpcrd file in the NetCDF format (full double-precision
+    coordinates, velocities, and unit cell parameters). Reads NetCDF restarts
+    written by LEaP and pmemd/sander. Requires scipy to parse NetCDF files.
+
+    Parameters
+    ----------
+    filename : str
+        Name of the restart file
+    asNumpy : bool (False)
+        Load the coordinates, velocities, and box as numpy ndarray objects
+
+    Attributes
+    ----------
+    coordinates : natom x 3 array, Quantity
+        Particle positions with units of length
+    velocities : natom x 3 array, Quantity
+        Particle velocities with units of length per time (None if velocities
+        are not present in the inpcrd file)
+    boxVectors : 3 x 3 array, Quantity
+        Box vectors with units of length (None if no box is present in the
+        inpcrd file)
+    time : float, Quantity
+        Simulation time (None if not present) with units of time
+    title : str
+        Title of the inpcrd file
+    filename : str
+        Name of the file we are parsing
+    natom : int
+        Number of atoms in the inpcrd file
+
+    Raises
+    ------
+        `IOError' if the file does not exist
+        `TypeError' if the file is not a NetCDF v3 file
+        `ImportError' if scipy is not available
+    Example
+    -------
+    >>> f = AmberNetcdfRestart('alanine-dipeptide.ncrst')
+    >>> coordinates = f.coordinates
+    """
+    def __init__(self, filename, asNumpy=False):
+        try:
+            from scipy.io.netcdf import NetCDFFile
+        except ImportError:
+            raise ImportError('scipy is necessary to parse NetCDF restarts')
+        
+        self.filename = filename
+        self.velocities = self.boxVectors = self.time = None
+
+        # Extract the information from the NetCDF file. We need to make copies
+        # here because the NetCDF variables are mem-mapped, but is only mapped
+        # to valid memory while the file handle is open. Since the context
+        # manager GCs the ncfile handle, the memory for the original variables
+        # is no longer valid. So copy those arrays while the handle is still
+        # open. This is unnecessary in scipy v.0.12 and lower because NetCDFFile
+        # accidentally leaks the file handle, but that was 'fixed' in 0.13. This
+        # fix taken from MDTraj
+        with NetCDFFile(filename, 'r') as ncfile:
+            self.coordinates = np.array(ncfile.variables['coordinates'][:])
+            if 'velocities' in ncfile.variables:
+                vels = ncfile.variables['velocities']
+                self.velocities = np.array(vels[:]) * vels.scale_factor
+            if ('cell_lengths' in ncfile.variables and
+                'cell_angles' in ncfile.variables):
+                self.boxVectors = np.zeros((3,3), np.float32)
+                _box_vectors_from_lengths_angles(
+                        ncfile.variables['cell_lengths'][:],
+                        ncfile.variables['cell_angles'][:],
+                        self.boxVectors,
+                )
+            if 'time' in ncfile.variables:
+                self.time = ncfile.variables['time'].getValue()
+
+        # They are already numpy -- convert to list if we don't want numpy
+        if not asNumpy:
+            self.coordinates = self.coordinates.tolist()
+            if self.velocities is not None:
+                self.velocities = self.velocities.tolist()
+            if self.boxVectors is not None:
+                self.boxVectors = self.boxVectors.tolist()
+
+        # Now add the units
+        self.coordinates = units.Quantity(self.coordinates, units.angstroms)
+        if self.velocities is not None:
+            self.velocities = units.Quantity(self.velocities,
+                                             units.angstroms/units.picoseconds)
+        if self.boxVectors is not None:
+            self.boxVectors = units.Quantity(self.boxVectors, units.angstroms)
+        self.time = units.Quantity(self.time, units.picosecond)
+
+def _box_vectors_from_lengths_angles(lengths, angles, boxVectors):
+    """
+    Converts lengths and angles into a series of box vectors and modifies
+    boxVectors in-place (it must be a mutable sequence)
+
+    Parameters
+    ----------
+    lengths : 3-element array of floats
+        Lengths of the 3 periodic box vectors
+    angles : 3-element array of floats
+        Angles (in degrees) between the 3 periodic box vectors
+    boxVectors : mutable 3x3 sequence
+    """
+    alpha = angles[0] * pi / 180.0
+    beta = angles[1] * pi / 180.0
+    gamma = angles[2] * pi / 180.0
+
+    boxVectors[0][0] = lengths[0]
+
+    boxVectors[1][0] = lengths[1] * cos(gamma)
+    boxVectors[1][1] = lengths[1] * sin(gamma)
+
+    boxVectors[2][0] = cx = lengths[2] * cos(beta)
+    boxVectors[2][1] = cy = lengths[2] * (cos(alpha) - cos(beta) * cos(gamma))
+    boxVectors[2][2] = sqrt(lengths[2]*lengths[2] - cx*cx - cy*cy)
+
+    boxVectors[0][1] = boxVectors[0][2] = boxVectors[1][2] = 0.0
+
+def readAmberCoordinates(filename, asNumpy=False):
     """
     Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.
 
     ARGUMENTS
 
     filename (string) - name of Amber coordinates file to be read in
-    system (simtk.openmm.System) - System object for which coordinates are to be read
 
     OPTIONAL ARGUMENTS
 
-    verbose (boolean) - if True, will print out verbose information about the file being read
     asNumpy (boolean) - if True, results will be returned as Numpy arrays instead of lists of Vec3s
 
+    RETURNS
+
+    coordinates, velocities, boxVectors
+        The velocities and boxVectors will be None if they are not found in the
+        restart file
+
     EXAMPLES
 
     Read coordinates in vacuum.
@@ -950,101 +1267,45 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
 
     """
 
-    # Open coordinate file for reading.
-    infile = open(filename, 'r')
-
-    # Read title
-    title = infile.readline().strip()
-    if verbose: print "title: '%s'" % title
-
-    # Read number of atoms
-    natoms = int(infile.readline().split()[0])
-    if verbose: print "%d atoms" % natoms
-
-    # Allocate storage for coordinates
-    coordinates = []
-
-    # Read coordinates
-    mm = simtk.openmm
-    natoms_read = 0
-    while (natoms_read < natoms):
-        line = infile.readline()
-        if len(line) == 0:
-            raise ValueError("Unexpected end of file while reading coordinates")
-        line = line.strip()
-        elements = line.split()
-        while (len(elements) > 0):
-            coordinates.append(mm.Vec3(float(elements.pop(0)), float(elements.pop(0)), float(elements.pop(0))))
-            natoms_read += 1
-    if asNumpy:
-        newcoords = numpy.zeros([natoms,3], numpy.float32)
-        for i in range(len(coordinates)):
-            for j in range(3):
-                newcoords[i,j] = coordinates[i][j]
-        coordinates = newcoords
-    # Assign units.
-    coordinates = units.Quantity(coordinates, units.angstroms)
-
-    # Read velocities if requested.
-    velocities = None
-    if (read_velocities):
-        # Read velocities
-        velocities = []
-        natoms_read = 0
-        while (natoms_read < natoms):
-            line = infile.readline()
-            if len(line) == 0:
-                raise ValueError("Unexpected end of file while reading velocities")
-            line = line.strip()
-            elements = line.split()
-            while (len(elements) > 0):
-                velocities.append(20.455*mm.Vec3(float(elements.pop(0)), float(elements.pop(0)), float(elements.pop(0))))
-                natoms_read += 1
-        if asNumpy:
-            newvel = numpy.zeros([natoms,3], numpy.float32)
-            for i in range(len(velocities)):
-                for j in range(3):
-                    newvel[i,j] = velocities[i][j]
-            velocities = newvel
-        # Assign units.
-        velocities = units.Quantity(velocities, units.angstroms/units.picoseconds)
-
-    # Read box size if present
-    box_vectors = None
-    if (read_box):
-        line = infile.readline()
-        if len(line) == 0:
-            raise ValueError("Unexpected end of file while reading box vectors")
-        line = line.strip()
-        elements = line.split()
-        nelements = len(elements)
-        box_dimensions = [0.0]*nelements
-        for i in range(nelements):
-            box_dimensions[i] = float(elements[i])
-        # TODO: Deal with non-standard box sizes.
-        if nelements == 6:
-            if asNumpy:
-                a = units.Quantity(numpy.array([box_dimensions[0], 0.0, 0.0]), units.angstroms)
-                b = units.Quantity(numpy.array([0.0, box_dimensions[1], 0.0]), units.angstroms)
-                c = units.Quantity(numpy.array([0.0, 0.0, box_dimensions[2]]), units.angstroms)
-            else:
-                a = units.Quantity(mm.Vec3(box_dimensions[0], 0.0, 0.0), units.angstroms)
-                b = units.Quantity(mm.Vec3(0.0, box_dimensions[1], 0.0), units.angstroms)
-                c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)
-            box_vectors = [a,b,c]
-        else:
-            raise Exception("Don't know what to do with box vectors: %s" % line)
-
-    # Close file
-    infile.close()
-
-    if box_vectors and velocities:
-        return (coordinates, box_vectors, velocities)
-    if box_vectors:
-        return (coordinates, box_vectors)
-    if velocities:
-        return (coordinates, velocities)
-    return coordinates
+    try:
+        crdfile = AmberNetcdfRestart(filename)
+    except ImportError:
+        # See if it's an ASCII file.  If so, no need to complain
+        try:
+            crdfile = AmberAsciiRestart(filename)
+            """
+        `TypeError' if the format of the file is not recognized
+        `ValueError' if not all fields are numbers (for example, if a field is
+                     filled with ****'s)
+        `IndexError' if the file is empty
+        `ImportError' if numpy is requested but could not be imported
+        """
+        except TypeError:
+            raise TypeError('Problem parsing %s as an ASCII Amber restart file '
+                            'and scipy could not be imported to try reading as '
+                            'a NetCDF restart file.' % filename)
+        except (IndexError, ValueError):
+            raise TypeError('Could not parse Amber ASCII restart file %s' %
+                            filename)
+        except ImportError:
+            raise ImportError('Could not find numpy; cannot use asNumpy=True')
+    except TypeError:
+        # We had scipy, but this is not a NetCDF v3 file. Try as ASCII now
+        try:
+            crdfile = AmberAsciiRestart(filename)
+        except TypeError:
+            raise TypeError('Problem parsing %s as an ASCII Amber restart file'
+                            % filename)
+        except (IndexError, ValueError):
+            raise TypeError('Could not parse Amber ASCII restart file %s' %
+                            filename)
+        # Import error cannot happen, since we had scipy which has numpy as a
+        # prereq. Do not catch that exception (only catch what you intend to
+        # catch...)
+
+    # We got here... one of the file types worked. Return the coordinates,
+    # velocities, and boxVectors
+    return crdfile.coordinates, crdfile.velocities, crdfile.boxVectors
 
 #=============================================================================================
 # MAIN AND TESTS