Make check_oh_consist a proper CLI

tests/test_check_oh_consist.py

@@ -0,0 +1,34 @@
+import importlib.util
+import os
+import subprocess
+import sys
+import unittest
+
+
+class TestCheckOhConsist(unittest.TestCase):
+    def setUp(self):
+        self.repository_root = os.path.abspath("..")
+        self.script = os.path.join(self.repository_root, "tools", "check_oh_consist.py")
+
+    def test_import_has_no_dump_file_side_effect(self):
+        """Loading the tool defines its API without reading dump.hti."""
+        spec = importlib.util.spec_from_file_location("check_oh_consist", self.script)
+        module = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(module)
+
+        self.assertTrue(callable(module.get_oh_distance_stats))
+
+    def test_help_runs_from_repository_root(self):
+        result = subprocess.run(
+            [sys.executable, self.script, "--help"],
+            cwd=self.repository_root,
+            check=True,
+            capture_output=True,
+            text=True,
+        )
+
+        self.assertIn("LAMMPS dump trajectory", result.stdout)
+
+
+if __name__ == "__main__":
+    unittest.main()

tools/check_oh_consist.py

@@ -1,47 +1,53 @@
 #!/usr/bin/env python3
 
-import lib.dump as dump
-import lib.water as water
+import argparse
+import sys
+from pathlib import Path
+
 import numpy as np
-from lib.dump import split_traj
-
-# def func (xx) :
-#     return 0.02*xx*xx+0.01*xx+0.03
-
-# x0 = np.arange(0,10.1)
-# x1 = np.arange(0,10.1, 0.5)
-# x2 = np.arange(0,10.1, 0.25)
-
-
-# i0 = integrate(x0, func(x0), np.zeros(x0.shape))
-# i1 = integrate(x1, func(x1), np.zeros(x1.shape))
-# i2 = integrate(x2, func(x2), np.zeros(x2.shape))
-# e0 = integrate_sys_err(x0, func(x0))
-# e1 = integrate_sys_err(x1, func(x1))
-# e2 = integrate_sys_err(x2, func(x2))
-
-# print(i0[0], e0)
-# print(i1[0], e1)
-# print(i2[0], e2)
-
-# get_thermo('log.lammps')
-
-lines = open("dump.hti").read().split("\n")
-ret = split_traj(lines)
-# print(get_posi(ret[0]))
-# print(get_posi(ret[0])[127:130])
-# print(get_atype(ret[0]))
-# print(get_atype(ret[0])[127:130])
-
-bd, tl = dump.get_dumpbox(ret[0])
-orig, box = dump.dumpbox2box(bd, tl)
-atype = dump.get_atype(ret[0])
-posi = dump.get_posi(ret[0])
-oh_list = water.min_oh_list(box, atype, posi)
-
-for idx, ii in enumerate(ret):
-    bd, tl = dump.get_dumpbox(ii)
-    orig, box = dump.dumpbox2box(bd, tl)
-    posi = dump.get_posi(ii)
-    dists = water.dist_via_oh_list(box, posi, oh_list)
-    print(idx, np.min(dists), np.max(dists), np.average(dists))
+
+# Direct execution sets sys.path to tools/, so add the repository for source-tree use.
+if __package__ in (None, ""):
+    sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+
+from dpti.lib import dump, water
+
+
+def get_oh_distance_stats(dump_path):
+    """Return per-frame O-H distance statistics from a LAMMPS dump trajectory."""
+    trajectories = dump.split_traj(Path(dump_path).read_text().splitlines())
+    if not trajectories:
+        raise ValueError(f"no LAMMPS trajectory frames found in {dump_path}")
+
+    bounds, tilt = dump.get_dumpbox(trajectories[0])
+    _, box = dump.dumpbox2box(bounds, tilt)
+    atom_types = dump.get_atype(trajectories[0])
+    positions = dump.get_posi(trajectories[0])
+    oh_list = water.min_oh_list(box, atom_types, positions)
+
+    stats = []
+    for index, trajectory in enumerate(trajectories):
+        bounds, tilt = dump.get_dumpbox(trajectory)
+        _, box = dump.dumpbox2box(bounds, tilt)
+        positions = dump.get_posi(trajectory)
+        distances = water.dist_via_oh_list(box, positions, oh_list)
+        stats.append(
+            (index, np.min(distances), np.max(distances), np.average(distances))
+        )
+    return stats
+
+
+def main(argv=None):
+    """Parse command-line arguments and print O-H distance statistics."""
+    parser = argparse.ArgumentParser(
+        description="Check O-H bond consistency across a LAMMPS dump trajectory"
+    )
+    parser.add_argument("DUMP", help="LAMMPS dump trajectory to inspect")
+    args = parser.parse_args(argv)
+
+    for values in get_oh_distance_stats(args.DUMP):
+        print(*values)
+
+
+if __name__ == "__main__":
+    main()