Fix lammps

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

@@ -192,6 +192,7 @@ class DCDWriter(base.Writer):
     .. _Issue 187: https://github.com/MDAnalysis/mdanalysis/issues/187
     """
     format = 'DCD'
+    flavor = 'CHARMM'
     units = {'time': 'AKMA', 'length': 'Angstrom'}
 
     def __init__(self, filename, n_atoms, start=0, step=1,
@@ -399,6 +400,7 @@ class DCDReader(base.Reader):
        Removed skip keyword and functionality
     """
     format = 'DCD'
+    flavor = 'CHARMM'
     units = {'time': 'AKMA', 'length': 'Angstrom'}
     _Timestep = Timestep
 
@@ -426,7 +428,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 +482,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 +494,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
 
@@ -592,7 +592,8 @@ def Writer(self, filename, **kwargs):
 
     @property
     def dt(self):
-        return self.skip_timestep * self.convert_time_from_native(self.delta)
+        """Time between two trajectory frames in picoseconds."""
+        return self.ts.dt
 
 
 DCDReader._read_dcd_header = types.MethodType(_dcdmodule.__read_dcd_header, None, DCDReader)

package/MDAnalysis/coordinates/LAMMPS.py

@@ -22,32 +22,52 @@
 trajectories. Trajectories can be read regardless of system-endianness
 as this is auto-detected.
 
-LAMMPS can write DCD_ trajectories but unlike a `CHARMM trajectory`_
+LAMMPS can `write DCD`_ trajectories but unlike a `CHARMM trajectory`_
 (which is often called a DCD even though CHARMM itself calls them
 "trj") the time unit is not fixed to be the AKMA_ time unit (20 AKMA
 is 0.978 picoseconds or 1 AKMA = 4.888821e-14 s) but can depend on
-settings in LAMMPS. The most common case appears to be that the time
-step is actually recorded in picoseconds. Other cases are unit-less
-Lennard-Jones time units.
+settings in LAMMPS. The most common case for biomolecular simulations
+appears to be that the time step is recorded in femtoseconds (command
+`units real`_ in the input file) and lengths in ångströms. Other cases
+are unit-less Lennard-Jones time units.
 
 This presents a problem for MDAnalysis because it cannot autodetect
 the unit from the file. By default we are assuming that the unit for
-length is the ångström and for the time is the picosecond. If this is
-not true then the user *should supply the appropriate units* in the as
+length is the ångström and for the time is the femtosecond. If this is
+not true then the user *should supply the appropriate units* in the
 keywords *timeunit* and/or *lengthunit* to :class:`DCDWriter` and
-:class:`DCDReader`.
+:class:`~MDAnalysis.core.AtomGroup.Universe` (which then calls
+:class:`DCDReader`).
+
+.. Rubric:: Example: Loading a LAMMPS simulation
+
+To load a LAMMPS simulation from a LAMMPS data file (using the
+:class:`~MDAnalysis.topology.LAMMPSParser.DATAParser`) together with a
+LAMMPS DCD with "*real*" provide the keyword *format="LAMMPS*"::
+
+  u = MDAnalysis.Universe("lammps.data", "lammps_real.dcd", format="LAMMPS")
+
+If the trajectory uses *units nano* then use ::
+
+  u = MDAnalysis.Universe("lammps.data", "lammps_nano.dcd", format="LAMMPS",
+                          lengthunit="nm", timeunit="ns")
 
 .. Note::
 
-   Lennard-Jones units are not implemented. See
-   :mod:`MDAnalysis.units` for other recognized values.
+   Lennard-Jones units are not implemented. See :mod:`MDAnalysis.units` for
+   other recognized values and the documentation for the LAMMPS `units
+   command`_.
 
-.. SeeAlso:: For further discussion follow the reports for `Issue 84`_ and `Issue 64`_.
+.. SeeAlso::
+
+   For further discussion follow the reports for `Issue 84`_ and `Issue 64`_.
 
 .. _LAMMPS: http://lammps.sandia.gov/
-.. _DCD: http://lammps.sandia.gov/doc/dump.html
+.. _write DCD: http://lammps.sandia.gov/doc/dump.html
 .. _CHARMM trajectory: http://www.charmm.org/documentation/c36b1/dynamc.html#%20Trajectory
 .. _AKMA: http://www.charmm.org/documentation/c36b1/usage.html#%20AKMA
+.. _units real: http://lammps.sandia.gov/doc/units.html
+.. _units command: http://lammps.sandia.gov/doc/units.html
 .. _`Issue 64`: https://github.com/MDAnalysis/mdanalysis/issues/64
 .. _`Issue 84`: https://github.com/MDAnalysis/mdanalysis/issues/84
 
@@ -79,10 +99,11 @@ class DCDWriter(DCD.DCDWriter):
     "Angstrom". See :mod:`MDAnalysis.units` for other recognized
     values.
     """
-    format = "DCD"
+    format = 'DCD'
+    flavor = 'LAMMPS'
 
     def __init__(self, *args, **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():
@@ -98,14 +119,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"
+    format = 'DCD'
+    flavor = 'LAMMPS'
 
     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,14 @@ 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']))
+    assert_almost_equal(u.trajectory.dt, dt,
+                        err_msg="trajectory.dt does not match set dt")
 
 class TestDCDWriter(TestCase):
     def setUp(self):