FIX: LAMMPSDUMP Convert forces from kcal to kJ units

package/MDAnalysis/coordinates/LAMMPS.py

@@ -615,6 +615,12 @@ class DumpReader(base.ReaderBase):
     **kwargs
        Other keyword arguments used in :class:`~MDAnalysis.coordinates.base.ReaderBase`
 
+    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.
 
     .. versionchanged:: 2.8.0
        Reading of arbitrary, additional columns is now supported.
@@ -631,6 +637,12 @@ class DumpReader(base.ReaderBase):
     """
 
     format = "LAMMPSDUMP"
+    units = {
+        'time': 'fs',
+        'length': 'Angstrom',
+        'velocity': 'Angstrom/fs',
+        'force': 'kcal/(mol*Angstrom)'
+    }
     _conventions = [
         "auto",
         "unscaled",
@@ -910,4 +922,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/coordinates/test_lammps.py

@@ -768,6 +768,100 @@ 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"""
+        from MDAnalysis.coordinates.LAMMPS import DumpReader
+
+        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"""
+        from MDAnalysis import units
+
+        # 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}")
+
 
 @pytest.mark.parametrize(
     "convention", ["unscaled", "unwrapped", "scaled_unwrapped"]