Merge pull request #294 from richardjgowers/feature-newtimestepconstr…

…uctor

Reworked coordinates.base.Timestep construction

package/CHANGELOG

@@ -17,11 +17,18 @@ The rules for this file:
   Enhancements
   * Added reading of DL_Poly format CONFIG and HISTORY files (Issue #298)
     These can both act as both Topology and Coordinate information
+  * Timestep objects now have __eq__ method
+  * coordinates.base.Timestep now can handle velocities and forces
   * Waterdynamics analysis module added, including five analysis classes: Hydrogen Bond
     Lifetimes, Water Orientational Relaxation, Angular Distribution, Mean Square 
     Displacement and Survival Probability. Documentation and test are included.
 
   Changes
+  * Timestep can now only init using an integer argument (which
+    represents the number of atoms)
+  * Added from_timestep and from_coordinates construction methods
+    to base.Timestep
+  * Removed FormatError, now replaced by ValueError
   * distances.contact_matrix(): changed keyword 'suppress_progmet' to
     'quiet'
 

package/MDAnalysis/__init__.py

@@ -206,10 +206,6 @@ class NoDataError(ValueError):
     """Raised when empty input is not allowed or required data are missing."""
 
 
-class FormatError(EnvironmentError):
-    """Raised when there appears to be a problem with format of input files."""
-
-
 class ApplicationError(OSError):
     """Raised when an external application failed.
 

package/MDAnalysis/coordinates/CRD.py

@@ -23,15 +23,17 @@
 
 """
 
-import numpy
+import numpy as np
 
 import MDAnalysis.core.util as util
 from . import base
-from MDAnalysis import FormatError
 
 
 class CRDReader(base.SingleFrameReader):
     """CRD reader that implements the standard and extended CRD coordinate formats
+
+    .. versionchanged:: 0.11.0
+       Now returns a ValueError instead of FormatError
     """
     format = 'CRD'
     units = {'time': None, 'length': 'Angstrom'}
@@ -61,22 +63,23 @@ def _read_first_frame(self):
                 # process coordinates
                 try:
                     if extended:
-                        coords_list.append(numpy.array(map(float, line[45:100].split()[0:3])))
+                        coords_list.append(np.array(map(float, line[45:100].split()[0:3])))
                     else:
-                        coords_list.append(numpy.array(map(float, line[20:50].split()[0:3])))
+                        coords_list.append(np.array(map(float, line[20:50].split()[0:3])))
                 except:
-                    raise FormatError("Check CRD format at line %d: %s" % (linenum, line.rstrip()))
+                    raise ValueError("Check CRD format at line {0}: {1}"
+                                     "".format(linenum, line.rstrip()))
 
         self.numatoms = len(coords_list)
 
-        self.ts = self._Timestep(numpy.array(coords_list))
+        self.ts = self._Timestep.from_coordinates(np.array(coords_list))
         self.ts.frame = 1  # 1-based frame number
         # if self.convert_units:
         #    self.convert_pos_from_native(self.ts._pos)             # in-place !
 
         # sanity check
         if self.numatoms != natoms:
-            raise FormatError("Found %d coordinates in %r but the header claims that there "
+            raise ValueError("Found %d coordinates in %r but the header claims that there "
                               "should be %d coordinates." % (self.numatoms, self.filename, natoms))
 
     def Writer(self, filename, **kwargs):

package/MDAnalysis/coordinates/DCD.py

@@ -68,7 +68,7 @@
 """
 import os
 import errno
-import numpy
+import numpy as np
 import struct
 
 from . import base
@@ -131,15 +131,15 @@ def dimensions(self):
         # override for other formats; this strange ordering is kept for historical reasons
         # (the user should not need concern themselves with this)
         ## orig MDAnalysis 0.8.1/dcd.c (~2004)
-        ##return numpy.take(self._unitcell, [0,2,5,1,3,4])
+        ##return np.take(self._unitcell, [0,2,5,1,3,4])
 
         # MDAnalysis 0.9.0 with recent dcd.c (based on 2013 molfile
         # DCD plugin, which implements the ordering of recent NAMD
         # (>2.5?)). See Issue 187.
-        uc = numpy.take(self._unitcell, self._ts_order)
+        uc = np.take(self._unitcell, self._ts_order)
         # heuristic sanity check: uc = A,B,C,alpha,beta,gamma
         # XXX: should we worry about these comparisons with floats?
-        if numpy.any(uc < 0.) or numpy.any(uc[3:] > 180.):
+        if np.any(uc < 0.) or np.any(uc[3:] > 180.):
             # might be new CHARMM: box matrix vectors
             H = self._unitcell
             e1, e2, e3 = H[[0,1,3]],  H[[1,2,4]], H[[3,4,5]]
@@ -153,7 +153,7 @@ def dimensions(self, box):
         .. versionadded:: 0.9.0
         """
         # note that we can re-use self._ts_order with put!
-        numpy.put(self._unitcell, self._ts_order, box)
+        np.put(self._unitcell, self._ts_order, box)
 
 
 class DCDWriter(base.Writer):
@@ -311,9 +311,9 @@ def write_next_timestep(self, ts=None):
                 ts = self.ts
         elif not ts.numatoms == self.numatoms:
             raise ValueError("DCDWriter: Timestep does not have the correct number of atoms")
-        unitcell = self.convert_dimensions_to_unitcell(ts).astype(numpy.float32)  # must be float32 (!)
+        unitcell = self.convert_dimensions_to_unitcell(ts).astype(np.float32)  # must be float32 (!)
         if not ts._pos.flags.f_contiguous:  # Not in fortran format
-            ts = Timestep(ts)  # wrap in a new fortran formatted Timestep
+            ts = Timestep.from_timestep(ts)  # wrap in a new fortran formatted Timestep
         if self.convert_units:
             pos = self.convert_pos_to_native(ts._pos,
                                              inplace=False)  # possibly make a copy to avoid changing the trajectory
@@ -337,17 +337,17 @@ def convert_dimensions_to_unitcell(self, ts, _ts_order=Timestep._ts_order):
         # unitcell is A,B,C,alpha,beta,gamma - convert to order expected by low level DCD routines
         lengths, angles = ts.dimensions[:3], ts.dimensions[3:]
         self.convert_pos_to_native(lengths)
-        unitcell = numpy.zeros_like(ts.dimensions)
+        unitcell = np.zeros_like(ts.dimensions)
         # __write_DCD_frame() wants uc_array = [A, gamma, B, beta, alpha, C] and
         # will write the lengths and the angle cosines to the DCD. NOTE: It
         # will NOT write CHARMM36+ box matrix vectors. When round-tripping a
         # C36+ DCD, we loose the box vector information. However, MDAnalysis
         # itself will not detect this because the DCDReader contains the
         # heuristic to deal with either lengths/angles or box matrix on the fly.
         #
-        # use numpy.put so that we can re-use ts_order (otherwise we
+        # use np.put so that we can re-use ts_order (otherwise we
         # would need a ts_reverse_order such as _ts_reverse_order = [0, 5, 1, 4, 3, 2])
-        numpy.put(unitcell, _ts_order, numpy.concatenate([lengths, angles]))
+        np.put(unitcell, _ts_order, np.concatenate([lengths, angles]))
         return unitcell
 
     def close(self):
@@ -485,7 +485,7 @@ def _read_next_timestep(self, ts=None):
     def _read_frame(self, frame):
         self._jump_to_frame(frame)
         ts = self.ts
-        ts.frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1)  # XXX required!!
+        ts.frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1)
         return ts
 
     def timeseries(self, asel, start=0, stop=-1, skip=1, format='afc'):

