fixed dt for LAMPPS DCD and setting of dt by user

- LAMMPS default time unit is fs (for unit real, see http://lammps.sandia.gov/doc/units.html)
  but we had it as ps
- Universe(..., dt=<dt>) is working now for DCDReader (and anything derived from it such as
  the LAMMPS DCDReader)

package/CHANGELOG

@@ -15,7 +15,7 @@ The rules for this file:
 ------------------------------------------------------------------------------
 10/??/15
 
-  * 0.12.1 kain88-de
+  * 0.12.1 kain88-de, orbeckst
 
 API Changes
 
@@ -25,6 +25,10 @@ Changes
 
 Fixes
   * Fixed OpenMP detection on Linux/OSX #459
+  * Fixed reading of LAMMPS trajectory times: default unit ought
+    to be fs and not ps
+  * Fixed setting of dt for DCDReader (and LAMMPS DCDReader) with
+    keyword argument Universe(..., dt=<dt>)
 
 10/04/15 kain88-de, richardjgowers, dotsdl, sseyler, orbeckst, jbarnoud
 

package/MDAnalysis/coordinates/DCD.py

@@ -426,7 +426,7 @@ def __init__(self, dcdfilename, **kwargs):
         # Convert delta to ps
         delta = mdaunits.convert(self.delta, self.units['time'], 'ps')
 
-        self._ts_kwargs.update({'dt':self.skip_timestep * delta})
+        self._ts_kwargs.setdefault('dt', self.skip_timestep * delta)
         self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs)
         # Read in the first timestep
         self._read_next_timestep()
@@ -480,7 +480,6 @@ 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, 1)
-
         ts.frame += 1
         return ts
 
@@ -493,7 +492,6 @@ def _read_frame(self, frame):
         self._jump_to_frame(frame)
         ts = self.ts
         ts._frame = self._read_next_frame(ts._x, ts._y, ts._z, ts._unitcell, 1)
-
         ts.frame = frame
         return ts
 

package/MDAnalysis/coordinates/LAMMPS.py

@@ -98,14 +98,17 @@ class DCDReader(DCD.DCDReader):
     """Read a LAMMPS_ DCD trajectory.
 
     The units can be set from the constructor with the keyword
-    arguments *timeunit* and *lengthunit*. The defaults are "ps" and
-    "Angstrom". See :mod:`MDAnalysis.units` for other recognized
-    values.
+    arguments *timeunit* and *lengthunit*. The defaults are "fs" and
+    "Angstrom", corresponding to LAMMPS `units style`_ "**real**". See
+    :mod:`MDAnalysis.units` for other recognized values.
+
+    .. _units style: http://lammps.sandia.gov/doc/units.html
+    ..
     """
     format = "DCD"
 
     def __init__(self, dcdfilename, **kwargs):
-        self.units = {'time': 'ps', 'length': 'Angstrom'}  # must be instance level
+        self.units = {'time': 'fs', 'length': 'Angstrom'}  # must be instance level
         self.units['time'] = kwargs.pop('timeunit', self.units['time'])
         self.units['length'] = kwargs.pop('lengthunit', self.units['length'])
         for unit_type, unit in self.units.items():

testsuite/MDAnalysisTests/coordinates/reference.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from MDAnalysisTests.datafiles import (PDB_small, PDB, PDB_full, LAMMPSdata,
+                                       LAMMPSdata2, LAMMPSdcd2,
                                        LAMMPSdata_mini, PSF_TRICLINIC,
                                        DCD_TRICLINIC, PSF_NAMD_TRICLINIC,
                                        DCD_NAMD_TRICLINIC)
@@ -192,6 +193,17 @@ class RefLAMMPSData(object):
                            ],
                           dtype=np.float32)
 
+class RefLAMMPSDataDCD(object):
+    format = "LAMMPS"
+    topology = LAMMPSdata2
+    trajectory = LAMMPSdcd2
+    n_atoms = 12421
+    n_frames = 5
+    dt = 0.5  # ps per frame
+    mean_dimensions = np.array(
+        [ 50.66186142,  47.18824387,  52.33762741,
+          90.        ,  90.        ,  90.        ], dtype=np.float32)
+
 
 class RefLAMMPSDataMini(object):
     filename = LAMMPSdata_mini
@@ -231,3 +243,4 @@ class RefNAMDtriclinicDCD(object):
     ref_dimensions = np.array([
         [1., 38.426594, 38.393101, 44.759800, 90.000000, 90.000000, 60.028915],
     ])
+

testsuite/MDAnalysisTests/coordinates/test_dcd.py

@@ -134,6 +134,12 @@ def test_volume(self):
                             err_msg="wrong volume for unitcell (no unitcell "
                             "in DCD so this should be 0)")
 
