From 9fe966536915e3c43d2af165583d2ab0f117eccb Mon Sep 17 00:00:00 2001 From: Roy Kid Date: Fri, 28 Aug 2026 10:54:37 +0200 Subject: [PATCH] feat: re-export molrs.md as molpy.md plus UnitPreset/VerletSkin molpy.md is a verbatim molrs.md surface. UnitPreset and VerletSkin join the molrs identity re-exports. Lint no longer installs unpublished molrs 0.14. --- .github/workflows/ci.yml | 2 +- .pre-commit-config.yaml | 2 +- src/molpy/__init__.py | 6 ++ src/molpy/md/__init__.py | 8 +++ tests/test_full_import.py | 39 ++++++++++++ tests/test_md/test_driver.py | 60 ++++++++++++++++++ tests/test_md/test_integrators.py | 101 ++++++++++++++++++++++++++++++ tests/test_md/test_maxwell.py | 27 ++++++++ tests/test_md/test_neighbors.py | 66 +++++++++++++++++++ tests/test_md/test_reexport.py | 44 +++++++++++++ tests/test_md/test_units.py | 15 +++++ tests/test_molrs_mirror.py | 16 +++++ tests/test_sink_policy.py | 35 +++++++++++ 13 files changed, 419 insertions(+), 2 deletions(-) create mode 100644 src/molpy/md/__init__.py create mode 100644 tests/test_full_import.py create mode 100644 tests/test_md/test_driver.py create mode 100644 tests/test_md/test_integrators.py create mode 100644 tests/test_md/test_maxwell.py create mode 100644 tests/test_md/test_neighbors.py create mode 100644 tests/test_md/test_reexport.py create mode 100644 tests/test_md/test_units.py create mode 100644 tests/test_molrs_mirror.py create mode 100644 tests/test_sink_policy.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f726d83d..2f92ff43 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,7 +17,7 @@ jobs: with: python-version: "3.14" - uses: astral-sh/setup-uv@v9.0.0 - - run: uv run --extra dev tox -e lint + - run: uv run --no-project --with 'tox>=4.23' --with ruff==0.16.1 --with ty==0.0.65 tox -e lint test: runs-on: ${{ matrix.os }} diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 486f99d7..744350d8 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -22,7 +22,7 @@ repos: hooks: - id: tox-lint name: tox -e lint - entry: uv run --extra dev tox -e lint + entry: uv run --no-project --with 'tox>=4.23' --with ruff==0.16.1 --with ty==0.0.65 tox -e lint language: system pass_filenames: false always_run: true diff --git a/src/molpy/__init__.py b/src/molpy/__init__.py index ed9d34d5..adce3d3e 100644 --- a/src/molpy/__init__.py +++ b/src/molpy/__init__.py @@ -24,6 +24,7 @@ data, engine, io, + md, optimize, pack, parser, @@ -42,6 +43,7 @@ "data", "engine", "io", + "md", "optimize", "pack", "parser", @@ -200,9 +202,11 @@ def __dir__() -> list[str]: ScalarObservable, Sphere, Unit, + UnitPreset, UnitRegistry, UnitsError, VectorObservable, + VerletSkin, keys, schema, signal, @@ -402,12 +406,14 @@ def __dir__() -> list[str]: "NeighborList", "NeighborQuery", "Neighbors", + "VerletSkin", "Block", "FRAME_SCHEMA_VERSION", "Frame", "MetaValue", "Quantity", "Unit", + "UnitPreset", "UnitRegistry", "ScalarObservable", "VectorObservable", diff --git a/src/molpy/md/__init__.py b/src/molpy/md/__init__.py new file mode 100644 index 00000000..31eab49c --- /dev/null +++ b/src/molpy/md/__init__.py @@ -0,0 +1,8 @@ +"""molpy.md — the user-facing MD namespace, a verbatim re-export of molrs.md. + +Users spell everything ``molpy.md.``; the objects are identical to +their ``molrs.md`` counterparts. +""" + +from molrs.md import * # noqa: F403 +from molrs.md import __all__ as __all__ # noqa: F401 diff --git a/tests/test_full_import.py b/tests/test_full_import.py new file mode 100644 index 00000000..7b140983 --- /dev/null +++ b/tests/test_full_import.py @@ -0,0 +1,39 @@ +"""Every molpy submodule imports. Enumerated, not a hand list.""" + +from __future__ import annotations + +import importlib +import pkgutil +import warnings + +import molpy + + +def test_walk_packages_imports() -> None: + failures: list[str] = [] + for mod in pkgutil.walk_packages(molpy.__path__, molpy.__name__ + "."): + if mod.name.endswith(".__main__"): + continue + try: + importlib.import_module(mod.name) + except ModuleNotFoundError as exc: + # Optional extras (rdkit, openbabel, …) are not a hard import. + if exc.name in {"rdkit", "openbabel"}: + continue + failures.append(f"{mod.name}: {type(exc).__name__}: {exc}") + except Exception as exc: # noqa: BLE001 — the point is to surface any miss + failures.append(f"{mod.name}: {type(exc).__name__}: {exc}") + assert not failures, "import failures:\n" + "\n".join(failures) + + +def test_import_molpy_and_md_are_silent() -> None: + import sys + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always", FutureWarning) + importlib.reload(molpy) + sys.modules.pop("molpy.md", None) + sys.modules.pop("molrs.md", None) + importlib.import_module("molpy.md") + fw = [w for w in caught if issubclass(w.category, FutureWarning)] + assert not fw diff --git a/tests/test_md/test_driver.py b/tests/test_md/test_driver.py new file mode 100644 index 00000000..93d3dac6 --- /dev/null +++ b/tests/test_md/test_driver.py @@ -0,0 +1,60 @@ +"""Tests for the ForceField + Frame MD driver (molpy.md.MD).""" + +import numpy as np +import pytest + +import molrs +from molpy.md import MD + + +def _bond_frame() -> tuple[object, object]: + ff = molrs.ff.ForceField("bond-only") + ff.def_bondtype("harmonic", "CT", "CT", {"k": 300.0, "r0": 1.5}) + + frame = molrs.Frame() + atoms = molrs.Block() + atoms.insert("x", np.array([0.0, 2.0])) + atoms.insert("y", np.array([0.0, 0.0])) + atoms.insert("z", np.array([0.0, 0.0])) + atoms.insert("mass", np.array([12.0, 12.0])) + frame["atoms"] = atoms + bonds = molrs.Block() + bonds.insert("atomi", np.array([0], dtype=np.uint64)) + bonds.insert("atomj", np.array([1], dtype=np.uint64)) + bonds.insert("type", np.array(["CT-CT"], dtype=str)) + frame["bonds"] = bonds + return ff, frame + + +def test_set_potential_runs_bonded_dimer(): + ff, frame = _bond_frame() + pots = ff.to_potentials(frame) + state = MD().set_potential(pots).run(frame, 20, dt=0.1) + assert state.pos.shape == (2, 3) + assert np.all(np.isfinite(state.pos)) + + +def test_set_forcefield_compiles_per_run(): + ff, frame = _bond_frame() + driver = MD().set_forcefield(ff) + first = driver.run(frame, 20, dt=0.1) + assert np.isfinite(first.energy) + second = driver.run(frame, 20, dt=0.1) + assert np.isfinite(second.energy) + + +def test_set_forcefield_returns_self(): + ff, _frame = _bond_frame() + driver = MD() + assert driver.set_forcefield(ff) is driver + + +def test_md_requires_forcefield(): + with pytest.raises(RuntimeError, match="set_forcefield"): + MD().run(molrs.Frame(), 1, dt=0.01) + + +def test_thermo_requires_kb(): + ff, frame = _bond_frame() + with pytest.raises(ValueError, match="kb="): + MD().set_forcefield(ff).run(frame, 1, dt=0.01, thermo=1) diff --git a/tests/test_md/test_integrators.py b/tests/test_md/test_integrators.py new file mode 100644 index 00000000..1128dc26 --- /dev/null +++ b/tests/test_md/test_integrators.py @@ -0,0 +1,101 @@ +"""Tests for VelocityVerlet / Langevin — constructor-owned LJCut + VerletSkin.""" + +import numpy as np +import pytest + +from molpy import Box, NeighborList, VerletSkin +from molpy.md import LJCut, Langevin, MD, VelocityVerlet + + +def _dimer(*, skin: float = 0.3, rc: float = 2.5): + """Two atoms inside cutoff; returns pos, lj, skin_nl, mass.""" + pos = np.array([[0.0, 0.0, 0.0], [1.1, 0.0, 0.0]], dtype=np.float64) + box = Box.cubic(20.0) + nl = VerletSkin(NeighborList(rc + skin), rc, pos, box, skin=skin) + lj = LJCut(1.0, 1.0, rc, shifted=True) + mass = np.ones(2, dtype=np.float64) + return pos, lj, nl, mass + + +def _nve(dt: float = 0.01) -> tuple[np.ndarray, VelocityVerlet]: + pos, lj, nl, mass = _dimer() + return pos, VelocityVerlet(dt, potential=lj, neighbors=nl, mass=mass) + + +def _langevin(**kw) -> tuple[np.ndarray, Langevin]: + dt = kw.pop("dt", 0.05) + gamma = kw.pop("gamma", 2.0) + kbt = kw.pop("kbt", 1.5) + seed = kw.pop("seed", 0) + mass = kw.pop("mass", None) + pos, lj, nl, default_mass = _dimer() + if mass is None: + mass = default_mass + else: + mass = np.atleast_1d(np.asarray(mass, dtype=np.float64)) + if mass.size == 1: + mass = np.full(2, float(mass[0])) + return pos, Langevin( + dt, gamma=gamma, kbt=kbt, potential=lj, neighbors=nl, mass=mass, seed=seed + ) + + +def test_langevin_constants_match_closed_form(): + dt, gamma, kbt, mass = 0.05, 2.0, 1.5, 2.0 + _, ig = _langevin(dt=dt, gamma=gamma, kbt=kbt, mass=mass) + assert ig.c1 == pytest.approx(np.exp(-gamma * dt)) + assert ig.c2 == pytest.approx(np.sqrt(1.0 - np.exp(-2.0 * gamma * dt))) + assert float(ig.sigma[0, 0]) == pytest.approx(np.sqrt(kbt / mass)) + assert float(ig.inv_mass[0, 0]) == pytest.approx(1.0 / mass) + + +def test_langevin_rejects_gamma_zero(): + pos, lj, nl, mass = _dimer() + with pytest.raises(ValueError, match="VelocityVerlet"): + Langevin(0.01, gamma=0.0, kbt=1.0, potential=lj, neighbors=nl, mass=mass) + + +def test_removed_dof_follows_the_scheme(): + _, nve = _nve() + _, lgv = _langevin(dt=0.01, gamma=2.0, kbt=1.0) + assert nve.removed_dof == 3 + assert lgv.removed_dof == 0 + + +def test_mass_must_be_positive(): + pos, lj, nl, _ = _dimer() + with pytest.raises(ValueError, match="strictly positive"): + VelocityVerlet(0.01, potential=lj, neighbors=nl, mass=-1.0) + pos, lj, nl, _ = _dimer() + with pytest.raises(ValueError, match="strictly positive"): + VelocityVerlet(0.01, potential=lj, neighbors=nl, mass=np.array([1.0, -2.0])) + + +def test_non_double_dtype_is_reserved_on_the_driver(): + with pytest.raises(ValueError, match="float64"): + MD(dtype=np.float32) + + +def test_advance_n_matches_manual_advance_loop(): + pos0 = np.array([[0.0, 0.0, 0.0], [1.1, 0.0, 0.0]], dtype=np.float64) + vel0 = np.array([[0.01, 0.0, 0.0], [-0.01, 0.0, 0.0]], dtype=np.float64) + _, a = _langevin(dt=0.05, gamma=3.0, kbt=1.0, mass=1.0, seed=11) + end_a = a.advance_n(a.initial(pos0.copy(), vel0.copy()), 5) + _, b = _langevin(dt=0.05, gamma=3.0, kbt=1.0, mass=1.0, seed=11) + state = b.initial(pos0.copy(), vel0.copy()) + for _ in range(5): + state = b.advance(state) + np.testing.assert_array_equal(end_a.pos, state.pos) + np.testing.assert_array_equal(end_a.vel, state.vel) + + +def test_nve_force_changes_when_atoms_move_with_skin(): + """Skin>0 must not freeze forces on live geometry (stale sorted_pos bug).""" + pos, lj, nl, mass = _dimer(skin=1.0) + ig = VelocityVerlet(0.01, potential=lj, neighbors=nl, mass=mass) + vel = np.zeros_like(pos) + s0 = ig.initial(pos, vel) + pos1 = pos.copy() + pos1[1, 0] += 0.05 + s1 = ig.initial(pos1, vel) + assert not np.allclose(s0.forces, s1.forces) diff --git a/tests/test_md/test_maxwell.py b/tests/test_md/test_maxwell.py new file mode 100644 index 00000000..32d0f1f1 --- /dev/null +++ b/tests/test_md/test_maxwell.py @@ -0,0 +1,27 @@ +"""Tests for MaxwellBoltzmann.""" + +import numpy as np +import pytest + +from molpy.md import MaxwellBoltzmann + + +def test_same_seed_is_reproducible(): + pos = np.zeros((8, 3)) + mass = np.ones(8) + a = MaxwellBoltzmann(300.0, seed=7).velocities(pos, mass) + b = MaxwellBoltzmann(300.0, seed=7).velocities(pos, mass) + np.testing.assert_array_equal(a, b) + + +def test_remove_com_leaves_zero_com(): + pos = np.zeros((6, 3)) + mass = np.full(6, 2.0) + vel = MaxwellBoltzmann(200.0, seed=1).velocities(pos, mass) + com = (mass.reshape(-1, 1) * vel).sum(0) / mass.sum() + np.testing.assert_allclose(com, 0.0, atol=1e-12) + + +def test_rejects_nonpositive_kbt(): + with pytest.raises(ValueError, match="strictly positive"): + MaxwellBoltzmann(0.0) diff --git a/tests/test_md/test_neighbors.py b/tests/test_md/test_neighbors.py new file mode 100644 index 00000000..4c88860e --- /dev/null +++ b/tests/test_md/test_neighbors.py @@ -0,0 +1,66 @@ +"""Tests for MD using the core :class:`molpy.NeighborList`.""" + +import numpy as np + +from molpy import Box, NeighborList +from molpy.md import LJCut + + +def test_md_neighborlist_is_the_core_engine(): + import molrs + + assert NeighborList is molrs.NeighborList + + +def test_pair_inside_cutoff_is_half_shell(): + nl = NeighborList(2.5) + nl.build(np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]), Box.cubic(20.0)) + neigh = nl.neighbors() + assert neigh.n_pairs == 1 + pairs = set(zip(neigh.query_point_indices(), neigh.point_indices(), strict=True)) + assert pairs == {(0, 1)} + + +def test_pair_outside_cutoff_is_absent(): + nl = NeighborList(1.0) + nl.build(np.array([[0.0, 0.0, 0.0], [1.5, 0.0, 0.0]]), Box.cubic(20.0)) + assert nl.neighbors().n_pairs == 0 + + +def test_lj_flags_bake_the_kernel(): + cut = LJCut(1.0, 1.0, 2.5, shifted=False, smeared=False) + assert cut.n == 12 and cut.m == 6 + assert not cut.shifted + assert not cut.smeared + shifted = LJCut(1.0, 1.0, 2.5, shifted=True) + assert shifted.shifted and not shifted.smeared + smeared = LJCut(1.0, 1.0, 2.5, smeared=True) + assert smeared.smeared and smeared.shifted + + +def test_lj_consumes_neighbors_table(): + pos = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) + nl = NeighborList(2.5) + nl.build(pos, Box.cubic(20.0)) + lj = LJCut(1.0, 1.0, 2.5, shifted=True) + energy, forces = lj.eval_table(2, nl.neighbors()) + assert forces.shape == (2, 3) + np.testing.assert_allclose(forces.sum(axis=0), 0.0, atol=1e-12) + e2, f2 = lj.eval_pairs( + 2, + nl.neighbors().query_point_indices(), + nl.neighbors().point_indices(), + nl.neighbors().disp(), + nl.neighbors().dist_sq(), + ) + np.testing.assert_allclose(f2, forces) + np.testing.assert_allclose(e2, energy) + + +def test_update_reindexes_moved_points(): + nl = NeighborList(2.5) + box = Box.cubic(20.0) + nl.build(np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]), box) + assert nl.neighbors().n_pairs == 1 + nl.update(np.array([[0.0, 0.0, 0.0], [3.0, 0.0, 0.0]])) + assert nl.neighbors().n_pairs == 0 diff --git a/tests/test_md/test_reexport.py b/tests/test_md/test_reexport.py new file mode 100644 index 00000000..f760f223 --- /dev/null +++ b/tests/test_md/test_reexport.py @@ -0,0 +1,44 @@ +"""molpy.md is THE user-facing MD namespace — a verbatim molrs.md re-export.""" + +import numpy as np +import pytest + +import molpy.md +import molrs.md + + +def test_every_molrs_md_public_name_is_identical_in_molpy_md(): + """Completeness: molpy.md re-exports the full molrs.md surface verbatim.""" + for name in molrs.md.__all__: + assert hasattr(molpy.md, name), f"molpy.md misses {name}" + assert getattr(molpy.md, name) is getattr(molrs.md, name), ( + f"molpy.md.{name} is not the molrs.md object" + ) + assert name in molpy.md.__all__ + + +def test_potential_subclass_runs_through_the_integrator(): + """``class MyPot(molpy.md.Potential)`` drives VelocityVerlet end to end.""" + + class Tether(molpy.md.Potential): + def __init__(self, k): + self.k = float(k) + self.calls = 0 + + def calc_energy_forces(self, pos): + self.calls += 1 + return 0.5 * self.k * float((pos * pos).sum()), -self.k * pos + + pos = np.array([[0.0, 0.0, 0.0], [1.5, 0.0, 0.0]]) + tether = Tether(1.0e-3) + vv = molpy.md.VelocityVerlet(0.5, potential=tether, mass=np.ones(2)) + state = vv.advance_n(vv.initial(pos, np.zeros_like(pos)), 3) + assert tether.calls > 0 + assert np.isfinite(state.energy) + # The tether pulls the off-origin atom inward: the override's forces acted. + assert abs(float(state.pos[1, 0])) < 1.5 + + +def test_unoverridden_potential_base_raises(): + with pytest.raises((TypeError, NotImplementedError)): + molpy.md.Potential() diff --git a/tests/test_md/test_units.py b/tests/test_md/test_units.py new file mode 100644 index 00000000..e5293faf --- /dev/null +++ b/tests/test_md/test_units.py @@ -0,0 +1,15 @@ +"""UnitPreset is the MD-facing constants view.""" + +from molpy import UnitPreset + + +def test_real_boltzmann_is_positive(): + kb = UnitPreset("real").boltzmann() + assert isinstance(kb, float) + assert kb > 0.0 + + +def test_real_energy_and_time_names(): + p = UnitPreset("real") + assert p.energy() == "kilocalorie_per_mole" + assert p.time() == "femtosecond" diff --git a/tests/test_molrs_mirror.py b/tests/test_molrs_mirror.py new file mode 100644 index 00000000..34346c1e --- /dev/null +++ b/tests/test_molrs_mirror.py @@ -0,0 +1,16 @@ +"""molpy.md is a verbatim re-export of molrs.md.""" + +from __future__ import annotations + +import molpy.md +import molrs.md + + +def test_every_molrs_md_name_is_the_same_object() -> None: + for name in molrs.md.__all__: + assert getattr(molpy.md, name) is getattr(molrs.md, name) + + +def test_molpy_md_adds_no_extra_public_names() -> None: + extra = set(molpy.md.__all__) - set(molrs.md.__all__) + assert not extra, extra diff --git a/tests/test_sink_policy.py b/tests/test_sink_policy.py new file mode 100644 index 00000000..be513e0f --- /dev/null +++ b/tests/test_sink_policy.py @@ -0,0 +1,35 @@ +"""Permanent anti-duplication: molrs-covered capabilities stay in molrs.""" + +from __future__ import annotations + +from pathlib import Path + +_IO = Path(__file__).resolve().parents[1] / "src" / "molpy" / "io" + +# Formats molrs already exposes on the Python surface. molpy readers must +# delegate (contain ``molrs.io`` / ``molrs.ff``), not grow a second parser. +_SUNK = { + "data/pdb.py": "molrs.io", + "data/xyz.py": "molrs.io", + "data/gro.py": "molrs.io", +} + +# molpy-native extensions: molrs does not provide these on the Python surface. +_NATIVE = { + "data/ac.py", + "forcefield/moltemplate.py", + "data/lammps_bond_react.py", + "emit/openmm.py", +} + + +def test_sunk_formats_delegate_to_molrs() -> None: + for rel, token in _SUNK.items(): + text = (_IO / rel).read_text(encoding="utf-8") + assert token in text, f"{rel} must delegate via {token}" + + +def test_native_extensions_are_not_required_to_delegate() -> None: + for rel in _NATIVE: + path = _IO / rel + assert path.exists(), f"native extension missing: {rel}"