Dcd use istart information

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.

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);

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,22 +224,22 @@ 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']))
     assert_almost_equal(u.trajectory.dt, 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)
 

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'))