diff --git a/tests/test_check_oh_consist.py b/tests/test_check_oh_consist.py new file mode 100644 index 00000000..6343d063 --- /dev/null +++ b/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() diff --git a/tools/check_oh_consist.py b/tools/check_oh_consist.py index 155c9c0c..c0b60ec5 100644 --- a/tools/check_oh_consist.py +++ b/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()