Merge pull request #441 from kain88-de/cython-xdrlib

Cython xdrlib

.gitignore

@@ -26,5 +26,6 @@ examples/output.txt
 .adk_oplsaa.trr_offsets.pkl
 .adk_oplsaa.xtc_offsets.pkl
 .gram_A_identical_frames.xtc_offsets.pkl
+*.npz
 # Ignore html build
 doc/html/

package/CHANGELOG

@@ -13,16 +13,22 @@ The rules for this file:
   * release numbers follow "Semantic Versioning" http://semver.org
 
 ------------------------------------------------------------------------------
-??/??/16
+??/??/16 kain88-de
 
   * 0.13.1
 
 API Changes
+  * Offsets files for Gromcas trajectory formats have been changed to a numpy
+    style format '.npz'. Offsets files will be regenerated when you load a
+    xtc/trr trajectory again.
 
 Enhancement
+  * Offsets reading for xtc/trr files has been sped up.
 
 Changes
 
+  * xdrlib has been ported to cython.
+
 Fixes
 
 01/16/16 tyler.je.reddy, kain88-de, richardjgowers, manuel.nuno.melo,

package/MDAnalysis/coordinates/TRR.py

@@ -1,9 +1,9 @@
 # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
-# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
 #
 # MDAnalysis --- http://www.MDAnalysis.org
-# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein
-# and contributors (see AUTHORS for the full list)
+# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver
+# Beckstein and contributors (see AUTHORS for the full list)
 #
 # Released under the GNU Public Licence, v2 or any higher version
 #
@@ -13,111 +13,131 @@
 # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 #