+def test_DCDReader_set_dt(dt=100., frame=3):
+    u = mda.Universe(PSF, DCD, dt=dt)
+    assert_almost_equal(u.trajectory[frame].time, frame*dt,
+                        err_msg="setting time step dt={0} failed: "
+                        "actually used dt={1}".format(
+            dt, u.trajectory._ts_kwargs['dt']))
 
 class TestDCDWriter(TestCase):
     def setUp(self):

testsuite/MDAnalysisTests/coordinates/test_lammps.py

@@ -1,11 +1,13 @@
-import MDAnalysis as mda
+import numpy as np
 
+import MDAnalysis as mda
 
-from numpy.testing import (assert_equal, assert_raises)
+from numpy.testing import (assert_equal, assert_almost_equal, assert_raises)
 from unittest import TestCase
 
 from MDAnalysisTests.coordinates.reference import (RefLAMMPSData,
-                                                   RefLAMMPSDataMini)
+                                                   RefLAMMPSDataMini,
+                                                   RefLAMMPSDataDCD)
 
 
 def test_datareader_VE():
@@ -58,3 +60,46 @@ class TestLammpsData_Coords(_TestLammpsData_Coords, RefLAMMPSData):
 
 class TestLammpsDataMini_Coords(_TestLammpsData_Coords, RefLAMMPSDataMini):
     pass
+
+# need more tests of the LAMMPS DCDReader
+
+class TestLAMPPSDCDReader(TestCase, RefLAMMPSDataDCD):
+    def setUp(self):
+        self.u = mda.Universe(self.topology, self.trajectory,
+                              format=self.format)
+
+    def get_frame_from_end(self, offset):
+        iframe = self.u.trajectory.n_frames - 1 - offset
+        iframe = iframe if iframe > 0 else 0
+        return iframe
+
+    def test_n_atoms(self):
+        assert_equal(self.u.atoms.n_atoms, self.n_atoms)
+
+    def test_n_frames(self):
+        assert_equal(self.u.trajectory.n_frames, self.n_frames)
+
+    def test_dimensions(self):
+        mean_dimensions = np.mean([ts.dimensions for ts in self.u.trajectory],
+                                  axis=0)
+        assert_almost_equal(mean_dimensions, self.mean_dimensions)
+
+    def test_dt(self):
+        assert_almost_equal(self.u.trajectory.dt, self.dt,
+                            err_msg="Time between frames dt is wrong.")
+
+    def test_Timestep_time(self):
+        iframe = self.get_frame_from_end(1)
+        assert_almost_equal(self.u.trajectory[iframe].time,
+                            iframe * self.dt,
+                            err_msg="Time for frame {0} (dt={1}) is wrong.".format(
+                iframe, self.dt))
+
+    def test_LAMMPSDCDReader_set_dt(self, dt=1500.):
+        u = mda.Universe(self.topology, self.trajectory, format=self.format,
+                         dt=dt)
+        iframe = self.get_frame_from_end(1)
+        assert_almost_equal(u.trajectory[iframe].time, iframe*dt,
+                            err_msg="setting time step dt={0} failed: "
+                            "actually used dt={1}".format(
+                dt, u.trajectory._ts_kwargs['dt']))

testsuite/MDAnalysisTests/datafiles.py

@@ -77,6 +77,7 @@
     "mol2_molecules", "mol2_molecule", "mol2_broken_molecule",
     "capping_input", "capping_output", "capping_ace", "capping_nma",
     "LAMMPSdata", "trz4data", "LAMMPSdata_mini",
+    "LAMMPSdata2", "LAMMPSdcd2",
     "unordered_res",  # pdb file with resids non sequential
     "GMS_ASYMOPT",  # GAMESS C1  optimization
     "GMS_SYMOPT",   # GAMESS D4h optimization
@@ -214,12 +215,11 @@
 capping_ace = resource_filename(__name__, "data/capping/ace.pdb")
 capping_nma = resource_filename(__name__, "data/capping/nma.pdb")
 
-trz4data = resource_filename(__name__, "data/datatest.trz")
-LAMMPSdata = resource_filename(__name__, "data/datatest.data")
-LAMMPSdata_mini = resource_filename(__name__, "data/mini.data")
-
-LAMMPSdata2 = resource_filename(__name__, "data/lammps2.data")
-LAMMPSdcd2 = resource_filename(__name__, "data/lammps2.dcd")
+trz4data = resource_filename(__name__, "data/lammps/datatest.trz")
+LAMMPSdata = resource_filename(__name__, "data/lammps/datatest.data")
+LAMMPSdata_mini = resource_filename(__name__, "data/lammps/mini.data")
+LAMMPSdata2 = resource_filename(__name__, "data/lammps/ifabp_apo_100mM.data.bz2")
+LAMMPSdcd2 = resource_filename(__name__, "data/lammps/ifabp_apo_100mM.dcd")
 
 unordered_res = resource_filename(__name__, "data/unordered_res.pdb")