diff --git a/package/CHANGELOG b/package/CHANGELOG index e1768284d03..370000f0aa2 100644 --- a/package/CHANGELOG +++ b/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 diff --git a/package/MDAnalysis/coordinates/LAMMPS.py b/package/MDAnalysis/coordinates/LAMMPS.py index 6e0043321d3..d4d175d6bc6 100644 --- a/package/MDAnalysis/coordinates/LAMMPS.py +++ b/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 `__) + .. 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 diff --git a/testsuite/MDAnalysisTests/analysis/conftest.py b/testsuite/MDAnalysisTests/analysis/conftest.py index 7918d25f137..61ff0614ce5 100644 --- a/testsuite/MDAnalysisTests/analysis/conftest.py +++ b/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 diff --git a/testsuite/MDAnalysisTests/coordinates/test_lammps.py b/testsuite/MDAnalysisTests/coordinates/test_lammps.py index 3b203f5bae2..e4e3024cf7b 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_lammps.py +++ b/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"]