-
-"""
-Gromacs TRR file IO --- :mod:`MDAnalysis.coordinates.TRR`
-=========================================================
-
-The Gromacs `TRR trajectory format`_ is a lossless format like
-e.g. the DCD format (see :mod:`~MDAnalysis.coordinates.DCD`) and
-unlike the :mod:`~MDAnalysis.coordinates.XTC` format, which stores
-reduced precision coordinates. Therefore, if one wants to convert
-*to* Gromacs trajectories without loss of precision then one should
-use the TRR format.
-
-The TRR format can store *velocities* and *forces* in addition to
-coordinates. It is also used by other Gromacs tools to store and
-process other data such as modes from a principal component analysis.
-
-The TRR I/O interface uses :mod:`~MDAnalysis.coordinates.xdrfile.libxdrfile2`
-to implement random access to frames. This works by initially building an
-internal index of all frames and then using this index for direct
-seeks. Building the index is triggered by
-:func:`~MDAnalysis.coordinates.xdrfile.libxdrfile2.read_trr_n_frames`, which
-typically happens when one accesses the :attr:`TRRReader.n_frames` attribute
-for the first time. Building the index may take many minutes for large
-trajectories but afterwards access is faster than with native Gromacs tools.
-
-
-.. _TRR trajectory format:
-   http://www.gromacs.org/Documentation/File_Formats/.trr_File
-
-.. versionchanged:: 0.8.0
-   The TRR I/O interface now uses
-   :mod:`~MDAnalysis.coordinates.xdrfile.libxdrfile2`, which has seeking and
-   indexing capabilities. Note that unlike
-   :mod:`~MDAnalysis.coordinates.xdrfile.libxdrfile` before it,
-   :mod:`~MDAnalysis.coordinates.xdrfile.libxdrfile2` is distributed under the
-   GNU GENERAL PUBLIC LICENSE, version 2 (or higher).
-   :class:`~MDAnalysis.coordinates.TRR.Timestep` now correctly deals
-   with presence/absence of coordinate/velocity/force information on a
-   per-frame basis.
-
-
-Tips and Tricks
----------------
-
-Filling a TRR with PCA modes
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The following `recipe by Ramon Crehuet`_ shows how to convert modes
-stored in a NumPy-like array (e.g. from a PCA analysis with MMTK_) to
-a TRR usable by Gromacs. The idea is to manually fill a
-:class:`~MDAnalysis.coordinates.xdrfile.TRR.Timestep` with the desired
-values and then write it to a file with the appropriate
-:class:`~MDAnalysis.coordinates.xdrfile.TRR.TRRWriter`. In order to
-respect the Gromacs format for modes in a TRR file, one must write the
-average coordinates in the first frame of the TRR and the modes into
-subsequent ones. The mode number is stored in the
-:attr:`~MDAnalysis.coordinates.xdrfile.TRR.Timestep.step` attribute
-and the mode coordinates are filling the
-:attr:`~MDAnalysis.coordinates.xdrfile.TRR.Timestep.positions` attribute of
-:class:`~MDAnalysis.coordinates.xdrfile.TRR.Timestep`::
-
-   # 'modes' is a mode object with M PCs, similar to a MxNx3 array
-   # 'xav' the average coordinates, a Nx3 array for N atoms
-
-   N = len(xav)   # number of atoms, i.e. number of coordinates
-
-   W = Writer('pca.trr', n_atoms=N)            # TRR writer
-   ts = MDAnalysis.coordinates.TRR.Timestep(N)  # TRR time step
-                                                # N of atoms is passed.
-   for frame,mode in enumerate(modes[4:16]):
-       ts.lmbda = -1
-
-       ts.frame = frame         # manually change the frame number
-       ts._frame = frame - 1
-
-       if frame<=1:
-          ts.positions = xav
-       else:
-          ts.positions = mode.scaledToNorm(1.).array*10   # nm to angstroms
-       if frame <= 1:
-          ts.time = frame-1
-       else:
-          ts.time = mode.frequency
-       W.write(ts)             # converts angstrom to nm for gmx
-
-    W.close()
-
-.. _MMTK: http://dirac.cnrs-orleans.fr/Manuals/MMTK/index.html
-
-.. _`recipe by Ramon Crehuet`: https://github.com/MDAnalysis/mdanalysis/issues/79
-
-Module reference
-----------------
-
-.. autoclass:: Timestep
-   :members:
-   :inherited-members:
-
-.. autoclass:: TRRReader
-   :members:
-   :inherited-members:
-
-.. autoclass:: TRRWriter
-   :members:
-   :inherited-members:
-"""
-
-from .xdrfile.TRR import TRRReader, TRRWriter, Timestep
+from .XDR import XDRBaseReader, XDRBaseWriter
+from ..lib.formats.xdrlib import TRRFile
+from ..lib.mdamath import triclinic_vectors
+
+
+class TRRWriter(XDRBaseWriter):
+    """The Gromacs TRR trajectory format is a lossless format. The TRR format can
+    store *velocoties* and *forces* in addition to the coordinates. It is also
+    used by other Gromacs tools to store and process other data such as modes
+    from a principal component analysis.
+
+    Parameter
+    ---------
+    filename : str
+        filename of the trajectory
+    n_atoms : int
+        number of atoms to write
+    convert_units : bool (optional)
+        convert into MDAnalysis units
+    precision : float (optional)
+        set precision of saved trjactory to this number of decimal places.
+    """
+
+    format = 'TRR'
+    units = {'time': 'ps', 'length': 'nm', 'velocity': 'nm/ps',
+             'force': 'kJ/(mol*nm)'}
+    _file = TRRFile
+
+    def write_next_timestep(self, ts):
+        """Write timestep object into trajectory.
+
+        Parameters
+        ----------
+        ts : TimeStep
+
+        See Also
+        --------
+        <FormatWriter>.write(AtomGroup/Universe/TimeStep)
+        The normal write() method takes a more general input
+        """
+        xyz = None
+        if ts.has_positions:
+            xyz = ts.positions.copy()
+            if self._convert_units:
+                self.convert_pos_to_native(xyz)
+
+        velo = None
+        if ts.has_velocities:
+            velo = ts.velocities.copy()
+            if self._convert_units:
+                self.convert_velocities_to_native(velo)
+
+        forces = None
+        if ts.has_forces:
+            forces = ts.forces.copy()
+            if self._convert_units:
+                self.convert_forces_to_native(forces)
+
+        time = ts.time
+        step = ts.frame
+
+        if self._convert_units:
+            dimensions = self.convert_dimensions_to_unitcell(ts, inplace=False)
+
+        box = triclinic_vectors(dimensions)
+
+        self._xdr.write(xyz, velo, forces, box, step, time, 1, self.n_atoms)
+
+
+class TRRReader(XDRBaseReader):
+    """The Gromacs TRR trajectory format is a lossless format. The TRR format can
+    store *velocoties* and *forces* in addition to the coordinates. It is also
+    used by other Gromacs tools to store and process other data such as modes
+    from a principal component analysis.
+
+    Parameter
+    ---------
+    filename : str
+        filename of the trajectory
+    convert_units : bool (optional)
+        convert into MDAnalysis units
+    sub : atomgroup (optional)
+        Yeah what is that exactly
+    refresh_offsets : bool (optional)
+        Recalculate offsets for random access from file. If ``False`` try to
+        retrieve offsets from hidden offsets file.
+    """
+    format = 'TRR'
+    units = {'time': 'ps', 'length': 'nm', 'velocity': 'nm/ps',
+             'force': 'kJ/(mol*nm)'}
+    _writer = TRRWriter
+    _file = TRRFile
+
+    def _frame_to_ts(self, frame, ts):
+        """convert a trr-frame to a mda TimeStep"""
+        ts.time = frame.time
+        ts.frame = self._frame
+        ts.data['step'] = frame.step
+
+        ts.has_positions = frame.hasx
+        ts.has_velocities = frame.hasv
+        ts.has_forces = frame.hasf
+
+        if ts.has_positions:
+            if self._sub is not None:
+                ts.positions = frame.x[self._sub]
+            else:
+                ts.positions = frame.x
+            if self.convert_units:
+                self.convert_pos_from_native(ts.positions)
+
+        if ts.has_velocities:
+            if self._sub is not None:
+                ts.velocities = frame.v[self._sub]
+            else:
+                ts.velocities = frame.v
+            if self.convert_units:
+                self.convert_velocities_from_native(ts.velocities)
+
+        if ts.has_forces:
+            if self._sub is not None:
+                ts.forces = frame.f[self._sub]
+            else:
+                ts.forces = frame.f
+            if self.convert_units:
+                self.convert_forces_from_native(ts.forces)
+
+        return ts