IOError: Error reading frame from DCD file

[hainm]

hi,

I had a wrapper for Universe object in pytraj (named TrajectoryMDAnalysisIterator). I've tried AMBER, GROMACS files and everythings works fine. But got above error with DCD file. Below is full error log. Any idea why? thanks

ERROR: test_1 (__main__.Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/pytraj/decorators.py", line 80, in _no_test
    return func(*args, **kwd)
  File "test_MDAnalysisIterator.py", line 76, in test_1
    aa_eq(t[:2].xyz, traj[:2].xyz)
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/pytraj/trajs/TrajectoryMDAnalysisIterator.py", line 53, in __getitem__
    with self:
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/pytraj/trajs/TrajectoryMDAnalysisIterator.py", line 77, in __enter__
    self._traj_holder.rewind()
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/coordinates/base.py", line 626, in rewind
    self.next()
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/coordinates/base.py", line 621, in next
    return self._read_next_timestep()
  File "/home/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/coordinates/DCD.py", line 471, in _read_next_timestep
    ts.frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, self.skip)
IOError: Error reading frame from DCD file

The test script is here
https://github.com/pytraj/pytraj/blob/master/tests/test_MDAnalysisIterator.py#L64

and this is wrapper code for Universer in pytraj (if anyone is interested)
https://github.com/pytraj/pytraj/blob/master/pytraj/trajs/TrajectoryMDAnalysisIterator.py#L36

cheers
Hai

[orbeckst]

Do the unit tests for MDAnalysis itself pass in your installation?

[hainm]

@orbeckst sadly No. This is a part of log file

======================================================================
ERROR: MDAnalysisTests.test_topology.TestDLPolyHistoryParser.test_usage
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/nose-1.3.6-py2.7.egg/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysisTests-0.11.0.dev0-py2.7.egg/MDAnalysisTests/test_topology.py", line 1161, in test_usage
    struc = parser.parse()
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/topology/DLPolyParser.py", line 105, in parse
    with openany(self.filename, 'r') as inf:
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/contextlib.py", line 17, in __enter__
    return self.gen.next()
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/core/util.py", line 252, in openany
    stream = anyopen(datasource, mode=mode, reset=reset)
  File "/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysis-0.11.0.dev0-py2.7-linux-x86_64.egg/MDAnalysis/core/util.py", line 317, in anyopen
    raise IOError(errno.EIO, "Cannot open file or stream in mode=%(mode)r." % vars(), repr(filename))
IOError: [Errno 5] Cannot open file or stream in mode=u'r'.: "'/mnt/raidc2/haichit/anaconda3/envs/python2/lib/python2.7/site-packages/MDAnalysisTests-0.11.0.dev0-py2.7.egg/MDAnalysisTests/data/dlpoly/HISTORY'"

----------------------------------------------------------------------
Ran 1688 tests in 236.351s

FAILED (SKIP=15, errors=59)

[hainm]

this is the full log file if you're interested (uploaded to by box account with txt file, so you can preview it without downloading for your safety :P)

https://app.box.com/s/2musxsfkjyyq5kz3gfionh1lwmz3ceze

btw, I like the test suite and test cases in mdanalysis. Will learn this for pytraj. :D

[richardjgowers]

The tests that are failing there are related to the new DL_Poly stuff, I think I've packaged the test files wrong and you're not getting them with the package for some reason. It's completely unrelated to the DCD module.

Your error seems to be when you try and start an __iter__ it's failing to read the first frame even... Which is strange because I think you're using the same DCD file used in our tests.

[hainm]

The tests that are failing there are related to the new DL_Poly stuff, I think I've packaged the test files wrong and you're not getting them with the package for some reason. It's completely unrelated to the DCD module.

I agree that those errors don't relate to DCD stuff.

Is is a proper way to get update all (package + testsuite)?

~/programs/mdanalysis_git/mdanalysis (develop) $ git pull
Already up-to-date.

Your error seems to be when you try and start an iter it's failing to read the first frame even... Which is strange because I think you're using the same DCD file used in our tests.

Yes, it's weird. Like I said things are fine with Amber and Gromacs files (from mdanalysis). Only happen with DCD. Can you confirm with me that mdanalysis does not support slicing for DCD file? This is only reason I can think of. thanks.

Hai

[hainm]

I think this is the source of problem with DCD file

                try:
                    # if MDAnalysis object support slicing
                    for _ in self._traj_holder[idx]:
                        frame = Frame(self.n_atoms)
                        frame.xyz[:] = atom_groups.positions
                        fa.append(frame, copy=False)
                    return fa
                except:
                    # old fashion way
                    count = start
                    while count < stop:
                        fa.append(self[count], copy=False)
                        count += stride
                    return fa

[hainm]

