diff --git a/.travis.yml b/.travis.yml index 1fdbab25f06..2bac9adee56 100644 --- a/.travis.yml +++ b/.travis.yml @@ -31,8 +31,8 @@ env: - MAIN_CMD="pytest ${PYTEST_LIST} ${PYTEST_FLAGS}; python ./testsuite/MDAnalysisTests/mda_nosetests ${NOSE_TEST_LIST} ${NOSE_FLAGS}" - SETUP_CMD="" - BUILD_CMD="pip install -v package/ && pip install testsuite/" - - CONDA_DEPENDENCIES="mmtf-python nose=1.3.7 mock six biopython networkx cython joblib nose-timer matplotlib scipy griddataformats pytest-cov pytest-xdist" - - CONDA_ALL_DEPENDENCIES="mmtf-python nose=1.3.7 mock six biopython networkx cython joblib nose-timer matplotlib netcdf4 scikit-learn scipy griddataformats seaborn coveralls clustalw=2.1 pytest-cov pytest-xdist" + - CONDA_DEPENDENCIES="mmtf-python nose=1.3.7 mock six biopython networkx cython joblib nose-timer matplotlib scipy griddataformats pytest-cov pytest-xdist hypothesis" + - CONDA_ALL_DEPENDENCIES="mmtf-python nose=1.3.7 mock six biopython networkx cython joblib nose-timer matplotlib netcdf4 scikit-learn scipy griddataformats seaborn coveralls clustalw=2.1 pytest-cov pytest-xdist hypothesis" # Install griddataformats from PIP so that scipy is only installed in the full build (#1147) - PIP_DEPENDENCIES='pytest-raises' - CONDA_CHANNELS='biobuilds conda-forge' diff --git a/package/CHANGELOG b/package/CHANGELOG index 170dbe84cbb..03289d0dc30 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -13,11 +13,13 @@ The rules for this file: * release numbers follow "Semantic Versioning" http://semver.org ------------------------------------------------------------------------------ -mm/dd/17 richardjgowers, rathann, orbeckst, tylerjereddy, mtiberti +mm/dd/17 richardjgowers, rathann, orbeckst, tylerjereddy, mtiberti, kain88-de * 0.17.0 Enhancements + * add low level lib.formats.libdcd module for reading/writing DCD (PR #1372) + * replace old DCD reader with a Python 3 ready DCD reader (Issue #659) Deprecations diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 3c769b2a762..3175c756c6a 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -19,8 +19,6 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # - - """DCD trajectory I/O --- :mod:`MDAnalysis.coordinates.DCD` ============================================================ @@ -29,27 +27,14 @@ system-endianness as this is auto-detected. Generally, DCD trajectories produced by any code can be read (with the -:class:`DCDReader`) although there can be issues with the unitcell -(simulation box) representation (see -:attr:`Timestep.dimensions`). DCDs can also be written but the -:class:`DCDWriter` follows recent NAMD/VMD convention for the unitcell -but still writes AKMA time. Reading and writing these trajectories -within MDAnalysis will work seamlessly but if you process those -trajectories with other tools you might need to watch out that time -and unitcell dimensions are correctly interpreted. - -Note ----- -The DCD file format is not well defined. In particular, NAMD and -CHARMM use it differently. Currently, MDAnalysis tries to guess the -correct **format for the unitcell representation** but it can be -wrong. **Check the unitcell dimensions**, especially for triclinic -unitcells (see `Issue 187`_ and :attr:`Timestep.dimensions`). A -second potential issue are the units of time which are AKMA for the -:class:`DCDReader` (following CHARMM) but ps for NAMD. As a -workaround one can employ the configurable -:class:`MDAnalysis.coordinates.LAMMPS.DCDReader` for NAMD -trajectories. +:class:`DCDReader`) although there can be issues with the unitcell (simulation +box) representation (see :attr:`DCDReader.dimensions`). DCDs can also be +written but the :class:`DCDWriter` follows recent NAMD/VMD convention for the +unitcell but still writes AKMA time. Reading and writing these trajectories +within MDAnalysis will work seamlessly but if you process those trajectories +with other tools you might need to watch out that time and unitcell dimensions +are correctly interpreted. + See Also -------- @@ -57,9 +42,6 @@ module provides a more flexible DCD reader/writer. -The classes in this module are the reference implementations for the -Trajectory API. - .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187 @@ -67,8 +49,6 @@ Classes ------- -.. autoclass:: Timestep - :inherited-members: .. autoclass:: DCDReader :inherited-members: .. autoclass:: DCDWriter @@ -89,438 +69,203 @@ from ..core import flags from .. import units as mdaunits # use mdaunits instead of units to avoid a clash from ..exceptions import NoDataError -from . import base -from . import core -# Add the c functions to their respective classes so they act as class methods -from . import _dcdmodule - - -class Timestep(base.Timestep): - #: Indices into :attr:`Timestep._unitcell` (``[A, gamma, B, beta, alpha, - #: C]``, provided by the :class:`DCDReader` C code) to pull out - #: ``[A, B, C, alpha, beta, gamma]``. - _ts_order = [0, 2, 5, 4, 3, 1] - - @property - def dimensions(self): - """unitcell dimensions (*A*, *B*, *C*, *alpha*, *beta*, *gamma*) - - lengths *A*, *B*, *C* are in the MDAnalysis length unit (Å), and - angles are in degrees. - - :attr:`dimensions` is read-only because it transforms the actual format - of the unitcell (which differs between different trajectory formats) to - the representation described here, which is used everywhere in - MDAnalysis. - - The ordering of the angles in the unitcell is the same as in recent - versions of VMD's DCDplugin_ (2013), namely the `X-PLOR DCD format`_: - The original unitcell is read as ``[A, gamma, B, beta, alpha, C]`` from - the DCD file (actually, the direction cosines are stored instead of the - angles but the underlying C code already does this conversion); if any - of these values are < 0 or if any of the angles are > 180 degrees then - it is assumed it is a new-style CHARMM unitcell (at least since c36b2) - in which box vectors were recorded. - - - .. warning:: The DCD format is not well defined. Check your unit cell - dimensions carefully, especially when using triclinic - boxes. Different software packages implement different conventions - and MDAnalysis is currently implementing the newer NAMD/VMD convention - and tries to guess the new CHARMM one. Old CHARMM trajectories might - give wrong unitcell values. For more details see `Issue 187`_. - - .. versionchanged:: 0.9.0 - Unitcell is now interpreted in the newer NAMD DCD format as ``[A, - gamma, B, beta, alpha, C]`` instead of the old MDAnalysis/CHARMM - ordering ``[A, alpha, B, beta, gamma, C]``. We attempt to detect the - new CHARMM DCD unitcell format (see `Issue 187`_ for a discussion). - - .. _`X-PLOR DCD format`: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html - .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187 - .. _DCDplugin: http://www.ks.uiuc.edu/Research/vmd/plugins/doxygen/dcdplugin_8c-source.html#l00947 - """ - - # Layout of unitcell is [A, alpha, B, beta, gamma, C] --- (originally CHARMM DCD) - # 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 np.take(self._unitcell, [0,2,5,1,3,4]) +from . import base, core +from ..lib.formats.libdcd import DCDFile +from ..lib.mdamath import triclinic_box - # 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 = 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 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]] - uc = core.triclinic_box(e1, e2, e3) - return uc - @dimensions.setter - def dimensions(self, box): - """Set unitcell with (*A*, *B*, *C*, *alpha*, *beta*, *gamma*) - - .. versionadded:: 0.9.0 - """ - # note that we can re-use self._ts_order with put! - np.put(self._unitcell, self._ts_order, box) - - -class DCDWriter(base.WriterBase): - """Write to a CHARMM/NAMD DCD trajectory file. - - Parameters - ---------- - filename : str - name of output file - n_atoms : int (optional) - number of atoms in dcd file - start : int (optional) - starting frame number - step : int (optional) - skip between subsequent timesteps (indicate that `step` **MD - integrator steps (!)** make up one trajectory frame); default is 1. - delta : float (optional) - timestep (**MD integrator time step (!)**, in AKMA units); default is - 20.45482949774598 (corresponding to 1 ps). - remarks : str (optional) - comments to annotate dcd file - dt : float (optional) - **Override** `step` and `delta` so that the DCD records that - `dt` ps lie between two frames. (It sets `step` = 1 and - `delta` ``= AKMA(dt)``.) The default is ``None``, in which case - `step` and `delta` are used. - convert_units : bool (optional) - units are converted to the MDAnalysis base format; ``None`` selects - the value of :data:`MDAnalysis.core.flags` ['convert_lengths']. - (see :ref:`flags-label`) - - Note - ---- - The keyword arguments set the **low-level attributes of the DCD** according - to the CHARMM format. The time between two frames would be `delta` * - `step`! For convenience, one can alternatively supply the `dt` keyword - (see above) to just tell the writer that it should record "There are `dt` - ps between each frame". - - The Writer will write the **unit cell information** to the DCD in a format - compatible with NAMD and older CHARMM versions, namely the unit cell - lengths in Angstrom and the angle cosines (see :class:`Timestep`). Newer - versions of CHARMM (at least c36b2) store the matrix of the box - vectors. Writing this matrix to a DCD is currently not supported (although - reading is supported with the :class:`DCDReader`); instead the angle - cosines are written, which *might make the DCD file unusable in CHARMM - itself*. See `Issue 187`_ for further information. - - The writing behavior of the :class:`DCDWriter` is identical to that of the - DCD molfile plugin of VMD with the exception that by default it will use - AKMA time units. - - Example - ------- - Typical usage:: - - with DCDWriter("new.dcd", u.atoms.n_atoms) as w: - for ts in u.trajectory - w.write_next_timestep(ts) - - Keywords are available to set some of the low-level attributes of the DCD. +class DCDReader(base.ReaderBase): + """Reader for the DCD format. + DCD is used by NAMD, CHARMM and LAMMPS as the default trajectory format. + The DCD file format is not well defined. In particular, NAMD and CHARMM use + it differently. Currently, MDAnalysis tries to guess the correct **format + for the unitcell representation** but it can be wrong. **Check the unitcell + dimensions**, especially for triclinic unitcells (see `Issue 187`_). DCD + trajectories produced by CHARMM and NAMD( >2.5) record time in AKMA units. + If other units have been recorded (e.g., ps) then employ the configurable + :class:MDAnalysis.coordinates.LAMMPS.DCDReader and set the time unit as a + optional argument. You can find a list of units used in the DCD formats on + the MDAnalysis `wiki`_. + + + MDAnalysis always uses ``(*A*, *B*, *C*, *alpha*, *beta*, *gamma*)`` to + represent the unit cell. Lengths *A*, *B*, *C* are in the MDAnalysis length + unit (Å), and angles are in degrees. + + The ordering of the angles in the unitcell is the same as in recent + versions of VMD's DCDplugin_ (2013), namely the `X-PLOR DCD format`_: The + original unitcell is read as ``[A, gamma, B, beta, alpha, C]`` from the DCD + file. If any of these values are < 0 or if any of the angles are > 180 + degrees then it is assumed it is a new-style CHARMM unitcell (at least + since c36b2) in which box vectors were recorded. + + .. warning:: + The DCD format is not well defined. Check your unit cell + dimensions carefully, especially when using triclinic boxes. + Different software packages implement different conventions and + MDAnalysis is currently implementing the newer NAMD/VMD convention + and tries to guess the new CHARMM one. Old CHARMM trajectories might + give wrong unitcell values. For more details see `Issue 187`_. + + .. _`X-PLOR DCD format`: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187 - + .. _DCDplugin: http://www.ks.uiuc.edu/Research/vmd/plugins/doxygen/dcdplugin_8c-source.html#l00947 + .. _wiki: https://github.com/MDAnalysis/mdanalysis/wiki/FileFormats#dcd """ format = 'DCD' - multiframe = True flavor = 'CHARMM' units = {'time': 'AKMA', 'length': 'Angstrom'} - def __init__(self, filename, n_atoms, start=0, step=1, - delta=mdaunits.convert(1., 'ps', 'AKMA'), dt=None, - remarks="Created by DCDWriter", convert_units=None): - if n_atoms == 0: - raise ValueError("DCDWriter: no atoms in output trajectory") - elif n_atoms is None: - # probably called from MDAnalysis.Writer() so need to give user a gentle heads up... - raise ValueError("DCDWriter: REQUIRES the number of atoms in the 'n_atoms' argument\n" + - " " * len("ValueError: ") + - "For example: n_atoms=universe.atoms.n_atoms") - self.filename = filename - # convert length and time to base units on the fly? - self.convert_units = flags['convert_lengths'] if convert_units is None \ - else convert_units - self.n_atoms = n_atoms - - self.frames_written = 0 - self.start = start - if dt is not None: - if dt > 0: - # ignore step and delta - self.step = 1 - self.delta = mdaunits.convert(dt, 'ps', 'AKMA') - else: - raise ValueError("DCDWriter: dt must be > 0, not {0}".format(dt)) - else: - self.step = step - self.delta = delta - self.dcdfile = open(self.filename, 'wb') - self.remarks = remarks - self._write_dcd_header(self.n_atoms, self.start, self.step, self.delta, self.remarks) - - def _dcd_header(self): - """Returns contents of the DCD header C structure:: - typedef struct { - fio_fd fd; // FILE * - fio_size_t header_size; // size_t == sizeof(int) - int natoms; - int nsets; - int setsread; - int istart; - int nsavc; - double delta; - int nfixed; - int *freeind; - float *fixedcoords; - int reverse; - int charmm; - int first; - int with_unitcell; - } dcdhandle; - - .. deprecated:: 0.7.5 - This function only exists for debugging purposes and might - be removed without notice. Do not rely on it. - + def __init__(self, filename, convert_units=True, dt=None, **kwargs): """ - # was broken (no idea why [orbeckst]), see Issue 27 - # 'PiiiiiidiPPiiii' should be the unpack string according to the struct. - # struct.unpack("LLiiiiidiPPiiii",self._dcd_C_str) - # seems to do the job on Mac OS X 10.6.4 ... but I have no idea why, - # given that the C code seems to define them as normal integers - desc = [ - 'file_desc', 'header_size', 'natoms', 'nsets', 'setsread', 'istart', - 'nsavc', 'delta', 'nfixed', 'freeind_ptr', 'fixedcoords_ptr', - 'reverse', 'charmm', 'first', 'with_unitcell'] - return dict(zip(desc, struct.unpack("LLiiiiidiPPiiii", self._dcd_C_str))) - - def write_next_timestep(self, ts=None): - """Write a new timestep to the DCD file. - Parameters ---------- - ts : `Timestep` (optional) - `Timestep` object containing coordinates to be written to DCD file; - by default it uses the current `Timestep` associated with the - Writer. - - Raises - ------ - :exc:`ValueError` - if wrong number of atoms supplied - :exc:`~MDAnalysis.NoDataError` - if no coordinates to be written. - + filename : str + trajectory filename + convert_units : bool (optional) + convert units to MDAnalysis units + dt : float (optional) + overwrite time delta stored in DCD + **kwargs : dict + General reader arguments. - .. versionchanged:: 0.7.5 - Raises :exc:`ValueError` instead of generic :exc:`Exception` - if wrong number of atoms supplied and :exc:`~MDAnalysis.NoDataError` - if no coordinates to be written. + .. versionchanged:: 0.17.0 + Changed to use libdcd.pyx library and removed the correl function """ - if ts is None: - try: - ts = self.ts - except AttributeError: - raise NoDataError("DCDWriter: no coordinate data to write to trajectory file") - if not ts.n_atoms == self.n_atoms: - raise ValueError("DCDWriter: Timestep does not have the correct number of atoms") - 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.from_timestep(ts) # wrap in a new fortran formatted Timestep + super(DCDReader, self).__init__( + filename, convert_units=convert_units, **kwargs) + self._file = DCDFile(self.filename) + self.n_atoms = self._file.header['natoms'] + + delta = mdaunits.convert(self._file.header['delta'], + self.units['time'], 'ps') + if dt is None: + dt = delta * self._file.header['nsavc'] + self.skip_timestep = self._file.header['nsavc'] + + self._ts_kwargs['dt'] = dt + self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) + frame = self._file.read() + # reset trajectory + if self._file.n_frames > 1: + self._file.seek(1) + else: + self._file.seek(0) + self._frame = 0 + self.ts = self._frame_to_ts(frame, self.ts) + # these should only be initialized once + self.ts.dt = dt if self.convert_units: - pos = self.convert_pos_to_native(ts._pos, - inplace=False) # possibly make a copy to avoid changing the trajectory - self._write_next_frame(pos[:, 0], pos[:, 1], pos[:, 2], unitcell) - self.frames_written += 1 - - def convert_dimensions_to_unitcell(self, ts, _ts_order=Timestep._ts_order): - """Read dimensions from timestep *ts* and return appropriate native unitcell. - - See Also - -------- - :attr:`Timestep.dimensions` - """ - # 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 = 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 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]) - np.put(unitcell, _ts_order, np.concatenate([lengths, angles])) - return unitcell + self.convert_pos_from_native(self.ts.dimensions[:3]) def close(self): - """Close trajectory and flush buffers.""" - if hasattr(self, 'dcdfile') and self.dcdfile is not None: - self._finish_dcd_write() - self.dcdfile.close() - self.dcdfile = None - - -class DCDReader(base.ReaderBase): - """Reads from a DCD file - - :Data: - ts - :class:`Timestep` object containing coordinates of current frame - - :Methods: - ``dcd = DCD(dcdfilename)`` - open dcd file and read header - ``len(dcd)`` - return number of frames in dcd - ``for ts in dcd:`` - iterate through trajectory - ``for ts in dcd[start:stop:skip]:`` - iterate through a trajectory - ``dcd[i]`` - random access into the trajectory (i corresponds to frame number) - ``data = dcd.timeseries(...)`` - retrieve a subset of coordinate information for a group of atoms - - Note - ---- - The DCD file format is not well defined. In particular, NAMD and CHARMM use - it differently. Currently, MDAnalysis tries to guess the correct format for - the unitcell representation but it can be wrong. **Check the unitcell - dimensions**, especially for triclinic unitcells (see `Issue 187`_ and - :attr:`Timestep.dimensions`). A second potential issue are the units of - time (TODO). - - - .. versionchanged:: 0.9.0 - The underlying DCD reader (written in C and derived from the - catdcd/molfile plugin code of VMD) is now reading the unitcell in - NAMD ordering: ``[A, B, C, sin(gamma), sin(beta), - sin(alpha)]``. See `Issue 187`_ for further details. + """close reader""" + self._file.close() - .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187 - - .. versionchanged:: 0.11.0 - Frames now 0-based instead of 1-based. - Native frame number read into `ts._frame`. - Removed `skip` keyword and functionality. - - """ - format = 'DCD' - flavor = 'CHARMM' - units = {'time': 'AKMA', 'length': 'Angstrom'} - _Timestep = Timestep - - def __init__(self, dcdfilename, **kwargs): - super(DCDReader, self).__init__(dcdfilename, **kwargs) - - self.dcdfilename = self.filename # dcdfilename is legacy - self.dcdfile = None # set right away because __del__ checks - - # Issue #32: segfault if dcd is 0-size - # Hack : test here... (but should be fixed in dcd.c) - stats = os.stat(self.filename) - if stats.st_size == 0: - raise IOError(errno.EIO, "DCD file is zero size", self.filename) - - self.dcdfile = open(self.filename, 'rb') - self.n_atoms = 0 - self.n_frames = 0 - self.fixed = 0 - self.periodic = False - - # This reads skip_timestep and delta from header - self._read_dcd_header() - - # Convert delta to ps - delta = mdaunits.convert(self.delta, self.units['time'], 'ps') - - self._ts_kwargs.setdefault('dt', self.skip_timestep * delta) - self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) - # Read in the first timestep - self._read_next_timestep() - - def _dcd_header(self): # pragma: no cover - """Returns contents of the DCD header C structure:: - typedef struct { - fio_fd fd; // FILE * - fio_size_t header_size; // size_t == sizeof(int) - int natoms; - int nsets; - int setsread; - int istart; - int nsavc; - double delta; - int nfixed; - int *freeind; - float *fixedcoords; - int reverse; - int charmm; - int first; - int with_unitcell; - } dcdhandle; - - .. deprecated:: 0.7.5 - This function only exists for debugging purposes and might - be removed without notice. Do not rely on it. - - """ - # was broken (no idea why [orbeckst]), see Issue 27 - # 'PiiiiiidiPPiiii' should be the unpack string according to the struct. - # struct.unpack("LLiiiiidiPPiiii",self._dcd_C_str) - # seems to do the job on Mac OS X 10.6.4 ... but I have no idea why, - # given that the C code seems to define them as normal integers - desc = [ - 'file_desc', 'header_size', 'natoms', 'nsets', 'setsread', 'istart', - 'nsavc', 'delta', 'nfixed', 'freeind_ptr', 'fixedcoords_ptr', 'reverse', - 'charmm', 'first', 'with_unitcell'] - return dict(zip(desc, struct.unpack("LLiiiiidiPPiiii", self._dcd_C_str))) + @property + def n_frames(self): + """number of frames in trajectory""" + return len(self._file) def _reopen(self): - self.ts.frame = -1 - self._reset_dcd_read() + """reopen trajectory""" + self.ts.frame = 0 + self._frame = -1 + self._file.close() + self._file.open('r') + + def _read_frame(self, i): + """read frame i""" + self._frame = i - 1 + self._file.seek(i) + return self._read_next_timestep() def _read_next_timestep(self, ts=None): - """Read the next frame - - .. versionchanged 0.11.0:: - Native frame read into ts._frame, ts.frame naively iterated - """ + """copy next frame into timestep""" + if self._frame == self.n_frames - 1: + raise IOError('trying to go over trajectory limit') if ts is None: - ts = self.ts - ts._frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1) - ts.frame += 1 + # use a copy to avoid that ts always points to the same reference + # removing this breaks lammps reader + ts = self.ts.copy() + frame = self._file.read() + self._frame += 1 + ts = self._frame_to_ts(frame, ts) + self.ts = ts return ts - def _read_frame(self, frame): - """Skip to frame and read + def Writer(self, filename, n_atoms=None, **kwargs): + """Return writer for trajectory format""" + if n_atoms is None: + n_atoms = self.n_atoms + return DCDWriter( + filename, + n_atoms=n_atoms, + dt=self.ts.dt, + convert_units=self.convert_units, + **kwargs) + + def _frame_to_ts(self, frame, ts): + """convert a dcd-frame to a :class:`TimeStep`""" + ts.frame = self._frame + ts.time = ts.frame * self.ts.dt + ts.data['step'] = self._file.tell() + + # The original unitcell is read as ``[A, gamma, B, beta, alpha, C]`` + _ts_order = [0, 2, 5, 4, 3, 1] + uc = np.take(frame.unitcell, _ts_order) + + pi_2 = np.pi / 2 + if (-1.0 <= uc[3] <= 1.0) and (-1.0 <= uc[4] <= 1.0) and ( + -1.0 <= uc[5] <= 1.0): + # This file was generated by Charmm, or by NAMD > 2.5, with the + # angle cosines of the periodic cell angles written to the DCD + # file. This formulation improves rounding behavior for orthogonal + # cells so that the angles end up at precisely 90 degrees, unlike + # acos(). (changed in MDAnalysis 0.9.0 to have NAMD ordering of the + # angles; see Issue 187) */ + uc[3] = 90.0 - np.arcsin(uc[3]) * 90.0 / pi_2 + uc[4] = 90.0 - np.arcsin(uc[4]) * 90.0 / pi_2 + uc[5] = 90.0 - np.arcsin(uc[5]) * 90.0 / pi_2 + # heuristic sanity check: uc = A,B,C,alpha,beta,gamma + elif np.any(uc < 0.) or np.any(uc[3:] > 180.): + # might be new CHARMM: box matrix vectors + H = frame.unitcell.copy() + e1, e2, e3 = H[[0, 1, 3]], H[[1, 2, 4]], H[[3, 4, 5]] + uc = triclinic_box(e1, e2, e3) + else: + # This file was likely generated by NAMD 2.5 and the periodic cell + # angles are specified in degrees rather than angle cosines. + pass + + ts.dimensions = uc + ts.positions = frame.xyz + + if self.convert_units: + self.convert_pos_from_native(ts.dimensions[:3]) + self.convert_pos_from_native(ts.positions) - .. versionchanged:: 0.11.0 - Native frame read into ts._frame, ts.frame naively set to frame - """ - self._jump_to_frame(frame) - ts = self.ts - ts._frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1) - ts.frame = frame return ts - def timeseries(self, asel=None, start=None, stop=None, step=None, skip=None, + @property + def dimensions(self): + """unitcell dimensions (*A*, *B*, *C*, *alpha*, *beta*, *gamma*) + """ + return self.ts.dimensions + + @property + def dt(self): + """timestep between frames""" + return self.ts.dt + + def timeseries(self, + asel=None, + start=None, + stop=None, + step=None, + skip=None, format='afc'): """Return a subset of coordinate data for an AtomGroup @@ -530,16 +275,17 @@ def timeseries(self, asel=None, start=None, stop=None, step=None, skip=None, The :class:`~MDAnalysis.core.groups.AtomGroup` to read the coordinates from. Defaults to None, in which case the full set of coordinate data is returned. - start : int (optional) - Begin reading the trajectory at frame index `start` (where 0 is the index - of the first frame in the trajectory); the default ``None`` starts - at the beginning. + start : int (optional) + Begin reading the trajectory at frame index `start` (where 0 is the + index of the first frame in the trajectory); the default ``None`` + starts at the beginning. stop : int (optional) - End reading the trajectory at frame index `stop`-1, i.e, `stop` is excluded. - The trajectory is read to the end with the default ``None``. + End reading the trajectory at frame index `stop`-1, i.e, `stop` is + excluded. The trajectory is read to the end with the default + ``None``. step : int (optional) - Step size for reading; the default ``None`` is equivalent to 1 and means to - read every frame. + Step size for reading; the default ``None`` is equivalent to 1 and + means to read every frame. format : str (optional) the order/shape of the return data array, corresponding to (a)tom, (f)rame, (c)oordinates all six combinations @@ -550,106 +296,122 @@ def timeseries(self, asel=None, start=None, stop=None, step=None, skip=None, .. deprecated:: 0.16.0 `skip` has been deprecated in favor of the standard keyword `step`. - """ if skip is not None: step = skip - warnings.warn("Skip is deprecated and will be removed in" - "in 1.0. Use step instead.", - category=DeprecationWarning) + warnings.warn( + "Skip is deprecated and will be removed in" + "in 1.0. Use step instead.", + category=DeprecationWarning) start, stop, step = self.check_slice_indices(start, stop, step) if asel is not None: if len(asel) == 0: - raise NoDataError("Timeseries requires at least one atom to analyze") + raise NoDataError( + "Timeseries requires at least one atom to analyze") atom_numbers = list(asel.indices) else: atom_numbers = list(range(self.n_atoms)) - if len(format) != 3 and format not in ['afc', 'acf', 'caf', 'cfa', 'fac', 'fca']: - raise ValueError("Invalid timeseries format") - # Check if the atom numbers can be grouped for efficiency, then we can read partial buffers - # from trajectory file instead of an entire timestep - # XXX needs to be implemented - return self._read_timeseries(atom_numbers, start, stop, step, format) + frames = self._file.readframes( + start, stop, step, order=format, indices=atom_numbers) + return frames.xyz - def close(self): - if self.dcdfile is not None: - self._finish_dcd_read() - self.dcdfile.close() - self.dcdfile = None - def Writer(self, filename, **kwargs): - """Returns a DCDWriter for `filename` with the same parameters as this DCD. +class DCDWriter(base.WriterBase): + """DCD Writer class - Defaults for all values are obtained from the `DCDReader` itself but - all values can be changed through keyword arguments. + The writer follows recent NAMD/VMD convention for the unitcell (box lengths + in Å and angle-cosines, ``[A, cos(gamma), B, cos(beta), cos(alpha), C]``) + and writes positions in Å and time in AKMA time units. - Parameters + """ + format = 'DCD' + multiframe = True + flavor = 'NAMD' + units = {'time': 'AKMA', 'length': 'Angstrom'} + + def __init__(self, + filename, + n_atoms, + convert_units=True, + step=1, + dt=1, + remarks='', + nsavc=1, + **kwargs): + """Parameters ---------- filename : str - filename of the output DCD trajectory - n_atoms : int (optional) - number of atoms - start : int (optional) - number of the first recorded MD step + filename of trajectory + n_atoms : int + number of atoms to be written + convert_units : bool (optional) + convert from MDAnalysis units to format specific units step : int (optional) - indicate that `step` **MD steps (!)** make up one trajectory frame - delta : float (optional) - **MD integrator time step (!)**, in AKMA units + number of steps between frames to be written dt : float (optional) - **Override** `step` and `delta` so that the DCD records that - `dt` ps lie between two frames. (It sets `step` `` = 1`` and - `delta` `` = AKMA(dt)``.) The default is ``None``, in which case - `step` and `delta` are used. - + time between two frames. If ``None`` guess from first written + :class:`TimeStep` remarks : str (optional) - string that is stored in the DCD header - - Returns - ------- - :class:`DCDWriter` + remarks to be stored in DCD. Shouldn't be more then 240 characters + nsavc : int (optional) + DCD files can also store the number of integrator time steps that + correspond to the interval between two frames as nsavc (i.e., every + how many MD steps is a frame saved to the DCD). By default, this + number is just set to one and this should be sufficient for almost + all cases but if required, nsavc can be changed. + **kwargs : dict + General writer arguments + """ + self.filename = filename + self._convert_units = convert_units + if n_atoms is None: + raise ValueError("n_atoms argument is required") + self.n_atoms = n_atoms + self._file = DCDFile(self.filename, 'w') + self.step = step + self.dt = dt + dt = mdaunits.convert(dt, 'ps', self.units['time']) + delta = float(dt) / nsavc + self._file.write_header( + remarks=remarks, + natoms=self.n_atoms, + nsavc=nsavc, + delta=delta, + is_periodic=1, + istart=0) + + def write_next_timestep(self, ts): + """Write timestep object into trajectory. - Note - ---- - The keyword arguments set the low-level attributes of the DCD according - to the CHARMM format. The time between two frames would be `delta` * - `step`! - - Here `step` is really the number of MD integrator time steps that - occured after this frame, including the frame itself that is the - coordinate snapshot and `delta` is the integrator stime step. The DCD - file format contains this information so it needs to be provided here. - + Parameters + ---------- + ts: TimeStep See Also -------- - :class:`DCDWriter` - + :meth:`DCDWriter.write` takes a more general input """ - n_atoms = kwargs.pop('n_atoms', self.n_atoms) - kwargs.setdefault('start', self.start_timestep) - kwargs.setdefault('step', self.skip_timestep) - kwargs.setdefault('delta', self.delta) - kwargs.setdefault('remarks', self.remarks) - # dt keyword is simply passed through if provided - return DCDWriter(filename, n_atoms, **kwargs) + xyz = ts.positions.copy() + dimensions = ts.dimensions.copy() - @property - def dt(self): - """Time between two trajectory frames in picoseconds.""" - return self.ts.dt + if self._convert_units: + xyz = self.convert_pos_to_native(xyz, inplace=True) + dimensions = self.convert_dimensions_to_unitcell(ts, inplace=True) + + # we only support writing charmm format unit cell info + # The DCD unitcell is written as ``[A, gamma, B, beta, alpha, C]`` + _ts_order = [0, 5, 1, 4, 3, 2] + box = np.take(dimensions, _ts_order) + self._file.write(xyz=xyz, box=box) -DCDReader._read_dcd_header = types.MethodType(_dcdmodule.__read_dcd_header, None, DCDReader) -DCDReader._read_next_frame = types.MethodType(_dcdmodule.__read_next_frame, None, DCDReader) -DCDReader._jump_to_frame = types.MethodType(_dcdmodule.__jump_to_frame, None, DCDReader) -DCDReader._reset_dcd_read = types.MethodType(_dcdmodule.__reset_dcd_read, None, DCDReader) -DCDReader._finish_dcd_read = types.MethodType(_dcdmodule.__finish_dcd_read, None, DCDReader) -DCDReader._read_timeseries = types.MethodType(_dcdmodule.__read_timeseries, None, DCDReader) + def close(self): + """close trajectory""" + self._file.close() -DCDWriter._write_dcd_header = types.MethodType(_dcdmodule.__write_dcd_header, None, DCDWriter) -DCDWriter._write_next_frame = types.MethodType(_dcdmodule.__write_next_frame, None, DCDWriter) -DCDWriter._finish_dcd_write = types.MethodType(_dcdmodule.__finish_dcd_write, None, DCDWriter) + def __del__(self): + self.close() diff --git a/package/MDAnalysis/coordinates/include/correl.h b/package/MDAnalysis/coordinates/include/correl.h deleted file mode 100644 index 7eb8cbae6d6..00000000000 --- a/package/MDAnalysis/coordinates/include/correl.h +++ /dev/null @@ -1,195 +0,0 @@ -/* -*- Mode: C; tab-width: 4; indent-tabs-mode:nil; -*- */ -/* vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 */ -/* - - MDAnalysis --- http://www.mdanalysis.org - Copyright (c) 2006-2016 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. - - 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 - -*/ - -#ifndef CORREL_H -#define CORREL_H - -#include -/* Python.h for 'typedef int Py_intptr_t;' (fixes Issue 19) */ -#include - -static void -copyseries(int frame, char *data, const Py_intptr_t *strides, - const float *tempX, const float *tempY, const float *tempZ, - const char* datacode, int numdata, const int* atomlist, const int* atomcounts, - int lowerb, double* aux) -{ - char code; - int index1 = 0, index2 = 0, index3 = 0, index4 = 0; - double x1, x2, y1, y2, z1, z2, x3, y3, z3, aux1, aux2; - int i = 0, j = 0, atomno = 0; - int dataset = 0; - int stride0, stride1; - - stride0 = strides[0]; - stride1 = strides[1]; - - /* If I eventually switch to using frame,property ordering for timeseries data - stride0 = strides[1]; - stride1 = strides[0]; - */ - - for (i=0;i 0) || (aux2 > 0 && aux1 < 0)) { - aux1 *= -1; - } - // Check if the dihedral has wrapped around 2 pi - aux2 = *(double*)(data + dataset*stride0 + (frame-1)*stride1); - if (fabs(aux1-aux2) > M_PI) { - if (aux1 > 0) { aux1 -= 2*M_PI; } - else { aux1 += 2*M_PI; } - } - *(double*)(data + dataset++*stride0 + frame*stride1) = aux1; - break; - case 'w': - /* dipole orientation of 3-site water: ^ d - index1 = oxygen, index2, index3 = hydrogen | - returns d ,O, - d = rO - (rH1 + rH2)/2 H | H - | - */ - index1 = atomlist[atomno++]-lowerb; // O - index2 = atomlist[atomno++]-lowerb; // H1 - index3 = atomlist[atomno++]-lowerb; // H2 - x1 = tempX[index1] - 0.5*(tempX[index2] + tempX[index3]); // dx - y1 = tempY[index1] - 0.5*(tempY[index2] + tempY[index3]); // dy - z1 = tempZ[index1] - 0.5*(tempZ[index2] + tempZ[index3]); // dz - *(double*)(data + dataset++*stride0 + frame*stride1) = x1; - *(double*)(data + dataset++*stride0 + frame*stride1) = y1; - *(double*)(data + dataset++*stride0 + frame*stride1) = z1; - break; - } - } -} - -// This accounts for periodic boundary conditions -// taken from MMTK -#define distance_vector_2(d, r1, r2, data) \ - { \ - double xh = 0.5*(data)[0]; \ - double yh = 0.5*(data)[1]; \ - double zh = 0.5*(data)[2]; \ - d[0] = r2[0]-r1[0]; \ - if (d[0] > xh) d[0] -= (data)[0]; \ - if (d[0] <= -xh) d[0] += (data)[0]; \ - d[1] = r2[1]-r1[1]; \ - if (d[1] > yh) d[1] -= (data)[1]; \ - if (d[1] <= -yh) d[1] += (data)[1]; \ - d[2] = r2[2]-r1[2]; \ - if (d[2] > zh) d[2] -= (data)[2]; \ - if (d[2] <= -zh) d[2] += (data)[2]; \ - } - -#endif diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c deleted file mode 100644 index 782c163b35a..00000000000 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ /dev/null @@ -1,949 +0,0 @@ -/* -*- Mode: C; tab-width: 4; indent-tabs-mode:nil; -*- */ -/* vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 */ -/* - - MDAnalysis --- http://www.mdanalysis.org - Copyright (c) 2006-2016 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. - - 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 - -*/ - - -/* I have to fix error handling in these functions */ - -#include -#include - -#include "readdcd.h" - -typedef struct { - fio_fd fd; - fio_size_t header_size; - int natoms; - int nsets; - int setsread; - int istart; - int nsavc; - double delta; - int nfixed; - int *freeind; - float *fixedcoords; - int reverse; - int charmm; - int first; - int with_unitcell; -} dcdhandle; - - -int jump_to_frame(dcdhandle *dcd, int frame); - -static PyObject * -__write_dcd_header(PyObject *self, PyObject *args) -{ - /* At this point we assume the file has been opened for writing */ - PyObject *temp = NULL; - dcdhandle *dcd = NULL; - fio_fd fd; - int rc = 0; - int natoms = 0; - const char *remarks = "DCD"; - int istart = 0; /* starting timestep of DCD file */ - int nsavc = 1; /* number of timesteps between written DCD frames */ - double delta = 1.0; /* length of a timestep */ - int with_unitcell = 1; /* contains unit cell information (charmm format) */ - int charmm = DCD_IS_CHARMM; /* charmm-formatted DCD file */ - if (with_unitcell) - charmm |= DCD_HAS_EXTRA_BLOCK; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "Oi|iids", &self, &natoms, &istart, &nsavc, &delta, &remarks) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "i|iids", &natoms, &istart, &nsavc, &delta, &remarks) ) - return NULL; - } - - // Get the file object from the class - if (!PyObject_HasAttrString(self, "dcdfile")) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "dcdfile")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); - return NULL; - } - - if (!PyFile_CheckExact(temp)) { - // Raise exception - PyErr_SetString(PyExc_TypeError, "dcdfile does not refer to a file object"); - Py_DECREF(temp); - return NULL; - } - fd = fileno(PyFile_AsFile(temp)); - // No longer need the reference to temp - Py_DECREF(temp); - - dcd = (dcdhandle *)malloc(sizeof(dcdhandle)); - memset(dcd, 0, sizeof(dcdhandle)); - dcd->fd = fd; - - if ((rc = write_dcdheader(dcd->fd, remarks, natoms, istart, nsavc, delta, with_unitcell, charmm)) < 0) { - PyErr_SetString(PyExc_IOError, "Cannot write header of DCD file"); - free(dcd); - return NULL; - } - - dcd->natoms = natoms; - dcd->nsets = 0; - dcd->istart = istart; - dcd->nsavc = nsavc; - dcd->delta = delta; - dcd->with_unitcell = with_unitcell; - dcd->charmm = charmm; - temp = PyCObject_FromVoidPtr(dcd, free); // Creates a New Reference - if (PyObject_SetAttrString(self, "_dcd_C_ptr", temp) == -1) { - // Raise exception - who knows what exception to raise?? - PyErr_SetString(PyExc_AttributeError, "Could not create attribute _dcd_C_ptr"); - Py_DECREF(temp); - return NULL; - } - Py_DECREF(temp); - - // For debugging purposes - temp = PyBuffer_FromMemory(dcd, sizeof(dcdhandle)); // Creates a New Reference - if (PyObject_SetAttrString(self, "_dcd_C_str", temp) == -1) { - // Raise exception - who knows what exception to raise?? - PyErr_SetString(PyExc_AttributeError, "Could not create attribute _dcd_C_str"); - Py_DECREF(temp); - return NULL; - } - Py_DECREF(temp); - - Py_INCREF(Py_None); - return Py_None; -} - -static PyObject * -__write_next_frame(PyObject *self, PyObject *args) -{ - dcdhandle *dcd; - PyObject *temp; - PyArrayObject *x, *y, *z, *uc; - int rc, curstep; - float* uc_array; - double unitcell[6]; - - if (!self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "OO!O!O!O!", &self, &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "O!O!O!O!", &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc) ) - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); - dcd->nsets++; - curstep = dcd->istart + dcd->nsets * dcd->nsavc; - - uc_array = (float*) uc->data; - unitcell[0] = uc_array[0]; /* A */ - unitcell[2] = uc_array[2]; /* B */ - unitcell[5] = uc_array[5]; /* C */ - /* write angle cosines with NAMD ordering [orbeckst] */ - /* (changed in MDAnalysis 0.9.0) */ - unitcell[4] = sin((M_PI_2 / 90.0 ) * (90.0 - uc_array[4])); /* cos(alpha) */ - unitcell[3] = sin((M_PI_2 / 90.0 ) * (90.0 - uc_array[3])); /* cos(beta) */ - unitcell[1] = sin((M_PI_2 / 90.0 ) * (90.0 - uc_array[1])); /* cos(gamma) */ - - if ((rc = write_dcdstep(dcd->fd, dcd->nsets, curstep, dcd->natoms, (float*)x->data, (float*)y->data, (float*)z->data, - dcd->with_unitcell ? unitcell: NULL, - dcd->charmm)) < 0) - { - PyErr_SetString(PyExc_IOError, "Could not write timestep to dcd file"); - return NULL; - } - - Py_INCREF(Py_None); - return Py_None; -} - -static PyObject * -__finish_dcd_write(PyObject *self, PyObject *args) -{ - //PyObject* temp; - //dcdhandle *dcd; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "") ) - return NULL; - } - - /*if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - free(dcd); - Py_DECREF(temp);*/ - - if ( PyObject_DelAttrString(self, "_dcd_C_ptr") == -1) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - Py_INCREF(Py_None); - return Py_None; -} - -static PyObject * -__read_dcd_header(PyObject *self, PyObject *args) -{ - - /* At this point we assume the file has been checked for existence and opened, and the DCD class - * contains a reference to the file (which can be used to retrieve the file descriptor) - */ - PyObject *temp; - fio_fd fd; - int rc; - char *remarks = NULL; - int len_remarks = 0; - dcdhandle *dcd = NULL; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "") ) - return NULL; - } - - // Get the file object from the class - if (!PyObject_HasAttrString(self, "dcdfile")) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "dcdfile")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); - return NULL; - } - - if (!PyFile_CheckExact(temp)) { - // Raise exception - PyErr_SetString(PyExc_TypeError, "dcdfile does not refer to a file object"); - goto error; - } - fd = fileno(PyFile_AsFile(temp)); - // No longer need the reference to temp - Py_DECREF(temp); - - dcd = (dcdhandle *)malloc(sizeof(dcdhandle)); - memset(dcd, 0, sizeof(dcdhandle)); - dcd->fd = fd; - - if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart, - &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind, - &dcd->fixedcoords, &dcd->reverse, &dcd->charmm, &remarks, &len_remarks))) { - // Raise exception - PyErr_SetString(PyExc_IOError, "Cannot read DCD header"); - goto error; - } - - /* - * Check that the file is big enough to really hold all the frames it claims to have - */ - { - off_t ndims, firstframesize, framesize, extrablocksize; - off_t filesize; - struct stat stbuf; - - extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; - ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3; - firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize; - framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) - + extrablocksize; - - // Get filesize from file descriptor - memset(&stbuf, 0, sizeof(struct stat)); - if (fstat(dcd->fd, &stbuf)) { - PyErr_SetString(PyExc_IOError, "Could not stat file"); - goto error; - } - - // Size of header - dcd->header_size = fio_ftell(dcd->fd); - - /* - * It's safe to use ftell, even though ftell returns a long, because the - * header size is < 4GB. - */ - filesize = stbuf.st_size - fio_ftell(dcd->fd) - firstframesize; - if (filesize < 0) { - PyErr_SetString(PyExc_IOError, "DCD file appears to contain no timesteps"); - goto error; - } - dcd->nsets = filesize / framesize + 1; - dcd->setsread = 0; - } - - // We are at the first frame - dcd->first = 1; - - temp = Py_BuildValue("s#", remarks, len_remarks); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "remarks", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute remarks"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("i", dcd->natoms); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "n_atoms", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute n_atoms"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("i", dcd->nsets); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "n_frames", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute n_frames"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("i", dcd->nfixed); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "fixed", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute fixed"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("i", dcd->istart); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "start_timestep", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute fixed"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("i", dcd->nsavc); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "skip_timestep", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute fixed"); - goto error; - } - Py_DECREF(temp); - temp = Py_BuildValue("d", dcd->delta); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "delta", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute fixed"); - goto error; - } - Py_DECREF(temp); - - temp = PyCObject_FromVoidPtr(dcd, NULL); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "_dcd_C_ptr", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute _dcd_C_ptr"); - goto error; - } - - if ((dcd->charmm & DCD_IS_CHARMM) && (dcd->charmm & DCD_HAS_EXTRA_BLOCK)) { - dcd->with_unitcell = 1; - temp = Py_True; - } else { - temp = Py_False; - } - Py_INCREF(temp); - if (PyObject_SetAttrString(self, "periodic", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute periodic"); - goto error; - } - Py_DECREF(temp); - - // For debugging purposes - temp = PyBuffer_FromMemory(dcd, sizeof(dcdhandle)); - if (temp == NULL) goto error; - if (PyObject_SetAttrString(self, "_dcd_C_str", temp) == -1) { - PyErr_SetString(PyExc_AttributeError, "Could not create attribute _dcd_C_str"); - goto error; - } - Py_DECREF(temp); - - Py_INCREF(Py_None); - return Py_None; - - error: - Py_XDECREF(temp); - if (dcd != NULL) free(dcd); - if (remarks != NULL) free(remarks); - return NULL; -} - -static PyObject * -__read_timeseries(PyObject *self, PyObject *args) -{ - PyObject *temp = NULL; - PyArrayObject *coord = NULL; - PyListObject *atoms = NULL; - int lowerb = 0, upperb = 0, range=0; - float *tempX = NULL, *tempY = NULL, *tempZ = NULL; - int rc; - int i, j, index; - int n_atoms = 0, n_frames = 0; - - /* Stop = -1 is incorrect and causes an error, look for a fix of this in line - 469 */ - int start = 0, stop = -1, step = 1, numskip = 0, remaining_frames=0; - dcdhandle *dcd = NULL; - int *atomlist = NULL; - npy_intp dimensions[3]; - float unitcell[6]; - const char* format = "afc"; - - if (!self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "OO!|iiis", &self, &PyList_Type, &atoms, &start, &stop, &step, &format) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "O!|iiis", &PyList_Type, &atoms, &start, &stop, &step, &format) ) - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); - - n_frames = ((stop-start) / step); - if ((stop-start) % step > 0) { n_frames++; } - //n_frames = dcd->nsets / skip; - n_atoms = PyList_Size((PyObject*)atoms); - if (n_atoms == 0) { - PyErr_SetString(PyExc_Exception, "No atoms passed into _read_timeseries function"); - return NULL; - } - atomlist = (int*)malloc(sizeof(int)*n_atoms); - memset(atomlist, 0, sizeof(int)*n_atoms); - - // Get the atom indexes - for (i=0;ifd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - - // Jump to starting frame - jump_to_frame(dcd, start); - - // Now read through trajectory and get the atom pos - tempX = (float*)malloc(sizeof(float)*range); - tempY = (float*)malloc(sizeof(float)*range); - tempZ = (float*)malloc(sizeof(float)*range); - if ((tempX == NULL) | (tempY == NULL) || (tempZ == NULL)) { - PyErr_SetString(PyExc_MemoryError, "Can't allocate temporary space for coordinate arrays"); - goto error; - } - - remaining_frames = stop-start; - for (i=0;i 1 && i>0) { - // Check if we have fixed atoms - // XXX not done - /* Figure out how many steps to step over, if step = n, np array - slicing treats this as skip over n-1, read the nth. */ - numskip = step -1; - /* If the number to skip is greater than the number of frames left - to be jumped over, just take one more step to reflect np slicing - if there is a remainder, guaranteed to have at least one more - frame. - */ - if(remaining_frames < numskip){ - numskip = 1; - } - rc = skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm, numskip); - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error skipping frame from DCD file"); - goto error; - } - } - // on first iteration, numskip == 0, first set is always read. - dcd->setsread += numskip; - //now read from subset - rc = read_dcdsubset(dcd->fd, dcd->natoms, lowerb, upperb, tempX, tempY, tempZ, - unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, - dcd->reverse, dcd->charmm); - dcd->first = 0; - dcd->setsread++; - remaining_frames = stop - dcd->setsread; - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file"); - goto error; - - } - - - // Copy into Numeric array only those atoms we are interested in - for (j=0;jdata + a*coord->strides[0] + b*coord->strides[1] + c*coord->strides[2]) - */ - if (strncasecmp(format, "afc", 3) == 0) { - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; - } else if (strncasecmp(format, "acf", 3) == 0) { - *(double*)(coord->data + j*coord->strides[0] + 0*coord->strides[1] + i*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + j*coord->strides[0] + 1*coord->strides[1] + i*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + j*coord->strides[0] + 2*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; - } else if (strncasecmp(format, "fac", 3) == 0) { - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; - } else if (strncasecmp(format, "fca", 3) == 0) { - *(double*)(coord->data + i*coord->strides[0] + 0*coord->strides[1] + j*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + i*coord->strides[0] + 1*coord->strides[1] + j*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + i*coord->strides[0] + 2*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; - } else if (strncasecmp(format, "caf", 3) == 0) { - *(double*)(coord->data + 0*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + 1*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + 2*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; - } else if (strncasecmp(format, "cfa", 3) == 0) { - *(double*)(coord->data + 0*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + 1*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + 2*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; - } - } - // Check if we've been interupted by the user - if (PyErr_CheckSignals() == 1) goto error; - } - - // Reset trajectory - rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - free(atomlist); - free(tempX); - free(tempY); - free(tempZ); - return PyArray_Return(coord); - - error: - // Reset trajectory - rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - Py_XDECREF(coord); - if (atomlist != NULL) free(atomlist); - if (tempX != NULL) free(tempX); - if (tempY != NULL) free(tempY); - if (tempZ != NULL) free(tempZ); - return NULL; -} - - -static PyObject * -__read_next_frame(PyObject *self, PyObject *args) -{ - dcdhandle *dcd; - PyObject *temp; - PyArrayObject *x, *y, *z, *uc; - int skip=1; - int rc,numskip; - float* unitcell; - float alpha, beta, gamma; - - if (!self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "OO!O!O!O!|i", &self, &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc, &skip) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "O!O!O!O!|i", &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc, &skip) ) - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); - - unitcell = (float*) uc->data; - unitcell[0] = unitcell[2] = unitcell[5] = 0.0f; - unitcell[1] = unitcell[3] = unitcell[4] = 90.0f; - - /* Check for EOF here; that way all EOF's encountered later must be errors */ - if (dcd->setsread == dcd->nsets) { - // Is this an exception? - PyErr_SetString(PyExc_IOError, "End of file reached for dcd file"); - return NULL; - //return MOLFILE_EOF; - } - if (skip > 1) { - if (dcd->first && dcd->nfixed) { - /* We can't just skip it because we need the fixed atom coordinates */ - rc = read_dcdstep(dcd->fd, dcd->natoms, (float*)x->data, (float*)y->data, (float*)z->data, - unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, - dcd->reverse, dcd->charmm); - dcd->first = 0; - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error reading first frame from DCD file"); - //fprintf(stderr, "read_dcdstep returned %d\n", rc); - return NULL; - //return MOLFILE_ERROR; - } - dcd->setsread++; - temp = Py_BuildValue("i", dcd->setsread); - return temp; - //return rc; /* XXX this needs to be updated */ - } - dcd->first = 0; - // Figure out how many steps to skip - numskip = skip - (dcd->setsread % skip) - 1; - /* XXX this needs to be changed */ - rc = skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm, numskip); - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error skipping frame from DCD file"); - //fprintf(stderr, "read_dcdstep returned %d\n", rc); - return NULL; - //return MOLFILE_ERROR; - } - dcd->setsread+=numskip; - } - rc = read_dcdstep(dcd->fd, dcd->natoms, (float*)x->data, (float*)y->data, (float*)z->data, unitcell, - dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, - dcd->reverse, dcd->charmm); - dcd->first = 0; - dcd->setsread++; - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file"); - //fprintf(stderr, "read_dcdstep returned %d\n", rc); - return NULL; - //return MOLFILE_ERROR; - } - - if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 && - unitcell[3] >= -1.0 && unitcell[3] <= 1.0 && - unitcell[4] >= -1.0 && unitcell[4] <= 1.0) { - /* This file was generated by Charmm, or by NAMD > 2.5, with the angle */ - /* cosines of the periodic cell angles written to the DCD file. */ - /* This formulation improves rounding behavior for orthogonal cells */ - /* so that the angles end up at precisely 90 degrees, unlike acos(). */ - /* (changed in MDAnalysis 0.9.0 to have NAMD ordering of the angles; */ - /* see Issue 187) */ - alpha = 90.0 - asin(unitcell[4]) * 90.0 / M_PI_2; - beta = 90.0 - asin(unitcell[3]) * 90.0 / M_PI_2; - gamma = 90.0 - asin(unitcell[1]) * 90.0 / M_PI_2; - } else { - /* This file was likely generated by NAMD 2.5 and the periodic cell */ - /* angles are specified in degrees rather than angle cosines. */ - alpha = unitcell[4]; - beta = unitcell[3]; - gamma = unitcell[1]; - } - unitcell[4] = alpha; - unitcell[3] = beta; - unitcell[1] = gamma; - - // Return the frame read - temp = Py_BuildValue("i", dcd->setsread); - return temp; -} - -static PyObject * -__finish_dcd_read(PyObject *self, PyObject *args) -{ - PyObject* temp; - dcdhandle *dcd; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "") ) - return NULL; - } - - if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - - close_dcd_read(dcd->freeind, dcd->fixedcoords); - free(dcd); - Py_DECREF(temp); - Py_INCREF(Py_None); - return Py_None; -} -int jump_to_frame(dcdhandle *dcd, int frame) -{ - int rc; - if (frame > dcd->nsets) { - return -1; - } - // Calculate file offset - { - off_t extrablocksize, ndims, firstframesize, framesize; - off_t pos; - extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; - ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3; - firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize; - framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) - + extrablocksize; - // Use zero indexing - if (frame == 0) { - pos = dcd->header_size; - dcd->first = 1; - } - else { - pos = dcd->header_size + firstframesize + framesize * (frame-1); - dcd->first = 0; - } - rc = fio_fseek(dcd->fd, pos, FIO_SEEK_SET); - } - dcd->setsread = frame; - return rc; -} - -static PyObject * -__jump_to_frame(PyObject *self, PyObject *args) -{ - PyObject* temp; - dcdhandle *dcd; - int frame; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "Oi", &self, &frame) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "i", &frame) ) - return NULL; - } - - if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); - - /*if (frame > dcd->nsets) { - PyErr_SetString(PyExc_IndexError, "Invalid frame"); - return NULL; - } - - // Calculate file offset - { - off_t extrablocksize, ndims, firstframesize, framesize; - off_t pos; - extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; - ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3; - firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize; - framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) - + extrablocksize; - // Use zero indexing - if (frame == 0) { - pos = dcd->header_size; - } - else { - pos = dcd->header_size + firstframesize + framesize * (frame-1); - } - rc = fio_fseek(dcd->fd, pos, FIO_SEEK_SET); - } - dcd->setsread = frame; - dcd->first = 0;*/ - jump_to_frame(dcd, frame); - - temp = Py_BuildValue("i", dcd->setsread); - return temp; -} - -static PyObject * -__reset_dcd_read(PyObject *self, PyObject *args) -{ - PyObject* temp; - dcdhandle *dcd; - int rc; - - if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "") ) - return NULL; - } - - if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } - - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - Py_DECREF(temp); - Py_INCREF(Py_None); - return Py_None; -} - -static PyMethodDef DCDMethods[] = { - {"__write_dcd_header", __write_dcd_header, METH_VARARGS, "Write a DCD header."}, - {"__write_next_frame", __write_next_frame, METH_VARARGS, "Write the next timestep."}, - {"__finish_dcd_write", __finish_dcd_write, METH_VARARGS, "Clean up and data for handling dcd file."}, - {"__read_dcd_header", __read_dcd_header, METH_VARARGS, "Read in a DCD header."}, - {"__read_next_frame", __read_next_frame, METH_VARARGS, "Read in the next timestep."}, - {"__jump_to_frame", __jump_to_frame, METH_VARARGS, "Jump to specified timestep."}, - {"__reset_dcd_read", __reset_dcd_read, METH_VARARGS, "Reset dcd file reading."}, - {"__finish_dcd_read", __finish_dcd_read, METH_VARARGS, "Clean up any data for handling dcd file."}, - {"__read_timeseries", __read_timeseries, METH_VARARGS, "Return a Numpy array of the coordinates for a set of atoms for the whole trajectory."}, - {NULL, NULL, 0, NULL} /* Sentinel */ -}; - -PyMODINIT_FUNC -init_dcdmodule(void) -{ - (void) Py_InitModule("_dcdmodule", DCDMethods); - import_array(); -} diff --git a/package/MDAnalysis/lib/formats/__init__.py b/package/MDAnalysis/lib/formats/__init__.py index 5a851008ea8..cf9eb30c6ed 100644 --- a/package/MDAnalysis/lib/formats/__init__.py +++ b/package/MDAnalysis/lib/formats/__init__.py @@ -21,5 +21,6 @@ # from __future__ import absolute_import from . import libmdaxdr +from . import libdcd -__all__ = ['libmdaxdr'] +__all__ = ['libmdaxdr', 'libdcd'] diff --git a/package/MDAnalysis/coordinates/include/endianswap.h b/package/MDAnalysis/lib/formats/include/endianswap.h similarity index 100% rename from package/MDAnalysis/coordinates/include/endianswap.h rename to package/MDAnalysis/lib/formats/include/endianswap.h diff --git a/package/MDAnalysis/coordinates/include/fastio.h b/package/MDAnalysis/lib/formats/include/fastio.h similarity index 92% rename from package/MDAnalysis/coordinates/include/fastio.h rename to package/MDAnalysis/lib/formats/include/fastio.h index 998b5ebd5fb..4eeda73e21f 100644 --- a/package/MDAnalysis/coordinates/include/fastio.h +++ b/package/MDAnalysis/lib/formats/include/fastio.h @@ -22,10 +22,13 @@ * ***************************************************************************/ +#ifndef FASTIO_H +#define FASTIO_H + /* Compiling on windows */ #if defined(_MSC_VER) -#if 1 +#if 1 /* use native Windows I/O calls */ #define FASTIO_NATIVEWIN32 1 @@ -51,10 +54,10 @@ typedef struct { static int fio_win32convertfilename(const char *filename, char *newfilename, int maxlen) { int i; int len=strlen(filename); - + if ((len + 1) >= maxlen) return -1; - + for (i=0; i @@ -352,10 +355,10 @@ static int fio_open(const char *filename, int mode, fio_fd *fd) { int nfd; int oflag = 0; - if (mode == FIO_READ) + if (mode == FIO_READ) oflag = O_RDONLY; - if (mode == FIO_WRITE) + if (mode == FIO_WRITE) oflag = O_WRONLY | O_CREAT | O_TRUNC; nfd = open(filename, oflag, 0666); @@ -371,10 +374,10 @@ static int fio_fclose(fio_fd fd) { return close(fd); } -static fio_size_t fio_fread(void *ptr, fio_size_t size, +static fio_size_t fio_fread(void *ptr, fio_size_t size, fio_size_t nitems, fio_fd fd) { int i; - fio_size_t len = 0; + fio_size_t len = 0; int cnt = 0; for (i=0; i= 0) return 0; /* success (emulate behavior of fseek) */ - else + else return -1; /* failure (emulate behavior of fseek) */ } @@ -452,3 +455,4 @@ static int fio_write_str(fio_fd fd, const char *str) { return (fio_fwrite((void *) str, len, 1, fd) != 1); } +#endif // FASTIO_H diff --git a/package/MDAnalysis/coordinates/include/readdcd.h b/package/MDAnalysis/lib/formats/include/readdcd.h similarity index 97% rename from package/MDAnalysis/coordinates/include/readdcd.h rename to package/MDAnalysis/lib/formats/include/readdcd.h index c9d0f9999c9..b248097f4db 100644 --- a/package/MDAnalysis/coordinates/include/readdcd.h +++ b/package/MDAnalysis/lib/formats/include/readdcd.h @@ -70,7 +70,7 @@ static int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, int * * unitcell holds unit cell data if present. */ static int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, - float *unitcell, int nfixed, int first, int *freeind, + double *unitcell, int nfixed, int first, int *freeind, float *fixedcoords, int reverse, int charmm); /* @@ -89,7 +89,7 @@ static int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, * unitcell holds unit cell data if present. */ static int read_dcdsubset(fio_fd fd, int natoms, int lowerb, int upperb, float *x, float *y, float *z, - float *unitcell, int nfixed, int first, int *freeind, + double *unitcell, int nfixed, int first, int *freeind, float *fixedcoords, int reverse, int charmm); /* @@ -136,7 +136,7 @@ static int write_dcdheader(fio_fd fd, const char *remarks, int natoms, * Output: 0 on success, negative error code on failure. * Side effects: coordinates are written to the dcd file. */ -static int write_dcdstep(fio_fd fd, int curstep, int curframe, +static int write_dcdstep(fio_fd fd, int curframe, int curstep, int natoms, const float *x, const float *y, const float *z, const double *unitcell, int charmm); @@ -357,7 +357,7 @@ static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART, } static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian, - float *unitcell) { + double *unitcell) { int i, input_integer; if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) { @@ -370,7 +370,7 @@ static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian, if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD; if (reverseEndian) swap8_aligned(tmp, 6); - for (i=0; i<6; i++) unitcell[i] = (float)tmp[i]; + for (i=0; i<6; i++) unitcell[i] = tmp[i]; } else { /* unrecognized block, just skip it */ if (fio_fseek(fd, input_integer, FIO_SEEK_CUR)) return DCD_BADREAD; @@ -426,7 +426,7 @@ static int read_charmm_4dim(fio_fd fd, int charmm, int reverseEndian) { /* XXX This is completely broken for fixed coordinates */ static int read_dcdsubset(fio_fd fd, int N, int lowerb, int upperb, float *X, float *Y, float *Z, - float *unitcell, int num_fixed, int first, int *indexes, float *fixedcoords, + double *unitcell, int num_fixed, int first, int *indexes, float *fixedcoords, int reverseEndian, int charmm) { //int ret_val; /* Return value from read */ fio_size_t seekpos; @@ -496,11 +496,12 @@ static int read_dcdsubset(fio_fd fd, int N, int lowerb, int upperb, float *X, fl } else { return DCD_BADFORMAT; } + return DCD_SUCCESS; } static int read_dcdstep(fio_fd fd, int N, float *X, float *Y, float *Z, - float *unitcell, int num_fixed, + double *unitcell, int num_fixed, int first, int *indexes, float *fixedcoords, int reverseEndian, int charmm) { int ret_val; /* Return value from read */ @@ -692,10 +693,9 @@ static int write_dcdheader(fio_fd fd, const char *remarks, int N, int charmm) { int out_integer; float out_float; - char title_string[200]; + char title_string[241]; time_t cur_time; struct tm *tmbuf; - char time_str[81]; out_integer = 84; WRITE(fd, (char *) & out_integer, sizeof(int)); @@ -736,16 +736,10 @@ static int write_dcdheader(fio_fd fd, const char *remarks, int N, } fio_write_int32(fd, 84); fio_write_int32(fd, 164); - fio_write_int32(fd, 2); - - strncpy(title_string, remarks, 80); - title_string[79] = '\0'; - WRITE(fd, title_string, 80); + fio_write_int32(fd, 3); /* the number of 80 character title strings */ - cur_time=time(NULL); - tmbuf=localtime(&cur_time); - strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf); - WRITE(fd, time_str, 80); + strncpy(title_string, remarks, 240); + WRITE(fd, title_string, 240); fio_write_int32(fd, 164); fio_write_int32(fd, 4); diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx new file mode 100644 index 00000000000..d080d8df6b6 --- /dev/null +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -0,0 +1,716 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- http://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. +# +# 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 +# +"""\ +Low level DCD trajectory reading - :mod:`MDAnalysis.lib.formats.libdcd` +------------------------------------------------------------------------ + +:mod:`libdcd` contains the class :class:`DCDFile` to read and write frames of a +DCD file. The class tries to behave similar to a normal file object. + +:mod:`libdcd` contains the classes :class:`DCDFile`, which can be used to read +and write frames from and to DCD files. These classes are used internally by +MDAnalysis in :mod:`MDAnalysis.coordinates.DCD`. They behave similar to normal +file objects. + +For example, one can use a :class:`DCDFile` to directly calculate mean +coordinates (where the coordinates are stored in `x` attribute of the +:class:`namedtuple` `frame`): + +.. code-block:: python + :emphasize-lines: 1,2,5 + + with DCDFile("trajectory.dcd") as dcd: + header = dcd.header + mean = np.zeros((header['natoms'], 3)) + # iterate over trajectory + for frame in dcd: + mean += frame.x + mean /= header['natoms'] + + +Besides iteration one can also seek to arbitrary frames using the +:meth:`~DCDFile.seek` method. Note that instead of seeking to a byte-offset as +for normal Python streams, the seek and tell method of DCDFile operate on +complete trajectory frames. + +.. rubric:: Acknowledgements + +:mod:`libdcd` contains and is originally based on DCD reading and writing code +from VMD's `molfile`_ plugin and `catdcd`_. + +.. _molfile: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/ +.. _catdcd: http://www.ks.uiuc.edu/Development/MDTools/catdcd/ + +""" +from six.moves import range + + +from os import path +import numpy as np +from collections import namedtuple +import six +import string +import sys + +cimport numpy as np +from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END +from libc.stdint cimport uintptr_t +from libc.stdlib cimport free + +_whence_vals = {"FIO_SEEK_SET": SEEK_SET, + "FIO_SEEK_CUR": SEEK_CUR, + "FIO_SEEK_END": SEEK_END} + +# Tell cython about the off_t type. It doesn't need to match exactly what is +# defined since we don't expose it to python but the cython compiler needs to +# know about it. +cdef extern from 'sys/types.h': + ctypedef int off_t + +ctypedef int fio_fd; +ctypedef off_t fio_size_t + +ctypedef np.float32_t FLOAT_T +ctypedef np.float64_t DOUBLE_T +FLOAT = np.float32 +DOUBLE = np.float64 + +cdef enum: + FIO_READ = 0x01 + FIO_WRITE = 0x02 + +DCD_IS_CHARMM = 0x01 +DCD_HAS_4DIMS = 0x02 +DCD_HAS_EXTRA_BLOCK = 0x04 + +DCD_ERRORS = { + 0: 'Success', + -1: 'Normal EOF', + -2: 'DCD file does not exist', + -3: 'Open of DCD file failed', + -4: 'read call on DCD file failed', + -5: 'premature EOF found in DCD file', + -6: 'format of DCD file is wrong', + -7: 'output file already exists', + -8: 'malloc failed' +} + +cdef extern from 'include/fastio.h': + int fio_open(const char *filename, int mode, fio_fd *fd) + int fio_fclose(fio_fd fd) + fio_size_t fio_ftell(fio_fd fd) + fio_size_t fio_fseek(fio_fd fd, fio_size_t offset, int whence) + +cdef extern from 'include/readdcd.h': + int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, + int *nsavc, double *delta, int *nfixed, int **freeind, + float **fixedcoords, int *reverse_endian, int *charmm, + char **remarks, int *len_remarks) + void close_dcd_read(int *freeind, float *fixedcoords) + int read_dcdstep(fio_fd fd, int natoms, float *X, float *Y, float *Z, + double *unitcell, int num_fixed, + int first, int *indexes, float *fixedcoords, + int reverse_endian, int charmm) + int read_dcdsubset(fio_fd fd, int natoms, int lowerb, int upperb, + float *X, float *Y, float *Z, + double *unitcell, int num_fixed, + int first, int *indexes, float *fixedcoords, + int reverse_endian, int charmm) + int write_dcdheader(fio_fd fd, const char *remarks, int natoms, + int istart, int nsavc, double delta, int with_unitcell, + int charmm); + int write_dcdstep(fio_fd fd, int curframe, int curstep, + int natoms, const float *x, const float *y, const float *z, + const double *unitcell, int charmm); + +DCDFrame = namedtuple('DCDFrame', 'xyz unitcell') + +cdef class DCDFile: + """DCDFile(fname, mode='r') + + File like wrapper for DCD files + + 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 protocol is supported as well. + + DCDFile can read typical DCD files created by e.g., CHARMM, NAMD, or LAMMPS. It + reads raw data from the trajectory and hence interpretation of, for instance, + different unitcell conventions or time and length units, has to be handled in + higher level code. Reading and writing does not support fixed atoms or 4D + coordinates. + + Parameters + ---------- + fname : str + The filename to open. + mode : ('r', 'w') + The mode in which to open the file, either 'r' read or 'w' write + + Examples + -------- + >>> from MDAnalysis.lib.formats.libdcd import DCDFile + >>> with DCDFile('foo.dcd') as f: + >>> for frame in f: + >>> print(frame.x) + + + Notes + ----- + DCD is not a well defined format. One consequence of this is that different + programs like CHARMM and NAMD are using different convention to store the + unitcell information. The :class:`DCDFile` will read the unitcell information + as is when available. Post processing depending on the program this DCD file + was written with is necessary. Have a look at the MDAnalysis DCD reader for + possible post processing into a common unitcell data structure. You can also + find more information how different programs store unitcell information in DCD + on the `mdawiki`_ . + + .. _mdawiki: https://github.com/MDAnalysis/mdanalysis/wiki/FileFormats#dcd + """ + cdef fio_fd fp + cdef readonly fname + cdef int istart + cdef int nsavc + cdef double delta + cdef int natoms + cdef int nfixed + cdef int *freeind + cdef float *fixedcoords + cdef int reverse_endian + cdef int charmm + cdef readonly is_periodic + cdef remarks + cdef str mode + cdef readonly int ndims + cdef readonly int n_frames + cdef bint b_read_header + cdef int current_frame + cdef readonly int _firstframesize + cdef readonly int _framesize + cdef readonly int _header_size + cdef int is_open + cdef int reached_eof + cdef int wrote_header + + def __cinit__(self, fname, mode='r'): + self.fname = fname.encode('utf-8') + self.natoms = 0 + self.is_open = False + self.wrote_header = False + self.open(mode) + + def __dealloc__(self): + self.close() + + def __enter__(self): + """Support context manager""" + if not self.is_open: + self.open(self.mode) + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + """Support context manager""" + self.close() + # always propagate exceptions forward + return False + + def __iter__(self): + self.close() + self.open(self.mode) + return self + + def __next__(self): + if self.reached_eof: + raise StopIteration + return self.read() + + def __len__(self): + if not self.is_open: + raise IOError('No file currently opened') + return self.n_frames + + def tell(self): + """ + Returns + ------- + current frame (0-based) + """ + return self.current_frame + + def open(self, mode='r'): + """open(mode='r') + + Open a DCD file + + If another DCD file is currently opened it will be closed + + Parameters + ---------- + mode : ('r', 'w') + The mode in which to open the file, either 'r' read or 'w' write + + """ + if self.is_open: + self.close() + + if mode == 'r': + fio_mode = FIO_READ + elif mode == 'w': + fio_mode = FIO_WRITE + else: + raise IOError("unkown mode '{}', use either r or w".format(mode)) + self.mode = str(mode) + + ok = fio_open(self.fname, fio_mode, &self.fp) + if ok != 0: + raise IOError("couldn't open file: {}\n" + "ErrorCode: {}".format(self.fname, DCD_ERRORS[ok])) + self.is_open = True + self.current_frame = 0 + self.reached_eof = False + self.wrote_header = False + # Has to come last since it checks the reached_eof flag + if self.mode == 'r': + self._read_header() + + def close(self): + """Close the open DCD file + + """ + if self.is_open: + # In case there are fixed atoms we should free the memory again. + # Both pointers are guaranted to be non NULL if either one is. + if self.freeind != NULL: + close_dcd_read(self.freeind, self.fixedcoords); + + ok = fio_fclose(self.fp) + + self.is_open = False + if ok != 0: + raise IOError("couldn't close file: {}\n" + "ErrorCode: {}".format(self.fname, DCD_ERRORS[ok])) + + + cdef void _read_header(self): + """read header and populate internal fields""" + if not self.is_open: + raise IOError("No file open") + + cdef char* c_remarks + cdef int len_remarks = 0 + cdef int nsets + + ok = read_dcdheader(self.fp, &self.natoms, &nsets, &self.istart, + &self.nsavc, &self.delta, &self.nfixed, &self.freeind, + &self.fixedcoords, &self.reverse_endian, + &self.charmm, &c_remarks, &len_remarks) + if ok != 0: + raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + + self.is_periodic = bool((self.charmm & DCD_IS_CHARMM) and + (self.charmm & DCD_HAS_EXTRA_BLOCK)) + + if c_remarks != NULL: + py_remarks = c_remarks[:len_remarks] + free(c_remarks) + else: + py_remarks = "" + self.ndims = 3 if not self.charmm & DCD_HAS_4DIMS else 4 + # This function assumes that the dcd header was already read and + # self.ndims is set. It will only work when called here !!! + self.n_frames = self._estimate_n_frames() + self.b_read_header = True + + # make sure fixed atoms have been read + try: + self.read() + self.seek(0) + except IOError: + if self.n_frames != 0: + raise IOError("DCD is corrupted") + + if sys.version_info[0] < 3: + py_remarks = unicode(py_remarks, 'ascii', "ignore") + py_remarks = str(py_remarks.encode('ascii', 'ignore')) + else: + if isinstance(py_remarks, bytes): + py_remarks = py_remarks.decode('ascii', 'ignore') + + py_remarks = "".join(s for s in py_remarks if s in string.printable) + + self.remarks = py_remarks + + cdef int _estimate_n_frames(self): + """ Only call this function in _read_header!!! + """ + extrablocksize = 48 + 8 if self.charmm & DCD_HAS_EXTRA_BLOCK else 0 + self._firstframesize = (self.natoms + 2) * self.ndims * sizeof(float) + extrablocksize + self._framesize = ((self.natoms - self.nfixed + 2) * self.ndims * sizeof(float) + + extrablocksize) + filesize = path.getsize(self.fname) + # It's safe to use ftell, even though ftell returns a long, because the + # header size is < 4GB. + self._header_size = fio_ftell(self.fp) + nframessize = filesize - self._header_size - self._firstframesize + return nframessize / self._framesize + 1 + + def seek(self, frame): + """seek(frame) + + Seek to Frame. + + Parameters + ---------- + frame : int + seek the file to given frame (0-based) + + """ + if frame >= self.n_frames: + raise EOFError('Trying to seek over max number of frames') + self.reached_eof = False + + cdef fio_size_t offset + if frame == 0: + offset = self._header_size + else: + offset = self._header_size + self._firstframesize + self._framesize * (frame - 1) + + ok = fio_fseek(self.fp, offset, _whence_vals["FIO_SEEK_SET"]) + if ok != 0: + raise IOError("DCD seek failed with DCD error={}".format(DCD_ERRORS[ok])) + self.current_frame = frame + + @property + def header(self): + """ + Returns + ------- + dict of header values needed to write new dcd. + natoms: number of atoms + istart: starting frame number + nsavc: number of frames between saves + delta: integrator time step. + charm: bitfield integer if file contains special CHARMM information + remarks: remark string, max 240 bytes. + + + """ + return {'natoms': self.natoms, + 'istart': self.istart, + 'nsavc': self.nsavc, + 'delta': self.delta, + 'is_periodic': self.is_periodic, + 'remarks': self.remarks} + + @property + def charmm_bitfield(self): + """This DCDFile reader can process files written by different MD simulation + programs. For files produced by CHARMM or other programs that follow + the same convention we are reading a special CHARMM bitfield that + stores different flags about additional information that is stored in + the dcd. The bit flags are: + + .. code:: + + DCD_IS_CHARMM = 0x01 + DCD_HAS_4DIMS = 0x02 + DCD_HAS_EXTRA_BLOCK = 0x04 + + Here `DCD_HAS_EXTRA_BLOCK` means that unitcell information is stored. + + """ + return self.charmm + + def write_header(self, remarks, natoms, istart, nsavc, delta, is_periodic): + """write_header(remarks, natoms, istart, nsavc, delta, is_periodic) + Write DCD header + + This function needs to be called before the first frame can be written. + + Parameters + ---------- + remarks : str + remarks of DCD file. Shouldn't be more then 240 characters (ASCII) + natoms : int + number of atoms to write + istart : int + starting frame number + nsavc : int + number of frames between saves + delta : float + integrator time step. The time for 1 frame is nsavc * delta + is_periodic : bool + write unitcell information. Also pretends that file was written by CHARMM 24 + """ + if not self.is_open: + raise IOError("No file open") + if not self.mode=='w': + raise IOError("Incorrect file mode for writing.") + if self.wrote_header: + raise IOError("Header already written") + + cdef int with_unitcell = is_periodic + if is_periodic: + self.charmm = DCD_HAS_EXTRA_BLOCK | DCD_IS_CHARMM + self.natoms = natoms + + if isinstance(remarks, six.string_types): + try: + remarks = bytearray(remarks, 'ascii') + except UnicodeDecodeError: + remarks = bytearray(remarks) + + ok = write_dcdheader(self.fp, remarks, self.natoms, istart, + nsavc, delta, with_unitcell, + self.charmm) + if ok != 0: + raise IOError("Writing DCD header failed: {}".format(DCD_ERRORS[ok])) + self.wrote_header = True + + def write(self, xyz, box=None): + """write(xyz, box=None) + write one frame into DCD file. + + Parameters + ---------- + xyz : array_like, shape=(natoms, 3) + cartesion coordinates + box : array_like, shape=(6) (optional) + Box vectors for this frame. Can be left to skip writing a unitcell + + """ + if not self.is_open: + raise IOError("No file open") + if self.mode != 'w': + raise IOError('File opened in mode: {}. Writing only allowed ' + 'in mode "w"'.format('self.mode')) + if (self.charmm & DCD_HAS_EXTRA_BLOCK): + if len(box) != 6: + raise ValueError("box size is wrong should be 6, got: {}".format(box.size)) + else: + # use a dummy box. It won't be written anyway in readdcd. + box = np.zeros(6) + + if not self.wrote_header: + raise IOError("write header first before frames can be written") + xyz = np.asarray(xyz, order='F', dtype=FLOAT) + if xyz.shape != (self.natoms, 3): + raise ValueError("xyz shape is wrong should be (natoms, 3), got:".format(xyz.shape)) + + cdef DOUBLE_T[::1] c_box = np.asarray(box, order='C', dtype=DOUBLE) + cdef FLOAT_T[::1] x = xyz[:, 0] + cdef FLOAT_T[::1] y = xyz[:, 1] + cdef FLOAT_T[::1] z = xyz[:, 2] + + step = self.istart + self.current_frame * self.nsavc + ok = write_dcdstep(self.fp, self.current_frame + 1, step, + self.natoms, &x[0], + &y[0], &z[0], + &c_box[0], self.charmm) + if ok != 0: + raise IOError("Couldn't write DCD frame: reason {}".format(DCD_ERRORS[ok])) + + self.current_frame += 1 + + def read(self): + """ + Read next dcd frame + + Returns + ------- + DCDFrame : namedtuple + positions are in ``x`` and unitcell in ``unitcell`` attribute of DCDFrame + + Notes + ----- + unitcell is read as is from DCD. Post processing depending on the program this + DCD file was written with is necessary. Have a look at the MDAnalysis DCD reader + for possible post processing into a common unitcell data structure. + + """ + if self.reached_eof: + raise EOFError('Reached last frame in DCD, seek to 0') + if not self.is_open: + raise IOError("No file open") + if self.mode != 'r': + raise IOError('File opened in mode: {}. Reading only allow ' + 'in mode "r"'.format('self.mode')) + if self.n_frames == 0: + raise IOError("opened empty file. No frames are saved") + + cdef np.ndarray xyz = np.empty((self.natoms, self.ndims), dtype=FLOAT, order='F') + cdef np.ndarray unitcell = np.empty(6, dtype=DOUBLE) + unitcell[0] = unitcell[2] = unitcell[5] = 0.0; + unitcell[4] = unitcell[3] = unitcell[1] = 90.0; + + first_frame = self.current_frame == 0 + ok = self.c_readframes_helper(xyz[:, 0], xyz[:, 1], xyz[:, 2], unitcell, first_frame) + if ok != 0 and ok != -4: + raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + + # we couldn't read any more frames. + if ok == -4: + self.reached_eof = True + raise StopIteration + + self.current_frame += 1 + return DCDFrame(xyz, unitcell) + + + def readframes(self, start=None, stop=None, step=None, order='fac', indices=None): + """readframes(start=None, stop=None, step=None, order='fac', indices=None) + read multiple frames at once + + Parameters + ---------- + start : int (optional) + starting frame, default to 0 + stop : int (optional) + stop frame, default to ``n_frames`` + step : int (optional) + step between frames read, defaults to 1 + order : str (optional) + give order of returned array with `f`:frames, `a`:atoms, `c`:coordinates + indices : array_like (optional) + only read selected atoms. In ``None`` read all. + + Returns + ------- + DCDFrame : namedtuple + positions are in ``x`` and unitcell in ``unitcell`` attribute of DCDFrame. + Here the attributes contain the positions for all frames in the given order + + Notes + ----- + unitcell is read as it from DCD. Post processing depending the program + this DCD file was written with is necessary. Have a look at the + MDAnalysis DCD reader for possible post processing into a common + unitcell data structure. + + """ + if self.reached_eof: + raise EOFError('Reached last frame in DCD, seek to 0') + if not self.is_open: + raise IOError("No file open") + if self.mode != 'r': + raise IOError('File opened in mode: {}. Reading only allow ' + 'in mode "r"'.format('self.mode')) + if self.n_frames == 0: + raise IOError("opened empty file. No frames are saved") + + self.seek(0) + # if we only want to iterate backwards flip start and end + if start is None and stop is None and step is not None and step < 0: + stop = -1 + start = self.n_frames - 1 + stop = stop if not stop is None else self.n_frames + start = start if not start is None else 0 + step = step if not step is None else 1 + + + cdef int n + n = len(range(start, stop, step)) + + cdef np.ndarray[np.int64_t, ndim=1] c_indices + if indices is None: + c_indices = np.arange(self.natoms) + natoms = self.natoms + else: + natoms = len(indices) + c_indices = np.asarray(indices, dtype=np.int64) + + cdef int hash_order = -1 + if order == 'fac': + shape = (n, natoms, self.ndims) + hash_order = 1 + elif order == 'fca': + shape = (n, self.ndims, natoms) + hash_order = 2 + elif order == 'afc': + shape = (natoms, n, self.ndims) + hash_order = 3 + elif order == 'acf': + shape = (natoms, self.ndims, n) + hash_order = 4 + elif order == 'caf': + shape = (self.ndims, natoms, n) + hash_order = 5 + elif order == 'cfa': + hash_order = 6 + shape = (self.ndims, n, natoms) + else: + raise ValueError("unkown order '{}'".format(order)) + + cdef np.ndarray[FLOAT_T, ndim=3] xyz = np.empty(shape, dtype=FLOAT) + cdef np.ndarray[DOUBLE_T, ndim=2] box = np.empty((n, 6)) + + + cdef np.ndarray xyz_tmp = np.empty((self.natoms, self.ndims), dtype=FLOAT, order='F') + cdef int ok, i + + if start == 0 and step == 1 and stop == self.n_frames: + for i in range(n): + ok = self.c_readframes_helper(xyz_tmp[:, 0], xyz_tmp[:, 1], xyz_tmp[:, 2], box[i], i==0) + if ok != 0 and ok != -4: + raise IOError("Reading DCD frames failed: {}".format(DCD_ERRORS[ok])) + copy_in_order(xyz_tmp[c_indices], xyz, hash_order, i) + else: + counter = 0 + for i in range(start, stop, step): + self.seek(i) + ok = self.c_readframes_helper(xyz_tmp[:, 0], xyz_tmp[:, 1], xyz_tmp[:, 2], box[counter], i==0) + if ok != 0 and ok != -4: + raise IOError("Reading DCD frames failed: {}".format(DCD_ERRORS[ok])) + copy_in_order(xyz_tmp[c_indices], xyz, hash_order, counter) + counter += 1 + + return DCDFrame(xyz, box) + + # Helper to read current DCD frame + cdef int c_readframes_helper(self, FLOAT_T[::1] x, + FLOAT_T[::1] y, FLOAT_T[::1] z, + DOUBLE_T[::1] unitcell, int first_frame): + cdef int ok + ok = read_dcdstep(self.fp, self.natoms, + &x[0], + &y[0], &z[0], + &unitcell[0], self.nfixed, first_frame, + self.freeind, self.fixedcoords, + self.reverse_endian, self.charmm) + return ok + + +# Helper in readframes to copy given a specific memory layout +cdef void copy_in_order(FLOAT_T[:, :] source, FLOAT_T[:, :, :] target, int order, int index): + if order == 1: # 'fac': + target[index] = source + elif order == 2: # 'fca': + target[index] = source.T + elif order == 3: # 'afc': + target[:, index] = source + elif order == 4: # 'acf': + target[:, :, index] = source + elif order == 5: # 'caf': + target[:, :, index] = source.T + elif order == 6: # 'cfa': + target[:, index] = source.T diff --git a/package/doc/sphinx/source/documentation_pages/lib/formats/libdcd.rst b/package/doc/sphinx/source/documentation_pages/lib/formats/libdcd.rst new file mode 100644 index 00000000000..aaf528c1d63 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/lib/formats/libdcd.rst @@ -0,0 +1,3 @@ +.. automodule:: MDAnalysis.lib.formats.libdcd + :members: + :inherited-members: diff --git a/package/doc/sphinx/source/documentation_pages/lib_modules.rst b/package/doc/sphinx/source/documentation_pages/lib_modules.rst index 97368315bc3..5bbe3d041ed 100644 --- a/package/doc/sphinx/source/documentation_pages/lib_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/lib_modules.rst @@ -65,3 +65,4 @@ Python-based projects. :maxdepth: 1 ./lib/formats/libmdaxdr + ./lib/formats/libdcd diff --git a/package/setup.py b/package/setup.py index 8f9226f030b..54c10977236 100755 --- a/package/setup.py +++ b/package/setup.py @@ -295,11 +295,11 @@ def extensions(config): include_dirs = [get_numpy_include] - dcd = MDAExtension('coordinates._dcdmodule', - ['MDAnalysis/coordinates/src/dcd.c'], - include_dirs=include_dirs + ['MDAnalysis/coordinates/include'], - define_macros=define_macros, - extra_compile_args=extra_compile_args) + libdcd = MDAExtension('lib.formats.libdcd', + ['MDAnalysis/lib/formats/libdcd' + source_suffix], + include_dirs=include_dirs + ['MDAnalysis/lib/formats/include'], + define_macros=define_macros, + extra_compile_args=extra_compile_args) distances = MDAExtension('lib.c_distances', ['MDAnalysis/lib/c_distances' + source_suffix], include_dirs=include_dirs + ['MDAnalysis/lib/include'], @@ -352,7 +352,7 @@ def extensions(config): include_dirs = include_dirs+['MDAnalysis/analysis/encore/dimensionality_reduction/include'], libraries=["m"], extra_compile_args=["-O3", "-ffast-math","-std=c99"]) - pre_exts = [dcd, distances, distances_omp, qcprot, + pre_exts = [libdcd, distances, distances_omp, qcprot, transformation, libmdaxdr, util, encore_utils, ap_clustering, spe_dimred] @@ -524,8 +524,11 @@ def dynamic_author_list(): }, test_suite="MDAnalysisTests", tests_require=[ - 'pytest>=3.1.2', 'nose>=1.3.7', + 'pytest>=3.1.2', + 'pytest-cov', + 'pytest-xdist', + 'hypothesis', 'MDAnalysisTests=={0}'.format(RELEASE), # same as this release! ], zip_safe=False, # as a zipped egg the *.so files are not found (at diff --git a/testsuite/MDAnalysisTests/coordinates/base.py b/testsuite/MDAnalysisTests/coordinates/base.py index 9f6838b577f..45bea93541a 100644 --- a/testsuite/MDAnalysisTests/coordinates/base.py +++ b/testsuite/MDAnalysisTests/coordinates/base.py @@ -243,13 +243,13 @@ def test_get_writer_2(self): assert_equal(W.n_atoms, 100) def test_dt(self): - assert_equal(self.reader.dt, self.ref.dt) + assert_almost_equal(self.reader.dt, self.ref.dt) def test_ts_dt_matches_reader(self): assert_equal(self.reader.ts.dt, self.reader.dt) def test_total_time(self): - assert_equal(self.reader.totaltime, self.ref.totaltime) + assert_almost_equal(self.reader.totaltime, self.ref.totaltime, decimal=5) def test_first_dimensions(self): self.reader.rewind() @@ -1068,7 +1068,7 @@ def assert_timestep_almost_equal(A, B, decimal=6, verbose=True): 'A.frame = {}, B.frame={}'.format( A.frame, B.frame)) - if A.time != B.time: + if not np.allclose(A.time, B.time): raise AssertionError('A and B refer to different times: ' 'A.time = {}, B.time={}'.format( A.time, B.time)) diff --git a/testsuite/MDAnalysisTests/coordinates/reference.py b/testsuite/MDAnalysisTests/coordinates/reference.py index 7c71d0f08bb..40d9e0084f4 100644 --- a/testsuite/MDAnalysisTests/coordinates/reference.py +++ b/testsuite/MDAnalysisTests/coordinates/reference.py @@ -253,34 +253,3 @@ class RefLAMMPSDataMini(object): vel_atom1 = np.array([-0.005667593, 0.00791380978, -0.00300779533], dtype=np.float32) dimensions = np.array([60., 50., 30., 90., 90., 90.], dtype=np.float32) - - -class RefCHARMMtriclinicDCD(object): - topology = PSF_TRICLINIC - trajectory = DCD_TRICLINIC - # time(ps) A B C alpha beta gamma (length in Angstrome, angles in degrees) - # dcd starts at t = 1ps - ref_dimensions = np.array([ - [1., 35.44604, 35.06156, 34.1585, 91.32802, 61.73521, 44.40703], - [2., 34.65957, 34.22689, 33.09897, 90.56206, 61.79192, 44.14549], - [3., 34.52772, 34.66422, 33.53881, 90.55859, 63.11228, 40.14044], - [4., 34.43749, 33.38432, 34.02133, 88.82457, 64.98057, 36.77397], - [5., 33.73129, 32.47752, 34.18961, 89.88102, 65.89032, 36.10921], - [6., 33.78703, 31.90317, 34.98833, 90.03092, 66.12877, 35.07141], - [7., 33.24708, 31.18271, 34.9654, 93.11122, 68.17743, 35.73643], - [8., 32.92599, 30.31393, 34.99197, 93.89051, 69.3799, 33.48945], - [9., 32.15295, 30.43056, 34.96157, 96.01416, 71.50115, 32.56111], - [10., 31.99748, 30.21518, 35.24292, 95.85821, 71.08429, 31.85939] - ]) - - -class RefNAMDtriclinicDCD(object): - topology = PSF_NAMD_TRICLINIC - trajectory = DCD_NAMD_TRICLINIC - # vmd topology trajectory - # molinfo 0 get {a b c alpha beta gamma} - # time(ps) A B C alpha beta gamma (length in Angstrome, angles in degrees) - ref_dimensions = np.array([ - [1., 38.426594, 38.393101, 44.759800, 90.000000, 90.000000, 60.028915], - ]) - diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 3c300985d07..550fa5a1e64 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -19,35 +19,52 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -from __future__ import absolute_import -import MDAnalysis as mda -import numpy as np -import os +from __future__ import absolute_import, print_function, division from six.moves import zip, range +import numpy as np + +import MDAnalysis as mda +from MDAnalysis.coordinates.DCD import DCDReader +from MDAnalysis.exceptions import NoDataError -from nose.plugins.attrib import attr from numpy.testing import (assert_equal, assert_array_equal, assert_raises, - assert_almost_equal, assert_array_almost_equal, - assert_allclose, dec) -from unittest import TestCase + assert_almost_equal, assert_array_almost_equal) + +from MDAnalysisTests.datafiles import (DCD, PSF, DCD_empty, PRMncdf, NCDF, + COORDINATES_TOPOLOGY, COORDINATES_DCD, + PSF_TRICLINIC, DCD_TRICLINIC, + PSF_NAMD_TRICLINIC, DCD_NAMD_TRICLINIC) +from MDAnalysisTests.coordinates.base import (MultiframeReaderTest, + BaseReference, + BaseWriterTest) + +import pytest -from MDAnalysisTests.datafiles import (DCD, PSF, DCD_empty, CRD, PRMncdf, NCDF) -from MDAnalysisTests.coordinates.reference import (RefCHARMMtriclinicDCD, - RefNAMDtriclinicDCD) -from MDAnalysisTests.coordinates.base import BaseTimestepTest -from MDAnalysisTests import module_not_found, tempdir -from MDAnalysisTests.plugins.knownfailure import knownfailure -@attr('issue') -def TestDCD_Issue32(): - """Test for Issue 32: 0-size dcds lead to a segfault: now caught with - IOError""" - assert_raises(IOError, mda.Universe, PSF, DCD_empty) +class DCDReference(BaseReference): + def __init__(self): + super(DCDReference, self).__init__() + self.trajectory = COORDINATES_DCD + self.topology = COORDINATES_TOPOLOGY + self.reader = mda.coordinates.DCD.DCDReader + self.writer = mda.coordinates.DCD.DCDWriter + self.ext = 'xtc' + self.prec = 3 + self.changing_dimensions = True -class TestDCDReaderClass(TestCase): - def test_with_statement(self): - from MDAnalysis.coordinates.DCD import DCDReader +class TestDCDReader(MultiframeReaderTest): + def __init__(self, reference=None): + if reference is None: + reference = DCDReference() + super(TestDCDReader, self).__init__(reference) + + @staticmethod + def test_empty_dcd(): + assert_raises(IOError, mda.Universe, PSF, DCD_empty) + + @staticmethod + def test_with_statement(): try: with DCDReader(DCD) as trj: @@ -65,420 +82,284 @@ def test_with_statement(self): err_msg="with_statement: DCDReader does not read all frames") -class TestDCDReader(TestCase): - def setUp(self): - self.universe = mda.Universe(PSF, DCD) - self.dcd = self.universe.trajectory - self.ts = self.universe.coord - - def tearDown(self): - del self.universe - del self.dcd - del self.ts - - def test_rewind_dcd(self): - self.dcd.rewind() - assert_equal(self.ts.frame, 0, "rewinding to frame 0") - - def test_next_dcd(self): - self.dcd.rewind() - self.dcd.next() - assert_equal(self.ts.frame, 1, "loading frame 1") - - def test_jump_dcd(self): - self.dcd[15] # index is 0-based and frames are 0-based - assert_equal(self.ts.frame, 15, "jumping to frame 15") - - def test_jump_lastframe_dcd(self): - self.dcd[-1] - assert_equal(self.ts.frame, 97, "indexing last frame with dcd[-1]") - - def test_slice_dcd(self): - frames = [ts.frame for ts in self.dcd[5:17:3]] - assert_equal(frames, [5, 8, 11, 14], "slicing dcd [5:17:3]") - - def test_list_trajectory(self): - frames = [ts.frame for ts in self.dcd[[0, 3, 4, 5]]] - assert_equal(frames, [0, 3, 4, 5]) - - def test_array_trajectory(self): - frames = [ts.frame for ts in self.dcd[np.array([0, 3, 4, 5])]] - assert_equal(frames, [0, 3, 4, 5]) - - def test_list_reverse_trajectory(self): - frames = [ts.frame for ts in self.dcd[[0, 4, 2, 3, 0, 1]]] - assert_equal(frames, [0, 4, 2, 3, 0, 1]) - - def test_list_repeated_trajectory(self): - frames = [ts.frame for ts in self.dcd[[0, 0, 1, 1, 2, 1, 1]]] - assert_equal(frames, [0, 0, 1, 1, 2, 1, 1]) - - def test_reverse_dcd(self): - frames = [ts.frame for ts in self.dcd[20:5:-1]] - assert_equal(frames, list(range(20, 5, -1)), - "reversing dcd [20:5:-1]") - - def test_n_atoms(self): - assert_equal(self.universe.trajectory.n_atoms, 3341, - "wrong number of atoms") - - def test_n_frames(self): - assert_equal(self.universe.trajectory.n_frames, 98, - "wrong number of frames in dcd") - - def test_dt(self): - assert_almost_equal(self.universe.trajectory.dt, - 1.0, - 4, - err_msg="wrong timestep dt") - - def test_totaltime(self): - # test_totaltime(): need to reduce precision because dt is only precise - # to ~4 decimals and accumulating the inaccuracy leads to even lower - # precision in the totaltime (consequence of fixing Issue 64) - assert_almost_equal(self.universe.trajectory.totaltime, - 97.0, - 3, - err_msg="wrong total length of AdK trajectory") - - def test_frame(self): - self.dcd[15] # index is 0-based and frames are 0-based - assert_equal(self.universe.trajectory.frame, 15, "wrong frame number") - - def test_time(self): - self.dcd[15] # index is 0-based and frames are 0-based - assert_almost_equal(self.universe.trajectory.time, - 15.0, - 5, - err_msg="wrong time of frame") - - def test_volume(self): - assert_almost_equal(self.ts.volume, 0.0, 3, - err_msg="wrong volume for unitcell (no unitcell " - "in DCD so this should be 0)") - - def test_timeseries_slicing(self): - # check that slicing behaves correctly - # should before issue #914 resolved - x = [(0, 1, 1), (1,1,1), (1, 2, 1), (1, 2, 2), (1, 4, 2), (1, 4, 4), - (0, 5, 5), (3, 5, 1), (None, None, None)] - for start, stop, step in x: - yield self._slice_generation_test, start, stop, step - - def test_backwards_stepping(self): - x = [(4, 0, -1), (5, 0, -2), (5, 0, -4)] - for start, stop, step in x: - yield self._failed_slices_test, start, stop, step - - def _slice_generation_test(self, start, stop, step): - self.u = mda.Universe(PSF, DCD) - ts = self.u.trajectory.timeseries(self.u.atoms) - ts_skip = self.u.trajectory.timeseries(self.u.atoms, start, stop, step) - assert_array_almost_equal(ts[:, start:stop:step,:], ts_skip, 5) - - @knownfailure - def _failed_slices_test(self, start, stop, step): - self.u = mda.Universe(PSF, DCD) - ts = self.u.trajectory.timeseries(self.u.atoms) - ts_skip = self.u.trajectory.timeseries(self.u.atoms, start, stop, step) - assert_array_almost_equal(ts[:, start:stop:step,:], ts_skip, 5) - - -def test_DCDReader_set_dt(dt=100., frame=3): +class TestDCDWriter(BaseWriterTest): + def __init__(self, reference=None): + if reference is None: + reference = DCDReference() + super(TestDCDWriter, self).__init__(reference) + + +def test_write_random_unitcell(tmpdir): + testname = str(tmpdir.join('test.dcd')) + testname = 'test.dcd' + rstate = np.random.RandomState(1178083) + random_unitcells = rstate.uniform( + high=80, size=(98, 6)).astype(np.float64) + + u = mda.Universe(PSF, DCD) + with mda.Writer(testname, n_atoms=u.atoms.n_atoms) as w: + for index, ts in enumerate(u.trajectory): + u.atoms.dimensions = random_unitcells[index] + w.write(u.atoms) + + u2 = mda.Universe(PSF, testname) + for index, ts in enumerate(u2.trajectory): + assert_array_almost_equal(u2.trajectory.dimensions, + random_unitcells[index], + decimal=5) + + +################ +# Legacy tests # +################ + +@pytest.fixture(scope='module') +def universe_dcd(): + return mda.Universe(PSF, DCD) + + +def test_rewind(universe_dcd): + universe_dcd.trajectory.rewind() + assert universe_dcd.trajectory.ts.frame == 0 + + +def test_next(universe_dcd): + universe_dcd.trajectory.rewind() + universe_dcd.trajectory.next() + assert universe_dcd.trajectory.ts.frame == 1 + + +def test_jump_last_frame(universe_dcd): + universe_dcd.trajectory[-1] + assert universe_dcd.trajectory.ts.frame == 97 + + +@pytest.mark.parametrize("start, stop, step", ((5, 17, 3), + (20, 5, -1))) +def test_slice(universe_dcd, start, stop, step): + frames = [ts.frame for ts in universe_dcd.trajectory[start:stop:step]] + assert_array_equal(frames, np.arange(start, stop, step)) + + +@pytest.mark.parametrize("array_like", [list, np.array]) +def test_array_like(universe_dcd, array_like): + ar = array_like([0, 3, 4, 5]) + frames = [ts.frame for ts in universe_dcd.trajectory[ar]] + assert_array_equal(frames, ar) + + +@pytest.mark.parametrize("indices", ([0, 4, 2, 3, 0, 1], + [0, 0, 1, 1, 2, 1, 1])) +def test_list_indices(universe_dcd, indices): + frames = [ts.frame for ts in universe_dcd.trajectory[indices]] + assert_array_equal(frames, indices) + + +@pytest.mark.parametrize( + "slice, length", + [([None, None, None], 98), ([0, None, None], 98), ([None, 98, None], 98), + ([None, None, 1], 98), ([None, None, -1], 98), ([2, 6, 2], 2), + ([0, 10, None], 10), ([2, 10, None], 8), ([0, 1, 1], 1), ([1, 1, 1], 0), + ([1, 2, 1], 1), ([1, 2, 2], 1), ([1, 4, 2], 2), ([1, 4, 4], 1), + ([0, 5, 5], 1), ([3, 5, 1], 2), ([4, 0, -1], 4), ([5, 0, -2], 3), + ([5, 0, -4], 2)]) +def test_timeseries_slices(slice, length, universe_dcd): + start, stop, step = slice + allframes = universe_dcd.trajectory.timeseries(format='fac') + xyz = universe_dcd.trajectory.timeseries(start=start, stop=stop, step=step, + format='fac') + assert len(xyz) == length + assert_array_almost_equal(xyz, allframes[start:stop:step]) + + +@pytest.mark.parametrize("order, shape", ( + ('fac', (98, 3341, 3)), + ('fca', (98, 3, 3341)), + ('afc', (3341, 98, 3)), + ('acf', (3341, 3, 98)), + ('caf', (3, 3341, 98)), + ('cfa', (3, 98, 3341)), )) +def test_timeseries_order(order, shape, universe_dcd): + x = universe_dcd.trajectory.timeseries(format=order) + assert x.shape == shape + + +@pytest.mark.parametrize("indices", [[1, 2, 3, 4], [5, 10, 15, 19], + [9, 4, 2, 0, 50]]) +def test_timeseries_atomindices(indices, universe_dcd): + allframes = universe_dcd.trajectory.timeseries(format='afc') + asel = universe_dcd.atoms[indices] + xyz = universe_dcd.trajectory.timeseries(asel=asel, format='afc') + assert len(xyz) == len(indices) + assert_array_almost_equal(xyz, allframes[indices]) + + +def test_timeseries_empty_selection(universe_dcd): + with pytest.raises(NoDataError): + asel = universe_dcd.select_atoms('name FOO') + universe_dcd.trajectory.timeseries(asel=asel) + + +def test_timeseries_skip(universe_dcd): + with pytest.warns(DeprecationWarning): + xyz = universe_dcd.trajectory.timeseries(skip=2, format='fac') + assert len(xyz) == universe_dcd.trajectory.n_frames / 2 + + +def test_reader_set_dt(): + dt = 100 + frame = 3 u = mda.Universe(PSF, DCD, dt=dt) assert_almost_equal(u.trajectory[frame].time, frame*dt, err_msg="setting time step dt={0} failed: " "actually used dt={1}".format( - dt, u.trajectory._ts_kwargs['dt'])) + dt, u.trajectory._ts_kwargs['dt'])) assert_almost_equal(u.trajectory.dt, dt, err_msg="trajectory.dt does not match set dt") -class TestDCDWriter(TestCase): - def setUp(self): - self.universe = mda.Universe(PSF, DCD) - ext = ".dcd" - self.tmpdir = tempdir.TempDir() - self.outfile = self.tmpdir.name + '/dcd-writer' + ext - self.Writer = mda.coordinates.DCD.DCDWriter - def tearDown(self): - try: - os.unlink(self.outfile) - except OSError: - pass - del self.universe - del self.Writer - del self.tmpdir - - @attr('issue') - def test_write_trajectory(self): - """Test writing DCD trajectories (Issue 50)""" - t = self.universe.trajectory - W = self.Writer(self.outfile, t.n_atoms, dt=t.dt, step=t.skip_timestep) - for ts in self.universe.trajectory: +@pytest.mark.parametrize("ext, decimal", (("dcd", 5), + ("xtc", 3))) +def test_writer_dt(tmpdir, ext, decimal): + dt = 5.0 # set time step to 5 ps + universe_dcd = mda.Universe(PSF, DCD, dt=dt) + t = universe_dcd.trajectory + outfile = str(tmpdir.join("test.{}".format(ext))) + with mda.Writer(outfile, n_atoms=t.n_atoms, dt=dt) as W: + for ts in universe_dcd.trajectory: + W.write(universe_dcd.atoms) + + uw = mda.Universe(PSF, outfile) + assert_almost_equal(uw.trajectory.totaltime, + (uw.trajectory.n_frames - 1) * dt, decimal) + times = np.array([uw.trajectory.time for ts in uw.trajectory]) + frames = np.array([ts.frame for ts in uw.trajectory]) + assert_array_almost_equal(times, frames * dt, decimal) + + +@pytest.mark.parametrize("ext, decimal", (("dcd", 5), + ("xtc", 2))) +def test_other_writer(universe_dcd, tmpdir, ext, decimal): + t = universe_dcd.trajectory + outfile = str(tmpdir.join("test.{}".format(ext))) + with t.OtherWriter(outfile) as W: + for ts in universe_dcd.trajectory: W.write_next_timestep(ts) - W.close() - - uw = mda.Universe(PSF, self.outfile) - - # check that the coordinates are identical for each time step - for orig_ts, written_ts in zip(self.universe.trajectory, - uw.trajectory): - assert_array_almost_equal(written_ts._pos, orig_ts._pos, 3, - err_msg="coordinate mismatch between " - "original and written trajectory at " - "frame %d (orig) vs %d (written)" % ( - orig_ts.frame, written_ts.frame)) - - def test_dt(self): - DT = 5.0 - t = self.universe.trajectory - with self.Writer(self.outfile, - t.n_atoms, - dt=DT) as W: # set time step to 5 ps - for ts in self.universe.trajectory: - W.write_next_timestep(ts) - - uw = mda.Universe(PSF, self.outfile) - assert_almost_equal(uw.trajectory.totaltime, - (uw.trajectory.n_frames - 1) * DT, 5) - times = np.array([uw.trajectory.time for ts in uw.trajectory]) - frames = np.array([ts.frame for ts in uw.trajectory]) - assert_array_almost_equal(times, frames * DT, 5) - - def test_OtherWriter(self): - t = self.universe.trajectory - W = t.OtherWriter(self.outfile) - for ts in self.universe.trajectory: - W.write_next_timestep(ts) - W.close() - - uw = mda.Universe(PSF, self.outfile) - - # check that the coordinates are identical for each time step - for orig_ts, written_ts in zip(self.universe.trajectory, - uw.trajectory): - assert_array_almost_equal(written_ts._pos, orig_ts._pos, 3, - err_msg="coordinate mismatch between " - "original and written trajectory at " - "frame %d (orig) vs %d (written)" % ( - orig_ts.frame, written_ts.frame)) - - def test_single_frame(self): - u = mda.Universe(PSF, CRD) - W = mda.Writer(self.outfile, u.atoms.n_atoms) + + uw = mda.Universe(PSF, outfile) + # check that the coordinates are identical for each time step + for orig_ts, written_ts in zip(universe_dcd.trajectory, + uw.trajectory): + assert_array_almost_equal(written_ts.positions, orig_ts.positions, + decimal, + err_msg="coordinate mismatch between " + "original and written trajectory at " + "frame {} (orig) vs {} (written)".format( + orig_ts.frame, written_ts.frame)) + + +def test_single_frame(universe_dcd, tmpdir): + u = universe_dcd + outfile = str(tmpdir.join("test.dcd")) + with mda.Writer(outfile, u.atoms.n_atoms) as W: W.write(u.atoms) - W.close() - w = mda.Universe(PSF, self.outfile) - assert_equal(w.trajectory.n_frames, 1, - "single frame trajectory has wrong number of frames") - assert_almost_equal(w.atoms.positions, - u.atoms.positions, - 3, - err_msg="coordinates do not match") - - def test_with_statement(self): - u = mda.Universe(PSF, CRD) - try: - with mda.Writer(self.outfile, u.atoms.n_atoms) as W: - W.write(u.atoms) - except AttributeError: # misses __exit__ - raise AssertionError("DCDWriter: does not support with statement") - w = mda.Universe(PSF, self.outfile) - assert_equal(w.trajectory.n_frames, 1, - "with_statement: single frame trajectory has wrong " - "number of frames") - assert_almost_equal(w.atoms.positions, - u.atoms.positions, - 3, - err_msg="with_statement: coordinates do not match") - - -class TestDCDWriter_Issue59(TestCase): - def setUp(self): - """Generate input xtc.""" - self.u = mda.Universe(PSF, DCD) - self.tmpdir = tempdir.TempDir() - self.xtc = self.tmpdir.name + '/dcd-writer-issue59-test.xtc' - wXTC = mda.Writer(self.xtc, self.u.atoms.n_atoms) - for ts in self.u.trajectory: - wXTC.write(ts) - wXTC.close() - - def tearDown(self): - try: - os.unlink(self.xtc) - except OSError: - pass - try: - os.unlink(self.dcd) - except (AttributeError, OSError): - pass - del self.u - del self.tmpdir - - @attr('issue') - def test_issue59(self): - """Test writing of XTC to DCD (Issue 59)""" - xtc = mda.Universe(PSF, self.xtc) - self.dcd = self.tmpdir.name + '/dcd-writer-issue59-test.dcd' - wDCD = mda.Writer(self.dcd, xtc.atoms.n_atoms) - for ts in xtc.trajectory: - wDCD.write(ts) - wDCD.close() - - dcd = mda.Universe(PSF, self.dcd) - - xtc.trajectory.rewind() - dcd.trajectory.rewind() - - assert_array_almost_equal( - xtc.atoms.positions, - dcd.atoms.positions, - 3, - err_msg="XTC -> DCD: DCD coordinates are messed up (Issue 59)") - - def test_OtherWriter(self): - dcd = self.u - wXTC = dcd.trajectory.OtherWriter(self.xtc) - for ts in dcd.trajectory: - wXTC.write(ts) - wXTC.close() - - xtc = mda.Universe(PSF, self.xtc) - xtc.trajectory.rewind() - dcd.trajectory.rewind() - - assert_array_almost_equal( - dcd.atoms.positions, - xtc.atoms.positions, - 2, - err_msg="DCD -> XTC: coordinates are messed up (frame {0:d})".format( - dcd.trajectory.frame)) - xtc.trajectory[3] - dcd.trajectory[3] - assert_array_almost_equal( - dcd.atoms.positions, - xtc.atoms.positions, - 2, - err_msg="DCD -> XTC: coordinates are messed up (frame {0:d})".format( - dcd.trajectory.frame)) - - -class _TestDCDReader_TriclinicUnitcell(TestCase): - - __test__ = False - - def setUp(self): - self.u = mda.Universe(self.topology, self.trajectory) - self.tempdir = tempdir.TempDir() - self.dcd = self.tempdir.name + '/dcd-reader-triclinic.dcd' - - def tearDown(self): - try: - os.unlink(self.dcd) - except (AttributeError, OSError): - pass - del self.u - del self.tempdir - - @attr('issue') - def test_read_triclinic(self): - """test reading of triclinic unitcell (Issue 187) for NAMD or new - CHARMM format (at least since c36b2)""" - for ts, box in zip(self.u.trajectory, - self.ref_dimensions[:, 1:]): - assert_array_almost_equal(ts.dimensions, box, 4, - err_msg="box dimensions A,B,C,alpha," - "beta,gamma not identical at frame " - "{}".format(ts.frame)) - - @attr('issue') - def test_write_triclinic(self): - """test writing of triclinic unitcell (Issue 187) for NAMD or new - CHARMM format (at least since c36b2)""" - with self.u.trajectory.OtherWriter(self.dcd) as w: - for ts in self.u.trajectory: + w = mda.Universe(PSF, outfile) + assert w.trajectory.n_frames == 1 + assert_almost_equal(w.atoms.positions, + u.atoms.positions, + 3, + err_msg="coordinates do not match") + + +def test_write_no_natoms(): + with pytest.raises(ValueError): + mda.Writer('foobar.dcd') + + +def test_writer_trajectory_no_natoms(tmpdir, universe_dcd): + with tmpdir.as_cwd(): + universe_dcd.trajectory.Writer("foo.dcd") + + +class RefCHARMMtriclinicDCD(object): + topology = PSF_TRICLINIC + trajectory = DCD_TRICLINIC + # time(ps) A B C alpha beta gamma (length in Angstrome, angles in degrees) + # dcd starts at t = 1ps + ref_dimensions = np.array([ + [1., 35.44604, 35.06156, 34.1585, 91.32802, 61.73521, 44.40703], + [2., 34.65957, 34.22689, 33.09897, 90.56206, 61.79192, 44.14549], + [3., 34.52772, 34.66422, 33.53881, 90.55859, 63.11228, 40.14044], + [4., 34.43749, 33.38432, 34.02133, 88.82457, 64.98057, 36.77397], + [5., 33.73129, 32.47752, 34.18961, 89.88102, 65.89032, 36.10921], + [6., 33.78703, 31.90317, 34.98833, 90.03092, 66.12877, 35.07141], + [7., 33.24708, 31.18271, 34.9654, 93.11122, 68.17743, 35.73643], + [8., 32.92599, 30.31393, 34.99197, 93.89051, 69.3799, 33.48945], + [9., 32.15295, 30.43056, 34.96157, 96.01416, 71.50115, 32.56111], + [10., 31.99748, 30.21518, 35.24292, 95.85821, 71.08429, 31.85939] + ]) + + +class RefNAMDtriclinicDCD(object): + topology = PSF_NAMD_TRICLINIC + trajectory = DCD_NAMD_TRICLINIC + # vmd topology trajectory + # molinfo 0 get {a b c alpha beta gamma} + # time(ps) A B C alpha beta gamma (length in Angstrome, angles in degrees) + ref_dimensions = np.array([ + [1., 38.426594, 38.393101, 44.759800, 90.000000, 90.000000, 60.028915], + ]) + + +@pytest.mark.parametrize("ref", (RefCHARMMtriclinicDCD, RefNAMDtriclinicDCD)) +def test_read_unitcell_triclinic(ref): + u = mda.Universe(ref.topology, ref.trajectory) + for ts, box in zip(u.trajectory, ref.ref_dimensions[:, 1:]): + assert_array_almost_equal(ts.dimensions, box, 4, + err_msg="box dimensions A,B,C,alpha," + "beta,gamma not identical at frame " + "{}".format(ts.frame)) + + +@pytest.mark.parametrize("ref", (RefCHARMMtriclinicDCD, RefNAMDtriclinicDCD)) +def test_write_unitcell_triclinic(ref, tmpdir): + u = mda.Universe(ref.topology, ref.trajectory) + outfile = 'triclinic.dcd' + with tmpdir.as_cwd(): + with u.trajectory.OtherWriter(outfile) as w: + for ts in u.trajectory: w.write(ts) - w = mda.Universe(self.topology, self.dcd) - for ts_orig, ts_copy in zip(self.u.trajectory, - w.trajectory): + + w = mda.Universe(ref.topology, outfile) + for ts_orig, ts_copy in zip(u.trajectory, w.trajectory): assert_almost_equal(ts_orig.dimensions, ts_copy.dimensions, 4, err_msg="DCD->DCD: unit cell dimensions wrong " "at frame {0}".format(ts_orig.frame)) - del w -class TestDCDReader_CHARMM_Unitcell(_TestDCDReader_TriclinicUnitcell, - RefCHARMMtriclinicDCD): - __test__ = True - - -class TestDCDReader_NAMD_Unitcell(_TestDCDReader_TriclinicUnitcell, - RefNAMDtriclinicDCD): - __test__ = True - - -class TestNCDF2DCD(TestCase): - @dec.skipif(module_not_found("netCDF4"), - "Test skipped because netCDF is not available.") - def setUp(self): - self.u = mda.Universe(PRMncdf, NCDF) - # create the DCD - self.tmpdir = tempdir.TempDir() - self.dcd = self.tmpdir.name + '/ncdf-2-dcd.dcd' - DCD = mda.Writer(self.dcd, n_atoms=self.u.atoms.n_atoms) - for ts in self.u.trajectory: - DCD.write(ts) - DCD.close() - self.w = mda.Universe(PRMncdf, self.dcd) - - def tearDown(self): - try: - os.unlink(self.dcd) - except (AttributeError, OSError): - pass - del self.u - del self.w - del self.tmpdir - - @attr('issue') - def test_unitcell(self): - """NCDFReader: Test that DCDWriter correctly writes the CHARMM - unit cell""" - for ts_orig, ts_copy in zip(self.u.trajectory, - self.w.trajectory): - assert_almost_equal( - ts_orig.dimensions, - ts_copy.dimensions, - 3, - err_msg="NCDF->DCD: unit cell dimensions wrong at frame {0:d}".format( - ts_orig.frame)) - - def test_coordinates(self): - for ts_orig, ts_copy in zip(self.u.trajectory, - self.w.trajectory): - assert_almost_equal( - self.u.atoms.positions, - self.w.atoms.positions, - 3, - err_msg="NCDF->DCD: coordinates wrong at frame {0:d}".format( - ts_orig.frame)) - - -class TestDCDTimestep(BaseTimestepTest): - Timestep = mda.coordinates.DCD.Timestep - name = "DCD" - has_box = True - set_box = True - unitcell = np.array([10., 90., 11., 90., 90., 12.]) - uni_args = (PSF, DCD) - - def test_ts_order_define(self): - """Check that users can hack in a custom unitcell order""" - old = self.Timestep._ts_order - self.ts._ts_order = [0, 2, 5, 1, 3, 4] - self.ts.dimensions = np.array([10, 11, 12, 80, 85, 90]) - assert_allclose(self.ts._unitcell, np.array([10, 80, 11, 85, 90, 12])) - self.ts._ts_order = old - self.ts.dimensions = np.zeros(6) +@pytest.fixture(scope='module') +def ncdf2dcd(tmpdir_factory): + pytest.importorskip("netCDF4") + testfile = tmpdir_factory.mktemp('dcd').join('ncdf2dcd.dcd') + testfile = str(testfile) + ncdf = mda.Universe(PRMncdf, NCDF) + with mda.Writer(testfile, n_atoms=ncdf.atoms.n_atoms) as w: + for ts in ncdf.trajectory: + w.write(ts) + return ncdf, mda.Universe(PRMncdf, testfile) + + +def test_ncdf2dcd_unitcell(ncdf2dcd): + ncdf, dcd = ncdf2dcd + for ts_ncdf, ts_dcd in zip(ncdf.trajectory, dcd.trajectory): + assert_almost_equal(ts_ncdf.dimensions, + ts_dcd.dimensions, + 3) + + +def test_ncdf2dcd_coords(ncdf2dcd): + ncdf, dcd = ncdf2dcd + for ts_ncdf, ts_dcd in zip(ncdf.trajectory, dcd.trajectory): + assert_almost_equal(ts_ncdf.positions, + ts_dcd.positions, + 3) diff --git a/testsuite/MDAnalysisTests/coordinates/test_memory.py b/testsuite/MDAnalysisTests/coordinates/test_memory.py index 672680f8d56..499fa61200e 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_memory.py +++ b/testsuite/MDAnalysisTests/coordinates/test_memory.py @@ -37,7 +37,7 @@ class MemoryReference(BaseReference): 'DCD parser not available. Are you using python 3?') def __init__(self): super(MemoryReference, self).__init__() - + self.topology = PSF self.trajectory = DCD self.universe = mda.Universe(PSF, DCD) @@ -103,7 +103,7 @@ def test_default_memory_layout(self): universe2 = mda.Universe(PSF, DCD, in_memory=True, order='fac') assert_equal(universe1.trajectory.get_array().shape, universe2.trajectory.get_array().shape) - + def test_iteration(self): frames = 0 for i, frame in enumerate(self.reader): diff --git a/testsuite/MDAnalysisTests/coordinates/test_timestep_api.py b/testsuite/MDAnalysisTests/coordinates/test_timestep_api.py index 25b608ebe19..eee666dbb84 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_timestep_api.py +++ b/testsuite/MDAnalysisTests/coordinates/test_timestep_api.py @@ -52,8 +52,7 @@ def test_other_timestep(self): # can't do TRR Timestep here as it always has vels and forces # so isn't actually equal to a position only timestep - for otherTS in [mda.coordinates.DCD.Timestep, - mda.coordinates.TRJ.Timestep, + for otherTS in [mda.coordinates.TRJ.Timestep, mda.coordinates.DMS.Timestep, mda.coordinates.GRO.Timestep, mda.coordinates.TRZ.Timestep, diff --git a/testsuite/MDAnalysisTests/data/coordinates/create_data.py b/testsuite/MDAnalysisTests/data/coordinates/create_data.py index d9f8790a502..9a2b3ede853 100644 --- a/testsuite/MDAnalysisTests/data/coordinates/create_data.py +++ b/testsuite/MDAnalysisTests/data/coordinates/create_data.py @@ -32,6 +32,7 @@ def main(): create_test_trj(u, 'test.xtc') create_test_trj(u, 'test.trr') create_test_trj(u, 'test.gro') + create_test_trj(u, 'test.dcd') if __name__ == '__main__': main() diff --git a/testsuite/MDAnalysisTests/data/coordinates/test.dcd b/testsuite/MDAnalysisTests/data/coordinates/test.dcd new file mode 100644 index 00000000000..1c625476dba Binary files /dev/null and b/testsuite/MDAnalysisTests/data/coordinates/test.dcd differ diff --git a/testsuite/MDAnalysisTests/data/legacy_DCD_NAMD_coords.npy b/testsuite/MDAnalysisTests/data/legacy_DCD_NAMD_coords.npy new file mode 100644 index 00000000000..4bcb7c0f70c Binary files /dev/null and b/testsuite/MDAnalysisTests/data/legacy_DCD_NAMD_coords.npy differ diff --git a/testsuite/MDAnalysisTests/data/legacy_DCD_adk_coords.npy b/testsuite/MDAnalysisTests/data/legacy_DCD_adk_coords.npy new file mode 100644 index 00000000000..6d1a498c158 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/legacy_DCD_adk_coords.npy differ diff --git a/testsuite/MDAnalysisTests/data/legacy_DCD_c36_coords.npy b/testsuite/MDAnalysisTests/data/legacy_DCD_c36_coords.npy new file mode 100644 index 00000000000..065a51cfeaa Binary files /dev/null and b/testsuite/MDAnalysisTests/data/legacy_DCD_c36_coords.npy differ diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index ad5bff96916..c47da37fcc4 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -136,6 +136,7 @@ "Martini_membrane_gro", # for testing the leaflet finder "COORDINATES_XTC", "COORDINATES_TRR", + "COORDINATES_DCD", "COORDINATES_TOPOLOGY", "NUCLsel", "GRO_empty_atom", "GRO_missing_atomname", # for testing GROParser exception raise @@ -147,10 +148,19 @@ "MMTF", "MMTF_gz", "ALIGN_BOUND", # two component bound system "ALIGN_UNBOUND", # two component unbound system + "legacy_DCD_ADK_coords", # frames 5 and 29 read in for adk_dims.dcd using legacy DCD reader + "legacy_DCD_NAMD_coords", # frame 0 read in for SiN_tric_namd.dcd using legacy DCD reader + "legacy_DCD_c36_coords", # frames 1 and 4 read in for tip125_tric_C36.dcd using legacy DCD reader ] from pkg_resources import resource_filename +legacy_DCD_NAMD_coords = resource_filename(__name__, +'data/legacy_DCD_NAMD_coords.npy') +legacy_DCD_ADK_coords = resource_filename(__name__, +'data/legacy_DCD_adk_coords.npy') +legacy_DCD_c36_coords = resource_filename(__name__, +'data/legacy_DCD_c36_coords.npy') AUX_XVG_LOWF = resource_filename(__name__, 'data/test_lowf.xvg') AUX_XVG_HIGHF = resource_filename(__name__, 'data/test_highf.xvg') XVG_BAD_NCOL = resource_filename(__name__, 'data/bad_num_col.xvg') @@ -167,6 +177,7 @@ __name__, 'data/coordinates/test.xyz.bz2') COORDINATES_XTC = resource_filename(__name__, 'data/coordinates/test.xtc') COORDINATES_TRR = resource_filename(__name__, 'data/coordinates/test.trr') +COORDINATES_DCD = resource_filename(__name__, 'data/coordinates/test.dcd') COORDINATES_TOPOLOGY = resource_filename(__name__, 'data/coordinates/test_topology.pdb') PSF = resource_filename(__name__, 'data/adk.psf') diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py new file mode 100644 index 00000000000..15d127871ff --- /dev/null +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -0,0 +1,511 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# 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) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# 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 +# +from __future__ import absolute_import, print_function +from six.moves import zip + +from collections import namedtuple +import os +import string + +import hypothesis.strategies as strategies +from hypothesis import example, given +import numpy as np + +from numpy.testing import (assert_array_almost_equal, assert_equal, + assert_array_equal) + +from MDAnalysis.lib.formats.libdcd import DCDFile, DCD_IS_CHARMM, DCD_HAS_EXTRA_BLOCK + +from MDAnalysisTests.datafiles import ( + DCD, DCD_NAMD_TRICLINIC, legacy_DCD_ADK_coords, legacy_DCD_NAMD_coords, + legacy_DCD_c36_coords, DCD_TRICLINIC) + +import pytest + + +@pytest.mark.parametrize("dcdfile, is_periodic", + [(DCD, False), (DCD_NAMD_TRICLINIC, True), + (DCD_TRICLINIC, True)]) +def test_is_periodic(dcdfile, is_periodic): + with DCDFile(dcdfile) as f: + assert f.is_periodic == is_periodic + + +@pytest.mark.parametrize("dcdfile, natoms", + [(DCD, 3341), (DCD_NAMD_TRICLINIC, 5545), + (DCD_TRICLINIC, 375)]) +def test_read_coordsshape(dcdfile, natoms): + # confirm shape of coordinate data against result from previous + # MDAnalysis implementation of DCD file handling + with DCDFile(dcdfile) as dcd: + dcd_frame = dcd.read() + xyz = dcd_frame[0] + assert xyz.shape == (natoms, 3) + + +@pytest.mark.parametrize( + "dcdfile, unit_cell", + [(DCD, [0., 90., 0., 90., 90., 0.]), + (DCD_NAMD_TRICLINIC, [38.42659378, 0.499563, 38.393102, 0., 0., 44.7598]), + (DCD_TRICLINIC, + [30.841836, 14.578635, 31.780088, 9.626323, -2.60815, 32.67009])]) +def test_read_unit_cell(dcdfile, unit_cell): + # confirm unit cell read against result from previous + # MDAnalysis implementation of DCD file handling + with DCDFile(dcdfile) as dcd: + dcd_frame = dcd.read() + assert_array_almost_equal(dcd_frame.unitcell, unit_cell) + + +def test_seek_over_max(): + with DCDFile(DCD) as dcd: + with pytest.raises(EOFError): + dcd.seek(102) + + +@pytest.fixture +def dcd(): + with DCDFile(DCD) as dcd: + yield dcd + + +@pytest.mark.parametrize("new_frame", (10, 42, 21)) +def test_seek_normal(new_frame, dcd): + # frame seek within range is tested + dcd.seek(new_frame) + assert dcd.tell() == new_frame + + +def test_seek_negative(dcd): + with pytest.raises(IOError): + dcd.seek(-78) + + +def test_iteration(dcd): + num_iters = 10 + for _ in range(num_iters): + dcd.__next__() + assert dcd.tell() == num_iters + + +def test_open_wrong_mode(): + with pytest.raises(IOError): + DCDFile(DCD, 'e') + + +def test_raise_not_existing(): + with pytest.raises(IOError): + DCDFile('foo') + + +def test_zero_based_frames_counting(dcd): + assert dcd.tell() == 0 + + +@pytest.mark.parametrize("dcdfile, natoms", + [(DCD, 3341), (DCD_NAMD_TRICLINIC, 5545), + (DCD_TRICLINIC, 375)]) +def test_natoms(dcdfile, natoms): + with DCDFile(dcdfile) as dcd: + assert dcd.header['natoms'] == natoms + + +def test_read_closed(dcd): + dcd.close() + with pytest.raises(IOError): + dcd.read() + + +@pytest.mark.parametrize("dcdfile, nframes", + [(DCD, 98), (DCD_NAMD_TRICLINIC, 1), + (DCD_TRICLINIC, 10)]) +def test_length_traj(dcdfile, nframes): + with DCDFile(dcdfile) as dcd: + assert len(dcd) == nframes + + +def test_read_write_mode_file(tmpdir): + fname = str(tmpdir.join('foo')) + with DCDFile(fname, 'w') as f: + with pytest.raises(IOError): + f.read() + + +def test_iterating_twice(dcd): + with dcd as f: + for i, _ in enumerate(f): + assert_equal(i + 1, f.tell()) + # second iteration should work from start again + for i, _ in enumerate(f): + assert_equal(i + 1, f.tell()) + + +DCD_HEADER = '''* DIMS ADK SEQUENCE FOR PORE PROGRAM * WRITTEN BY LIZ DENNING (6.2008) * DATE: 6/ 6/ 8 17:23:56 CREATED BY USER: denniej0 ''' +DCD_NAMD_TRICLINIC_HEADER = 'Created by DCD pluginREMARKS Created 06 July, 2014 at 17:29Y5~CORD,' +DCD_TRICLINIC_HEADER = '* CHARMM TRICLINIC BOX TESTING * (OLIVER BECKSTEIN 2014) * BASED ON NPTDYN.INP : SCOTT FELLER, NIH, 7/15/95 * TEST EXTENDED SYSTEM CONSTANT PRESSURE AND TEMPERATURE * DYNAMICS WITH WATER BOX. * DATE: 7/ 7/14 13:59:46 CREATED BY USER: oliver ' + + +@pytest.mark.parametrize("dcdfile, remarks", ( + (DCD, DCD_HEADER), (DCD_NAMD_TRICLINIC, DCD_NAMD_TRICLINIC_HEADER), + (DCD_TRICLINIC, DCD_TRICLINIC_HEADER))) +def test_header_remarks(dcdfile, remarks): + # confirm correct header remarks section reading + with DCDFile(dcdfile) as f: + assert len(f.header['remarks']) == len(remarks) + + +@pytest.mark.parametrize("dcdfile, legacy_data, frames", + ((DCD, legacy_DCD_ADK_coords, [5, 29]), + (DCD_NAMD_TRICLINIC, legacy_DCD_NAMD_coords, [0]), + (DCD_TRICLINIC, legacy_DCD_c36_coords, [1, 4]))) +def test_read_coord_values(dcdfile, legacy_data, frames): + # test the actual values of coordinates read in versus + # stored values read in by the legacy DCD handling framework + # to reduce repo storage burden, we only compare for a few + # randomly selected frames + legacy = np.load(legacy_data) + with DCDFile(dcdfile) as dcd: + for index, frame_num in enumerate(frames): + dcd.seek(frame_num) + actual_coords = dcd.read()[0] + desired_coords = legacy[index] + assert_array_equal(actual_coords, desired_coords) + + +@pytest.mark.parametrize("dcdfile, legacy_data, frame_idx", + ((DCD, legacy_DCD_ADK_coords, [5, 29]), + (DCD_NAMD_TRICLINIC, legacy_DCD_NAMD_coords, [0]), + (DCD_TRICLINIC, legacy_DCD_c36_coords, [1, 4]))) +def test_readframes(dcdfile, legacy_data, frame_idx): + legacy = np.load(legacy_data) + with DCDFile(dcdfile) as dcd: + frames = dcd.readframes() + xyz = frames.xyz + assert_equal(len(xyz), len(dcd)) + for index, frame_num in enumerate(frame_idx): + assert_array_almost_equal(xyz[frame_num], legacy[index]) + + +def test_write_header(tmpdir): + # test that _write_header() can produce a very crude + # header for a new / empty file + testfile = str(tmpdir.join('test.dcd')) + with DCDFile(testfile, 'w') as dcd: + dcd.write_header( + remarks='Crazy!', + natoms=22, + istart=12, + nsavc=10, + delta=0.02, + is_periodic=1) + + with DCDFile(testfile) as dcd: + header = dcd.header + assert header['remarks'] == 'Crazy!' + assert header['natoms'] == 22 + assert header['istart'] == 12 + assert header['is_periodic'] == 1 + assert header['nsavc'] == 10 + assert np.allclose(header['delta'], .02) + + +def test_write_no_header(tmpdir): + fname = str(tmpdir.join('test.dcd')) + with DCDFile(fname, 'w') as dcd: + with pytest.raises(IOError): + dcd.write(np.ones(3), np.ones(6)) + + +def test_write_header_twice(tmpdir): + # an IOError should be raised if a duplicate + # header writing is attempted + + header = { + "remarks": 'Crazy!', + "natoms": 22, + "istart": 12, + "nsavc": 10, + "delta": 0.02, + "is_periodic": 1 + } + + fname = str(tmpdir.join('test.dcd')) + with DCDFile(fname, 'w') as dcd: + dcd.write_header(**header) + with pytest.raises(IOError): + dcd.write_header(**header) + + +def test_write_header_wrong_mode(dcd): + # an exception should be raised on any attempt to use + # write_header with a DCDFile object in 'r' mode + with pytest.raises(IOError): + dcd.write_header( + remarks='Crazy!', + natoms=22, + istart=12, + nsavc=10, + delta=0.02, + is_periodic=1) + + +def test_write_mode(dcd): + # ensure that writing of DCD files only occurs with properly + # opened files + with pytest.raises(IOError): + dcd.write(xyz=np.zeros((3, 3)), box=np.zeros(6, dtype=np.float64)) + + +def write_dcd(in_name, out_name, remarks='testing', header=None): + with DCDFile(in_name) as f_in, DCDFile(out_name, 'w') as f_out: + if header is None: + header = f_in.header + f_out.write_header(**header) + for frame in f_in: + f_out.write(xyz=frame.xyz, box=frame.unitcell) + + +@given(remarks=strategies.text( + alphabet=string.printable, min_size=0, + max_size=240)) # handle the printable ASCII strings +@example(remarks='') +def test_written_remarks_property(remarks, tmpdir, dcd): + # property based testing for writing of a wide range of string + # values to REMARKS field + testfile = str(tmpdir.join('test.dcd')) + header = dcd.header + header['remarks'] = remarks + write_dcd(DCD, testfile, header=header) + expected_remarks = remarks[:240] + with DCDFile(testfile) as f: + assert f.header['remarks'] == expected_remarks + + +@pytest.fixture(scope='session') +def written_dcd(tmpdir_factory): + testfile = 'test.dcd' + with DCDFile(DCD) as dcd: + header = dcd.header + testfile = tmpdir_factory.mktemp('dcd').join('test.dcd') + testfile = str(testfile) + write_dcd(DCD, testfile) + Result = namedtuple("Result", "testfile, header, orgfile") + return Result(testfile, header, DCD) + + +def test_written_header(written_dcd): + header = written_dcd.header + with DCDFile(written_dcd.testfile) as dcd: + dcdheader = dcd.header + assert dcdheader == header + + +def test_written_num_frames(written_dcd): + with DCDFile(written_dcd.testfile) as dcd, DCDFile( + written_dcd.orgfile) as other: + assert len(dcd) == len(other) + + +def test_written_dcd_coordinate_data_shape(written_dcd): + with DCDFile(written_dcd.testfile) as dcd, DCDFile( + written_dcd.orgfile) as other: + for frame, other_frame in zip(dcd, other): + assert frame.xyz.shape == other_frame.xyz.shape + + +def test_written_seek(written_dcd): + # ensure that we can seek properly on written DCD file + with DCDFile(written_dcd.testfile) as f: + f.seek(40) + assert_equal(f.tell(), 40) + + +def test_written_coord_match(written_dcd): + with DCDFile(written_dcd.testfile) as test, DCDFile( + written_dcd.orgfile) as ref: + for frame, o_frame in zip(test, ref): + assert_array_almost_equal(frame.xyz, o_frame.xyz) + + +def test_written_unit_cell(written_dcd): + with DCDFile(written_dcd.testfile) as test, DCDFile( + written_dcd.orgfile) as ref: + for frame, o_frame in zip(test, ref): + assert_array_almost_equal(frame.unitcell, o_frame.unitcell) + + +@pytest.mark.parametrize( + "dtype", (np.int32, np.int64, np.float32, np.float64, int, float)) +def test_write_all_dtypes(tmpdir, dtype): + fname = str(tmpdir.join('foo.dcd')) + with DCDFile(fname, 'w') as out: + natoms = 10 + xyz = np.ones((natoms, 3), dtype=dtype) + box = np.ones(6, dtype=dtype) + out.write_header( + remarks='test', + natoms=natoms, + is_periodic=1, + delta=1, + nsavc=1, + istart=1) + out.write(xyz=xyz, box=box) + + +@pytest.mark.parametrize("array_like", (np.array, list)) +def test_write_array_like(tmpdir, array_like): + fname = str(tmpdir.join('foo.dcd')) + with DCDFile(fname, 'w') as out: + natoms = 10 + xyz = array_like([[1, 1, 1] for i in range(natoms)]) + box = array_like([i for i in range(6)]) + out.write_header( + remarks='test', + natoms=natoms, + is_periodic=1, + delta=1, + nsavc=1, + istart=1) + out.write(xyz=xyz, box=box) + + +def test_write_wrong_shape_xyz(tmpdir): + fname = str(tmpdir.join('foo.dcd')) + with DCDFile(fname, 'w') as out: + natoms = 10 + xyz = np.ones((natoms + 1, 3)) + box = np.ones(6) + out.write_header( + remarks='test', + natoms=natoms, + is_periodic=1, + delta=1, + nsavc=1, + istart=1) + with pytest.raises(ValueError): + out.write(xyz=xyz, box=box) + + +def test_write_wrong_shape_box(tmpdir): + fname = str(tmpdir.join('foo.dcd')) + with DCDFile(fname, 'w') as out: + natoms = 10 + xyz = np.ones((natoms, 3)) + box = np.ones(7) + out.write_header( + remarks='test', + natoms=natoms, + is_periodic=1, + delta=1, + nsavc=1, + istart=1) + with pytest.raises(ValueError): + out.write(xyz=xyz, box=box) + + +@pytest.mark.parametrize("dcdfile", (DCD, DCD_TRICLINIC, DCD_NAMD_TRICLINIC)) +def test_relative_frame_sizes(dcdfile): + # the first frame of a DCD file should always be >= in size + # to subsequent frames, as the first frame contains the same + # atoms + (optional) fixed atoms + with DCDFile(dcdfile) as dcd: + first_frame_size = dcd._firstframesize + general_frame_size = dcd._framesize + + assert first_frame_size >= general_frame_size + + +@pytest.mark.parametrize("dcdfile", (DCD, DCD_TRICLINIC, DCD_NAMD_TRICLINIC)) +def test_file_size_breakdown(dcdfile): + # the size of a DCD file is equivalent to the sum of the header + # size, first frame size, and (N - 1 frames) * size per general + # frame + + expected = os.path.getsize(dcdfile) + with DCDFile(dcdfile) as dcd: + actual = dcd._header_size + dcd._firstframesize + ( + (dcd.n_frames - 1) * dcd._framesize) + assert actual == expected + + +@pytest.mark.parametrize("dcdfile", (DCD, DCD_TRICLINIC, DCD_NAMD_TRICLINIC)) +def test_nframessize_int(dcdfile): + # require that the (nframessize / framesize) value used by DCDFile + # is an integer (because nframessize / framesize + 1 = total frames, + # which must also be an int) + filesize = os.path.getsize(dcdfile) + with DCDFile(dcdfile) as dcd: + nframessize = filesize - dcd._header_size - dcd._firstframesize + assert float(nframessize) % float(dcd._framesize) == 0 + + +@pytest.mark.parametrize( + "slice, length", + [([None, None, None], 98), ([0, None, None], 98), ([None, 98, None], 98), + ([None, None, 1], 98), ([None, None, -1], 98), ([2, 6, 2], 2), + ([0, 10, None], 10), ([2, 10, None], 8), ([0, 1, 1], 1), ([1, 1, 1], 0), + ([1, 2, 1], 1), ([1, 2, 2], 1), ([1, 4, 2], 2), ([1, 4, 4], 1), + ([0, 5, 5], 1), ([3, 5, 1], 2), ([4, 0, -1], 4), ([5, 0, -2], 3), + ([5, 0, -4], 2)]) +def test_readframes_slices(slice, length, dcd): + start, stop, step = slice + allframes = dcd.readframes().xyz + frames = dcd.readframes(start=start, stop=stop, step=step) + xyz = frames.xyz + assert len(xyz) == length + assert_array_almost_equal(xyz, allframes[start:stop:step]) + + +@pytest.mark.parametrize("order, shape", ( + ('fac', (98, 3341, 3)), + ('fca', (98, 3, 3341)), + ('afc', (3341, 98, 3)), + ('acf', (3341, 3, 98)), + ('caf', (3, 3341, 98)), + ('cfa', (3, 98, 3341)), )) +def test_readframes_order(order, shape, dcd): + x = dcd.readframes(order=order).xyz + assert x.shape == shape + + +@pytest.mark.parametrize("indices", [[1, 2, 3, 4], [5, 10, 15, 19], + [9, 4, 2, 0, 50]]) +def test_readframes_atomindices(indices, dcd): + allframes = dcd.readframes(order='afc').xyz + frames = dcd.readframes(indices=indices, order='afc') + xyz = frames.xyz + assert len(xyz) == len(indices) + assert_array_almost_equal(xyz, allframes[indices]) + + +def test_write_random_unitcell(tmpdir): + testname = str(tmpdir.join('test.dcd')) + testname = 'test.dcd' + rstate = np.random.RandomState(1178083) + random_unitcells = rstate.uniform( + high=80, size=(98, 6)).astype(np.float64) + + with DCDFile(DCD) as f_in, DCDFile(testname, 'w') as f_out: + header = f_in.header + header['is_periodic'] = True + f_out.write_header(**header) + for index, frame in enumerate(f_in): + f_out.write(xyz=frame.xyz, box=random_unitcells[index]) + + with DCDFile(testname) as test: + for index, frame in enumerate(test): + assert_array_almost_equal(frame.unitcell, + random_unitcells[index]) diff --git a/testsuite/setup.py b/testsuite/setup.py index 006e0a4c3f3..a61410bc470 100755 --- a/testsuite/setup.py +++ b/testsuite/setup.py @@ -214,8 +214,11 @@ def dynamic_author_list(): install_requires=[ 'MDAnalysis=={0!s}'.format(RELEASE), # same as this release! 'numpy>=1.10.4', - 'pytest>=3.1.2', 'nose>=1.3.7', + 'pytest>=3.1.2', + 'pytest-cov', + 'pytest-xdist', + 'hypothesis', 'psutil>=4.0.2', 'mock>=2.0.0', ],