package/MDAnalysis/coordinates/DLPoly.py

@@ -49,7 +49,7 @@ def dimensions(self, new):
 class ConfigReader(base.SingleFrameReader):
     """DLPoly Config file Reader
 
-    .. versionadded:: 0.10.1
+    .. versionadded:: 0.11.0
     """
     format = 'CONFIG'
     units = _DLPOLY_UNITS
@@ -116,8 +116,9 @@ def _read_first_frame(self):
             if has_forces:
                 forces = forces[order]
 
-        # TODO: Merge with "new" style Timestep when finished
-        ts = self.ts = self._Timestep(self.numatoms)
+        ts = self.ts = self._Timestep(self.numatoms,
+                                      velocities=has_vels,
+                                      forces=has_forces)
         ts._pos = coords
         if has_vels:
             ts._velocities = velocities
@@ -133,7 +134,7 @@ def _read_first_frame(self):
 class HistoryReader(base.Reader):
     """Reads DLPoly format HISTORY files
 
-    .. versionadded:: 0.10.1
+    .. versionadded:: 0.11.0
     """
     format = 'HISTORY'
     units = _DLPOLY_UNITS
@@ -157,13 +158,12 @@ def __init__(self, filename, convert_units=None, **kwargs):
         self._file = open(self.filename, 'r')
         self.title = self._file.readline().strip()
         self._levcfg, self._imcon, self.numatoms = map(int, self._file.readline().split()[:3])
+        self._has_vels = True if self._levcfg > 0 else False
+        self._has_forces = True if self._levcfg == 2 else False
 
-        # TODO: Replace with new style Timestep
-        self.ts = self._Timestep(self.numatoms)
-        if self._levcfg > 0:
-            self.ts._velocities = np.zeros((self.numatoms, 3), dtype=np.float32, order='F')
-        if self._levcfg == 2:
-            self.ts._forces = np.zeros((self.numatoms, 3), dtype=np.float32, order='F')
+        self.ts = self._Timestep(self.numatoms,
+                                 velocities=self._has_vels,
+                                 forces=self._has_forces)
         self._read_next_timestep()
 
     def _read_next_timestep(self, ts=None):
@@ -193,9 +193,9 @@ def _read_next_timestep(self, ts=None):
 
             # Read in this order for now, then later reorder in place
             ts._pos[i] = map(float, self._file.readline().split())
