From 0e364610daaac4d9efc4d2b3e0bae9443cbf0a54 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 19 Jul 2015 16:40:18 +0100 Subject: [PATCH] Updated Reader API fixed, skip_timestep, periodic, skip and delta no longer required for Readers (Issue #350) All writers now refer to time between steps as "dt" (was previously delta in XTC, TRR and NCDF Writers) (Issue #206) Timesteps now have a default dt of 1.0 ps (and issue warning when using this) Timestep attribute time now tries to return data['time'] first, then dt * frame, using the default value of 1.0 ps if necessary. (Issue #308) --- package/CHANGELOG | 8 ++ package/MDAnalysis/coordinates/DCD.py | 12 ++- package/MDAnalysis/coordinates/DLPoly.py | 1 - package/MDAnalysis/coordinates/GMS.py | 5 +- package/MDAnalysis/coordinates/INPCRD.py | 9 ++- package/MDAnalysis/coordinates/MOL2.py | 3 - package/MDAnalysis/coordinates/PDB.py | 12 --- package/MDAnalysis/coordinates/PDBQT.py | 6 +- package/MDAnalysis/coordinates/TRJ.py | 32 ++++---- package/MDAnalysis/coordinates/TRZ.py | 48 +++++------- package/MDAnalysis/coordinates/XYZ.py | 23 ++---- package/MDAnalysis/coordinates/__init__.py | 24 +++--- package/MDAnalysis/coordinates/base.py | 75 +++++++++++-------- .../MDAnalysis/coordinates/xdrfile/core.py | 75 +++++++++---------- testsuite/MDAnalysisTests/test_coordinates.py | 36 ++++----- testsuite/MDAnalysisTests/test_reader_api.py | 8 +- 16 files changed, 176 insertions(+), 201 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 56171dd4b24..1a0e43d43f9 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -74,8 +74,16 @@ The rules for this file: should not affect any code because MDAnalysis.selections.get_writer() already existed. * Moved norm, normal, angle, stp and dihedral from lib.util to lib.mdamath + * All Readers have a default "dt" of 1.0 ps. This applies also to single + frame readers, these would previously raise an error on accessing dt. + (Issue #350) + * NCDF Reader no longer has a default skip_timestep of 1 (Issue #350) + * All Writers now refer to time between written Timesteps as "dt", + was previously "delta" in some Writers. (Issue #206) Fixes + * Amber TRJ and NCDF Reader & Writers now use 'dt' instead of 'delta' + to refer to time elapsed between timesteps. (Issue #350) * Fixed TPRParser considering LJ 1..4 exclusions as bonds (Issue #351) * ChainReaders no longer leak (Issue #312) * analysis.hbonds.HydrogenBondAnalysis performs a sanity check for diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index a83b6b0406b..d61df209a31 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -419,7 +419,13 @@ def __init__(self, dcdfilename, **kwargs): self.fixed = 0 self.periodic = False + # This reads skip_timestep and delta from header self._read_dcd_header() + + # Convert delta to ps + delta = mdaunits.convert(self.delta, 'AKMA', 'ps') + + self._ts_kwargs.update({'dt':self.skip_timestep * delta}) self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) # Read in the first timestep self._read_next_timestep() @@ -472,7 +478,7 @@ def _read_next_timestep(self, ts=None): """ if ts is None: ts = self.ts - ts._frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, self.skip) + ts._frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1) ts.frame += 1 return ts @@ -583,6 +589,10 @@ def Writer(self, filename, **kwargs): # dt keyword is simply passed through if provided return DCDWriter(filename, numatoms, **kwargs) + @property + def dt(self): + return self.skip_timestep * self.convert_time_from_native(self.delta) + DCDReader._read_dcd_header = new.instancemethod(_dcdmodule.__read_dcd_header, None, DCDReader) DCDReader._read_next_frame = new.instancemethod(_dcdmodule.__read_next_frame, None, DCDReader) diff --git a/package/MDAnalysis/coordinates/DLPoly.py b/package/MDAnalysis/coordinates/DLPoly.py index 3c7140fe0d8..e6a9dc10a77 100644 --- a/package/MDAnalysis/coordinates/DLPoly.py +++ b/package/MDAnalysis/coordinates/DLPoly.py @@ -26,7 +26,6 @@ import numpy as np -import MDAnalysis from . import base from . import core diff --git a/package/MDAnalysis/coordinates/GMS.py b/package/MDAnalysis/coordinates/GMS.py index 7247623a327..fdbe4895d77 100644 --- a/package/MDAnalysis/coordinates/GMS.py +++ b/package/MDAnalysis/coordinates/GMS.py @@ -76,13 +76,11 @@ def __init__(self, outfilename, **kwargs): self._runtyp = None self.ts = self._Timestep(0) # need for properties initial calculations - self.periodic = False - self.delta = kwargs.pop("delta", 1.0) # there is no time, so let delta be 1 + # update runtyp property self.runtyp if not self.runtyp in ['optimize', 'surface']: raise AttributeError('Wrong RUNTYP= '+self.runtyp) - self.skip_timestep = 1 self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) # update numframes property @@ -248,7 +246,6 @@ def open_trajectory(self): # reset ts ts = self.ts ts.frame = -1 - ts.time = 0 return self.outfile def close(self): diff --git a/package/MDAnalysis/coordinates/INPCRD.py b/package/MDAnalysis/coordinates/INPCRD.py index 5dee31c52fe..a12623f72f5 100644 --- a/package/MDAnalysis/coordinates/INPCRD.py +++ b/package/MDAnalysis/coordinates/INPCRD.py @@ -37,13 +37,14 @@ def _read_first_frame(self): self.title = inf.readline().strip() line = inf.readline().split() self.numatoms = int(line[0]) + + self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) try: time = float(line[1]) except IndexError: - time = 0.0 - - self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) - self.ts.time = time + pass + else: + self.ts.time = time self.ts.frame = 0 for p in xrange(self.numatoms // 2): diff --git a/package/MDAnalysis/coordinates/MOL2.py b/package/MDAnalysis/coordinates/MOL2.py index f7068daba48..f7670b192fc 100644 --- a/package/MDAnalysis/coordinates/MOL2.py +++ b/package/MDAnalysis/coordinates/MOL2.py @@ -85,9 +85,6 @@ def __init__(self, filename, **kwargs): self.substructure = {} self.frames = blocks self.numframes = len(blocks) - self.periodic = False - self.delta = 0 - self.skip_timestep = 1 def parse_block(self, block): sections = {} diff --git a/package/MDAnalysis/coordinates/PDB.py b/package/MDAnalysis/coordinates/PDB.py index 16a5ace3b78..8cd76bb4f7b 100644 --- a/package/MDAnalysis/coordinates/PDB.py +++ b/package/MDAnalysis/coordinates/PDB.py @@ -476,15 +476,6 @@ def __init__(self, filename, **kwargs): frames = {} - ts_kwargs = {} - for att in ('dt', 'time_offset'): - try: - val = kwargs[att] - except KeyError: - pass - else: - ts_kwargs[att] = val - self.ts = self._Timestep(self._numatoms, **self._ts_kwargs) pos = 0 # atom position for filling coordinates array @@ -556,9 +547,6 @@ def __init__(self, filename, **kwargs): self.frames = frames self.numframes = len(frames) if frames else 1 - self.periodic = False - self.delta = 0 - self.skip_timestep = 1 def get_occupancy(self): """Return an array of occupancies in atom order.""" diff --git a/package/MDAnalysis/coordinates/PDBQT.py b/package/MDAnalysis/coordinates/PDBQT.py index e80ddc68a02..c3956c9f773 100644 --- a/package/MDAnalysis/coordinates/PDBQT.py +++ b/package/MDAnalysis/coordinates/PDBQT.py @@ -176,11 +176,7 @@ def _c(start, stop, typeclass=float): self.convert_pos_from_native(self.ts._pos) # in-place ! self.convert_pos_from_native(self.ts._unitcell[:3]) # in-place ! (only lengths) self.numframes = 1 - self.fixed = 0 - self.skip = 1 - self.periodic = False - self.delta = 0 - self.skip_timestep = 1 + # hack for PDBQTParser: self._atoms = np.rec.fromrecords(atoms, names="serial,name,resName,chainID,resSeq,occupancy,tempFactor," diff --git a/package/MDAnalysis/coordinates/TRJ.py b/package/MDAnalysis/coordinates/TRJ.py index b3ec267b229..a88d3f68e26 100644 --- a/package/MDAnalysis/coordinates/TRJ.py +++ b/package/MDAnalysis/coordinates/TRJ.py @@ -76,7 +76,7 @@ boxes. * The trajectory does not contain time information so we simply set - the time step to 1 ps (or the user could provide it as kwarg *delta*) + the time step to 1 ps (or the user could provide it as kwarg *dt*) * No direct access of frames is implemented, only iteration through the trajectory. @@ -177,7 +177,7 @@ class TRJReader(base.Reader): `numatoms` is left at its default value of ``None``. The length of a timestep is not stored in the trajectory itself but can - be set by passing the `delta` keyword argument to the constructor; it + be set by passing the `dt` keyword argument to the constructor; it is assumed to be in ps. The default value is 1 ps. Functionality is currently limited to simple iteration over the @@ -187,6 +187,7 @@ class TRJReader(base.Reader): .. versionchanged:: 0.11.0 Frames now 0-based instead of 1-based + kwarg 'delta' renamed to 'dt', for uniformity with other Readers """ format = 'TRJ' units = {'time': 'ps', 'length': 'Angstrom'} @@ -201,8 +202,6 @@ def __init__(self, filename, numatoms=None, **kwargs): self._numframes = None self.trjfile = None # have _read_next_timestep() open it properly! - self.skip_timestep = 1 # always 1 for trj at the moment - self.delta = kwargs.pop("delta", 1.0) # can set delta manually, default is 1ps self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) # FORMAT(10F8.3) (X(i), Y(i), Z(i), i=1,NATOM) @@ -351,7 +350,7 @@ def open_trajectory(self): # reset ts ts = self.ts ts.frame = -1 - ts.time = 0 + return self.trjfile def close(self): @@ -403,6 +402,7 @@ class NCDFReader(base.Reader): Added ability to read Forces .. versionchanged:: 0.11.0 Frame labels now 0-based instead of 1-based + kwarg 'delta' renamed to 'dt', for uniformity with other Readers """ format = 'NCDF' @@ -477,8 +477,6 @@ def __init__(self, filename, numatoms=None, **kwargs): self.has_velocities = 'velocities' in self.trjfile.variables self.has_forces = 'forces' in self.trjfile.variables - self.skip_timestep = 1 # always 1 for trj at the moment ? CHECK DOCS?? - self.delta = kwargs.pop("delta", 1.0) # SHOULD GET FROM NCDF (as mean(diff(time)))?? self.periodic = 'cell_lengths' in self.trjfile.variables self._current_frame = 0 @@ -549,7 +547,7 @@ def Writer(self, filename, **kwargs): :Keywords: *numatoms* number of atoms - *delta* + *dt* length of one timestep in picoseconds *remarks* string that is stored in the title field @@ -558,7 +556,7 @@ def Writer(self, filename, **kwargs): """ numatoms = kwargs.pop('numatoms', self.numatoms) kwargs.setdefault('remarks', self.remarks) - kwargs.setdefault('delta', self.delta) + kwargs.setdefault('dt', self.dt) return NCDFWriter(filename, numatoms, **kwargs) @@ -583,6 +581,8 @@ class NCDFWriter(base.Writer): .. versionchanged:: 0.10.0 Added ability to write velocities and forces + .. versionchanged:: 0.11.0 + kwarg 'delta' renamed to 'dt', for uniformity with other Readers """ format = 'NCDF' @@ -590,7 +590,7 @@ class NCDFWriter(base.Writer): units = {'time': 'ps', 'length': 'Angstrom', 'velocity': 'Angstrom/ps', 'force': 'kcal/(mol*Angstrom)'} - def __init__(self, filename, numatoms, start=0, step=1, delta=1.0, remarks=None, + def __init__(self, filename, numatoms, start=0, step=1, dt=1.0, remarks=None, convert_units=None, zlib=False, cmplevel=1, **kwargs): """Create a new NCDFWriter @@ -605,7 +605,7 @@ def __init__(self, filename, numatoms, start=0, step=1, delta=1.0, remarks=None, starting timestep *step* skip between subsequent timesteps - *delta* + *dt* timestep *convert_units* ``True``: units are converted to the AMBER base format; ``None`` selects @@ -630,7 +630,7 @@ def __init__(self, filename, numatoms, start=0, step=1, delta=1.0, remarks=None, self.start = start # do we use those? self.step = step # do we use those? - self.delta = delta + self.dt = dt self.remarks = remarks or "AMBER NetCDF format (MDAnalysis.coordinates.trj.NCDFWriter)" self.zlib = zlib @@ -771,14 +771,12 @@ def _write_next_timestep(self, ts): try: time = self.convert_time_to_native(ts.time, inplace=False) except AttributeError: - time = ts.frame * self.convert_time_to_native(self.delta, inplace=False) + time = ts.frame * self.convert_time_to_native(self.dt, inplace=False) unitcell = self.convert_dimensions_to_unitcell(ts) else: pos = ts._pos - try: - time = ts.time - except AttributeError: - time = ts.frame * self.delta + time = ts.time + unitcell = ts.dimensions # write step diff --git a/package/MDAnalysis/coordinates/TRZ.py b/package/MDAnalysis/coordinates/TRZ.py index 9b979da842b..f3a7975883b 100644 --- a/package/MDAnalysis/coordinates/TRZ.py +++ b/package/MDAnalysis/coordinates/TRZ.py @@ -81,6 +81,7 @@ from . import base from ..core import flags from ..lib import util +from ..lib.util import cached from .core import triclinic_box, triclinic_vectors @@ -154,12 +155,8 @@ def __init__(self, trzfilename, numatoms=None, **kwargs): raise ValueError('TRZReader requires the numatoms keyword') self.trzfile = util.anyopen(self.filename, 'rb') - + self._cache = dict() self._numatoms = numatoms - self._numframes = None - self._delta = None - self._dt = None - self._skip_timestep = None self._read_trz_header() self.ts = Timestep(self.numatoms, velocities=True, forces=self.has_force, @@ -270,16 +267,13 @@ def numatoms(self): return self._numatoms @property + @cached('numframes') def numframes(self): """Total number of frames in a trajectory""" - if not self._numframes is None: - return self._numframes try: - self._numframes = self._read_trz_numframes(self.trzfile) + return self._read_trz_numframes(self.trzfile) except IOError: return 0 - else: - return self._numframes def _read_trz_numframes(self, trzfile): """Uses size of file and dtype information to determine how many frames exist @@ -297,10 +291,7 @@ def _read_trz_numframes(self, trzfile): return nframes @property - def time(self): - return self.ts.time - - @property + @cached('dt') def dt(self): """The amount of time between frames in ps @@ -308,42 +299,41 @@ def dt(self): stitched together) Returns 0 in case of IOError """ - if not self._dt is None: - return self._dt + curr_frame = self.ts.frame try: t0 = self.ts.time self.next() t1 = self.ts.time - self._dt = t1 - t0 + dt = t1 - t0 except IOError: return 0 + else: + return dt finally: - self.rewind() - return self._dt + self._read_frame(curr_frame) @property + @cached('delta') def delta(self): """MD integration timestep""" - if not self._delta is None: - return self._delta - self._delta = self.dt / self.skip_timestep - return self._delta + return self.dt / self.skip_timestep @property + @cached('skip_timestep') def skip_timestep(self): """Timesteps between trajectory frames""" - if not self._skip_timestep is None: - return self._skip_timestep + curr_frame = self.ts.frame try: t0 = self.ts._frame self.next() t1 = self.ts._frame - self._skip_timestep = t1 - t0 + skip_timestep = t1 - t0 except IOError: return 0 + else: + return skip_timestep finally: - self.rewind() - return self._skip_timestep + self._read_frame(curr_frame) def _read_frame(self, frame): """Move to *frame* and fill timestep with data. @@ -402,7 +392,7 @@ def open_trajectory(self): #Reset ts ts = self.ts ts.frame = -1 - ts.time = 0 + return self.trzfile def Writer(self, filename, numatoms=None): diff --git a/package/MDAnalysis/coordinates/XYZ.py b/package/MDAnalysis/coordinates/XYZ.py index 74e6b154538..72eaf5310c3 100644 --- a/package/MDAnalysis/coordinates/XYZ.py +++ b/package/MDAnalysis/coordinates/XYZ.py @@ -257,14 +257,10 @@ def __init__(self, filename, **kwargs): # note that, like for xtc and trr files, _numatoms and _numframes are used quasi-private variables # to prevent the properties being recalculated # this is because there is no indexing so the way it measures the number of frames is to read the whole file! - self._numatoms = None + self._numatoms = self._read_xyz_natoms(self.filename) self._numframes = None - self.periodic = False - self.delta = kwargs.pop("delta", 1.0) # can set delta manually, default is 1ps (taken from TRJReader) - self.skip_timestep = 1 - - self.ts = self._Timestep(self.numatoms, **self._ts_kwargs) + self.ts = self._Timestep(self._numatoms, **self._ts_kwargs) # Haven't quite figured out where to start with all the self._reopen() etc. # (Also cannot just use seek() or reset() because that would break with urllib2.urlopen() streams) @@ -273,20 +269,13 @@ def __init__(self, filename, **kwargs): @property def numatoms(self): """number of atoms in a frame""" - if not self._numatoms is None: # return cached value - return self._numatoms - try: - self._numatoms = self._read_xyz_natoms(self.filename) - except IOError: - return 0 - else: - return self._numatoms + return self._numatoms def _read_xyz_natoms(self, filename): # this assumes that this is only called once at startup and that the filestream is already open # (FRAGILE) - n = self.xyzfile.readline() - self.close() + with util.anyopen(self.filename, 'r') as f: + n = f.readline() # need to check type of n return int(n) @@ -367,7 +356,7 @@ def open_trajectory(self): # reset ts ts = self.ts ts.frame = -1 - ts.time = 0 + return self.xyzfile def Writer(self, filename, **kwargs): diff --git a/package/MDAnalysis/coordinates/__init__.py b/package/MDAnalysis/coordinates/__init__.py index b5ee3cb05f9..d443de4ce8a 100644 --- a/package/MDAnalysis/coordinates/__init__.py +++ b/package/MDAnalysis/coordinates/__init__.py @@ -531,18 +531,6 @@ number of atoms (coordinate sets) in a frame (constant) ``numframes`` total number of frames (if known) -- ``None`` if not known - ``fixed`` - bool, saying if there are fixed atoms (e.g. dcds) - ``skip`` - step size for iterating through the trajectory [1] - ``skip_timestep`` - number of integrator steps between frames + 1 (i.e. - the stride at which the MD simulation was sampled) - ``delta`` - integrator time step (in native units); hence the "length" - of a trajctory frame is ``skip_timestep*delta`` time units - ``periodic`` - contains box information for periodic boundary conditions ``ts`` the :class:`~base.Timestep` object; typically customized for each trajectory format and derived from :class:`base.Timestep`. @@ -567,8 +555,20 @@ **Optional attributes** + ``delta`` + integrator time step (in native units); hence the "length" + of a trajctory frame is ``skip_timestep*delta`` time units ``compressed`` string that identifies the compression (e.g. "gz" or "bz2") or ``None``. + ``fixed`` + bool, saying if there are fixed atoms (e.g. dcds) + ``periodic`` + contains box information for periodic boundary conditions + ``skip`` + step size for iterating through the trajectory [1] + ``skip_timestep`` + number of integrator steps between frames + 1 (i.e. + the stride at which the MD simulation was sampled) Trajectory Writer class diff --git a/package/MDAnalysis/coordinates/base.py b/package/MDAnalysis/coordinates/base.py index 6526da13cef..81f9905a847 100644 --- a/package/MDAnalysis/coordinates/base.py +++ b/package/MDAnalysis/coordinates/base.py @@ -173,7 +173,6 @@ def __init__(self, numatoms, **kwargs): # self.frame to 0 self.frame = -1 self._numatoms = numatoms - self.time = 0.0 self.data = {} @@ -592,6 +591,40 @@ def volume(self): """volume of the unitcell""" return core.box_volume(self.dimensions) + @property + def dt(self): + """The time difference in ps between timesteps + + .. Note:: + This defaults to 1.0 ps in the absence of time data + + .. versionadded:: 0.11.0 + """ + try: + return self.data['dt'] + except KeyError: + warnings.warn("Reader has no dt information, set to 1.0 ps") + return 1.0 + + @dt.setter + def dt(self, new): + self.data['dt'] = new + + @property + def time(self): + """The time in ps of this timestep + + .. versionadded:: 0.11.0 + """ + try: + return self.data['time'] + except KeyError: + return self.dt * self.frame + + @time.setter + def time(self, new): + self.data['time'] = new + class IObase(object): """Base class bundling common functionality for trajectory I/O. @@ -802,7 +835,7 @@ def rewind(self): @property def dt(self): """Time between two trajectory frames in picoseconds.""" - return self.skip_timestep * self.convert_time_from_native(self.delta) + return self.ts.dt @property def totaltime(self): @@ -821,15 +854,10 @@ def frame(self): def time(self): """Time of the current frame in MDAnalysis time units (typically ps). - time = :attr:`Timestep.frame` * :attr:`Reader.dt` + This is either read straight from the Timestep, or calculated as + time = :attr:`Timestep.frame` * :attr:`Timestep.dt` """ - try: - return self.ts.frame * self.dt - except KeyError: - # single frame formats fail with KeyError because they do not define - # a unit for time so we just always return 0 because time is a required - # attribute of the Reader - return 0.0 + return self.ts.time def Writer(self, filename, **kwargs): """Returns a trajectory writer with the same properties as this trajectory.""" @@ -839,16 +867,15 @@ def Writer(self, filename, **kwargs): def OtherWriter(self, filename, **kwargs): """Returns a writer appropriate for *filename*. - Sets the default keywords *start*, *step* and *delta* (if + Sets the default keywords *start*, *step* and *dt* (if available). *numatoms* is always set from :attr:`Reader.numatoms`. .. SeeAlso:: :meth:`Reader.Writer` and :func:`MDAnalysis.Writer` """ kwargs['numatoms'] = self.numatoms # essential kwargs.setdefault('start', self.frame) - kwargs.setdefault('step', self.skip_timestep) try: - kwargs.setdefault('delta', self.dt) + kwargs.setdefault('dt', self.dt) except KeyError: pass return core.writer(filename, **kwargs) @@ -970,8 +997,8 @@ def _check_slice_indices(self, start, stop, step): return start, stop, step def __repr__(self): - return "< %s %r with %d frames of %d atoms (%d fixed) >" % \ - (self.__class__.__name__, self.filename, self.numframes, self.numatoms, self.fixed) + return "< %s %r with %d frames of %d atoms>" % \ + (self.__class__.__name__, self.filename, self.numframes, self.numatoms) class Reader(ProtoReader): """Base class for trajectory readers that extends :class:`ProtoReader` with a :meth:`__del__` method. @@ -998,13 +1025,6 @@ def __init__(self, filename, convert_units=None, **kwargs): convert_units = flags['convert_lengths'] self.convert_units = convert_units - self.fixed = False - self.periodic = True - self.skip = 1 - self._delta = None - self._dt = None - self._skip_timestep = None - ts_kwargs = {} for att in ('dt', 'time_offset'): try: @@ -1079,7 +1099,7 @@ def __init__(self, filenames, **kwargs): self.skip = kwargs.get('skip', 1) self._default_delta = kwargs.pop('delta', None) self.numatoms = self._get_same('numatoms') - self.fixed = self._get_same('fixed') + #self.fixed = self._get_same('fixed') # Translation between virtual frames and frames in individual # trajectories. @@ -1293,10 +1313,10 @@ def get_flname(self, filename): # retrieve the actual filename of the list elem return filename[0] if isinstance(filename, tuple) else filename def __repr__(self): - return "< %s %r with %d frames of %d atoms (%d fixed) >" % \ + return "< %s %r with %d frames of %d atoms>" % \ (self.__class__.__name__, [os.path.basename(self.get_flname(fn)) for fn in self.filenames], - self.numframes, self.numatoms, self.fixed) + self.numframes, self.numatoms) class Writer(IObase): @@ -1392,11 +1412,6 @@ def __init__(self, filename, convert_units=None, **kwargs): self.convert_units = convert_units self.numframes = 1 - self.fixed = 0 - self.skip = 1 - self.periodic = False - self.delta = 0 - self.skip_timestep = 1 ts_kwargs = {} for att in ('dt', 'time_offset'): diff --git a/package/MDAnalysis/coordinates/xdrfile/core.py b/package/MDAnalysis/coordinates/xdrfile/core.py index 155a20e489b..9869294e6e2 100644 --- a/package/MDAnalysis/coordinates/xdrfile/core.py +++ b/package/MDAnalysis/coordinates/xdrfile/core.py @@ -78,6 +78,7 @@ from MDAnalysis.coordinates import base from MDAnalysis.coordinates.core import triclinic_box, triclinic_vectors import MDAnalysis.core +from ...lib.util import cached # This is the XTC class. The TRR overrides with it's own. @@ -129,7 +130,7 @@ class TrjWriter(base.Writer): #: override to define trajectory format of the reader (XTC or TRR) format = None - def __init__(self, filename, numatoms, start=0, step=1, delta=None, precision=1000.0, remarks=None, + def __init__(self, filename, numatoms, start=0, step=1, dt=None, precision=1000.0, remarks=None, convert_units=None): """ Create a new TrjWriter @@ -141,14 +142,16 @@ def __init__(self, filename, numatoms, start=0, step=1, delta=None, precision=10 :Keywords: *start* - starting timestep; only used when *delta* is set. + starting timestep frame; only used when *dt* is set. *step* - skip between subsequent timesteps; only used when *delta* is set. - *delta* - timestep to use. If set will override any time information contained in the - passed :class:`Timestep` objects; otherwise that will be used. If in the - latter case :attr:`~Timestep.time` is unavailable the TrjWriter will default - to setting the trajectory time at 1 MDAnalysis unit (typically 1ps) per step. + skip in frames between subsequent timesteps; only used when *dt* is set. + *dt* + time between frames to use. If set will override any time information + contained in thee + passed :class:`Timestep` objects, which will otherwise be used. + The :attr:`~Timestep.time` attribute defaults to a timestep of + to setting the trajectory time at 1 ps per step if there is no + time information. *precision* accuracy for lossy XTC format as a power of 10 (ignored for TRR) [1000.0] @@ -161,9 +164,9 @@ def __init__(self, filename, numatoms, start=0, step=1, delta=None, precision=10 The TRR writer is now able to write TRRs without coordinates/velocities/forces, depending on the properties available in the :class:`Timestep` objects passed to :meth:`~TRRWriter.write`. + .. versionchanged:: 0.11.0 + Keyword "delta" renamed to "dt" """ - assert self.format in ('XTC', 'TRR') - if numatoms == 0: raise ValueError("TrjWriter: no atoms in output trajectory") self.filename = filename @@ -182,7 +185,7 @@ def __init__(self, filename, numatoms, start=0, step=1, delta=None, precision=10 self.frames_written = 0 self.start = start self.step = step - self.delta = delta + self.dt = dt self.remarks = remarks self.precision = precision # only for XTC self.xdrfile = libxdrfile2.xdrfile_open(self.filename, 'w') @@ -221,17 +224,13 @@ def _write_next_timestep(self, ts): # (1) data common to XTC and TRR - # Time-writing logic: if the writer was created with a delta parameter, - # use delta*(start+step*frames_written) - # otherwise use the provided Timestep obj time attribute. Fallback to 1 if not present. - if self.delta is None: - try: - time = ts.time - except AttributeError: - # Default to 1 time unit. - time = self.start + self.step * self.frames_written + # Time-writing logic: if the writer was created with a dt parameter, + # use dt*(start+step*frames_written) + # otherwise use the provided Timestep obj time attribute + if self.dt is None: + time = ts.time else: - time = (self.start + self.step * self.frames_written) * self.delta + time = (self.start + self.step * self.frames_written) * self.dt if self.convert_units: time = self.convert_time_to_native(time, inplace=False) @@ -239,7 +238,7 @@ def _write_next_timestep(self, ts): try: step = int(ts._frame) except AttributeError: - # bogus, should be actual MD step number, i.e. frame * delta/dt + # bogus, should be actual MD step number, i.e. frame * dt/delta step = ts.frame unitcell = self.convert_dimensions_to_unitcell(ts).astype(np.float32) # must be float32 (!) @@ -353,10 +352,11 @@ def __init__(self, filename, sub=None, **kwargs): .. versionchanged:: 0.9.0 New keyword *refresh_offsets* - + .. versionchanged:: 0.11.0 + Renamed "delta" attribute to "dt" """ super(TrjReader, self).__init__(filename, **kwargs) - + self._cache = dict() # Convert filename to ascii because of SWIG bug. # See: http://sourceforge.net/p/swig/feature-requests/75 # Only needed for Python < 3 @@ -367,11 +367,8 @@ def __init__(self, filename, sub=None, **kwargs): self.xdrfile = None self._numframes = None # takes a long time, avoid accessing self.numframes - self.skip_timestep = 1 # always 1 for xdr files - self._delta = None # compute from time in first two frames! + self._dt = None # compute from time in first two frames! self._offsets = None # storage of offsets in the file - self.skip = 1 - self.periodic = False # actual number of atoms in the trr file # first time file is opened, exception should be thrown if bad file @@ -468,25 +465,24 @@ def offsets(self): return self._offsets @property - def delta(self): + @cached('dt') + def dt(self): """Time step length in ps. The result is computed from the trajectory and cached. If for any reason the trajectory cannot be read then 0 is returned. """ # no need for conversion: it's alread in our base unit ps - if not self._delta is None: # return cached value - return self._delta try: t0 = self.ts.time self.next() t1 = self.ts.time - self._delta = t1 - t0 + dt = t1 - t0 + return dt except IOError: return 0 finally: self.rewind() - return self._delta def _offset_filename(self): head, tail = os.path.split(self.filename) @@ -656,7 +652,7 @@ def open_trajectory(self): ts.data['status'] = libxdrfile2.exdrOK ts.frame = -1 ts._frame = 0 - ts.time = 0 + # additional data for XTC ts.data['prec'] = 0 # additional data for TRR @@ -682,19 +678,22 @@ def Writer(self, filename, **kwargs): :Keywords: *numatoms* number of atoms - *delta* + *dt* Time interval between frames. *precision* accuracy for lossy XTC format as a power of 10 (ignored for TRR) [1000.0] :Returns: appropriate :class:`TrjWriter` + + .. versionchanged:: 0.11.0 + Changed "delta" keyword to "dt" """ numatoms = kwargs.pop('numatoms', self.numatoms) - kwargs['step'] = self.skip_timestep - kwargs.setdefault('delta', self.delta) + + kwargs.setdefault('dt', self.dt) try: - kwargs['start'] = self[0].time / kwargs['delta'] + kwargs['start'] = self[0].time / kwargs['dt'] except (AttributeError, ZeroDivisionError): kwargs['start'] = 0 try: diff --git a/testsuite/MDAnalysisTests/test_coordinates.py b/testsuite/MDAnalysisTests/test_coordinates.py index 30b9d664b9c..b510397efec 100644 --- a/testsuite/MDAnalysisTests/test_coordinates.py +++ b/testsuite/MDAnalysisTests/test_coordinates.py @@ -47,7 +47,7 @@ def atom_distance(a, b): """Calculate the distance between two atoms a and b.""" - r = a.pos - b.pos + r = a.position - b.position return np.sqrt(np.sum(r ** 2)) @@ -549,7 +549,7 @@ def tearDown(self): def test_write_trajectory(self): t = self.universe.trajectory - W = self.Writer(self.outfile, t.numatoms, delta=t.delta) + W = self.Writer(self.outfile, t.numatoms, dt=t.dt) self._copy_traj(W) def test_OtherWriter(self): @@ -769,8 +769,8 @@ def go_to_2(traj=self.universe.trajectory): assert_raises(IndexError, go_to_2) def test_dt(self): - """testing that accessing universe.trajectory.dt raises a KeyError for single frame readers""" - assert_raises(KeyError, self.universe.trajectory.__getattribute__, "dt") + """testing that accessing universe.trajectory.dt gives 1.0 (the default)""" + assert_equal(1.0, self.universe.trajectory.dt) def test_coordinates(self): A10CA = self.universe.atoms.CA[10] @@ -1410,8 +1410,8 @@ def test_frame(self): assert_equal(self.universe.trajectory.frame, 0, "wrong frame number (should be 0 for 0-based ts.frame)") def test_dt(self): - """testing that accessing universe.trajectory.dt raises a KeyError for single frame readers""" - assert_raises(KeyError, self.universe.trajectory.__getattribute__, "dt") + """testing that accessing universe.trajectory.dt gives default of 1.0""" + assert_equal(self.universe.trajectory.dt, 1.0) def test_coordinates(self): A10CA = self.universe.SYSTEM.CA[10] @@ -1640,8 +1640,8 @@ def test_frame(self): @dec.slow def test_dt(self): - """testing that accessing universe.trajectory.dt raises a KeyError for single frame readers""" - assert_raises(KeyError, self.universe.trajectory.__getattribute__, "dt") + """testing that accessing universe.trajectory.dt returns the default of 1.0 ps""" + assert_equal(self.universe.trajectory.dt, 1.0) @dec.slow def test_coordinates(self): @@ -1783,7 +1783,7 @@ def tearDown(self): def test_write_trajectory(self): """Test writing DCD trajectories (Issue 50)""" t = self.universe.trajectory - W = self.Writer(self.outfile, t.numatoms, delta=t.delta, step=t.skip_timestep) + W = self.Writer(self.outfile, t.numatoms, dt=t.dt, step=t.skip_timestep) for ts in self.universe.trajectory: W.write_next_timestep(ts) W.close() @@ -2378,7 +2378,7 @@ def test_frame(self): @dec.slow def test_time(self): self.trajectory[4] - assert_almost_equal(self.universe.trajectory.time, 400.0, 5, + assert_almost_equal(self.universe.trajectory.time, 400.0, 3, err_msg="wrong time of frame") @dec.slow @@ -2525,7 +2525,7 @@ def tearDown(self): def test_write_trajectory(self): """Test writing Gromacs trajectories (Issue 38)""" t = self.universe.trajectory - W = self.Writer(self.outfile, t.numatoms, delta=t.delta, step=t.skip_timestep) + W = self.Writer(self.outfile, t.numatoms, dt=t.dt) for ts in self.universe.trajectory: W.write_next_timestep(ts) W.close() @@ -2548,7 +2548,7 @@ def test_timestep_not_modified_by_writer(self): x = ts._pos.copy() time = ts.time - W = self.Writer(self.outfile, trj.numatoms, delta=trj.delta, step=trj.skip_timestep) + W = self.Writer(self.outfile, trj.numatoms, dt=trj.dt) trj[-1] # last timestep (so that time != 0) (say it again, just in case...) W.write_next_timestep(ts) W.close() @@ -2566,7 +2566,7 @@ class TestTRRWriter(_GromacsWriter): def test_velocities(self): t = self.universe.trajectory - W = self.Writer(self.outfile, t.numatoms, delta=t.delta, step=t.skip_timestep) + W = self.Writer(self.outfile, t.numatoms, dt=t.dt) for ts in self.universe.trajectory: W.write_next_timestep(ts) W.close() @@ -2583,7 +2583,7 @@ def test_velocities(self): def test_gaps(self): """Tests the writing and reading back of TRRs with gaps in any of the coordinates/velocities properties.""" t = self.universe.trajectory - W = self.Writer(self.outfile, t.numatoms, delta=t.delta, step=t.skip_timestep) + W = self.Writer(self.outfile, t.numatoms, dt=t.dt) for ts in self.universe.trajectory: # Inset some gaps in the properties: coords every 4 steps, vels every 2. if not ts.frame % 4: @@ -2677,9 +2677,7 @@ def setUp(self): self.universe = mda.Universe(PRMncdf, NCDF) fd, self.outfile = tempfile.mkstemp(suffix=self.ext) os.close(fd) - self.Writer = MDAnalysis.Writer(self.outfile, numatoms=self.universe.atoms.numberOfAtoms(), - delta=self.universe.trajectory.delta, - step=self.universe.trajectory.skip_timestep) + self.Writer = MDAnalysis.Writer(self.outfile, numatoms=self.universe.atoms.numberOfAtoms()) def tearDown(self): try: @@ -2807,10 +2805,6 @@ def test_velocities(self): def test_delta(self): assert_almost_equal(self.trz.delta, self.ref_delta, self.prec, "wrong time delta in trz") - def test_delta_2(self): - self.trz._delta = newref = 4.00 - assert_almost_equal(self.trz.delta, newref, self.prec, "couldn't inject to delta") - def test_time(self): assert_almost_equal(self.trz.time, self.ref_time, self.prec, "wrong time value in trz") diff --git a/testsuite/MDAnalysisTests/test_reader_api.py b/testsuite/MDAnalysisTests/test_reader_api.py index 0a5cbe75184..efbc53ae2ed 100644 --- a/testsuite/MDAnalysisTests/test_reader_api.py +++ b/testsuite/MDAnalysisTests/test_reader_api.py @@ -29,11 +29,6 @@ def __init__(self, filename, **kwargs): self.filename = filename self.numframes = 10 self.numatoms = 10 - self.skip = 1 - self.fixed = 0 - self.periodic = True - self.delta = 1.0 - self.skip_timestep = 1 # ts isn't a real timestep, but just an integer # whose value represents the frame number (0 based) self.ts = Timestep(self.numatoms) @@ -81,8 +76,7 @@ def setUp(self): def test_required_attributes(self): """Test that Reader has the required attributes""" - for attr in ['filename', 'numatoms', 'numframes', 'fixed', 'skip', - 'skip_timestep', 'delta', 'periodic', 'ts', + for attr in ['filename', 'numatoms', 'numframes', 'ts', 'units', 'format']: assert_equal(hasattr(self.reader, attr), True, "Missing attr: {0}".format(attr))