so I tested with bot DCD and TRR file

    @test_if_having("MDAnalysis")
    def test_0(self):
        # PSF, DCD
        from MDAnalysis import Universe
        from MDAnalysisTests.datafiles import PSF, DCD
        u = Universe(PSF, DCD)

        t = mdio.load_MDAnalysisIterator(u)

        for f in t:
            print (f)

        for f in t:
            print (f)

    @test_if_having("MDAnalysis")
    def test_1(self):
        # GRO, TRR
        from MDAnalysis import Universe
        from MDAnalysisTests.datafiles import GRO, TRR
        u = Universe(GRO, TRR)

        t = mdio.load_MDAnalysisIterator(u)

        for f in t:
            print (f)

        for f in t:
            print (f)

for the case of DCD file, the iteration is hanged out if calling for f in t twice. For TRR file, the result is what I am expecting.

Below is output from TRR file testing

.<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
<Frame with 47681 atoms>
.
----------------------------------------------------------------------
Ran 2 tests in 5.646s

OK

[orbeckst]

The DCDReader supports full-featured slicing and this is also explicitly tested, e.g.

nosetests -v test_coordinates:TestDCDReader.test_slice_dcd

Just using MDAnalysis and the very latest development branch I don't manage to break DCD reading:

import MDAnalysis; from MDAnalysis.tests.datafiles import PSF, DCD
u = MDAnalysis.Universe(PSF,DCD)

Just standard iteration:

frames = [ts.frame for ts in u.trajectory]
print(frames)
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98]

Slicing:

frames2 = [ts.frame for ts in u.trajectory[10:60:5]]
frames3 = [ts.frame for ts in u.trajectory[slice(10,60,5)]]

print(frames2)
# [11, 16, 21, 26, 31, 36, 41, 46, 51, 56]
print(frames3)
# [11, 16, 21, 26, 31, 36, 41, 46, 51, 56]

Reverse slicing

frames4 = [ts.frame for ts in u.trajectory[60:10:-5]]
print(frames4)
# [61, 56, 51, 46, 41, 36, 31, 26, 21, 16]

Can you reproduce the problem using code that only involves MDAnalysis?

[orbeckst]

Is is a proper way to get update all (package + testsuite)?

The git pull should do the job but depending on how you install you might have to rerun the installation of the testsuite as well. If in doubt use maintainer/install_all.sh (see -h for usage or just look at the source ;-) ).

Btw, the reason for having a separate testsuite is that we chose to use real-world trajectories for testing and even though we try to keep them small, all the test data weigh about 20 MiB and we really don't want to add this to the installation of the library itself.

[hainm]

Is is a proper way to get update all (package + testsuite)?
The git pull should do the job but depending on how you install you might have to rerun the installation of the testsuite as well. If in doubt use maintainer/install_all.sh (see -h for usage or just look at the source ;-) ).

thanks. I will try ./install_all.sh to make sure.

Btw, the reason for having a separate testsuite is that we chose to use real-world trajectories for testing and even though we try to keep them small, all the test data weigh about 20 MiB and we really don't want to add this to the installation of the library itself.

I agree with having test-suite in different branch.

[hainm]

@orbeckst I think I know what's going on here. I originally thought that I need to explicitly open/close file in mdanalysis to do iteration or slicing.

This is non-working code (with self):

            with self:
                fa = Trajectory()
                fa.top = self.top
                atom_groups = self._ext_holder.atoms
                start, stop, stride = idx.indices(self.n_frames)

                try:
                    # if MDAnalysis object support slicing
                    for _ in self._traj_holder[idx]:
                        frame = Frame(self.n_atoms)
                        frame.xyz[:] = atom_groups.positions
                        fa.append(frame, copy=False)
                    return fa
                except:
                    # old fashion way
                    count = start
                    while count < stop:
                        fa.append(self[count], copy=False)
                        count += stride
                    return fa

    def __enter__(self):
        self._traj_holder.rewind()
        return self

    def __exit__(self, *args):
        self._traj_holder.close()

if removing with self, there is not error anymore. Thanks.

Did you guys just update this feature in recent months? I remember that I did always get segmentation fault if not using with self context.

Hai

[hainm]

so now we can have more 'talk' between pytraj and mdanalysis

from pytraj import io as mdio
from MDAnalysisTests.datafiles import PSF, DCD
from pytraj import TrajectoryIterator
from pytraj.trajs.TrajectoryMDAnalysisIterator import TrajectoryMDAnalysisIterator
traj = mdio.iterload(DCD, PSF, engine='pytraj')
traj2 = mdio.iterload(DCD, PSF, engine='mdanalysis')

[richardjgowers]

Readers should work as context managers, so you can do things like this

with DCDReader('adk_dims.dcd') as r:
    for ts in r[::10]:
         print ts.frame

But some readers will need to be given a numatoms kwarg to know the number of atoms in the trajectory.

The trajectories get opened in the __init__ of the Reader, so they should start with the first frame read, so you shouldn't have to do extra work there.

[richardjgowers]

I think I've managed to recreate this error:

r = DCDReader(DCD)
with r:
    print r.ts
r.rewind()  # IOError
r.next()  # IOError

So it looks like after using it as a context, it closes the file handle, then further use of it doesn't work.