diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index e8e44ebf572..2768c734e43 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -64,8 +64,6 @@ Classes ------- -.. autoclass:: Timestep - :inherited-members: .. autoclass:: DCDReader :inherited-members: .. autoclass:: DCDWriter @@ -87,563 +85,182 @@ 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 +from ..lib.formats.libdcd import DCDFile # dcdtimeseries is implemented with Pyrex - hopefully all dcd reading functionality can move to pyrex -from . import dcdtimeseries - - - -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]) - - # 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) +# from . import dcdtimeseries -class DCDWriter(base.WriterBase): - """Writes to a DCD file - - 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. - - :Methods: - ``d = DCDWriter(dcdfilename, n_atoms, start, step, delta, remarks)`` - - .. Note:: - - 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. +class DCDReader(base.ReaderBase): + """DCD Reader - .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187 """ 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): - """Create a new DCDWriter - - :Arguments: - *filename* - name of output file - *n_atoms* - number of atoms in dcd file - *start* - starting timestep - *step* - skip between subsequent timesteps (indicate that *step* MD - integrator steps (!) make up one trajectory frame); default is 1. - *delta* - timestep (MD integrator time step (!), in AKMA units); default is - 20.45482949774598 (corresponding to 1 ps). - *remarks* - comments to annotate dcd file - *dt* - **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* - 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". + def __init__(self, filename, convert_units=True, dt=None, **kwargs): + """Parameters + ---------- + 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. """ - 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 + super(DCDReader, self).__init__(filename, + convert_units=convert_units, + **kwargs) + self._file = DCDFile(self.filename) + self.n_atoms = self._file.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. - """ - # 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 - - *ts* - timestep object containing coordinates to be written to dcd file - - .. 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. - ''' - 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 - 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 + delta = mdaunits.convert(self._file.delta, self.units['time'], 'ps') + if dt is None: + dt = delta * self._file.nsavc + self.skip_timestep = self._file.nsavc - def convert_dimensions_to_unitcell(self, ts, _ts_order=Timestep._ts_order): - """Read dimensions from timestep *ts* and return appropriate native unitcell. - - .. SeeAlso:: :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._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._frame_to_ts(frame, self.ts) + # these should only be initialized once + self.ts.dt = dt + self.ts.dimensions = frame.unitcell + if self.convert_units: + 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 + """close reader""" + self._file.close() + @property + def n_frames(self): + """number of frames in trajectory""" + return len(self._file) -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 - ``data = dcd.correl(...)`` - populate a :class:`MDAnalysis.core.Timeseries.Collection` object with computed timeseries - - .. 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. - - .. _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 _reopen(self): + """reopen trajectory""" + self.ts.frame = 0 + self._frame = -1 + self._file.close() + self._file.open(self.filename.encode('utf-8'), 'r') + + def _read_frame(self, i): + """read frame i""" + self._frame = i - 1 + self._file.seek(i) + return self._read_next_timestep() - def __init__(self, dcdfilename, **kwargs): - super(DCDReader, self).__init__(dcdfilename, **kwargs) + def _read_next_timestep(self, ts=None): + """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 + frame = self._file.read() + self._frame += 1 + ts = self._frame_to_ts(frame, ts) + self.ts = ts + return ts - self.dcdfilename = self.filename # dcdfilename is legacy - self.dcdfile = None # set right away because __del__ checks + 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, **kwargs) - # 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) + def _frame_to_ts(self, frame, ts): + """convert a trr-frame to a mda TimeStep""" + ts.frame = self._frame + ts.time = ts.frame * self.ts.dt + ts.data['step'] = self._file.tell() - self.dcdfile = open(self.filename, 'rb') - self.n_atoms = 0 - self.n_frames = 0 - self.fixed = 0 - self.periodic = False + ts.dimensions = frame.unitcell + ts.positions = frame.x - # This reads skip_timestep and delta from header - self._read_dcd_header() + if self.convert_units: + self.convert_pos_from_native(ts.dimensions[:3]) + self.convert_pos_from_native(ts.positions) - # Convert delta to ps - delta = mdaunits.convert(self.delta, self.units['time'], 'ps') + return ts - 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. + @property + def dt(self): + return self.ts.dt - """ - # 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 _reopen(self): - self.ts.frame = -1 - self._reset_dcd_read() +class DCDWriter(base.WriterBase): + """Base class for libmdaxdr file formats xtc and trr""" - def _read_next_timestep(self, ts=None): - """Read the next frame + format = 'DCD' + multiframe = True + flavor = 'CHARMM' + units = {'time': 'AKMA', 'length': 'Angstrom'} - .. versionchanged 0.11.0:: - Native frame read into ts._frame, ts.frame naively iterated + def __init__(self, filename, n_atoms, convert_units=True, step=1, dt=None, **kwargs): """ - 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 - return ts - - def _read_frame(self, frame): - """Skip to frame and read - - .. versionchanged:: 0.11.0 - Native frame read into ts._frame, ts.frame naively set to frame + Parameters + ---------- + filename : str + 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) + number of steps between frames to be written + dt : float (optional) + use this time step in DCD. If ``None`` guess from written TimeStep + **kwargs : dict + General writer arguments """ - 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 + self.filename = filename + self._convert_units = convert_units + self.n_atoms = n_atoms + self._file = DCDFile(self.filename, 'w') + self.step = step + self.dt = 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 - - :Arguments: - *asel* - :class:`~MDAnalysis.core.groups.AtomGroup` object - Defaults to None, in which case the full set of coordinate data - is returned. - *start*, *stop*, *step* - A range of the trajectory to access, with start being inclusive - and stop being exclusive. - *format* - the order/shape of the return data array, corresponding - to (a)tom, (f)rame, (c)oordinates all six combinations - of 'a', 'f', 'c' are allowed ie "fac" - return array - where the shape is (frame, number of atoms, - coordinates) - :Deprecated: - *skip* - 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) - - 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") - 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) - - def correl(self, timeseries, start=None, stop=None, step=None, skip=None): - """Populate a TimeseriesCollection object with timeseries computed from the trajectory - - :Arguments: - *timeseries* - :class:`MDAnalysis.core.Timeseries.TimeseriesCollection` - *start, stop, step* - A subset of the trajectory to use, with start being inclusive - and stop being exclusive. - :Deprecated: - *skip* - 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) - - start, stop, step = self.check_slice_indices(start, stop, step) - atomlist = timeseries._getAtomList() - format = timeseries._getFormat() - lowerb, upperb = timeseries._getBounds() - sizedata = timeseries._getDataSize() - atomcounts = timeseries._getAtomCounts() - auxdata = timeseries._getAuxData() - return self._read_timecorrel(atomlist, atomcounts, format, auxdata, - sizedata, lowerb, upperb, start, stop, step) + def write_next_timestep(self, ts): + """Write timestep object into trajectory. - 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. - - All values can be changed through keyword arguments. - - :Arguments: - *filename* - filename of the output DCD trajectory - :Keywords: - *n_atoms* - number of atoms - *start* - number of the first recorded MD step - *step* - indicate that *step* MD steps (!) make up one trajectory frame - *delta* - MD integrator time step (!), in AKMA units - *dt* - **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. - *remarks* - string that is stored in the DCD header [XXX -- max length?] - - :Returns: :class:`DCDWriter` - - .. 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* ! - - .. SeeAlso:: :class:`DCDWriter` has detailed argument description + Parameters + ---------- + ts: TimeStep + + See Also + -------- + .write(AtomGroup/Universe/TimeStep) + The normal write() method 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() + time = ts.time + step = ts.frame + dimensions = ts.dimensions - @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=False) + dimensions = self.convert_dimensions_to_unitcell(ts, inplace=False) + dt = self.dt if self.dt is not None else ts.dt + dt = mdaunits.convert(dt, 'ps', self.units['time']) -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) + self._file.write(xyz=xyz, box=dimensions, step=step, natoms=xyz.shape[0], + charmm=1, time_step=dt, ts_between_saves=self.step, remarks='test') -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 close(self): + """close trajectory""" + self._file.close() -#DCDReader._read_timeseries = types.MethodType(dcdtimeseries.__read_timeseries, None, DCDReader) -DCDReader._read_timecorrel = types.MethodType(dcdtimeseries.__read_timecorrel, None, DCDReader) + def __del__(self): + self.close() 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 5a0fe576ab8..3e395d6a6f1 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/lib/formats/include/correl.h b/package/MDAnalysis/lib/formats/include/correl.h new file mode 100644 index 00000000000..f3dacb042e7 --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/correl.h @@ -0,0 +1,188 @@ +/* -*- Mode: C; tab-width: 4; indent-tabs-mode:nil; -*- */ +/* vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 */ +/* + MDAnalysis --- http://mdanalysis.googlecode.com + Copyright (c) 2006-2014 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, + in press. +*/ + +#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/lib/formats/include/endianswap.h b/package/MDAnalysis/lib/formats/include/endianswap.h new file mode 100644 index 00000000000..6f1ced37483 --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/endianswap.h @@ -0,0 +1,173 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2003 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: endianswap.h,v $ + * $Author: eamon $ $Locker: $ $State: Exp $ + * $Revision: 1.3 $ $Date: 2004/04/16 15:37:00 $ + * + *************************************************************************** + * DESCRIPTION: + * Byte swapping routines used in various plugins + * There are two versions of each routine, one that's safe to use in + * all cases (but is slow) and one that is only safe to use on memory + * addresses that are aligned to the word size that's being byte-swapped + * but are much much much faster. Use the aligned versions of these + * routines whenever possible. The 'ndata' length count parameters and + * internal loops should be safe to use on huge memory arrays on 64-bit + * machines. + * + ***************************************************************************/ + +#ifndef ENDIAN_SWAP_H +#define ENDIAN_SWAP_H + +/* works on unaligned 2-byte quantities */ +static void swap2_unaligned(void *v, long ndata) { + long i; + char * dataptr = (char *) v; + char tmp; + + for (i = 0; i < ndata-1; i += 2) { + tmp = dataptr[i]; + dataptr[i] = dataptr[i+1]; + dataptr[i+1] = tmp; + } +} + + +/* works on unaligned 4-byte quantities */ +static void swap4_unaligned(void *v, long ndata) { + long i; + char *dataptr; + char tmp; + + dataptr = (char *) v; + for (i=0; i>8)&0xff) | ((*N&0xff)<<8)); + } +} + + +/* Only works with aligned 4-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +static void swap4_aligned(void *v, long ndata) { + int *data = (int *) v; + long i; + int *N; + for (i=0; i>24)&0xff) | ((*N&0xff)<<24) | + ((*N>>8)&0xff00) | ((*N&0xff00)<<8)); + } +} + + +/* Only works with aligned 8-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +static void swap8_aligned(void *v, long ndata) { + /* Use int* internally to prevent bugs caused by some compilers */ + /* and hardware that would potentially load data into an FP reg */ + /* and hose everything, such as the old "jmemcpy()" bug in NAMD */ + int *data = (int *) v; + long i; + int *N; + int t0, t1; + + for (i=0; i>24)&0xff) | ((t0&0xff)<<24) | + ((t0>>8)&0xff00) | ((t0&0xff00)<<8)); + + t1 = N[1]; + t1=(((t1>>24)&0xff) | ((t1&0xff)<<24) | + ((t1>>8)&0xff00) | ((t1&0xff00)<<8)); + + N[0] = t1; + N[1] = t0; + } +} + +#if 0 +/* Other implementations that might be faster in some cases */ + +/* swaps the endianism of an eight byte word. */ +void mdio_swap8(double *i) { + char c; + char *n; + n = (char *) i; + c = n[0]; + n[0] = n[7]; + n[7] = c; + c = n[1]; + n[1] = n[6]; + n[6] = c; + c = n[2]; + n[2] = n[5]; + n[5] = c; + c = n[3]; + n[3] = n[4]; + n[4] = c; +} + +#endif + +#endif + diff --git a/package/MDAnalysis/lib/formats/include/fastio.h b/package/MDAnalysis/lib/formats/include/fastio.h new file mode 100644 index 00000000000..4eeda73e21f --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/fastio.h @@ -0,0 +1,458 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2009 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: fastio.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.20 $ $Date: 2009/04/29 15:45:29 $ + * + *************************************************************************** + * DESCRIPTION: + * This is a simple abstraction layer for system-dependent I/O calls + * that allow plugins to do binary I/O using the fastest possible method. + * + * This code is intended for use by binary trajectory reader plugins that + * work with multi-gigabyte data sets, reading only binary data. + * + ***************************************************************************/ + +#ifndef FASTIO_H +#define FASTIO_H + +/* Compiling on windows */ +#if defined(_MSC_VER) + +#if 1 +/* use native Windows I/O calls */ +#define FASTIO_NATIVEWIN32 1 + +#include +#include + +typedef HANDLE fio_fd; +typedef LONGLONG fio_size_t; +typedef void * fio_caddr_t; + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR FILE_CURRENT +#define FIO_SEEK_SET FILE_BEGIN +#define FIO_SEEK_END FILE_END + +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 + +typedef FILE * fio_fd; +typedef size_t fio_size_t; /* MSVC doesn't uinversally support ssize_t */ +typedef void * fio_caddr_t; /* MSVC doesn't universally support caddr_t */ + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR SEEK_CUR +#define FIO_SEEK_SET SEEK_SET +#define FIO_SEEK_END SEEK_END + +static int fio_open(const char *filename, int mode, fio_fd *fd) { + char * modestr; + FILE *fp; + + if (mode == FIO_READ) + modestr = "rb"; + + if (mode == FIO_WRITE) + modestr = "wb"; + + fp = fopen(filename, modestr); + if (fp == NULL) { + return -1; + } else { + *fd = fp; + return 0; + } +} + +static int fio_fclose(fio_fd fd) { + return fclose(fd); +} + +static fio_size_t fio_fread(void *ptr, fio_size_t size, + fio_size_t nitems, fio_fd fd) { + return fread(ptr, size, nitems, fd); +} + +static fio_size_t fio_readv(fio_fd fd, const fio_iovec * iov, int iovcnt) { + int i; + fio_size_t len = 0; + + for (i=0; i +#include +#include +#include + +typedef int fio_fd; +typedef off_t fio_size_t; /* off_t is 64-bits with LFS builds */ + +/* enable use of kernel readv() if available */ +#if defined(__sun) || defined(__APPLE_CC__) || defined(__linux) +#define USE_KERNEL_READV 1 +#endif + +typedef void * fio_caddr_t; + +#if defined(USE_KERNEL_READV) +#include +typedef struct iovec fio_iovec; +#else + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; +#endif + + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR SEEK_CUR +#define FIO_SEEK_SET SEEK_SET +#define FIO_SEEK_END SEEK_END + +static int fio_open(const char *filename, int mode, fio_fd *fd) { + int nfd; + int oflag = 0; + + if (mode == FIO_READ) + oflag = O_RDONLY; + + if (mode == FIO_WRITE) + oflag = O_WRONLY | O_CREAT | O_TRUNC; + + nfd = open(filename, oflag, 0666); + if (nfd < 0) { + return -1; + } else { + *fd = nfd; + return 0; + } +} + +static int fio_fclose(fio_fd fd) { + return close(fd); +} + +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; + int cnt = 0; + + for (i=0; i= 0) + return 0; /* success (emulate behavior of fseek) */ + else + return -1; /* failure (emulate behavior of fseek) */ +} + +static fio_size_t fio_ftell(fio_fd fd) { + return lseek(fd, 0, SEEK_CUR); +} + +#endif + + +/* higher level routines that are OS independent */ + +static int fio_write_int32(fio_fd fd, int i) { + return (fio_fwrite(&i, 4, 1, fd) != 1); +} + +static int fio_read_int32(fio_fd fd, int *i) { + return (fio_fread(i, 4, 1, fd) != 1); +} + +static int fio_write_str(fio_fd fd, const char *str) { + int len = strlen(str); + return (fio_fwrite((void *) str, len, 1, fd) != 1); +} + +#endif // FASTIO_H diff --git a/package/MDAnalysis/lib/formats/include/readdcd.h b/package/MDAnalysis/lib/formats/include/readdcd.h new file mode 100644 index 00000000000..b248097f4db --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/readdcd.h @@ -0,0 +1,758 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2003 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: readdcd.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.32 $ $Date: 2004/09/21 20:52:37 $ + * + ***************************************************************************/ + +#ifndef READ_DCD_H +#define READ_DCD_H + +#include +#include +#include +#include +#include "endianswap.h" +#include "fastio.h" + +/* DEFINE ERROR CODES THAT MAY BE RETURNED BY DCD ROUTINES */ +#define DCD_SUCCESS 0 /* No problems */ +#define DCD_EOF -1 /* Normal EOF */ +#define DCD_DNE -2 /* DCD file does not exist */ +#define DCD_OPENFAILED -3 /* Open of DCD file failed */ +#define DCD_BADREAD -4 /* read call on DCD file failed */ +#define DCD_BADEOF -5 /* premature EOF found in DCD file */ +#define DCD_BADFORMAT -6 /* format of DCD file is wrong */ +#define DCD_FILEEXISTS -7 /* output file already exists */ +#define DCD_BADMALLOC -8 /* malloc failed */ + +/* + * Read the header information from a dcd file. + * Input: fd - a file struct opened for binary reading. + * Output: 0 on success, negative error code on failure. + * Side effects: *natoms set to number of atoms per frame + * *nsets set to number of frames in dcd file + * *istart set to starting timestep of dcd file + * *nsavc set to timesteps between dcd saves + * *delta set to value of trajectory timestep + * *nfixed set to number of fixed atoms + * *freeind may be set to heap-allocated space + * *reverse set to one if reverse-endian, zero if not. + * *charmm set to internal code for handling charmm data. + */ +static int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, int *nsavc, + double *delta, int *nfixed, int **freeind, + float **fixedcoords, int *reverse, int *charmm, + char **remarks, int *len_remarks); + +/* + * Read a dcd timestep from a dcd file + * Input: fd - a file struct opened for binary reading, from which the + * header information has already been read. + * natoms, nfixed, first, *freeind, reverse, charmm - the corresponding + * items as set by read_dcdheader + * first - true if this is the first frame we are reading. + * x, y, z: space for natoms each of floats. + * unitcell - space for six floats to hold the unit cell data. + * Not set if no unit cell data is present. + * Output: 0 on success, negative error code on failure. + * Side effects: x, y, z contain the coordinates for the timestep read. + * unitcell holds unit cell data if present. + */ +static int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, + double *unitcell, int nfixed, int first, int *freeind, + float *fixedcoords, int reverse, int charmm); + +/* + * Read a subset of a timestep from a dcd file + * Input: fd - a file struct opened for binary reading, from which the + * header information has already been read + * natoms, nfixed, first, *freeind, reverse, charmm - the corresponding + * items as set by read_dcdheader + * first - true if this is the first frame we are reading. + * lowerb, upperb - the range of atoms to read data for + * x, y, z: space for upperb-lowerb+1 each of floats + * unitcell - space for six floats to hold the unit cell data. + * Not set if no unit cell data is present. + * Ouput: 0 on success, negative error code on failure. + * Side effects: x, y, z contain coordinates for the range of atoms + * 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, + double *unitcell, int nfixed, int first, int *freeind, + float *fixedcoords, int reverse, int charmm); + +/* + * Skip past a timestep. If there are fixed atoms, this cannot be used with + * the first timestep. + * Input: fd - a file struct from which the header has already been read + * natoms - number of atoms per timestep + * nfixed - number of fixed atoms + * charmm - charmm flags as returned by read_dcdheader + * Output: 0 on success, negative error code on failure. + * Side effects: One timestep will be skipped; fd will be positioned at the + * next timestep. + */ +static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numstep); + +/* + * clean up dcd data + * Input: nfixed, freeind - elements as returned by read_dcdheader + * Output: None + * Side effects: Space pointed to by freeind is freed if necessary. + */ +static void close_dcd_read(int *freeind, float *fixedcoords); + +/* + * Write a header for a new dcd file + * Input: fd - file struct opened for binary writing + * remarks - string to be put in the remarks section of the header. + * The string will be truncated to 70 characters. + * natoms, istart, nsavc, delta - see comments in read_dcdheader + * Output: 0 on success, negative error code on failure. + * Side effects: Header information is written to the dcd file. + */ +static int write_dcdheader(fio_fd fd, const char *remarks, int natoms, + int istart, int nsavc, double delta, int with_unitcell, + int charmm); + +/* + * Write a timestep to a dcd file + * Input: fd - a file struct for which a dcd header has already been written + * curframe: Count of frames written to this file, starting with 1. + * curstep: Count of timesteps elapsed = istart + curframe * nsavc. + * natoms - number of elements in x, y, z arrays + * x, y, z: pointers to atom coordinates + * 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 curframe, int curstep, + int natoms, const float *x, const float *y, const float *z, + const double *unitcell, int charmm); + + + +#define DCD_IS_CHARMM 0x01 +#define DCD_HAS_4DIMS 0x02 +#define DCD_HAS_EXTRA_BLOCK 0x04 + +/* READ Macro to make porting easier */ +#define READ(fd, buf, size) \ + fio_fread(((void *) buf), (size), 1, (fd)) + + +/* WRITE Macro to make porting easier */ +#define WRITE(fd, buf, size) \ + fio_fwrite(((void *) buf), (size), 1, (fd)) + +/* XXX This is broken - fread never returns -1 */ +#define CHECK_FREAD(X, msg) if (X==-1) \ + { \ + return(DCD_BADREAD); \ + } + +#define CHECK_FEOF(X, msg) if (X==0) \ + { \ + return(DCD_BADEOF); \ + } + +static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART, + int *NSAVC, double *DELTA, int *NAMNF, + int **FREEINDEXES, float **fixedcoords, int *reverseEndian, + int *charmm, char **remarks, int *len_remarks) +{ + int input_integer; /* buffer space */ + int ret_val; + char hdrbuf[84]; /* char buffer used to store header */ + int NTITLE; + + /* First thing in the file should be an 84 */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading first int from dcd file"); + CHECK_FEOF(ret_val, "reading first int from dcd file"); + + /* Check magic number in file header and determine byte order*/ + if (input_integer != 84) { + /* check to see if its merely reversed endianism */ + /* rather than a totally incorrect file magic number */ + swap4_aligned(&input_integer, 1); + + if (input_integer == 84) { + *reverseEndian=1; + } else { + /* not simply reversed endianism, but something rather more evil */ + return DCD_BADFORMAT; + } + } else { + *reverseEndian=0; + } + + /* Buffer the entire header for random access */ + ret_val = READ(fd, hdrbuf, 84); + CHECK_FREAD(ret_val, "buffering header"); + CHECK_FEOF(ret_val, "buffering header"); + + /* Check for the ID string "COORD" */ + if (hdrbuf[0] != 'C' || hdrbuf[1] != 'O' || + hdrbuf[2] != 'R' || hdrbuf[3] != 'D') { + return DCD_BADFORMAT; + } + + /* CHARMm-genereate DCD files set the last integer in the */ + /* header, which is unused by X-PLOR, to its version number. */ + /* Checking if this is nonzero tells us this is a CHARMm file */ + /* and to look for other CHARMm flags. */ + if (*((int *) (hdrbuf + 80)) != 0) { + (*charmm) = DCD_IS_CHARMM; + if (*((int *) (hdrbuf + 44)) != 0) + (*charmm) |= DCD_HAS_EXTRA_BLOCK; + + if (*((int *) (hdrbuf + 48)) == 1) + (*charmm) |= DCD_HAS_4DIMS; + } else { + (*charmm) = 0; + } + + /* Store the number of sets of coordinates (NSET) */ + (*NSET) = *((int *) (hdrbuf + 4)); + if (*reverseEndian) swap4_unaligned(NSET, 1); + + /* Store ISTART, the starting timestep */ + (*ISTART) = *((int *) (hdrbuf + 8)); + if (*reverseEndian) swap4_unaligned(ISTART, 1); + + /* Store NSAVC, the number of timesteps between dcd saves */ + (*NSAVC) = *((int *) (hdrbuf + 12)); + if (*reverseEndian) swap4_unaligned(NSAVC, 1); + + /* Store NAMNF, the number of fixed atoms */ + (*NAMNF) = *((int *) (hdrbuf + 36)); + if (*reverseEndian) swap4_unaligned(NAMNF, 1); + + /* Read in the timestep, DELTA */ + /* Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm */ + if ((*charmm) & DCD_IS_CHARMM) { + float ftmp; + ftmp = *((float *)(hdrbuf+40)); /* is this safe on Alpha? */ + if (*reverseEndian) + swap4_aligned(&ftmp, 1); + + *DELTA = (double)ftmp; + } else { + (*DELTA) = *((double *)(hdrbuf + 40)); + if (*reverseEndian) swap8_unaligned(DELTA, 1); + } + + /* Get the end size of the first block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading second 84 from dcd file"); + CHECK_FEOF(ret_val, "reading second 84 from dcd file"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 84) { + return DCD_BADFORMAT; + } + + /* Read in the size of the next block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of title block"); + CHECK_FEOF(ret_val, "reading size of title block"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (((input_integer-4) % 80) == 0) { + /* Read NTITLE, the number of 80 character title strings there are */ + ret_val = READ(fd, &NTITLE, sizeof(int)); + CHECK_FREAD(ret_val, "reading NTITLE"); + CHECK_FEOF(ret_val, "reading NTITLE"); + if (*reverseEndian) swap4_aligned(&NTITLE, 1); + *len_remarks = NTITLE*80; + *remarks = (char*)malloc(*len_remarks); + ret_val = fio_fread(*remarks, *len_remarks, 1, fd); + CHECK_FEOF(ret_val, "reading TITLE"); + + /* Get the ending size for this block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + + CHECK_FREAD(ret_val, "reading size of title block"); + CHECK_FEOF(ret_val, "reading size of title block"); + } else { + return DCD_BADFORMAT; + } + + /* Read in an integer '4' */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading a '4'"); + CHECK_FEOF(ret_val, "reading a '4'"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 4) { + return DCD_BADFORMAT; + } + + /* Read in the number of atoms */ + ret_val = READ(fd, N, sizeof(int)); + CHECK_FREAD(ret_val, "reading number of atoms"); + CHECK_FEOF(ret_val, "reading number of atoms"); + if (*reverseEndian) swap4_aligned(N, 1); + + /* Read in an integer '4' */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading a '4'"); + CHECK_FEOF(ret_val, "reading a '4'"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 4) { + return DCD_BADFORMAT; + } + + *FREEINDEXES = NULL; + *fixedcoords = NULL; + if (*NAMNF != 0) { + (*FREEINDEXES) = (int *) calloc(((*N)-(*NAMNF)), sizeof(int)); + if (*FREEINDEXES == NULL) + return DCD_BADMALLOC; + + *fixedcoords = (float *) calloc((*N)*4 - (*NAMNF), sizeof(float)); + if (*fixedcoords == NULL) + return DCD_BADMALLOC; + + /* Read in index array size */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != ((*N)-(*NAMNF))*4) { + return DCD_BADFORMAT; + } + + ret_val = READ(fd, (*FREEINDEXES), ((*N)-(*NAMNF))*sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + + if (*reverseEndian) + swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF))); + + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != ((*N)-(*NAMNF))*4) { + return DCD_BADFORMAT; + } + } + + return DCD_SUCCESS; +} + +static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian, + double *unitcell) { + int i, input_integer; + + if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) { + /* Leading integer must be 48 */ + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) + return DCD_BADREAD; + if (reverseEndian) swap4_aligned(&input_integer, 1); + if (input_integer == 48) { + double tmp[6]; + 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] = tmp[i]; + } else { + /* unrecognized block, just skip it */ + if (fio_fseek(fd, input_integer, FIO_SEEK_CUR)) return DCD_BADREAD; + } + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD; + } + + return DCD_SUCCESS; +} + +static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes, + int reverseEndian, const float *fixedcoords, + float *freeatoms, float *pos) { + int i, input_integer; + + /* Read leading integer */ + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD; + if (reverseEndian) swap4_aligned(&input_integer, 1); + if (input_integer != 4*num_free) return DCD_BADFORMAT; + + /* Read free atom coordinates */ + if (fio_fread(freeatoms, 4*num_free, 1, fd) != 1) return DCD_BADREAD; + if (reverseEndian) + swap4_aligned(freeatoms, num_free); + + /* Copy fixed and free atom coordinates into position buffer */ + memcpy(pos, fixedcoords, 4*N); + for (i=0; i 1) { + seekoffset *= numsteps; + } + + rc = fio_fseek(fd, seekoffset, FIO_SEEK_CUR); + if (rc == -1) return DCD_BADEOF; + + return DCD_SUCCESS; +} + +static int jump_to_dcdstep(fio_fd fd, int natoms, int nsets, int nfixed, int charmm, int header_size, int step) { + int rc; + if (step > nsets) { + return DCD_BADEOF; + } + // Calculate file offset + off_t extrablocksize, ndims, firstframesize, framesize; + off_t pos; + extrablocksize = charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; + ndims = charmm & DCD_HAS_4DIMS ? 4 : 3; + firstframesize = (natoms+2) * ndims * sizeof(float) + extrablocksize; + framesize = (natoms-nfixed+2) * ndims * sizeof(float) + extrablocksize; + // Use zero indexing + if (step == 0) { + pos = header_size; + } + else { + pos = header_size + firstframesize + framesize * (step-1); + } + rc = fio_fseek(fd, pos, FIO_SEEK_SET); + if (rc == -1) return DCD_BADEOF; + return DCD_SUCCESS; +} + +#define NFILE_POS 8L +#define NSTEP_POS 20L + +static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N, + const float *X, const float *Y, const float *Z, + const double *unitcell, int charmm) { + int out_integer; + + if (charmm) { + /* write out optional unit cell */ + if (unitcell != NULL) { + out_integer = 48; /* 48 bytes (6 doubles) */ + fio_write_int32(fd, out_integer); + WRITE(fd, unitcell, out_integer); + fio_write_int32(fd, out_integer); + } + } + + /* write out coordinates */ + out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */ + fio_write_int32(fd, out_integer); + WRITE(fd, X, out_integer); + fio_write_int32(fd, out_integer); + fio_write_int32(fd, out_integer); + WRITE(fd, Y, out_integer); + fio_write_int32(fd, out_integer); + fio_write_int32(fd, out_integer); + WRITE(fd, Z, out_integer); + fio_write_int32(fd, out_integer); + + /* update the DCD header information */ + fio_fseek(fd, NFILE_POS, FIO_SEEK_SET); + fio_write_int32(fd, curframe); + fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET); + fio_write_int32(fd, curstep); + fio_fseek(fd, 0, FIO_SEEK_END); + + return DCD_SUCCESS; +} + +static int write_dcdheader(fio_fd fd, const char *remarks, int N, + int ISTART, int NSAVC, double DELTA, int with_unitcell, + int charmm) { + int out_integer; + float out_float; + char title_string[241]; + time_t cur_time; + struct tm *tmbuf; + + out_integer = 84; + WRITE(fd, (char *) & out_integer, sizeof(int)); + strcpy(title_string, "CORD"); + WRITE(fd, title_string, 4); + fio_write_int32(fd, 0); /* Number of frames in file, none written yet */ + fio_write_int32(fd, ISTART); /* Starting timestep */ + fio_write_int32(fd, NSAVC); /* Timesteps between frames written to the file */ + fio_write_int32(fd, 0); /* Number of timesteps in simulation */ + fio_write_int32(fd, 0); /* NAMD writes NSTEP or ISTART - NSAVC here? */ + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + if (charmm) { + out_float = DELTA; + WRITE(fd, (char *) &out_float, sizeof(float)); + if (with_unitcell) { + fio_write_int32(fd, 1); + } else { + fio_write_int32(fd, 0); + } + } else { + WRITE(fd, (char *) &DELTA, sizeof(double)); + } + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + if (charmm) { + fio_write_int32(fd, 24); /* Pretend to be Charmm version 24 */ + } else { + fio_write_int32(fd, 0); + } + fio_write_int32(fd, 84); + fio_write_int32(fd, 164); + fio_write_int32(fd, 3); /* the number of 80 character title strings */ + + strncpy(title_string, remarks, 240); + WRITE(fd, title_string, 240); + + fio_write_int32(fd, 164); + fio_write_int32(fd, 4); + fio_write_int32(fd, N); + fio_write_int32(fd, 4); + + return DCD_SUCCESS; +} + +static void close_dcd_read(int *indexes, float *fixedcoords) { + free(indexes); + free(fixedcoords); +} + +#endif + diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx new file mode 100644 index 00000000000..32e753a85fb --- /dev/null +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -0,0 +1,485 @@ +# -*- 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 os import path +import numpy as np +from collections import namedtuple +from MDAnalysis.lib.mdamath import triclinic_box +import six +import string +import sys + +cimport numpy as np + + +from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END +from libc.math cimport M_PI_2, asin + +_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 + +from libc.stdint cimport uintptr_t +from libc.stdlib cimport free + +cdef enum: + FIO_READ = 0x01 + FIO_WRITE = 0x02 + +cdef enum: + DCD_IS_CHARMM = 0x01 + DCD_HAS_4DIMS = 0x02 + DCD_HAS_EXTRA_BLOCK = 0x04 + +DCD_ERRORS = { + 0: 'No Problem', + -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 exiss', + -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 *n_atoms, 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 n_atoms, 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 n_atoms, 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', 'x unitcell') + +cdef class DCDFile: + cdef fio_fd fp + cdef readonly fname + cdef int is_open + cdef readonly int n_atoms + cdef int nsets + cdef readonly int istart + cdef readonly int nsavc + cdef readonly double delta + cdef readonly int nfixed + cdef int *freeind + cdef float *fixedcoords + cdef int reverse_endian + cdef int charmm + cdef str mode + cdef readonly int n_dims + cdef readonly int n_frames + cdef bint b_read_header + cdef int current_frame + cdef readonly remarks + cdef int reached_eof + cdef readonly int _firstframesize + cdef readonly int _framesize + cdef readonly int _header_size + + def __cinit__(self, fname, mode='r'): + self.fname = fname.encode('utf-8') + self.n_atoms = 0 + self.is_open = False + self.open(self.fname, mode) + + def __dealloc__(self): + self.close() + if self.fixedcoords != NULL: + free(self.fixedcoords) + if self.freeind != NULL: + free(self.freeind) + + def __enter__(self): + """Support context manager""" + if not self.is_open: + self.open(self.fname, 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.fname, 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): + return self.current_frame + + def open(self, filename, mode='r'): + 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, ok)) + self.is_open = True + self.current_frame = 0 + self.reached_eof = False + # Has to come last since it checks the reached_eof flag + if self.mode == 'r': + self.remarks = self._read_header() + + def close(self): + 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, ok)) + + + def _read_header(self): + if not self.is_open: + raise IOError("No file open") + + cdef char* c_remarks + cdef int len_remarks = 0 + + ok = read_dcdheader(self.fp, &self.n_atoms, &self.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])) + + if c_remarks != NULL: + py_remarks = c_remarks[:len_remarks] + free(c_remarks) + else: + py_remarks = "" + + self.n_dims = 3 if not self.charmm & DCD_HAS_4DIMS else 4 + 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: + # if this fails the file is empty. Set flag and warn using during read + 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) + + return py_remarks + + @property + def header(self): + return {'n_atoms': self.n_atoms, + 'nsets': self.nsets, + 'istart': self.istart, + 'nsavc': self.nsacv, + 'delta': self.delta, + 'nfixed': self.nfixed, + 'reverse_endian': self.reverse_endian, + 'charmm': self.charmm, + 'remarks': self.remarks} + + def _estimate_n_frames(self): + extrablocksize = 48 + 8 if self.charmm & DCD_HAS_EXTRA_BLOCK else 0 + self._firstframesize = (self.n_atoms + 2) * self.n_dims * sizeof(float) + extrablocksize + self._framesize = ((self.n_atoms - self.nfixed + 2) * self.n_dims * 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) + # TODO: check that nframessize is larger then 0, the c-implementation + # used to do that. + nframessize = filesize - self._header_size - self._firstframesize + return nframessize / self._framesize + 1 + + @property + def periodic(self): + return bool((self.charmm & DCD_IS_CHARMM) and + (self.charmm & DCD_HAS_EXTRA_BLOCK)) + + + + def seek(self, frame): + """Seek to Frame. + + Please note that this function will generate internal file offsets if + they haven't been set before. For large file this means the first seek + can be very slow. Later seeks will be very fast. + + Parameters + ---------- + frame : int + seek the file to given frame + + Raises + ------ + IOError + If you seek for more frames than are available or if the + seek fails (the low-level system error is reported). + """ + if frame >= self.n_frames: + raise IOError('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 system errno={}".format(ok)) + self.current_frame = frame + + def _write_header(self, remarks, int n_atoms, int starting_step, + int ts_between_saves, double time_step, + int charmm): + + if not self.is_open: + raise IOError("No file open") + + if not self.mode=='w': + raise IOError("Incorrect file mode for writing.") + + #cdef char c_remarks + cdef int len_remarks = 0 + cdef int with_unitcell = 1 + + + if isinstance(remarks, six.string_types): + try: + remarks = bytearray(remarks, 'ascii') + except UnicodeDecodeError: + remarks = bytearray(remarks) + + ok = write_dcdheader(self.fp, remarks, n_atoms, starting_step, + ts_between_saves, time_step, with_unitcell, + charmm) + if ok != 0: + raise IOError("Writing DCD header failed: {}".format(DCD_ERRORS[ok])) + + def write(self, xyz, box, int step, int natoms, + int ts_between_saves, int charmm, double time_step, remarks): + """write one frame into DCD file. + + Parameters + ---------- + xyz : array_like, shape=(n_atoms, 3) + cartesion coordinates + box : array_like, shape=(6) + Box vectors for this frame + step : int + current step number, 1 indexed + time : float + current time + natoms : int + number of atoms in frame + + Raises + ------ + IOError + """ + 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 len(box) != 6: + raise ValueError("box size is wrong should be 6, got: {}".format(box.size)) + + # print("box to write ", box) + + #cdef double [:,:] unitcell = box + xyz = np.asarray(xyz, order='F', dtype=np.float32) + + # 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(box, _ts_order) + # print("writting box ", box) + cdef DOUBLE_T[::1] c_box = np.asarray(box, order='C', dtype=np.float64) + + if xyz.shape != (natoms, 3): + raise ValueError("xyz shape is wrong should be (natoms, 3), got:".format(xyz.shape)) + + cdef FLOAT_T[::1] x = xyz[:, 0] + cdef FLOAT_T[::1] y = xyz[:, 1] + cdef FLOAT_T[::1] z = xyz[:, 2] + cdef float alpha, beta, gamma, a, b, c; + + if self.current_frame == 0: + self._write_header(remarks=remarks, n_atoms=xyz.shape[0], starting_step=step, + ts_between_saves=ts_between_saves, + time_step=time_step, + charmm=charmm) + self.n_atoms = xyz.shape[0] + + # looks like self.nsavc is just 0 all the time + step = self.current_frame * self.nsavc + ok = write_dcdstep(self.fp, self.current_frame, step, + self.n_atoms, &x[0], + &y[0], &z[0], + &c_box[0], charmm) + + self.current_frame += 1 + + def read(self): + if self.reached_eof: + raise IOError('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.n_atoms, 3), 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; + + cdef FLOAT_T[::1] x = xyz[:, 0] + cdef FLOAT_T[::1] y = xyz[:, 1] + cdef FLOAT_T[::1] z = xyz[:, 2] + + first_frame = self.current_frame == 0 + ok = read_dcdstep(self.fp, self.n_atoms, + &x[0], + &y[0], &z[0], + unitcell.data, self.nfixed, first_frame, + self.freeind, self.fixedcoords, + self.reverse_endian, self.charmm) + 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 + + if (unitcell[1] >= -1.0 and unitcell[1] <= 1.0 and + unitcell[3] >= -1.0 and unitcell[3] <= 1.0 and + unitcell[4] >= -1.0 and 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; + + # The original unitcell is read as ``[A, gamma, B, beta, alpha, C]`` + _ts_order = [0, 2, 5, 4, 3, 1] + uc = np.take(unitcell, _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 = unitcell + e1, e2, e3 = H[[0,1,3]], H[[1,2,4]], H[[3,4,5]] + uc = triclinic_box(e1, e2, e3) + + unitcell = uc + + return DCDFrame(xyz, unitcell) diff --git a/package/setup.py b/package/setup.py index 961b82004a8..fe79aef5398 100755 --- a/package/setup.py +++ b/package/setup.py @@ -294,11 +294,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) dcd_time = MDAExtension('coordinates.dcdtimeseries', ['MDAnalysis/coordinates/dcdtimeseries' + source_suffix], include_dirs=include_dirs + ['MDAnalysis/coordinates/include'], @@ -356,10 +356,9 @@ 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, dcd_time, distances, distances_omp, qcprot, - transformation, libmdaxdr, util, encore_utils, - ap_clustering, spe_dimred] - + pre_exts = [libdcd, dcd_time, distances, distances_omp, qcprot, + transformation, libmdaxdr, util, encore_utils, + ap_clustering, spe_dimred] cython_generated = [] if use_cython: extensions = cythonize(pre_exts) diff --git a/testsuite/MDAnalysisTests/coordinates/base.py b/testsuite/MDAnalysisTests/coordinates/base.py index 2adc80231df..b482dd19407 100644 --- a/testsuite/MDAnalysisTests/coordinates/base.py +++ b/testsuite/MDAnalysisTests/coordinates/base.py @@ -242,13 +242,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) def test_first_dimensions(self): self.reader.rewind() @@ -1067,7 +1067,8 @@ 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 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/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index baf9d5989fc..e00b559ebad 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -21,6 +21,7 @@ # from __future__ import absolute_import import MDAnalysis as mda +from MDAnalysis.coordinates.DCD import DCDReader import numpy as np import os from six.moves import zip, range @@ -31,23 +32,44 @@ assert_allclose, dec) from unittest import TestCase -from MDAnalysisTests.datafiles import (DCD, PSF, DCD_empty, CRD, PRMncdf, NCDF) +from MDAnalysisTests.datafiles import (DCD, PSF, DCD_empty, CRD, PRMncdf, NCDF, + COORDINATES_TOPOLOGY, COORDINATES_DCD) from MDAnalysisTests.coordinates.reference import (RefCHARMMtriclinicDCD, RefNAMDtriclinicDCD) -from MDAnalysisTests.coordinates.base import BaseTimestepTest +# from MDAnalysisTests.coordinates.base import BaseTimestepTest + +from MDAnalysisTests.coordinates.base import (MultiframeReaderTest, BaseReference, + BaseWriterTest, + assert_timestep_almost_equal) + from MDAnalysisTests import module_not_found, tempdir -from MDAnalysisTests.plugins.knownfailure import knownfailure +# 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,7 +87,23 @@ def test_with_statement(self): err_msg="with_statement: DCDReader does not read all frames") -class TestDCDReader(object): +class TestDCDWriter(BaseWriterTest): + def __init__(self, reference=None): + if reference is None: + reference = DCDReference() + super(TestDCDWriter, self).__init__(reference) + + +################ +# Legacy tests # +################ + + +class TestDCDReaderClass(TestCase): + pass + + +class TestDCDReader_old(object): def setUp(self): self.universe = mda.Universe(PSF, DCD) self.dcd = self.universe.trajectory @@ -85,10 +123,6 @@ def test_next_dcd(self): 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]") @@ -118,70 +152,31 @@ def test_reverse_dcd(self): 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_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): @@ -193,7 +188,7 @@ def test_DCDReader_set_dt(dt=100., frame=3): assert_almost_equal(u.trajectory.dt, dt, err_msg="trajectory.dt does not match set dt") -class TestDCDWriter(TestCase): +class TestDCDWriter_old(object): def setUp(self): self.universe = mda.Universe(PSF, DCD) ext = ".dcd" @@ -210,26 +205,6 @@ def tearDown(self): 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: - 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 @@ -277,24 +252,8 @@ def test_single_frame(self): 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): +class TestDCDWriter_Issue59(object): def setUp(self): """Generate input xtc.""" self.u = mda.Universe(PSF, DCD) @@ -365,7 +324,7 @@ def test_OtherWriter(self): dcd.trajectory.frame)) -class _TestDCDReader_TriclinicUnitcell(TestCase): +class _TestDCDReader_TriclinicUnitcell(object): def setUp(self): self.u = mda.Universe(self.topology, self.trajectory) self.tempdir = tempdir.TempDir() @@ -394,10 +353,12 @@ def test_read_triclinic(self): def test_write_triclinic(self): """test writing of triclinic unitcell (Issue 187) for NAMD or new CHARMM format (at least since c36b2)""" + print("writing") with self.u.trajectory.OtherWriter(self.dcd) as w: for ts in self.u.trajectory: w.write(ts) w = mda.Universe(self.topology, self.dcd) + print("reading\n") for ts_orig, ts_copy in zip(self.u.trajectory, w.trajectory): assert_almost_equal(ts_orig.dimensions, ts_copy.dimensions, 4, @@ -416,7 +377,7 @@ class TestDCDReader_NAMD_Unitcell(_TestDCDReader_TriclinicUnitcell, pass -class TestNCDF2DCD(TestCase): +class TestNCDF2DCD(object): @dec.skipif(module_not_found("netCDF4"), "Test skipped because netCDF is not available.") def setUp(self): @@ -463,167 +424,172 @@ def test_coordinates(self): ts_orig.frame)) -class TestDCDCorrel(TestCase): - def setUp(self): - # Note: setUp is executed for *every* test ! - import MDAnalysis.core.Timeseries as TS - self.universe = mda.Universe(PSF, DCD) - self.dcd = self.universe.trajectory - self.ts = self.universe.coord - self.collection = TS.TimeseriesCollection() - self.collection_slicing = TS.TimeseriesCollection() - C = self.collection - C_step = self.collection_slicing - all = self.universe.atoms - ca = self.universe.s4AKE.atoms.CA - ca_termini = ca[[0, -1]] - # note that this is not quite phi... HN should be C of prec. residue - phi151 = self.universe.select_atoms('resid 151').select_atoms( - 'name HN', 'name N', 'name CA', 'name CB') - C.addTimeseries(TS.Atom('v', ca_termini)) # 0 - C.addTimeseries(TS.Bond(ca_termini)) # 1 - C.addTimeseries(TS.Bond([ca[0], ca[-1]])) # 2 - C.addTimeseries(TS.Angle(phi151[1:4])) # 3 - C.addTimeseries(TS.Dihedral(phi151)) # 4 - C.addTimeseries(TS.Distance('r', ca_termini)) # 5 - C.addTimeseries(TS.CenterOfMass(ca)) # 6 - C.addTimeseries(TS.CenterOfGeometry(ca)) # 7 - C.addTimeseries(TS.CenterOfMass(all)) # 8 - C.addTimeseries(TS.CenterOfGeometry(all)) # 9 - - C_step.addTimeseries(TS.Atom('v', ca_termini)) - - # cannot test WaterDipole because there's no water in the test dcd - C.compute(self.dcd) - C_step.compute(self.dcd, step=10) - - def tearDown(self): - del self.collection - del self.collection_slicing - del self.universe - del self.dcd - del self.ts - - def test_correl(self): - assert_equal(len(self.collection), 10, "Correl: len(collection)") - - def test_Atom(self): - assert_equal(self.collection[0].shape, (2, 3, 98), - "Correl: Atom positions") - - def test_Bonds(self): - C = self.collection - assert_array_equal(C[1].__data__, C[2].__data__, - "Correl: Bonds with lists and AtomGroup") - - def test_Angle(self): - C = self.collection - avg_angle = 1.9111695972912988 - assert_almost_equal(C[3].__data__.mean(), - avg_angle, - err_msg="Correl: average Angle") - - def test_Dihedral(self): - C = self.collection - avg_phi151 = 0.0088003870749735619 - assert_almost_equal(C[4].__data__.mean(), - avg_phi151, - err_msg="Correl: average Dihedral") - - def test_scalarDistance(self): - C = self.collection - avg_dist = 9.7960210987736236 - assert_almost_equal(C[5].__data__.mean(), - avg_dist, - err_msg="Correl: average scalar Distance") - - def test_CenterOfMass(self): - C = self.collection - avg_com_ca = np.array([0.0043688, -0.27812258, 0.0284051]) - avg_com_all = np.array([-0.10086529, -0.16357276, 0.12724672]) - assert_array_almost_equal(C[6].__data__.mean(axis=1), - avg_com_ca, - err_msg="Correl: average CA CenterOfMass") - assert_almost_equal(C[8].__data__.mean(axis=1), - avg_com_all, - err_msg="Correl: average all CenterOfMass") - - def test_CenterOfGeometry(self): - C = self.collection - avg_cog_all = np.array([-0.13554797, -0.20521885, 0.2118998]) - assert_almost_equal(C[9].__data__.mean(axis=1), - avg_cog_all, - err_msg="Correl: average all CenterOfGeometry") - - def test_CA_COMeqCOG(self): - C = self.collection - assert_array_almost_equal( - C[6].__data__, - C[7].__data__, - err_msg="Correl: CA CentreOfMass == CenterOfGeometry") - - def test_clear(self): - C = self.collection - C.clear() - assert_equal(len(C), 0, "Correl: clear()") - - def test_Atom_slicing(self): - assert_equal(self.collection_slicing[0].shape, (2, 3, 10), - "Correl: Atom positions") - assert_array_almost_equal(self.collection[0][:, :, ::10], - self.collection_slicing[0]) - -# notes: -def compute_correl_references(): - universe = mda.Universe(PSF, DCD) - - all = universe.atoms - ca = universe.s4AKE.CA - ca_termini = mda.core.AtomGroup.AtomGroup([ca[0], ca[-1]]) - phi151 = universe.select_atoms('resid 151').select_atoms( - 'name HN', 'name N', 'name CA', 'name CB') - - C = mda.collection - C.clear() - - C.addTimeseries(TS.Atom('v', ca_termini)) # 0 - C.addTimeseries(TS.Bond(ca_termini)) # 1 - C.addTimeseries(TS.Bond([ca[0], ca[-1]])) # 2 - C.addTimeseries(TS.Angle(phi151[1:4])) # 3 - C.addTimeseries(TS.Dihedral(phi151)) # 4 - C.addTimeseries(TS.Distance('r', ca_termini)) # 5 - C.addTimeseries(TS.CenterOfMass(ca)) # 6 - C.addTimeseries(TS.CenterOfGeometry(ca)) # 7 - C.addTimeseries(TS.CenterOfMass(all)) # 8 - C.addTimeseries(TS.CenterOfGeometry(all)) # 9 - - C.compute(universe.dcd) - - results = { - "avg_angle": C[3].__data__.mean(), - "avg_phi151": C[4].__data__.mean(), - "avg_dist": C[5].__data__.mean(), - "avg_com_ca": C[6].__data__.mean(axis=1), - "avg_com_all": C[8].__data__.mean(axis=1), - "avg_cog_all": C[9].__data__.mean(axis=1), - } - C.clear() - return results - - -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) +### +# Correl / Timeseries tests +### + + +# class TestDCDCorrel(TestCase): +# def setUp(self): +# # Note: setUp is executed for *every* test ! +# import MDAnalysis.core.Timeseries as TS +# self.universe = mda.Universe(PSF, DCD) +# self.dcd = self.universe.trajectory +# self.ts = self.universe.coord +# self.collection = TS.TimeseriesCollection() +# self.collection_slicing = TS.TimeseriesCollection() +# C = self.collection +# C_step = self.collection_slicing +# all = self.universe.atoms +# ca = self.universe.s4AKE.atoms.CA +# ca_termini = ca[[0, -1]] +# # note that this is not quite phi... HN should be C of prec. residue +# phi151 = self.universe.select_atoms('resid 151').select_atoms( +# 'name HN', 'name N', 'name CA', 'name CB') +# C.addTimeseries(TS.Atom('v', ca_termini)) # 0 +# C.addTimeseries(TS.Bond(ca_termini)) # 1 +# C.addTimeseries(TS.Bond([ca[0], ca[-1]])) # 2 +# C.addTimeseries(TS.Angle(phi151[1:4])) # 3 +# C.addTimeseries(TS.Dihedral(phi151)) # 4 +# C.addTimeseries(TS.Distance('r', ca_termini)) # 5 +# C.addTimeseries(TS.CenterOfMass(ca)) # 6 +# C.addTimeseries(TS.CenterOfGeometry(ca)) # 7 +# C.addTimeseries(TS.CenterOfMass(all)) # 8 +# C.addTimeseries(TS.CenterOfGeometry(all)) # 9 + +# C_step.addTimeseries(TS.Atom('v', ca_termini)) + +# # cannot test WaterDipole because there's no water in the test dcd +# C.compute(self.dcd) +# C_step.compute(self.dcd, step=10) + +# def tearDown(self): +# del self.collection +# del self.collection_slicing +# del self.universe +# del self.dcd +# del self.ts + +# def test_correl(self): +# assert_equal(len(self.collection), 10, "Correl: len(collection)") + +# def test_Atom(self): +# assert_equal(self.collection[0].shape, (2, 3, 98), +# "Correl: Atom positions") + +# def test_Bonds(self): +# C = self.collection +# assert_array_equal(C[1].__data__, C[2].__data__, +# "Correl: Bonds with lists and AtomGroup") + +# def test_Angle(self): +# C = self.collection +# avg_angle = 1.9111695972912988 +# assert_almost_equal(C[3].__data__.mean(), +# avg_angle, +# err_msg="Correl: average Angle") + +# def test_Dihedral(self): +# C = self.collection +# avg_phi151 = 0.0088003870749735619 +# assert_almost_equal(C[4].__data__.mean(), +# avg_phi151, +# err_msg="Correl: average Dihedral") + +# def test_scalarDistance(self): +# C = self.collection +# avg_dist = 9.7960210987736236 +# assert_almost_equal(C[5].__data__.mean(), +# avg_dist, +# err_msg="Correl: average scalar Distance") + +# def test_CenterOfMass(self): +# C = self.collection +# avg_com_ca = np.array([0.0043688, -0.27812258, 0.0284051]) +# avg_com_all = np.array([-0.10086529, -0.16357276, 0.12724672]) +# assert_array_almost_equal(C[6].__data__.mean(axis=1), +# avg_com_ca, +# err_msg="Correl: average CA CenterOfMass") +# assert_almost_equal(C[8].__data__.mean(axis=1), +# avg_com_all, +# err_msg="Correl: average all CenterOfMass") + +# def test_CenterOfGeometry(self): +# C = self.collection +# avg_cog_all = np.array([-0.13554797, -0.20521885, 0.2118998]) +# assert_almost_equal(C[9].__data__.mean(axis=1), +# avg_cog_all, +# err_msg="Correl: average all CenterOfGeometry") + +# def test_CA_COMeqCOG(self): +# C = self.collection +# assert_array_almost_equal( +# C[6].__data__, +# C[7].__data__, +# err_msg="Correl: CA CentreOfMass == CenterOfGeometry") + +# def test_clear(self): +# C = self.collection +# C.clear() +# assert_equal(len(C), 0, "Correl: clear()") + +# def test_Atom_slicing(self): +# assert_equal(self.collection_slicing[0].shape, (2, 3, 10), +# "Correl: Atom positions") +# assert_array_almost_equal(self.collection[0][:, :, ::10], +# self.collection_slicing[0]) + +# # notes: +# def compute_correl_references(): +# universe = mda.Universe(PSF, DCD) + +# all = universe.atoms +# ca = universe.s4AKE.CA +# ca_termini = mda.core.AtomGroup.AtomGroup([ca[0], ca[-1]]) +# phi151 = universe.select_atoms('resid 151').select_atoms( +# 'name HN', 'name N', 'name CA', 'name CB') + +# C = mda.collection +# C.clear() + +# C.addTimeseries(TS.Atom('v', ca_termini)) # 0 +# C.addTimeseries(TS.Bond(ca_termini)) # 1 +# C.addTimeseries(TS.Bond([ca[0], ca[-1]])) # 2 +# C.addTimeseries(TS.Angle(phi151[1:4])) # 3 +# C.addTimeseries(TS.Dihedral(phi151)) # 4 +# C.addTimeseries(TS.Distance('r', ca_termini)) # 5 +# C.addTimeseries(TS.CenterOfMass(ca)) # 6 +# C.addTimeseries(TS.CenterOfGeometry(ca)) # 7 +# C.addTimeseries(TS.CenterOfMass(all)) # 8 +# C.addTimeseries(TS.CenterOfGeometry(all)) # 9 + +# C.compute(universe.dcd) + +# results = { +# "avg_angle": C[3].__data__.mean(), +# "avg_phi151": C[4].__data__.mean(), +# "avg_dist": C[5].__data__.mean(), +# "avg_com_ca": C[6].__data__.mean(axis=1), +# "avg_com_all": C[8].__data__.mean(axis=1), +# "avg_cog_all": C[9].__data__.mean(axis=1), +# } +# C.clear() +# return results + + +# 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) 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 627b3d0d336..19eeec1c146 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -133,6 +133,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 @@ -142,12 +143,21 @@ "AUX_XVG", "XVG_BAD_NCOL", #for testing .xvg auxiliary reader "AUX_XVG_LOWF", "AUX_XVG_HIGHF", "MMTF", "MMTF_gz", + "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 "ALIGN_BOUND", # two component bound system "ALIGN_UNBOUND", # two component unbound system ] 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') @@ -164,6 +174,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..698966617b1 --- /dev/null +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -0,0 +1,577 @@ +from __future__ import print_function +from __future__ import absolute_import + +from nose.tools import raises +from numpy.testing import assert_equal +from numpy.testing import assert_allclose + +from MDAnalysis.lib.formats.libdcd import DCDFile +from MDAnalysisTests.datafiles import (PSF, DCD, DCD_NAMD_TRICLINIC, + PSF_NAMD_TRICLINIC, + legacy_DCD_ADK_coords, + legacy_DCD_NAMD_coords, + legacy_DCD_c36_coords, + PSF_TRICLINIC, + DCD_TRICLINIC) + +from MDAnalysisTests.tempdir import run_in_tempdir +from MDAnalysisTests import tempdir +import numpy as np +import os +import math +from hypothesis import given, example +import hypothesis.strategies as st +import string + + +class TestDCDReadFrame(): + + def setUp(self): + self.dcdfile = DCD + self.natoms = 3341 + self.traj_length = 98 + self.new_frame = 91 + self.context_frame = 22 + self.num_iters = 3 + self.selected_legacy_frames = [5, 29] + self.legacy_data = legacy_DCD_ADK_coords + self.expected_remarks = '''* DIMS ADK SEQUENCE FOR PORE PROGRAM * WRITTEN BY LIZ DENNING (6.2008) * DATE: 6/ 6/ 8 17:23:56 CREATED BY USER: denniej0 ''' + self.expected_unit_cell = np.array([ 0., 0., 0., 90., 90., 90.], + dtype=np.float32) + + def test_header_remarks(self): + # confirm correct header remarks section reading + with DCDFile(self.dcdfile) as f: + assert_equal(len(f.remarks), len(self.expected_remarks)) + + def test_read_coords(self): + # confirm shape of coordinate data against result from previous + # MDAnalysis implementation of DCD file handling + with DCDFile(self.dcdfile) as dcd: + dcd_frame = dcd.read() + xyz = dcd_frame[0] + assert_equal(xyz.shape, (self.natoms, 3)) + + def test_read_coord_values(self): + # 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_DCD_frame_data = np.load(self.legacy_data) + + with DCDFile(self.dcdfile) as dcd: + for index, frame_num in enumerate(self.selected_legacy_frames): + dcd.seek(frame_num) + actual_coords = dcd.read()[0] + desired_coords = legacy_DCD_frame_data[index] + assert_equal(actual_coords, + desired_coords) + + + def test_read_unit_cell(self): + # confirm unit cell read against result from previous + # MDAnalysis implementation of DCD file handling + with DCDFile(self.dcdfile) as dcd: + dcd_frame = dcd.read() + unitcell = dcd_frame[1] + assert_allclose(unitcell, self.expected_unit_cell, rtol=1e-05) + + @raises(IOError) + def test_seek_over_max(self): + # should raise IOError if beyond 98th frame + with DCDFile(DCD) as dcd: + dcd.seek(102) + + def test_seek_normal(self): + # frame seek within range is tested + with DCDFile(self.dcdfile) as dcd: + dcd.seek(self.new_frame) + assert_equal(dcd.tell(), self.new_frame) + + @raises(IOError) + def test_seek_negative(self): + # frame seek with negative number + with DCDFile(self.dcdfile) as dcd: + dcd.seek(-78) + + def test_iteration(self): + with DCDFile(self.dcdfile) as dcd: + for _ in range(self.num_iters): + dcd.__next__() + assert_equal(dcd.tell(), self.num_iters) + + def test_zero_based_frames(self): + expected_frame = 0 + with DCDFile(self.dcdfile) as dcd: + assert_equal(dcd.tell(), expected_frame) + + def test_length_traj(self): + expected = self.traj_length + with DCDFile(self.dcdfile) as dcd: + assert_equal(len(dcd), expected) + + @raises(IOError) + def test_open_wrong_mode(self): + DCDFile('foo', 'e') + + @raises(IOError) + def test_raise_not_existing(self): + DCDFile('foo') + + def test_n_atoms(self): + with DCDFile(self.dcdfile) as dcd: + assert_equal(dcd.n_atoms, self.natoms) + + @raises(IOError) + @run_in_tempdir() + def test_read_write_mode_file(self): + with DCDFile('foo', 'w') as f: + f.read() + + @raises(IOError) + def test_read_closed(self): + with DCDFile(self.dcdfile) as dcd: + dcd.close() + dcd.read() + + def test_iteration_2(self): + with DCDFile(self.dcdfile) as 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()) + + +class TestDCDWriteHeader(): + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.dcdfile = DCD + + def tearDown(self): + try: + os.unlink(self.testfile) + except OSError: + pass + del self.tmpdir + + def test_write_header_crude(self): + # test that _write_header() can produce a very crude + # header for a new / empty file + with DCDFile(self.testfile, 'w') as dcd: + dcd._write_header(remarks='Crazy!', n_atoms=22, + starting_step=12, ts_between_saves=10, + time_step=0.02, + charmm=1) + + # we're not actually asserting anything, yet + # run with: nosetests test_libdcd.py --nocapture + # to see printed output from nose + with DCDFile(self.testfile) as dcd: + assert_equal(dcd.remarks, 'Crazy!') + assert_equal(dcd.n_atoms, 22) + + @raises(IOError) + def test_write_header_only(self): + # test that _write_header() can produce a very crude + # header for a new / empty file + with DCDFile(self.testfile, 'w') as dcd: + dcd._write_header(remarks='Crazy!', n_atoms=22, + starting_step=12, ts_between_saves=10, + time_step=0.02, + charmm=1) + + # we're not actually asserting anything, yet + # run with: nosetests test_libdcd.py --nocapture + # to see printed output from nose + with DCDFile(self.testfile) as dcd: + dcd.read() + + @raises(IOError) + def test_write_header_mode_sensitivy(self): + # an exception should be raised on any attempt to use + # _write_header with a DCDFile object in 'r' mode + with DCDFile(self.dcdfile) as dcd: + dcd._write_header(remarks='Crazy!', n_atoms=22, + starting_step=12, ts_between_saves=10, + time_step=0.02, + charmm=1) + + + +class TestDCDWrite(): + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.testfile2 = self.tmpdir.name + '/test2.dcd' + self.readfile = DCD + self.natoms = 3341 + self.expected_frames = 98 + self.seek_frame = 91 + self.expected_remarks = '''* DIMS ADK SEQUENCE FOR PORE PROGRAM * WRITTEN BY LIZ DENNING (6.2008) * DATE: 6/ 6/ 8 17:23:56 CREATED BY USER: denniej0 ''' + self._write_files(testfile=self.testfile, + remarks_setting='input') + + def _write_files(self, testfile, remarks_setting): + + with DCDFile(self.readfile) as f_in, DCDFile(testfile, 'w') as f_out: + for frame in f_in: + if remarks_setting == 'input': + remarks = f_in.remarks + else: # accept the random remarks strings from hypothesis + remarks = remarks_setting + box=frame.unitcell.astype(np.float64) + f_out.write(xyz=frame.x, + box=box, + step=f_in.istart, + natoms=frame.x.shape[0], + charmm=1, # DCD should be CHARMM + time_step=f_in.delta, + ts_between_saves=f_in.nsavc, + remarks=remarks) + + def tearDown(self): + try: + os.unlink(self.testfile) + except OSError: + pass + del self.tmpdir + + @raises(IOError) + def test_write_mode(self): + # ensure that writing of DCD files only occurs with properly + # opened files + with DCDFile(self.readfile) as dcd: + dcd.write(xyz=np.zeros((3,3)), + box=np.zeros(6, dtype=np.float64), + step=0, + natoms=330, + charmm=0, + time_step=22.2, + ts_between_saves=3, + remarks='') + + def test_written_dcd_coordinate_data_shape(self): + # written coord shape should match for all frames + expected = (self.natoms, 3) + with DCDFile(self.testfile) as f: + if f.n_frames > 1: + for frame in f: + xyz = f.read()[0] + assert_equal(xyz.shape, expected) + else: + xyz = f.read()[0] + assert_equal(xyz.shape, expected) + + def test_written_unit_cell(self): + # written unit cell dimensions should match for all frames + with DCDFile(self.testfile) as test, DCDFile(self.readfile) as ref: + curr_frame = 0 + while curr_frame < test.n_frames: + written_unitcell = test.read().unitcell + ref_unitcell = ref.read().unitcell + curr_frame += 1 + assert_equal(written_unitcell, ref_unitcell) + + def test_written_num_frames(self): + with DCDFile(self.testfile) as f: + assert_equal(len(f), self.expected_frames) + + def test_written_seek(self): + # ensure that we can seek properly on written DCD file + with DCDFile(self.testfile) as f: + f.seek(self.seek_frame) + assert_equal(f.tell(), self.seek_frame) + + def test_written_zero_based_frames(self): + # ensure that the first written DCD frame is 0 + expected_frame = 0 + with DCDFile(self.testfile) as f: + assert_equal(f.tell(), expected_frame) + + def test_written_remarks(self): + # ensure that the REMARKS field *can be* preserved exactly + # in the written DCD file + with DCDFile(self.testfile) as f: + assert_equal(f.remarks, self.expected_remarks) + + @given(st.text(alphabet=string.printable, + min_size=0, + max_size=500)) # handle the printable ASCII strings + @example('') + def test_written_remarks_property(self, remarks_str): + # property based testing for writing of a wide range of string + # values to REMARKS field + self._write_files(testfile=self.testfile2, + remarks_setting=remarks_str) + expected_remarks = remarks_str[:240] + with DCDFile(self.testfile2) as f: + assert_equal(f.remarks, expected_remarks) + + def test_written_nsavc(self): + # ensure that nsavc, the timesteps between frames written + # to file, is preserved in the written DCD file + with DCDFile(self.readfile) as dcd_r, DCDFile(self.testfile) as dcd: + assert_equal(dcd.nsavc, dcd_r.nsavc) + + def test_written_istart(self): + # ensure that istart, the starting timestep, is preserved + # in the written DCD file + with DCDFile(self.readfile) as dcd_r, DCDFile(self.testfile) as dcd: + assert_equal(dcd.istart, dcd_r.istart) + + def test_written_delta(self): + # ensure that delta, the trajectory timestep, is preserved in + # the written DCD file + with DCDFile(self.readfile) as dcd_r, DCDFile(self.testfile) as dcd: + assert_equal(dcd.delta, dcd_r.delta) + + def test_coord_match(self): + # ensure that all coordinates match in each frame for the + # written DCD file relative to original + with DCDFile(self.testfile) as test, DCDFile(self.readfile) as ref: + curr_frame = 0 + while curr_frame < test.n_frames: + written_coords = test.read()[0] + ref_coords = ref.read()[0] + curr_frame += 1 + assert_equal(written_coords, ref_coords) + + def test_write_wrong_dtype(self): + """we should allow passing a range of dtypes""" + for dtype in (np.int32, np.int64, np.float32, np.float64): + with DCDFile(self.testfile, 'w') as out: + natoms = 10 + xyz = np.ones((natoms, 3), dtype=dtype) + box = np.ones(6, dtype=dtype) + out.write(xyz=xyz, box=box, step=1, natoms=natoms, charmm=1, time_step=0, + ts_between_saves=1, remarks='test') + + def test_write_array_like(self): + """we should allow passing a range of dtypes""" + for array_like in (np.array, list): + with DCDFile(self.testfile, '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(xyz=xyz, box=box, step=1, natoms=natoms, charmm=1, time_step=0, + ts_between_saves=1, remarks='test') + + @raises(ValueError) + def test_write_wrong_shape_xyz(self): + with DCDFile(self.testfile, 'w') as out: + natoms = 10 + xyz = np.ones((natoms+1, 3)) + box = np.ones(6) + out.write(xyz=xyz, box=box, step=1, natoms=natoms, charmm=1, time_step=0, + ts_between_saves=1, remarks='test') + + @raises(ValueError) + def test_write_wrong_shape_box(self): + with DCDFile(self.testfile, 'w') as out: + natoms = 10 + xyz = np.ones((natoms, 3)) + box = np.ones(8) + out.write(xyz=xyz, box=box, step=1, natoms=natoms, charmm=1, time_step=0, + ts_between_saves=1, remarks='test') + + +class TestDCDWriteNAMD(TestDCDWrite): + # repeat writing tests for NAMD format DCD + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.testfile2 = self.tmpdir.name + '/test2.dcd' + self.readfile = DCD_NAMD_TRICLINIC + self.natoms = 5545 + self.expected_frames = 1 + self.seek_frame = 0 + self.expected_remarks = '''Created by DCD pluginREMARKS Created 06 July, 2014 at 17:29Y5~CORD,''' + self._write_files(testfile=self.testfile, + remarks_setting='input') + + def test_written_unit_cell(self): + # there's no expectation that we can write unit cell + # data in NAMD format at the moment + pass + + +class TestDCDWriteCharmm36(TestDCDWrite): + # repeat writing tests for Charmm36 format DCD + # no expectation that we can write unit cell info though (yet) + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.testfile2 = self.tmpdir.name + '/test2.dcd' + self.readfile = DCD_TRICLINIC + self.natoms = 375 + self.expected_frames = 10 + self.seek_frame = 7 + self.expected_remarks = '* CHARMM TRICLINIC BOX TESTING * (OLIVER BECKSTEIN 2014) * BASED ON NPTDYN.INP : SCOTT FELLER, NIH, 7/15/95 ' + self._write_files(testfile=self.testfile, + remarks_setting='input') + + def test_written_unit_cell(self): + # there's no expectation that we can write unit cell + # data in NAMD format at the moment + pass + + +class TestDCDWriteHeaderNAMD(TestDCDWriteHeader): + # repeat header writing tests for NAMD format DCD + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.dcdfile = DCD_NAMD_TRICLINIC + + +class TestDCDWriteHeaderCharmm36(TestDCDWriteHeader): + # repeat header writing tests for Charmm36 format DCD + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.dcdfile = DCD_TRICLINIC + + +class TestDCDReadFrameTestNAMD(TestDCDReadFrame): + # repeat frame reading tests for NAMD format DCD + + def setUp(self): + self.dcdfile = DCD_NAMD_TRICLINIC + self.natoms = 5545 + self.traj_length = 1 + self.new_frame = 0 + self.context_frame = 0 + self.num_iters = 0 + self.selected_legacy_frames = [0] + self.legacy_data = legacy_DCD_NAMD_coords + self.expected_remarks = 'Created by DCD pluginREMARKS Created 06 July, 2014 at 17:29Y5~CORD,' + # expected unit cell based on previous DCD framework read in: + self.expected_unit_cell = np.array([ 38.42659378, 38.39310074, 44.75979996, + 90. , 90. , 60.02891541], + dtype=np.float32) + + +class TestDCDReadFrameTestCharmm36(TestDCDReadFrame): + # repeat frame reading tests for Charmm36 format DCD + + def setUp(self): + self.dcdfile = DCD_TRICLINIC + self.natoms = 375 + self.traj_length = 10 + self.new_frame = 2 + self.context_frame = 5 + self.num_iters = 7 + self.selected_legacy_frames = [1, 4] + self.legacy_data = legacy_DCD_c36_coords + self.expected_remarks = '* 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 ' + # expected unit cell based on previous DCD framework read in: + self.expected_unit_cell = np.array([ 35.44603729, 35.06156158, 34.15850067, + 91.32801819, 61.73519516, 44.4070282], + dtype=np.float32) + + +class TestDCDWriteRandom(object): + # should only be supported for Charmm24 format writing (for now) + + def setUp(self): + self.tmpdir = tempdir.TempDir() + self.testfile = self.tmpdir.name + '/test.dcd' + self.readfile = DCD + self.natoms = 3341 + self.expected_frames = 98 + self.seek_frame = 91 + self.expected_remarks = '''* DIMS ADK SEQUENCE FOR PORE PROGRAM * WRITTEN BY LIZ DENNING (6.2008) * DATE: 6/ 6/ 8 17:23:56 CREATED BY USER: denniej0 ''' + + np.random.seed(1178083) + self.random_unitcells = np.random.uniform(high=80,size=(self.expected_frames, 6)).astype(np.float64) + + with DCDFile(self.readfile) as f_in, DCDFile(self.testfile, 'w') as f_out: + for index, frame in enumerate(f_in): + box=frame.unitcell.astype(np.float64) + f_out.write(xyz=frame.x, + box=self.random_unitcells[index], + step=f_in.istart, + natoms=frame.x.shape[0], + charmm=1, # DCD should be CHARMM + time_step=f_in.delta, + ts_between_saves=f_in.nsavc, + remarks=f_in.remarks) + + def tearDown(self): + try: + os.unlink(self.testfile) + except OSError: + pass + del self.tmpdir + + def test_written_unit_cell_random(self): + with DCDFile(self.testfile) as test: + curr_frame = 0 + while curr_frame < test.n_frames: + written_unitcell = test.read()[1] + ref_unitcell = self.random_unitcells[curr_frame] + + curr_frame += 1 + assert_allclose(written_unitcell, ref_unitcell, + rtol=1e-05) + + +class TestDCDByteArithmetic(object): + + def setUp(self): + self.dcdfile = DCD + self._filesize = os.path.getsize(DCD) + + def test_relative_frame_sizes(self): + # 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(self.dcdfile) as dcd: + first_frame_size = dcd._firstframesize + general_frame_size = dcd._framesize + + assert_equal(first_frame_size >= general_frame_size, True) + + def test_file_size_breakdown(self): + # 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 = self._filesize + with DCDFile(self.dcdfile) as dcd: + actual = dcd._header_size + dcd._firstframesize + ((dcd.n_frames - 1) * dcd._framesize) + assert_equal(actual, expected) + + def test_nframessize_int(self): + # require that the (nframessize / framesize) value used by DCDFile + # is an integer (because nframessize / framesize + 1 = total frames, + # which must also be an int) + with DCDFile(self.dcdfile) as dcd: + nframessize = self._filesize - dcd._header_size - dcd._firstframesize + assert_equal(float(nframessize) % float(dcd._framesize), 0) + + + +class TestDCDByteArithmeticNAMD(TestDCDByteArithmetic): + # repeat byte arithmetic tests for NAMD format DCD + + def setUp(self): + self.dcdfile = DCD_NAMD_TRICLINIC + self._filesize = os.path.getsize(DCD_NAMD_TRICLINIC) + + +class TestDCDByteArithmeticCharmm36(TestDCDByteArithmetic): + # repeat byte arithmetic tests for Charmm36 format DCD + + def setUp(self): + self.dcdfile = DCD_TRICLINIC + self._filesize = os.path.getsize(DCD_TRICLINIC) diff --git a/testsuite/MDAnalysisTests/util.py b/testsuite/MDAnalysisTests/util.py index a4fbddd08b3..5311f262a32 100644 --- a/testsuite/MDAnalysisTests/util.py +++ b/testsuite/MDAnalysisTests/util.py @@ -37,10 +37,12 @@ from contextlib import contextmanager from functools import wraps import importlib -import mock +try: + import mock +except ImportError: # python 3 + from unittest import mock import os - def block_import(package): """Block import of a given package diff --git a/testsuite/setup.py b/testsuite/setup.py index 184052d5260..1480a88fb8a 100755 --- a/testsuite/setup.py +++ b/testsuite/setup.py @@ -206,7 +206,7 @@ def dynamic_author_list(): 'data/*.xml', 'data/coordinates/*', 'data/*xvg', - 'data/*.mmtf', 'data/*.mmtf.gz', 'data/analysis/*' + 'data/*.mmtf', 'data/*.mmtf.gz', 'data/analysis/*', ], }, classifiers=CLASSIFIERS,