MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.29k stars 646 forks source link

iterating with TRR reader ignores IOError from incomplete/damaged frame #2267

Open ianmkenney opened 5 years ago

ianmkenney commented 5 years ago

Expected behavior

I have a TRR "000001-00043359.trr" whose last frame is damaged (as shown by gmx check -f 000001-00043359.trr.

For a TRR trajectory with a damaged frame, iteration

for i in U1.trajectory:
     pass

should fail with an error such as IOError.

It should be the same error that comes up when I access the damaged frame directly:

In [34]: U1 = mda.Universe("input_structure.dms", "000001-00043359.trr")
~/.virtualenvs/napaanton2/local/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.py:195: UserWarning: Reload offsets from trajectory
 ctime or size or n_atoms did not match
  warnings.warn("Reload offsets from trajectory\n "

In [35]: U1.trajectory[-1]
---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-35-34a3a18136e5> in <module>()
----> 1 U1.trajectory[-1]

~/.virtualenvs/napaanton2/local/lib/python2.7/site-packages/MDAnalysis/coordinates/base.pyc in __getitem__(self, frame)
   1536         if isinstance(frame, numbers.Integral):
   1537             frame = self._apply_limits(frame)
-> 1538             return self._read_frame_with_aux(frame)
   1539         elif isinstance(frame, (list, np.ndarray)):
   1540             if len(frame) != 0 and isinstance(frame[0], (bool, np.bool_)):

~/.virtualenvs/napaanton2/local/lib/python2.7/site-packages/MDAnalysis/coordinates/base.pyc in _read_frame_with_aux(self, frame)
   1568     def _read_frame_with_aux(self, frame):
   1569         """Move to *frame*, updating ts with trajectory, transformations and auxiliary data."""
-> 1570         ts = self._read_frame(frame)  # pylint: disable=assignment-from-no-return
   1571         for aux in self.aux_list:
   1572             ts = self._auxs[aux].update_ts(ts)

~/.virtualenvs/napaanton2/local/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.pyc in _read_frame(self, i)
    240             self._read_offsets(store=True)
    241             self._xdr.seek(i)
--> 242             timestep = self._read_next_timestep()
    243         return timestep
    244 

~/.virtualenvs/napaanton2/local/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.pyc in _read_next_timestep(self, ts)
    249         if ts is None:
    250             ts = self.ts
--> 251         frame = self._xdr.read()
    252         self._frame += 1
    253         self._frame_to_ts(frame, ts)

MDAnalysis/lib/formats/libmdaxdr.pyx in MDAnalysis.lib.formats.libmdaxdr.TRRFile.read()

IOError: TRR read error = float

Actual behavior

Iterating seems to work just fine:

In [36]: for i in U1.trajectory:
    ...:     pass
    ...: 

In [37]: 

Actually, the last (damaged) frame is ignored and i only iterates up to the penultimate frame.

EDITS: @orbeckst added a few more details from conversation with @ianmkenney

orbeckst commented 5 years ago

@ianmkenney I clarified the description a bit, after our discussion.

orbeckst commented 5 years ago

Note that this was found because the ChainReader throws an error while the TRRReader on its own seems to take the IOError as the sign to just end iteration. This is bad behavior because one might never learn that one's trajectory is actually damaged and incomplete.

orbeckst commented 5 years ago

It seems that detecting the end of a TRR is pretty haphazard: https://github.com/MDAnalysis/mdanalysis/blob/f83498ec6b5cc546a0d75598a4d1d6ede6fe2782/package/MDAnalysis/lib/formats/libmdaxdr.pyx#L516-L532

However, the reported error here is float and not integer. Not sure why iteration finishes apparently cleanly.

orbeckst commented 5 years ago

I think the reason why "end of file" always shows up as a "EINTEGER" is because the reading of the magic integer fails https://github.com/MDAnalysis/mdanalysis/blob/f83498ec6b5cc546a0d75598a4d1d6ede6fe2782/package/MDAnalysis/lib/formats/src/xdrfile_trr.c#L72-L73

We could make a special exception that catches that particular case.

hmacdope commented 1 year ago

Possibly related to work in #3892