feat: molpy.md re-export plus UnitPreset/VerletSkin

.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 }}

.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

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",

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.<Name>``; the objects are identical to
+their ``molrs.md`` counterparts.
+"""
+
+from molrs.md import *  # noqa: F403
+from molrs.md import __all__ as __all__  # noqa: F401

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

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)

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)

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)

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