diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 8ba7f577350..6d6ad957381 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -215,7 +215,7 @@ def Writer(self, filename, n_atoms=None, **kwargs): def _frame_to_ts(self, frame, ts): """convert a dcd-frame to a :class:`TimeStep`""" - ts.frame = self._frame + ts.frame = self._frame + self._file.header['istart'] ts.time = ts.frame * self.ts.dt ts.data['step'] = self._file.tell() @@ -358,6 +358,7 @@ def __init__(self, dt=1, remarks='', nsavc=1, + istart=0, **kwargs): """Parameters ---------- @@ -380,6 +381,8 @@ def __init__(self, how many MD steps is a frame saved to the DCD). By default, this number is just set to one and this should be sufficient for almost all cases but if required, nsavc can be changed. + istart : int (optional) + starting frame number. CHARMM defaults to 1. **kwargs : dict General writer arguments @@ -400,7 +403,7 @@ def __init__(self, nsavc=nsavc, delta=delta, is_periodic=1, - istart=0) + istart=istart) def write_next_timestep(self, ts): """Write timestep object into trajectory. diff --git a/package/MDAnalysis/lib/formats/include/readdcd.h b/package/MDAnalysis/lib/formats/include/readdcd.h index 78544b0ca7f..153ac4723a5 100644 --- a/package/MDAnalysis/lib/formats/include/readdcd.h +++ b/package/MDAnalysis/lib/formats/include/readdcd.h @@ -735,13 +735,13 @@ int write_dcdheader(fio_fd fd, const char *remarks, int N, fio_write_int32(fd, 0); } fio_write_int32(fd, 84); - fio_write_int32(fd, 164); + fio_write_int32(fd, 244); 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, 244); fio_write_int32(fd, 4); fio_write_int32(fd, N); fio_write_int32(fd, 4); diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index ddda9bf5b24..03b9792f131 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -77,10 +77,20 @@ def test_with_statement(self): err_msg="with_statement: DCDReader reads wrong number of frames") assert_array_equal( frames, - np.arange(0, N), + np.arange(0, N) + 1000, err_msg="with_statement: DCDReader does not read all frames") +@pytest.mark.parametrize('istart', (0, 1, 2, 3)) +def test_write_istart(universe_dcd, tmpdir, istart): + u = universe_dcd + outfile = str(tmpdir.join('test.dcd')) + with mda.Writer(outfile, u.atoms.n_atoms, istart=istart) as w: + w.write(u.atoms) + u = mda.Universe(PSF, outfile) + assert u.trajectory._file.header['istart'] == istart + + class TestDCDWriter(BaseWriterTest): @staticmethod @pytest.fixture() @@ -107,6 +117,7 @@ def test_write_random_unitcell(tmpdir): decimal=5) + ################ # Legacy tests # ################ @@ -118,39 +129,39 @@ def universe_dcd(): def test_rewind(universe_dcd): universe_dcd.trajectory.rewind() - assert universe_dcd.trajectory.ts.frame == 0 + assert universe_dcd.trajectory.ts.frame == 1000 def test_next(universe_dcd): universe_dcd.trajectory.rewind() universe_dcd.trajectory.next() - assert universe_dcd.trajectory.ts.frame == 1 + assert universe_dcd.trajectory.ts.frame == 1001 def test_jump_last_frame(universe_dcd): universe_dcd.trajectory[-1] - assert universe_dcd.trajectory.ts.frame == 97 + assert universe_dcd.trajectory.ts.frame == 1097 @pytest.mark.parametrize("start, stop, step", ((5, 17, 3), (20, 5, -1))) def test_slice(universe_dcd, start, stop, step): frames = [ts.frame for ts in universe_dcd.trajectory[start:stop:step]] - assert_array_equal(frames, np.arange(start, stop, step)) + assert_array_equal(frames, np.arange(start, stop, step) + 1000) @pytest.mark.parametrize("array_like", [list, np.array]) def test_array_like(universe_dcd, array_like): ar = array_like([0, 3, 4, 5]) frames = [ts.frame for ts in universe_dcd.trajectory[ar]] - assert_array_equal(frames, ar) + assert_array_equal(frames, np.array(ar) + 1000) @pytest.mark.parametrize("indices", ([0, 4, 2, 3, 0, 1], [0, 0, 1, 1, 2, 1, 1])) def test_list_indices(universe_dcd, indices): frames = [ts.frame for ts in universe_dcd.trajectory[indices]] - assert_array_equal(frames, indices) + assert_array_equal(frames, np.array(indices) + 1000) @pytest.mark.parametrize( @@ -213,7 +224,7 @@ def test_reader_set_dt(): dt = 100 frame = 3 u = mda.Universe(PSF, DCD, dt=dt) - assert_almost_equal(u.trajectory[frame].time, frame*dt, + assert_almost_equal(u.trajectory[frame].time, (frame + 1000)*dt, err_msg="setting time step dt={0} failed: " "actually used dt={1}".format( dt, u.trajectory._ts_kwargs['dt'])) @@ -221,14 +232,14 @@ def test_reader_set_dt(): err_msg="trajectory.dt does not match set dt") -@pytest.mark.parametrize("ext, decimal", (("dcd", 5), +@pytest.mark.parametrize("ext, decimal", (("dcd", 3), ("xtc", 3))) def test_writer_dt(tmpdir, ext, decimal): dt = 5.0 # set time step to 5 ps universe_dcd = mda.Universe(PSF, DCD, dt=dt) t = universe_dcd.trajectory outfile = str(tmpdir.join("test.{}".format(ext))) - with mda.Writer(outfile, n_atoms=t.n_atoms, dt=dt) as W: + with mda.Writer(outfile, n_atoms=t.n_atoms, dt=dt, istart=t._file.header['istart']) as W: for ts in universe_dcd.trajectory: W.write(universe_dcd.atoms) diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py index 537ca36a03b..3b13e2ffa6a 100644 --- a/testsuite/MDAnalysisTests/formats/test_libdcd.py +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -20,6 +20,7 @@ from collections import namedtuple import os import string +import struct import hypothesis.strategies as strategies from hypothesis import example, given @@ -242,6 +243,22 @@ def test_write_header(tmpdir): assert header['nsavc'] == 10 assert np.allclose(header['delta'], .02) + # we also check the bytes written directly. + with open(testfile, 'rb') as fh: + header_bytes = fh.read() + # check for magic number + assert struct.unpack('i', header_bytes[:4])[0] == 84 + # magic number should be written again before remark section + assert struct.unpack('i', header_bytes[88:92])[0] == 84 + # length of remark section. We hard code this to 244 right now + assert struct.unpack('i', header_bytes[92:96])[0] == 244 + # say we have 3 block of length 80 + assert struct.unpack('i', header_bytes[96:100])[0] == 3 + # after the remark section the length should be reported again + assert struct.unpack('i', header_bytes[340:344])[0] == 244 + # this is a magic number as far as I see + assert struct.unpack('i', header_bytes[344:348])[0] == 4 + def test_write_no_header(tmpdir): fname = str(tmpdir.join('test.dcd'))