diff --git a/package/CHANGELOG b/package/CHANGELOG index c7b5c1e3a2f..7efe1751e54 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,6 +19,8 @@ The rules for this file: * 2.4.0 Fixes + * XTC and TRR readers now fail with IOError when a status except EOK (=0) is + detected on reading a frame instead of silently continuing * Consolidate license files across MDAnalysis library packages (PR #3939) * Kwargs passed through to `MemoryReader` when using `transfer_to_memory()` and `in_memory=True`. @@ -37,6 +39,10 @@ Fixes (e.g. bonds, angles) (PR #3779). Enhancements + * Improve C content of libxdr Cython, add `read_direct` methods to read + coordinates, velocities and forces directly into memoryviews of `Timestep` + attributes, make `TRR` timestep have positions, velocities and forces on + `__init__`. (Issue #3891 PR #3892). * Improve C content of libdcd Cython (Issue #3882, PR #3888) * The timeseries method for exporting coordinates from multiple frames to a NumPy array was added to ProtoReader (PR #3890) diff --git a/package/MDAnalysis/coordinates/TRR.py b/package/MDAnalysis/coordinates/TRR.py index fb05af54039..1eca1bdd40c 100644 --- a/package/MDAnalysis/coordinates/TRR.py +++ b/package/MDAnalysis/coordinates/TRR.py @@ -31,6 +31,7 @@ MDAnalysis.coordinates.XTC: Read and write GROMACS XTC trajectory files. MDAnalysis.coordinates.XDR: BaseReader/Writer for XDR based formats """ +import errno from . import base from .XDR import XDRBaseReader, XDRBaseWriter from ..lib.formats.libmdaxdr import TRRFile @@ -151,6 +152,28 @@ class TRRReader(XDRBaseReader): _writer = TRRWriter _file = TRRFile + def _read_next_timestep(self, ts=None): + """copy next frame into timestep + + versionadded:: 2.4.0 + TRRReader implements this method so that it can use + read_direct_xvf to read the data directly into the timestep + rather than copying it from a temporary array. + """ + if self._frame == self.n_frames - 1: + raise IOError(errno.EIO, 'trying to go over trajectory limit') + if ts is None: + ts = self.ts + # allocate arrays to read into, will set to proper values + # in _frame_to_ts + ts.has_positions = True + ts.has_velocities = True + ts.has_forces = True + frame = self._xdr.read_direct_xvf(ts.positions, ts.velocities, ts.forces) + self._frame += 1 + self._frame_to_ts(frame, ts) + return ts + def _frame_to_ts(self, frame, ts): """convert a trr-frame to a mda TimeStep""" ts.time = frame.time diff --git a/package/MDAnalysis/coordinates/XDR.py b/package/MDAnalysis/coordinates/XDR.py index c374c310c26..a3e15b2f35f 100644 --- a/package/MDAnalysis/coordinates/XDR.py +++ b/package/MDAnalysis/coordinates/XDR.py @@ -120,7 +120,8 @@ class XDRBaseReader(base.ReaderBase): XDR offsets read from trajectory if offsets file read-in fails .. versionchanged:: 2.0.0 Add a InterProcessLock when generating offsets - + .. versionchanged:: 2.4.0 + Use a direct read into ts attributes """ @store_init_arguments def __init__(self, filename, convert_units=True, sub=None, @@ -288,17 +289,6 @@ def _read_frame(self, i): timestep = self._read_next_timestep() return timestep - def _read_next_timestep(self, ts=None): - """copy next frame into timestep""" - if self._frame == self.n_frames - 1: - raise IOError(errno.EIO, 'trying to go over trajectory limit') - if ts is None: - ts = self.ts - frame = self._xdr.read() - self._frame += 1 - self._frame_to_ts(frame, ts) - return ts - def Writer(self, filename, n_atoms=None, **kwargs): """Return writer for trajectory format""" if n_atoms is None: diff --git a/package/MDAnalysis/coordinates/XTC.py b/package/MDAnalysis/coordinates/XTC.py index 3a1b6b439cb..be473669347 100644 --- a/package/MDAnalysis/coordinates/XTC.py +++ b/package/MDAnalysis/coordinates/XTC.py @@ -31,6 +31,8 @@ MDAnalysis.coordinates.TRR: Read and write GROMACS TRR trajectory files. MDAnalysis.coordinates.XDR: BaseReader/Writer for XDR based formats """ + +import errno from . import base from .XDR import XDRBaseReader, XDRBaseWriter from ..lib.formats.libmdaxdr import XTCFile @@ -130,12 +132,35 @@ class XTCReader(XDRBaseReader): See :ref:`Notes on offsets ` for more information about offsets. + + """ format = 'XTC' units = {'time': 'ps', 'length': 'nm'} _writer = XTCWriter _file = XTCFile + def _read_next_timestep(self, ts=None): + """ + copy next frame into timestep + + versionadded:: 2.4.0 + XTCReader implements this method so that it can use + read_direct_x method of XTCFile to read the data directly + into the timestep rather than copying it from a temporary array. + """ + if self._frame == self.n_frames - 1: + raise IOError(errno.EIO, 'trying to go over trajectory limit') + if ts is None: + ts = self.ts + if ts.has_positions: + frame = self._xdr.read_direct_x(ts.positions) + else: + frame = self._xdr.read() + self._frame += 1 + self._frame_to_ts(frame, ts) + return ts + def _frame_to_ts(self, frame, ts): """convert a xtc-frame to a mda TimeStep""" ts.frame = self._frame diff --git a/package/MDAnalysis/lib/formats/libmdaxdr.pxd b/package/MDAnalysis/lib/formats/libmdaxdr.pxd new file mode 100644 index 00000000000..999ddd69334 --- /dev/null +++ b/package/MDAnalysis/lib/formats/libmdaxdr.pxd @@ -0,0 +1,110 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +import numpy as np +cimport numpy as np + +np.import_array() +from libc.stdint cimport int64_t + +from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END + +cdef extern from 'include/xdrfile.h': + ctypedef struct XDRFILE: + pass + + XDRFILE* xdrfile_open (char * path, char * mode) + int xdrfile_close (XDRFILE * xfp) + int xdr_seek(XDRFILE *xfp, int64_t pos, int whence) + int64_t xdr_tell(XDRFILE *xfp) + ctypedef float matrix[3][3] + ctypedef float rvec[3] + + +cdef extern from 'include/xdrfile_xtc.h': + int read_xtc_natoms(char * fname, int * natoms) + int read_xtc(XDRFILE * xfp, int natoms, int * step, float * time, matrix box, + rvec * x, float * prec) + int write_xtc(XDRFILE * xfp, int natoms, int step, float time, matrix box, + rvec * x, float prec) + + + +cdef extern from 'include/xdrfile_trr.h': + int read_trr_natoms(char *fname, int *natoms) + int read_trr(XDRFILE *xfp, int natoms, int *step, float *time, float *_lambda, + matrix box, rvec *x, rvec *v, rvec *f, int *has_prop) + int write_trr(XDRFILE *xfp, int natoms, int step, float time, float _lambda, + matrix box, rvec *x, rvec *v, rvec *f) + + +cdef extern from 'include/xtc_seek.h': + int read_xtc_n_frames(char *fn, int *n_frames, int *est_nframes, int64_t **offsets) + + +cdef extern from 'include/trr_seek.h': + int read_trr_n_frames(char *fn, int *n_frames, int *est_nframes, int64_t **offsets) + + +cdef class _XDRFile: + """Base python wrapper for gromacs-xdr formats + + This class can be similar to the normal file objects in python. The read() + function will return a frame and all information in it instead of a single + line. Additionally the context-manager protocoll is supported as well. + + Note + ---- + This class can't be initialized; use one of the subclasses :class:`XTCFile` or :class:`TRRFile`. + """ + # number of atoms + cdef readonly int n_atoms + # whether the file is open + cdef int is_open + # whether we have reached the end of the file + cdef int reached_eof + # the XDR file pointer + cdef XDRFILE *xfp + # the name of the xdr file + cdef readonly fname + # the current frame in the file + cdef int current_frame + # the file mode + cdef str mode + # the simulation box + cdef np.ndarray box + # numpy array of offsets into the fle + cdef np.ndarray _offsets + # whether we have the offsets + cdef readonly int _has_offsets + + +cdef class XTCFile(_XDRFile): + """File-like wrapper for gromacs XTC files.""" + # precision of the XTC file + cdef float precision + + + +cdef class TRRFile(_XDRFile): + """File-like wrapper for gromacs TRR files""" diff --git a/package/MDAnalysis/lib/formats/libmdaxdr.pyx b/package/MDAnalysis/lib/formats/libmdaxdr.pyx index 81aaddebcb6..147b7be630c 100644 --- a/package/MDAnalysis/lib/formats/libmdaxdr.pyx +++ b/package/MDAnalysis/lib/formats/libmdaxdr.pyx @@ -20,6 +20,9 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # + +#cython: boundscheck=False, wraparound=False + """\ Low-level Gromacs XDR trajectory reading — :mod:`MDAnalysis.lib.formats.libmdaxdr` ---------------------------------------------------------------------------------- @@ -68,43 +71,6 @@ from libc.stdint cimport int64_t from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END _whence_vals = {"SEEK_SET": SEEK_SET, "SEEK_CUR": SEEK_CUR, "SEEK_END": SEEK_END} -cdef extern from 'include/xdrfile.h': - ctypedef struct XDRFILE: - pass - - XDRFILE* xdrfile_open (char * path, char * mode) - int xdrfile_close (XDRFILE * xfp) - int xdr_seek(XDRFILE *xfp, int64_t pos, int whence) - int64_t xdr_tell(XDRFILE *xfp) - ctypedef float matrix[3][3] - ctypedef float rvec[3] - - -cdef extern from 'include/xdrfile_xtc.h': - int read_xtc_natoms(char * fname, int * natoms) - int read_xtc(XDRFILE * xfp, int natoms, int * step, float * time, matrix box, - rvec * x, float * prec) - int write_xtc(XDRFILE * xfp, int natoms, int step, float time, matrix box, - rvec * x, float prec) - - - -cdef extern from 'include/xdrfile_trr.h': - int read_trr_natoms(char *fname, int *natoms) - int read_trr(XDRFILE *xfp, int natoms, int *step, float *time, float *_lambda, - matrix box, rvec *x, rvec *v, rvec *f, int *has_prop) - int write_trr(XDRFILE *xfp, int natoms, int step, float time, float _lambda, - matrix box, rvec *x, rvec *v, rvec *f) - - -cdef extern from 'include/xtc_seek.h': - int read_xtc_n_frames(char *fn, int *n_frames, int *est_nframes, int64_t **offsets) - - -cdef extern from 'include/trr_seek.h': - int read_trr_n_frames(char *fn, int *n_frames, int *est_nframes, int64_t **offsets) - - cdef enum: EOK = 0 EHEADER = 1 @@ -132,7 +98,7 @@ from collections import namedtuple np.import_array() -ctypedef np.float32_t DTYPE_T +ctypedef float DTYPE_T DTYPE = np.float32 cdef int DIMS = 3 cdef int HASX = 1 @@ -162,16 +128,6 @@ cdef class _XDRFile: ---- This class can't be initialized use one of the subclasses XTCFile, TRRFile """ - cdef readonly int n_atoms - cdef int is_open - cdef int reached_eof - cdef XDRFILE *xfp - cdef readonly fname - cdef int current_frame - cdef str mode - cdef np.ndarray box - cdef np.ndarray _offsets - cdef readonly int _has_offsets def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -222,6 +178,8 @@ cdef class _XDRFile: self.reached_eof = False self.current_frame = 0 + cdef int return_code + if mode == 'r': opening_mode = b'r' elif mode == 'w': @@ -318,10 +276,9 @@ cdef class _XDRFile: else: # pragma: no cover raise RuntimeError("Invalid frame number {} > {} -- this should" "not happen.".format(current_frame, - self.offsets.size) - ) + self.offsets.size)) - def seek(self, frame): + def seek(self, int frame): """Seek to Frame. Please note that this function will generate internal file offsets if @@ -357,7 +314,7 @@ cdef class _XDRFile: raise IOError("XDR seek failed with system errno={}".format(ok)) self.current_frame = frame - def _bytes_seek(self, offset, whence="SEEK_SET"): + def _bytes_seek(self, int64_t offset, whence="SEEK_SET"): """Low-level access to the stream repositioning xdr_seek call. Beware that this function will not update :attr:`current_frame`, @@ -412,7 +369,7 @@ cdef class _XDRFile: self._has_offsets = True return self._offsets - def set_offsets(self, offsets): + def set_offsets(self, np.ndarray offsets): """set frame offsets""" self._offsets = offsets self._has_offsets = True @@ -425,7 +382,6 @@ cdef class _XDRFile: """Low-level call to xdr_tell to get current byte offset.""" return xdr_tell(self.xfp) - TRRFrame = namedtuple('TRRFrame', 'x v f box step time lmbda hasx hasv hasf') @@ -454,14 +410,17 @@ cdef class TRRFile(_XDRFile): ----- This class can be pickled. The pickle will store filename, mode, current frame and offsets + + + .. versionchanged:: 2.4.0 + Added read_direct_xvf method to read into an existing positions array """ def _calc_natoms(self, fname): cdef int n_atoms - return_code = read_trr_natoms(fname, &n_atoms) + cdef int return_code = read_trr_natoms(fname, &n_atoms) return return_code, n_atoms - def calc_offsets(self): """read byte offsets from TRR file directly""" if not self.is_open: @@ -469,7 +428,7 @@ cdef class TRRFile(_XDRFile): cdef int n_frames = 0 cdef int est_nframes = 0 cdef int64_t* offsets = NULL - ok = read_trr_n_frames(self.fname, &n_frames, &est_nframes, &offsets) + cdef int ok = read_trr_n_frames(self.fname, &n_frames, &est_nframes, &offsets) if ok != EOK: raise IOError("TRR couldn't calculate offsets. " "XDR error = {}".format(error_message[ok])) @@ -477,7 +436,10 @@ cdef class TRRFile(_XDRFile): # overestimation. This number is saved in est_nframes and we need to # tell the new numpy array about the whole allocated memory to avoid # memory leaks. - cdef np.ndarray dims = np.array([est_nframes], dtype=np.int64) + cdef np.npy_intp[1] dim + dim[0] = 1 + cdef np.ndarray[np.int64_t, ndim=1] dims = np.PyArray_EMPTY(1, dim, np.NPY_INT64, 0) + dims[0] = est_nframes # this handles freeing the allocated memory correctly. cdef np.ndarray nd_offsets = ptr_to_ndarray( offsets, dims, np.NPY_INT64) return nd_offsets[:n_frames] @@ -513,12 +475,18 @@ cdef class TRRFile(_XDRFile): cdef float time = 0 cdef float lmbda = 0 - # Use this instead of memviews here to make sure that references are - # counted correctly - cdef np.ndarray xyz = np.empty((self.n_atoms, DIMS), dtype=DTYPE) - cdef np.ndarray velocity = np.empty((self.n_atoms, DIMS), dtype=DTYPE) - cdef np.ndarray forces = np.empty((self.n_atoms, DIMS), dtype=DTYPE) - cdef np.ndarray box = np.empty((DIMS, DIMS), dtype=DTYPE) + cdef np.npy_intp[2] dim + dim[0] = self.n_atoms + dim[1] = DIMS + + cdef np.npy_intp[2] unitcell_dim + unitcell_dim[0] = DIMS + unitcell_dim[1] = DIMS + + cdef np.ndarray[np.float32_t, ndim=2] xyz = np.PyArray_EMPTY(2, dim, np.NPY_FLOAT32, 0) + cdef np.ndarray[np.float32_t, ndim=2] velocity = np.PyArray_EMPTY(2, dim, np.NPY_FLOAT32, 0) + cdef np.ndarray[np.float32_t, ndim=2] forces = np.PyArray_EMPTY(2, dim, np.NPY_FLOAT32, 0) + cdef np.ndarray[np.float32_t, ndim=2] box = np.PyArray_EMPTY(2, unitcell_dim, np.NPY_FLOAT32, 0) return_code = read_trr(self.xfp, self.n_atoms, &step, &time, &lmbda, box.data, @@ -529,11 +497,84 @@ cdef class TRRFile(_XDRFile): # trr are a bit weird. Reading after the last frame always always # results in an integer error while reading. I tried it also with trr # produced by different codes (Gromacs, ...). - if return_code != EOK and return_code != EENDOFFILE \ - and return_code != EINTEGER: + # In a trr the integer error seems to indicate that the file is ending. + # There might be corrupted files where this is a legitimate error. But + # then we just can't read it and stop there which is not too bad. + if return_code == EENDOFFILE or return_code == EINTEGER: + self.reached_eof = True + raise StopIteration + + if return_code != EOK: raise IOError('TRR read error = {}'.format( error_message[return_code])) + self.current_frame += 1 + + has_x = bool(has_prop & HASX) + has_v = bool(has_prop & HASV) + has_f = bool(has_prop & HASF) + return TRRFrame(xyz, velocity, forces, box, step, time, lmbda, + has_x, has_v, has_f) + + def read_direct_xvf(self, np.float32_t[:, ::1] positions, + np.float32_t[:, ::1] velocities, + np.float32_t[:, ::1] forces,): + """ + Read next frame in the TRR file with positions read directly into + a pre-existing array. + + Parameters + ---------- + positions : np.ndarray + positions array to read positions into + + Returns + ------- + frame : libmdaxdr.TRRFrame + namedtuple with frame information + + See Also + -------- + TRRFrame + XTCFile + + Raises + ------ + IOError + + + .. versionadded:: 2.4.0 + """ + if self.reached_eof: + raise EOFError('Reached last frame in TRR, seek to 0') + if not self.is_open: + raise IOError('No file opened') + if self.mode != 'r': + raise IOError('File opened in mode: {}. Reading only allow ' + 'in mode "r"'.format('self.mode')) + + cdef int return_code = 1 + cdef int step = 0 + cdef int has_prop = 0 + cdef float time = 0 + cdef float lmbda = 0 + + cdef np.npy_intp[2] unitcell_dim + unitcell_dim[0] = DIMS + unitcell_dim[1] = DIMS + + + cdef np.ndarray[np.float32_t, ndim=2] box = np.PyArray_EMPTY(2, unitcell_dim, np.NPY_FLOAT32, 0) + + return_code = read_trr(self.xfp, self.n_atoms, &step, + &time, &lmbda, box.data, + &positions[0,0], + &velocities[0,0], + &forces[0,0], + &has_prop) + # trr are a bit weird. Reading after the last frame always always + # results in an integer error while reading. I tried it also with trr + # produced by different codes (Gromacs, ...). # In a trr the integer error seems to indicate that the file is ending. # There might be corrupted files where this is a legitimate error. But # then we just can't read it and stop there which is not too bad. @@ -541,13 +582,15 @@ cdef class TRRFile(_XDRFile): self.reached_eof = True raise StopIteration - if return_code == EOK: - self.current_frame += 1 + if return_code != EOK: + raise IOError('TRR read error = {}'.format(error_message[return_code])) + + self.current_frame += 1 has_x = bool(has_prop & HASX) has_v = bool(has_prop & HASV) has_f = bool(has_prop & HASF) - return TRRFrame(xyz, velocity, forces, box, step, time, lmbda, + return TRRFrame(positions, velocities, forces, box, step, time, lmbda, has_x, has_v, has_f) def write(self, xyz, velocity, forces, box, int step, float time, @@ -671,15 +714,16 @@ cdef class XTCFile(_XDRFile): ----- This class can be pickled. The pickle will store filename, mode, current frame and offsets + + .. versionchanged:: 2.4.0 + Added read_direct_x method to read into an existing positions array """ - cdef float precision def _calc_natoms(self, fname): cdef int n_atoms - return_code = read_xtc_natoms(fname, &n_atoms) + cdef int return_code = read_xtc_natoms(fname, &n_atoms) return return_code, n_atoms - def calc_offsets(self): """Calculate offsets from XTC file directly""" if not self.is_open: @@ -687,7 +731,7 @@ cdef class XTCFile(_XDRFile): cdef int n_frames = 0 cdef int est_nframes = 0 cdef int64_t* offsets = NULL - ok = read_xtc_n_frames(self.fname, &n_frames, &est_nframes, &offsets) + cdef int ok = read_xtc_n_frames(self.fname, &n_frames, &est_nframes, &offsets) if ok != EOK: raise IOError("XTC couldn't calculate offsets. " "XDR error = {}".format(error_message[ok])) @@ -695,7 +739,10 @@ cdef class XTCFile(_XDRFile): # overestimation. This number is saved in est_nframes and we need to # tell the new numpy array about the whole allocated memory to avoid # memory leaks. - cdef np.ndarray dims = np.array([est_nframes], dtype=np.int64) + cdef np.npy_intp[1] dim + dim[0] = 1 + cdef np.ndarray[np.int64_t, ndim=1] dims = np.PyArray_EMPTY(1, dim, np.NPY_INT64, 0) + dims[0] = est_nframes # this handles freeing the allocated memory correctly. cdef np.ndarray nd_offsets = ptr_to_ndarray( offsets, dims, np.NPY_INT64) return nd_offsets[:n_frames] @@ -729,24 +776,91 @@ cdef class XTCFile(_XDRFile): cdef int step cdef float time, prec - cdef np.ndarray xyz = np.empty((self.n_atoms, DIMS), dtype=DTYPE) - cdef np.ndarray box = np.empty((DIMS, DIMS), dtype=DTYPE) + cdef np.npy_intp[2] dim + dim[0] = self.n_atoms + dim[1] = DIMS + + cdef np.npy_intp[2] unitcell_dim + unitcell_dim[0] = DIMS + unitcell_dim[1] = DIMS + + cdef np.ndarray[np.float32_t, ndim=2] xyz = np.PyArray_EMPTY(2, dim, np.NPY_FLOAT32, 0) + cdef np.ndarray[np.float32_t, ndim=2] box = np.PyArray_EMPTY(2, unitcell_dim, np.NPY_FLOAT32, 0) return_code = read_xtc(self.xfp, self.n_atoms, &step, &time, box.data, xyz.data, &prec) - if return_code != EOK and return_code != EENDOFFILE: - raise IOError('XTC read error = {}'.format( - error_message[return_code])) if return_code == EENDOFFILE: self.reached_eof = True raise StopIteration - if return_code == EOK: - self.current_frame += 1 + if return_code != EOK: + raise IOError('XTC read error = {}'.format(error_message[return_code])) + self.current_frame += 1 + return XTCFrame(xyz, box, step, time, prec) + def read_direct_x(self, np.float32_t[:, ::1] positions): + """ + Read next frame in the XTC file with positions read directly into + a pre-existing array. + + Parameters + ---------- + positions : np.ndarray + positions array to read positions into + + Returns + ------- + frame : libmdaxdr.XTCFrame + namedtuple with frame information + + See Also + -------- + XTCFrame + TRRFile + + Raises + ------ + IOError + + + .. versionadded:: 2.4.0 + """ + if self.reached_eof: + raise EOFError('Reached last frame in XTC, seek to 0') + if not self.is_open: + raise IOError('No file opened') + if self.mode != 'r': + raise IOError('File opened in mode: {}. Reading only allow ' + 'in mode "r"'.format('self.mode')) + + return_code = 1 + cdef int step + cdef float time, prec + cdef np.npy_intp[2] unitcell_dim + unitcell_dim[0] = DIMS + unitcell_dim[1] = DIMS + + cdef np.ndarray[np.float32_t, ndim=2] box = np.PyArray_EMPTY(2, unitcell_dim, np.NPY_FLOAT32, 0) + + + return_code = read_xtc(self.xfp, self.n_atoms, &step, + &time, box.data, + &positions[0,0], &prec) + + if return_code == EENDOFFILE: + self.reached_eof = True + raise StopIteration + + if return_code != EOK: + raise IOError('XTC read error = {}'.format(error_message[return_code])) + self.current_frame += 1 + + return XTCFrame(positions, box, step, time, prec) + + def write(self, xyz, box, int step, float time, float precision=1000): """write one frame to the XTC file @@ -779,11 +893,11 @@ cdef class XTCFile(_XDRFile): raise IOError('File opened in mode: {}. Writing only allow ' 'in mode "w"'.format('self.mode')) - xyz = np.asarray(xyz) - box = np.asarray(box) + xyz = np.asarray(xyz, dtype=np.float32) + box = np.asarray(box, dtype=np.float32) - cdef DTYPE_T[:, ::1] xyz_view = np.ascontiguousarray(xyz, dtype=DTYPE) - cdef DTYPE_T[:, ::1] box_view = np.ascontiguousarray(box, dtype=DTYPE) + cdef DTYPE_T[:, ::1] xyz_view = np.PyArray_GETCONTIGUOUS(xyz) + cdef DTYPE_T[:, ::1] box_view = np.PyArray_GETCONTIGUOUS(box) if self.current_frame == 0: self.n_atoms = xyz.shape[0] diff --git a/package/MDAnalysis/lib/libmdanalysis/__init__.pxd b/package/MDAnalysis/lib/libmdanalysis/__init__.pxd index 318c411f55b..8f24ccac51e 100644 --- a/package/MDAnalysis/lib/libmdanalysis/__init__.pxd +++ b/package/MDAnalysis/lib/libmdanalysis/__init__.pxd @@ -2,4 +2,5 @@ # single place. from ..coordinates cimport timestep -from .formats cimport libdcd \ No newline at end of file +from .formats cimport libmdaxdr +from .formats cimport libdcd diff --git a/testsuite/MDAnalysisTests/coordinates/test_xdr.py b/testsuite/MDAnalysisTests/coordinates/test_xdr.py index 2c9cd22e8f0..1ca133acc6b 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_xdr.py +++ b/testsuite/MDAnalysisTests/coordinates/test_xdr.py @@ -213,7 +213,15 @@ def go_beyond_EOF(): with pytest.raises(StopIteration): go_beyond_EOF() - + + def test_read_next_timestep_ts_no_positions(self, universe): + # primarily tests branching on ts.has_positions in _read_next_timestep + ts = universe.trajectory[0] + ts.has_positions=False + ts_passed_in = universe.trajectory._read_next_timestep(ts=ts).copy() + universe.trajectory.rewind() + ts_returned = universe.trajectory._read_next_timestep(ts=None).copy() + assert(ts_passed_in == ts_returned) class TestXTCReader(_GromacsReader): filename = XTC