FIX: LAMMPSDUMP Convert forces from kcal to kJ units

package/CHANGELOG

@@ -48,6 +48,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)
  * Add conversion factor for speed unit A/ns and additional representations (Å/fs, A/µs,
    Å/μs, Å/ms) (PR #5053)
  * Improved performance of `analysis.rdf.InterRDF_s` > 20-fold

package/MDAnalysis/coordinates/LAMMPS.py

@@ -613,10 +613,41 @@ class DumpReader(base.ReaderBase):
        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
+    ----
+    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.
@@ -629,6 +660,12 @@ class DumpReader(base.ReaderBase):
     """
 
     format = "LAMMPSDUMP"
+    units = {
+        "time": "fs",
+        "length": "Angstrom",
+        "velocity": "Angstrom/fs",
+        "force": "kcal/(mol*Angstrom)",
+    }
     _conventions = [
         "auto",
         "unscaled",
@@ -657,6 +694,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)
@@ -671,6 +709,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 (
@@ -907,4 +1021,12 @@ def _read_next_timestep(self):
         # Transform to origin after transformation of scaled variables
         ts.positions -= np.array([xlo, ylo, zlo])[None, :]
 
+        # Convert units if requested
+        if self.convert_units:
+            self.convert_pos_from_native(ts.positions)
+            if self._has_vels:
+                self.convert_velocities_from_native(ts.velocities)
+            if self._has_forces:
+                self.convert_forces_from_native(ts.forces)
+
         return ts

testsuite/MDAnalysisTests/analysis/conftest.py

@@ -181,19 +181,13 @@ def client_DensityAnalysis(request):
     return request.param
 
 
-# MDAnalysis.analysis.lineardensity
-
-
-@pytest.fixture(scope="module", params=params_for_cls(LinearDensity))
-def client_LinearDensity(request):
+@pytest.fixture(scope="module", params=params_for_cls(InterRDF))
+def client_InterRDF(request):
     return request.param
 
 
-# MDAnalysis.analysis.polymer
-
-
-@pytest.fixture(scope="module", params=params_for_cls(PersistenceLength))
-def client_PersistenceLength(request):
+@pytest.fixture(scope="module", params=params_for_cls(InterRDF_s))
+def client_InterRDF_s(request):
     return request.param
 
 

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():
@@ -768,6 +771,186 @@ def test_warning(self, system, request):
         with pytest.warns(match="Some of the additional"):
             request.getfixturevalue(system)
 
+    def test_dump_reader_units_attribute(self):
+        """Test that DumpReader has proper units defined"""
+
+        expected_units = {
+            "time": "fs",
+            "length": "Angstrom",
+            "velocity": "Angstrom/fs",
+            "force": "kcal/(mol*Angstrom)",
+        }
+
+        assert DumpReader.units == expected_units
+
+    def test_force_unit_conversion_factor(self):
+        """Test that the force conversion factor is correct"""
+
+        # Get conversion factor from kcal/(mol*Angstrom) to kJ/(mol*Angstrom)
+        factor = units.get_conversion_factor(
+            "force",
+            "kcal/(mol*Angstrom)",  # from (LAMMPS native)
+            "kJ/(mol*Angstrom)",  # to (MDAnalysis base)
+        )
+
+        expected_factor = 4.184  # 1 kcal = 4.184 kJ
+        assert_allclose(factor, expected_factor, rtol=1e-6)
+
+    def test_force_conversion_with_image_vf_file(self):
+        """Test force unit conversion using existing test file with forces"""
+        # Test with convert_units=True (default)
+        u_converted = mda.Universe(
+            LAMMPS_image_vf, LAMMPSDUMP_image_vf, format="LAMMPSDUMP"
+        )
+
+        # Test with convert_units=False
+        u_native = mda.Universe(
+            LAMMPS_image_vf,
+            LAMMPSDUMP_image_vf,
+            format="LAMMPSDUMP",
+            convert_units=False,
+        )
+
+        # Both should have forces
+        assert hasattr(u_converted.atoms, "forces")
+        assert hasattr(u_native.atoms, "forces")
+
+        # Go to last frame where we know forces exist
+        u_converted.trajectory[-1]
+        u_native.trajectory[-1]
+
+        forces_converted = u_converted.atoms.forces
+        forces_native = u_native.atoms.forces
+
+        # Check that forces are different (converted vs native)
+        # The conversion factor should be 4.184
+        expected_factor = 4.184
+        forces_expected = forces_native * expected_factor
+
+        # Test that converted forces match expected values
+        assert_allclose(forces_converted, forces_expected, rtol=1e-6)
+
+        # Test that native forces are unchanged when convert_units=False
+        # Just check they are reasonable values (not zero everywhere)
+        assert not np.allclose(forces_native, 0.0)
+
+    def test_force_conversion_consistency_across_frames(self):
+        """Test that force conversion works consistently across all frames"""
+        u_converted = mda.Universe(
+            LAMMPS_image_vf, LAMMPSDUMP_image_vf, format="LAMMPSDUMP"
+        )
+
+        u_native = mda.Universe(
+            LAMMPS_image_vf,
+            LAMMPSDUMP_image_vf,
+            format="LAMMPSDUMP",
+            convert_units=False,
+        )
+
+        conversion_factor = 4.184
+
+        # Test conversion in all frames
+        for ts_conv, ts_native in zip(
+            u_converted.trajectory, u_native.trajectory
+        ):
+            if ts_conv.has_forces:
+                forces_converted = ts_conv.forces
+                forces_native = ts_native.forces
+                forces_expected = forces_native * conversion_factor
+
+                assert_allclose(
+                    forces_converted,
+                    forces_expected,
+                    rtol=1e-6,
+                    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"]