LAMMPS DumpReader: unit overrides (time/length/energy), default force conversion, and native-force preservation test

Author: Rupesh-Singh-Karki

package/CHANGELOG

@@ -42,6 +42,10 @@ Fixes
    directly passing them. (Issue #3520, PR #5006)
 
 Enhancements
+ * LAMMPS DumpReader now converts forces to MDAnalysis base units by default
+   and supports overriding native units via `timeunit`, `lengthunit`, and
+   `energyunit` kwargs to accommodate different LAMMPS unit styles such as
+   "real" and "metal" (Issue #5115, PR #5117)
  * Added support for reading and processing streamed data in `coordinates.base`
    with new `StreamFrameIteratorSliced` and `StreamReaderBase` (Issue #4827, PR #4923)
  * New coordinate reader: Added `IMDReader` for reading real-time streamed

package/MDAnalysis/coordinates/LAMMPS.py

@@ -615,16 +615,41 @@ class DumpReader(base.ReaderBase):
     **kwargs
        Other keyword arguments used in :class:`~MDAnalysis.coordinates.base.ReaderBase`
 
+    Additional keyword arguments
+    ----------------------------
+    timeunit : str, optional
+        Native time unit of the dump file (default: ``"fs"`` for LAMMPS
+        "real" units). Must be a valid MDAnalysis time unit.
+    lengthunit : str, optional
+        Native length unit of the dump file (default: ``"Angstrom"`` for
+        LAMMPS "real" units). Must be a valid MDAnalysis length unit.
+    energyunit : str, optional
+        Native energy unit per mole of the dump file (default:
+        ``"kcal/mol"`` for LAMMPS "real" units). Used together with
+        `lengthunit` to derive the native force unit (energy/length).
+
     Note
     ----
-    This reader assumes LAMMPS "real" units where time is in femtoseconds,
-    length is in Angstroms, velocities in Angstrom/femtosecond, and forces
-    in kcal/(mol*Angstrom). Forces are automatically converted to MDAnalysis
-    base units (kJ/(mol*Angstrom)) for consistency with other trajectory formats.
+    By default, this reader assumes LAMMPS "real" units where time is in
+    femtoseconds (``fs``), length is in Angstroms (``Angstrom``), velocities
+    in ``Angstrom/fs``, and forces in ``kcal/(mol*Angstrom)``. You can
+    override the native units with `timeunit`, `lengthunit`, and `energyunit`.
+    The velocity unit is derived as ``lengthunit/timeunit``; the force unit
+    is derived as ``energyunit/lengthunit`` (preserving ``/mol`` if present).
+
+    When ``convert_units=True`` (default), positions, velocities, and forces
+    are converted from the native units to MDAnalysis base units, i.e.,
+    positions to ``Angstrom``, time to ``ps`` (affecting velocity), and
+    forces to ``kJ/(mol*Angstrom)``. This ensures consistency with other
+    trajectory formats.
 
     .. versionchanged:: 2.8.0
        Reading of arbitrary, additional columns is now supported.
        (Issue `#3504 <https://github.com/MDAnalysis/mdanalysis/issues/3504>`__)
+     .. versionchanged:: 2.10.0
+         Forces are converted to MDAnalysis base units by default and native
+         units can be overridden via `timeunit`, `lengthunit`, and `energyunit`
+         to accommodate different LAMMPS unit styles.
     .. versionchanged:: 2.4.0
        Now imports velocities and forces, translates the box to the origin,
        and optionally unwraps trajectories with image flags upon loading.
@@ -671,6 +696,7 @@ def __init__(
         additional_columns=None,
         **kwargs,
     ):
+        # Initialize first to set convert_units and ts kwargs
         super(DumpReader, self).__init__(filename, **kwargs)
 
         root, ext = os.path.splitext(self.filename)
@@ -685,6 +711,82 @@ def __init__(
                 f"Please choose one of {option_string}"
             )
 
+        # Allow overriding native units per-instance to support LAMMPS unit styles
+        # Defaults correspond to "real": time=fs, length=Angstrom, energy=kcal/mol
+        timeunit = kwargs.pop("timeunit", None)
+        lengthunit = kwargs.pop("lengthunit", None)
+        energyunit = kwargs.pop("energyunit", None)
+
+        # Start from class defaults
+        self.units = self.units.copy()
+
+        # Helper to (re)compute velocity and force units from time/length/energy
+        def _compute_units(_time, _length, _energy):
+            # velocity as length/time
+            vel = None
+            if _time is not None and _length is not None:
+                vel = f"{_length}/{_time}"
+
+            # force as energy/length (add per-mol if provided in energy)
+            # Examples:
+            #  - energy="kcal/mol", length="Angstrom" -> "kcal/(mol*Angstrom)"
+            #  - energy="eV", length="Angstrom" -> "eV/Angstrom"
+            frc = None
+            if _energy is not None and _length is not None:
+                if "/mol" in _energy:
+                    base_energy = _energy.replace("/mol", "")
+                    frc = f"{base_energy}/(mol*{_length})"
+                else:
+                    frc = f"{_energy}/{_length}"
+            return vel, frc
+
+        # Apply overrides if provided
+        if timeunit is not None:
+            # Validate unit type if known
+            try:
+                if units.unit_types[timeunit] != 'time':
+                    raise TypeError(
+                        f"LAMMPS DumpReader: wrong unit {timeunit!r} for unit type 'time'"
+                    )
+            except KeyError:
+                raise ValueError(
+                    f"LAMMPS DumpReader: unknown time unit {timeunit!r}"
+                ) from None
+            self.units['time'] = timeunit
+
+        if lengthunit is not None:
+            try:
+                if units.unit_types[lengthunit] != 'length':
+                    raise TypeError(
+                        f"LAMMPS DumpReader: wrong unit {lengthunit!r} for unit type 'length'"
+                    )
+            except KeyError:
+                raise ValueError(
+                    f"LAMMPS DumpReader: unknown length unit {lengthunit!r}"
+                ) from None
+            self.units['length'] = lengthunit
+
+        # default energy for "real"
+        default_energy = 'kcal/mol'
+        if energyunit is None:
+            energyunit = default_energy
+        else:
+            try:
+                if units.unit_types[energyunit] != 'energy':
+                    raise TypeError(
+                        f"LAMMPS DumpReader: wrong unit {energyunit!r} for unit type 'energy'"
+                    )
+            except KeyError:
+                # Some compound forms like 'kcal/mol' may not be in unit_types; allow pass-through
+                pass
+
+        # Derive velocity and force units based on final time/length/energy
+        vunit, funit = _compute_units(self.units['time'], self.units['length'], energyunit)
+        if vunit is not None:
+            self.units['velocity'] = vunit
+        if funit is not None:
+            self.units['force'] = funit
+
         self._unwrap = unwrap_images
 
         if (

testsuite/MDAnalysisTests/coordinates/test_lammps.py

@@ -28,6 +28,7 @@
 
 import MDAnalysis as mda
 from MDAnalysis import NoDataError
+from MDAnalysis.lib.util import anyopen
 
 from numpy.testing import assert_equal, assert_allclose
 
@@ -53,6 +54,8 @@
     LAMMPSdata_additional_columns,
     LAMMPSDUMP_additional_columns,
 )
+from MDAnalysis.coordinates.LAMMPS import DumpReader
+from MDAnalysis import units
 
 
 def test_datareader_ValueError():
@@ -770,7 +773,6 @@ def test_warning(self, system, request):
 
     def test_dump_reader_units_attribute(self):
         """Test that DumpReader has proper units defined"""
-        from MDAnalysis.coordinates.LAMMPS import DumpReader
 
         expected_units = {
             "time": "fs",
@@ -783,7 +785,6 @@ def test_dump_reader_units_attribute(self):
 
     def test_force_unit_conversion_factor(self):
         """Test that the force conversion factor is correct"""
-        from MDAnalysis import units
 
         # Get conversion factor from kcal/(mol*Angstrom) to kJ/(mol*Angstrom)
         factor = units.get_conversion_factor(
@@ -864,6 +865,92 @@ def test_force_conversion_consistency_across_frames(self):
                     err_msg=f"Force conversion failed at frame {ts_conv.frame}",
                 )
 
+    def test_native_forces_preserved_first_last_atom(self):
+        """Check that native forces (convert_units=False) match raw dump values
+        for first and last atom (by LAMMPS id) on the last frame.
+
+        This tightens the previous loose check that merely ensured native forces
+        were non-zero, by validating exact preservation of raw values.
+        """
+
+        # Parse last frame forces from the raw dump (keyed by LAMMPS atom id)
+        def parse_last_frame_forces(path):
+            last_forces = None
+            n_atoms = None
+            id_idx = fx_idx = fy_idx = fz_idx = None
+            with anyopen(path) as f:
+                while True:
+                    line = f.readline()
+                    if not line:
+                        break
+                    if line.startswith("ITEM: TIMESTEP"):
+                        # timestep line and number
+                        _ = f.readline()
+                        # number of atoms header and value
+                        assert f.readline().startswith("ITEM: NUMBER OF ATOMS")
+                        n_atoms = int(f.readline().strip())
+                        # box bounds header + 3 lines
+                        assert f.readline().startswith("ITEM: BOX BOUNDS")
+                        f.readline(); f.readline(); f.readline()
+                        # atoms header with columns
+                        atoms_header = f.readline().strip()
+                        assert atoms_header.startswith("ITEM: ATOMS ")
+                        cols = atoms_header.split()[2:]  # after 'ITEM: ATOMS'
+                        # Identify indices for id and fx fy fz
+                        try:
+                            id_idx = cols.index("id")
+                            fx_idx = cols.index("fx")
+                            fy_idx = cols.index("fy")
+                            fz_idx = cols.index("fz")
+                        except ValueError as e:
+                            raise AssertionError(
+                                "Required columns 'id fx fy fz' not found in dump header"
+                            ) from e
+                        # Read this frame's atoms
+                        frame_forces = {}
+                        for _ in range(n_atoms):
+                            parts = f.readline().split()
+                            aid = int(parts[id_idx])
+                            fx = float(parts[fx_idx])
+                            fy = float(parts[fy_idx])
+                            fz = float(parts[fz_idx])
+                            frame_forces[aid] = np.array([fx, fy, fz], dtype=float)
+                        # Keep updating last_forces; at EOF it will be the last frame
+                        last_forces = frame_forces
+            assert last_forces is not None and n_atoms is not None
+            return last_forces, n_atoms
+
+        raw_forces_by_id, n_atoms = parse_last_frame_forces(LAMMPSDUMP_image_vf)
+
+        # Universe with native units preserved
+        u_native = mda.Universe(
+            LAMMPS_image_vf,
+            LAMMPSDUMP_image_vf,
+            format="LAMMPSDUMP",
+            convert_units=False,
+        )
+
+        u_native.trajectory[-1]
+        forces_native = u_native.atoms.forces
+
+        # Determine smallest and largest atom ids present in the frame
+        min_id = min(raw_forces_by_id.keys())
+        max_id = max(raw_forces_by_id.keys())
+
+        # Universe sorts by id, so index 0 corresponds to min_id, and -1 to max_id
+        expected_first = raw_forces_by_id[min_id]
+        expected_last = raw_forces_by_id[max_id]
+
+        # Allow tiny numerical differences due to float32 storage in trajectory
+        assert_allclose(
+            forces_native[0], expected_first, rtol=0, atol=1e-6,
+            err_msg="Native first-atom force does not match raw dump value",
+        )
+        assert_allclose(
+            forces_native[-1], expected_last, rtol=0, atol=1e-6,
+            err_msg="Native last-atom force does not match raw dump value",
+        )
+
 
 @pytest.mark.parametrize(
     "convention", ["unscaled", "unwrapped", "scaled_unwrapped"]