-            if self._levcfg > 0:
+            if self._has_vels:
                 ts._velocities[i] = map(float, self._file.readline().split())
-            if self._levcfg == 2:
+            if self._has_forces:
                 ts._forces[i] = map(float, self._file.readline().split())
 
         if ids:
@@ -204,9 +204,9 @@ def _read_next_timestep(self, ts=None):
             if not all(ids == (np.arange(self.numatoms) + 1)):
                 order = np.argsort(ids)
                 ts._pos[:] = ts._pos[order]
-                if self._levcfg > 0:
+                if self._has_vels:
                     ts._velocities[:] = ts._velocities[order]
-                if self._levcfg == 2:
+                if self._has_forces:
                     ts._forces[:] = ts._forces[order]
 
         ts.frame += 1
@@ -251,9 +251,9 @@ def _read_numframes(self):
                 for _ in range(self.numatoms):
                     f.readline()
                     f.readline()
-                    if self._levcfg > 0:  # vels
+                    if self._has_vels:
                         f.readline()
-                    if self._levcfg == 2:  # forces
+                    if self._has_forces:
                         f.readline()
                 position = f.tell()
                 line = f.readline()

package/MDAnalysis/coordinates/DMS.py

@@ -26,18 +26,18 @@
 .. _DMS: http://www.deshawresearch.com/Desmond_Users_Guide-0.7.pdf
 """
 
-import numpy
+import numpy as np
 import sqlite3
 
 from . import base
-from MDAnalysis.coordinates.core import triclinic_box, triclinic_vectors
+from .core import triclinic_box, triclinic_vectors
 
 
 class Timestep(base.Timestep):
     def _init_unitcell(self):
-        return {'x': numpy.zeros(3),
-                'y': numpy.zeros(3),
-                'z': numpy.zeros(3)}
+        return {'x': np.zeros(3),
+                'y': np.zeros(3),
+                'z': np.zeros(3)}
 
     @property
     def dimensions(self):
@@ -84,34 +84,37 @@ def get_global_cell(self, cur):
     def _read_first_frame(self):
         coords_list = None
         velocities_list = None
-        con = sqlite3.connect(self.filename)
 
         def dict_factory(cursor, row):
             d = {}
             for idx, col in enumerate(cursor.description):
                 d[col[0]] = row[idx]
             return d
 
-        with con:
+        with sqlite3.connect(self.filename) as con:
             # This will return dictionaries instead of tuples, when calling cur.fetch() or fetchall()
             con.row_factory = dict_factory
             cur = con.cursor()
             coords_list = self.get_coordinates(cur)
             velocities_list = self.get_particle_by_columns(cur, columns=['vx', 'vy', 'vz'])
             unitcell = self.get_global_cell(cur)
-        assert coords_list
+
+        if not coords_list:
+            raise IOError("Found no coordinates")
         self.numatoms = len(coords_list)
-        coords_list = numpy.array(coords_list)
-        self.ts = self._Timestep(coords_list)
+
+        velocities = np.array(velocities_list, dtype=np.float32)
+        if not velocities.any():
+            velocities = None
+
+        self.ts = self._Timestep.from_coordinates(
+            np.array(coords_list, dtype=np.float32), velocities=velocities)
         self.ts.frame = 1  # 1-based frame number
-        if velocities_list:
-            # TODO: use a Timestep that knows about velocities such as TRR.Timestep or better, TRJ.Timestep
-            velocities_arr = numpy.array(velocities_list, dtype=numpy.float32)
-            if numpy.any(velocities_arr):
-                self.ts._velocities = velocities_arr
-                self.convert_velocities_from_native(self.ts._velocities)  # converts nm/ps to A/ps units
-        # ts._unitcell layout is format dependent; Timestep.dimensions does the conversion
+
         self.ts._unitcell = unitcell
         if self.convert_units:
             self.convert_pos_from_native(self.ts._pos)  # in-place !
             self.convert_pos_from_native(self.ts._unitcell)  # in-place ! (all are lengths)
+            if self.ts.has_velocities:
+                # converts nm/ps to A/ps units
+                self.convert_velocities_from_native(self.ts._velocities)

package/MDAnalysis/coordinates/GMS.py

@@ -30,7 +30,6 @@
 import os
 import errno
 import re
-import numpy
 
 from . import base
 import MDAnalysis.core.util as util
@@ -54,7 +53,6 @@ class GMSReader(base.Reader):
           on the fly
     """
 
-    # this will be overidden when an instance is created and the file extension checked
     format = "GMS"
 
     # these are assumed!
@@ -228,7 +226,6 @@ def _read_next_timestep(self, ts=None):
 
             # stop when the cursor has reached the end of that block
             if counter == self._numatoms:
-                ts._unitcell = numpy.zeros((6), numpy.float)
                 ts._x[:] = x # more efficient to do it this way to avoid re-creating the numpy arrays
                 ts._y[:] = y
                 ts._z[:] = z