Merge pull request #461 from swails/amber_restarts

Rewrite Amber inpcrd file parsers

wrappers/python/CMakeLists.txt

@@ -55,8 +55,12 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.dms"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.top"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.par"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*psf"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/charmm22.*"
     )

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,87 +31,97 @@ 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
+        self._numpyVelocities = None
+        self._numpyBoxVectors = None
 
+    @numpy_protector
     def getPositions(self, asNumpy=False):
         """Get the atomic positions.
 
         Parameters:
          - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
-         """
+        """
         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/tests/TestAmberInpcrdFile.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
+
+def compareByElement(array1, array2, cmp):
+    for x, y in zip(array1, array2):
+        cmp(x, y)
+
+class TestAmberInpcrdFile(unittest.TestCase):
+    """Test the Amber inpcrd file parser"""
+ 
+    def test_CrdVelBox(self):
+        """ Test parsing ASCII restarts with crds, vels, and box """
+        cmp = self.assertAlmostEqual
+        inpcrd = AmberInpcrdFile('systems/crds_vels_box.rst7')
+        self.assertEqual(len(inpcrd.positions), 2101)
+        compareByElement(inpcrd.positions[-1].value_in_unit(angstroms),
+                         [3.5958082, 8.4176792, -8.2954064], cmp)
+        compareByElement(inpcrd.velocities[-1].value_in_unit(angstroms/picoseconds),
+                         [0.3091332*20.455, 0.7355925*20.455, -0.1206961*20.455], cmp)
+        compareByElement(inpcrd.boxVectors[0].value_in_unit(angstroms),
+                         [30.2642725, 0.0, 0.0], cmp)
+
+    def test_NetCDF(self):
+        """ Test NetCDF restart file parsing """
+        cmp = self.assertAlmostEqual
+        try:
+            from scipy.io import netcdf
+        except ImportError:
+            print('Not testing NetCDF file parser... scipy cannot be found')
+        else:
+            inpcrd = AmberInpcrdFile('systems/amber.ncrst')
+            self.assertEqual(len(inpcrd.positions), 2101)
+            compareByElement(inpcrd.positions[0].value_in_unit(angstroms),
+                             [6.82122492718229, 6.6276250662042, -8.51668999892245],
+                             cmp)
+            compareByElement(inpcrd.velocities[-1].value_in_unit(angstroms/picosecond),
+                             [0.349702202733541*20.455, 0.391525333168534*20.455,
+                              0.417941679767662*20.455], cmp)
+            self.assertAlmostEqual(inpcrd.boxVectors[0][0].value_in_unit(angstroms),
+                                   30.2642725, places=6)
+
+    def test_CrdBox(self):
+        """ Test parsing ASCII restarts with only crds and box """
+        inpcrd = AmberInpcrdFile('systems/crds_box.rst7')
+        self.assertEqual(len(inpcrd.positions), 18660)
+        self.assertTrue(inpcrd.velocities is None)
+        self.assertTrue(inpcrd.boxVectors is not None)
+
+    def test_CrdVel(self):
+        inpcrd = AmberInpcrdFile('systems/crds_vels.rst7')
+        self.assertTrue(inpcrd.boxVectors is None)
+        self.assertTrue(inpcrd.velocities is not None)
+
+    def test_CrdOnly(self):
+        inpcrd = AmberInpcrdFile('systems/crdsonly.rst7')
+        self.assertTrue(inpcrd.boxVectors is None)
+        self.assertTrue(inpcrd.velocities is None)
+
+if __name__ == '__main__':
+    unittest.main()

wrappers/python/tests/systems/crds_vels.rst7

+ACE                                                                             
+   28  0.1100000E+05
+   6.9962325   7.1345361  -3.9713983   7.4756482   6.1688310  -3.8111566
+   8.5024408   6.3930528  -3.5221598   7.4708173   5.7449547  -4.8153505
+   6.8002949   5.4664470  -2.6846071   7.1210336   4.3039353  -2.4099316
+   5.8138412   6.2015113  -2.0369948   5.6437332   7.0956281  -2.4748521
+   4.9765730   5.7180343  -0.9552686   5.5952335   5.0377318  -0.3699927
+   4.5105579   6.8291112   0.0795226   5.4196673   7.2692444   0.4892821
+   3.8405864   7.5537576  -0.3832260   3.7600399   6.2095975   1.2450279
+   3.0614461   7.0866468   2.0205499   3.9472641   4.9020081   1.4914431
+   3.4219922   4.6877864   2.2659097   3.7335795   4.9588789  -1.4954454
+   2.9265073   5.7168640  -2.1242842   4.4787994   4.4294475   0.8466484
+   2.5980237   6.7854333   2.8054767   3.2160253   7.9534080   1.6378901
+   3.5172878   3.6766505  -1.2008776   4.1472250   3.2061006  -0.5669477
+   2.4677851   2.8846866  -1.7794578   2.8829398   2.5139385  -2.7166333
+   2.3438054   2.0729996  -1.0625940   1.5440133   3.4553999  -1.8744943
+   0.3799223   0.0175576   1.1999574   0.3523292  -0.2043193   0.0353460
+   1.1942188  -1.4410315  -1.7982354  -0.5852708  -0.9087781   0.3223661
+   0.1459748   0.0526223   0.2682160  -0.0465667  -0.1870169   0.1459932
+  -0.2892203   0.0015456  -0.0112940   0.1854616  -0.2737819  -0.7785389
+  -0.2392158   0.1749611  -0.1239180  -0.3400171   0.0531231  -0.1587713
+  -0.4704835  -0.1493877  -0.2567553  -0.6101138  -0.0392770  -0.0635170
+   0.5715057   0.1942655  -1.2772686  -0.1237689  -0.3449472   0.2714899
+  -0.0185239  -0.0150607   0.0391274  -0.1280282  -0.0046097  -0.1994271
+  -0.3155035   0.3668607  -0.2215361   0.1294219  -0.1747917  -0.0917708
+   0.0742998   0.1509489  -0.2851266   0.4923665  -1.0386246   1.0229978
+  -0.8227836   1.0624132   0.0013427  -0.0778814  -0.1859122  -0.3774084
+   0.2098003   0.1213628  -0.2466628  -0.1936209  -0.2478820  -0.1148687
+   0.0315637   0.2866246  -0.0746640  -0.5854400  -0.8935872   0.0990994
+   0.6063127  -0.1110186  -0.4168450  -0.8390447  -1.1240213  -0.3891693

wrappers/python/tests/systems/crdsonly.rst7

+ACE 
+   28
+   1.8662825   1.2455949   0.6672110   1.9626006   2.1938026   0.1396779
+   1.2497428   2.9104392   0.5454310   1.7619492   2.0427757  -0.9201889
+   3.3659937   2.7323333   0.3106630   4.1907627   2.1039510   0.9734519
+   3.6273756   3.8949435  -0.2849918   2.8777368   4.3392797  -0.7950028
+   4.9163808   4.5976032  -0.2534247   5.4756521   4.2709258   0.6250439
+   5.7163632   4.2049342  -1.5101182   5.8371260   3.1200826  -1.5317512
+   5.1471196   4.4947601  -2.3959905   7.1055802   4.8539661  -1.5508595
+   7.9498622   4.6375603  -0.4939507   7.4316586   5.6201081  -2.6386852
+   8.3150784   5.9901727  -2.5863945   4.7269796   6.1275619  -0.1377803
+   3.7336230   6.6808666  -0.6238468   6.7262793   5.6692474  -3.2868494
+   8.7967113   5.0771061  -0.5933179   7.5650193   4.0788360   0.1842201
+   5.6831501   6.8163892   0.4999151   6.4763527   6.2959896   0.8529456
+   5.6745238   8.2651257   0.7063732   4.7905547   8.5573927   1.2769756
+   6.5663801   8.5705282   1.2568739   5.6602885   8.7808127  -0.2560452

wrappers/python/tests/systems/water-box-216.prmtop

@@ -1054,7 +1054,7 @@ BLA BLA BLA BLA BLA BLA BLA BLA
        3       3       3       3       3       3
 %FLAG BOX_DIMENSIONS                                                            
 %FORMAT(5E16.8)                                                                 
-  5.15661795E+03  1.87743490E+01  1.87743490E+01  1.87743490E+01
+  9.00000000E+01  1.87743490E+01  1.87743490E+01  1.87743490E+01
 %FLAG RADIUS_SET                                                                
 %FORMAT(1a80)                                                                   
 modified Bondi radii (mbondi)