From 067deb8bbe0188e2bb10d65829eb1fefe35feffc Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:35:13 -0500 Subject: [PATCH 01/25] Report lost saved frames on every run, and name the file when a segment cannot be read A deleted or empty segment went unnoticed: the run ended completed with exit 0 and the frames gone. A cut tail stopped every resume with a raw JSONDecodeError, and a cut earlier segment failed the run for good at its next refinement. Now run warns and its last line names the missing frames, a cut current segment is refused with the file and the way out, and refinement skips a segment it cannot read. MD continues, since lost frames cannot change the exact state in the checkpoint. --- src/reactionflow/cli.py | 5 ++- src/reactionflow/run.py | 88 +++++++++++++++++++++++++++++++++++++---- 2 files changed, 85 insertions(+), 8 deletions(-) diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 35f47e7..09ccce1 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -25,7 +25,7 @@ from .mlip import load_mlip_adapter from .model_setup import install_and_prepare, model_catalog, prepare_adapter, prepare_campaign from .pathway import check_saddle -from .run import ReactionRun, ReactionRunConfig, RunSummary +from .run import ReactionRun, ReactionRunConfig, RunSummary, missing_frames from .scan import ScanResult, file_frames, recorded_occurrences, run_frames, scan_frames from .status import campaign_status, format_status from .wizard import create_campaign @@ -362,6 +362,9 @@ def _run_and_record(campaign: CampaignConfig, arguments: argparse.Namespace) -> finally: temporary.unlink(missing_ok=True) raise + lost = missing_frames(campaign.output_root / trajectory_id, summary.global_frame) + if lost: + payload["missing_frames"] = lost print(json.dumps({**payload, **asdict(summary)}, sort_keys=True)) diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index 1a36d02..e53f26a 100644 --- a/src/reactionflow/run.py +++ b/src/reactionflow/run.py @@ -108,6 +108,50 @@ def _last_frame(path: Path) -> int: return int(trajectory[-1].info["reactionflow_global_frame"]) +def _readable_span(path: Path) -> tuple[int, int]: + """The first global frame of a segment and the last one left to read in it.""" + + with Trajectory(path) as trajectory: + first = int(trajectory[0].info["reactionflow_global_frame"]) + for index in range(len(trajectory) - 1, 0, -1): + try: + return first, int(trajectory[index].info["reactionflow_global_frame"]) + except Exception: # a cut tail + continue + return first, first + + +def missing_frames(root: Path, last_frame: int) -> list[list[int]]: + """Ranges of global frames up to ``last_frame`` that no segment of this run holds. + + Each frame is on disk before the checkpoint that counts it, except the very first frame of a + run, so a healthy run reports nothing however it was stopped. This reads each segment's ends, + so it finds a segment that is gone or cut; damage inside one shows only when ``scan`` reads it. + """ + + if last_frame < 1: + return [] + segments = [ + path for path in root.glob("segments/*/trajectory.traj") if path.parent.name.isdigit() + ] + missing, expected = [], 0 + for path in sorted(segments, key=lambda path: int(path.parent.name)): + try: + span = _readable_span(path) + except Exception: # an empty or cut file raises one of several types + continue + if span[0] > expected: + missing.append([expected, span[0] - 1]) + expected = max(expected, span[1] + 1) + if expected <= last_frame: + missing.append([expected, last_frame]) + return missing + + +def frame_ranges(ranges: list[list[int]]) -> str: + return ", ".join(str(a) if a == b else f"{a}-{b}" for a, b in ranges) + + def _same_atomic_state(first: Atoms, second: Atoms) -> bool: return ( set(first.arrays) == set(second.arrays) @@ -409,11 +453,25 @@ def _run_segment(self, total_steps: int, runtime_provider: ExactRuntimeProvider) ensure_directory(path.parent) # Each frame is written before its checkpoint, so an interrupted attempt may already # hold this boundary's frame or the next one, which the exact replay reproduces. - written = _last_frame(path) + way_out = ( + "Restore it from a copy, or move it aside to continue in a new file without its " + "frames, which every later run then reports as missing" + ) + try: + written = _last_frame(path) + except Exception as error: # a cut file raises one of several types + raise ValueError( + f"{path} cannot be read, so this run cannot append to it " + f"({type(error).__name__}: {error}). {way_out}" + ) from error if written > self._global_frame + 1: raise ValueError("trajectory is ahead of its exact runtime checkpoint") if 0 <= written < self._global_frame: - raise ValueError("trajectory is missing frames before its exact runtime checkpoint") + raise ValueError( + f"trajectory is missing frames before its exact runtime checkpoint: {path} " + f"ends at frame {written}, and the checkpoint is at {self._global_frame}. " + f"{way_out}" + ) with Trajectory(path, "w" if written < 0 else "a") as trajectory: self._write_frame(trajectory, runtime.atoms, written) sync_directory(path.parent) @@ -463,13 +521,21 @@ def _saved_frames(self, first: int, last: int) -> dict[int, Atoms]: for path in sorted(self.root.glob("segments/*/trajectory.traj")): if path.stat().st_size == 0: continue - with Trajectory(path) as trajectory: - if len(trajectory) == 0: + try: + trajectory = Trajectory(path) + except Exception: # cut, which run() reports: refinement can do without it + continue + with trajectory: + try: + # A generation holds one frame per observation, from the one it resumed at. + start = int(trajectory[0].info["reactionflow_global_frame"]) + except Exception: continue - # A generation holds one frame per observation, from the one it resumed at. - start = int(trajectory[0].info["reactionflow_global_frame"]) for frame in range(max(first, start), min(last, start + len(trajectory) - 1) + 1): - atoms = trajectory[frame - start] + try: + atoms = trajectory[frame - start] + except Exception: + break if int(atoms.info["reactionflow_global_frame"]) != frame: raise ValueError(f"{path} does not hold one frame per observation") frames.setdefault(frame, canonical_copy(atoms)) @@ -641,6 +707,14 @@ def run( elif atoms is not None: raise ValueError("initial atoms may only be supplied to a new run") self._refine_after_md = refine_after_md + lost = missing_frames(self.root, self._global_frame) + if lost: + warnings.warn( + f"saved frames {frame_ranges(lost)} of {self.root.name} are missing or unreadable. " + "The MD state comes from the checkpoint and is not affected, but scan cannot " + "replay those frames and refinement cannot use them", + stacklevel=2, + ) if self._phase == "refining" or self._refine_due(): # A reopened adapter has not checked its model against the saved MD state yet. From 3d16a344f47f496a71e3757298c8e500ecb02f8f Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:36:58 -0500 Subject: [PATCH 02/25] Scan past a lost segment and name the missing frames A deleted first segment gave frames 0 with exit 0 and no message, and a missing middle segment hid every later frame. A cut offset table raised a raw ValueError, a trajectory with results and no segments said not started, and segments without state.json raised FileNotFoundError. Now the scan reads what is left and lists what it could not read in the summary and on stderr. --- docs/exact-restart.md | 7 ++++++- docs/offline-scan.md | 11 +++++++---- src/reactionflow/cli.py | 26 ++++++++++++++++++++++---- src/reactionflow/scan.py | 27 ++++++++++++++++++++++----- 4 files changed, 57 insertions(+), 14 deletions(-) diff --git a/docs/exact-restart.md b/docs/exact-restart.md index 1346c80..10cd1e3 100644 --- a/docs/exact-restart.md +++ b/docs/exact-restart.md @@ -87,7 +87,12 @@ run keeps a single runtime checkpoint. `state.json` records the phase, generation, global counters, pending pathway occurrence IDs, detector state, configuration, any durable-state failure, and the active checkpoint. Each generation writes one frame per observation to `segments//trajectory.traj`; a frame is -written before its checkpoint, so an interrupted attempt never leaves a gap. An MD state with a +written before its checkpoint, so an interrupted attempt never leaves a gap. A segment lost later, +to a node crash or a cut copy, cannot change the state in the checkpoint, so MD continues, and every +`run`, `status` and `scan` then names the missing frames. `run` and `status` read each segment's +ends, so damage inside a segment shows only in `scan`, which reads every frame. The one refusal is a +current segment that cannot be appended to: restore it, or move it aside to continue in a new file. +An MD state with a non-finite position, momentum, cell, or integrator array gets neither: the job stops with a non-zero exit, and the next job repeats from the last checkpoint, or starts over if there is none. A frame or checkpoint an earlier version wrote from such a state is not repaired. diff --git a/docs/offline-scan.md b/docs/offline-scan.md index c4bdcb6..834c503 100644 --- a/docs/offline-scan.md +++ b/docs/offline-scan.md @@ -13,9 +13,11 @@ reactionflow scan XDATCAR --pair C-O 1.6 2.0 --frame-ps 0.05 ## Campaign runs A run saves one frame per observation in `segments//trajectory.traj`. A new generation -starts with the frame it resumed from, and the scan reads that frame once. It stops at the first -frame it cannot read, such as the one a running job is writing, and it prints `"scan": null` for a -trajectory with no saved frames yet. +starts with the frame it resumed from, and the scan reads that frame once. It reads a segment up to +the first frame it cannot read, such as the one a running job is writing, and passes over a segment +that is missing, empty or cut. Frames the run has checkpointed that the scan could not read are +listed as `missing_frames`, in the summary and on standard error, and the timing of a change across +such a gap is not known. It prints `"scan": null` for a trajectory with no saved frames. With no options the scan uses the run's own settings from `state.json`, so it finds the candidates the run recorded. `candidates.csv` gives each one's occurrence ID, and in `summary.json` @@ -62,7 +64,8 @@ directory. unresolved; - `net_changes_in_resolved` and `classes_in_resolved`, the bond changes and local classes in the resolved candidates; - - `recorded_by_run` and `matching_recorded` for a campaign run. + - `recorded_by_run` and `matching_recorded` for a campaign run; + - `missing_frames`, only when a campaign run's saved frames have gaps. These counts separate the stages between a threshold crossing and a candidate: what persistence filters, what is confirmed, how confirmed changes group into candidates, and how many never come diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 09ccce1..ac11d5a 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -8,6 +8,7 @@ import math import os import socket +import sys from collections import Counter from collections.abc import Mapping, Sequence from contextlib import suppress @@ -25,7 +26,7 @@ from .mlip import load_mlip_adapter from .model_setup import install_and_prepare, model_catalog, prepare_adapter, prepare_campaign from .pathway import check_saddle -from .run import ReactionRun, ReactionRunConfig, RunSummary, missing_frames +from .run import ReactionRun, ReactionRunConfig, RunSummary, frame_ranges, missing_frames from .scan import ScanResult, file_frames, recorded_occurrences, run_frames, scan_frames from .status import campaign_status, format_status from .wizard import create_campaign @@ -427,20 +428,37 @@ def main(argv: Sequence[str] | None = None) -> int: for index in indices: trajectory = campaign.trajectory(index) root = campaign.output_root / trajectory.id + state_path = root / "state.json" if not (root / "segments").is_dir(): - missing = {"scan": None, "trajectory_id": trajectory.id, "message": "not started"} + message = "no saved frames" if state_path.is_file() else "not started" + missing = {"scan": None, "trajectory_id": trajectory.id, "message": message} print(json.dumps(missing)) continue - state = json.loads((root / "state.json").read_text(encoding="utf-8")) - base = ReactionRunConfig.from_dict(state["config"]) + # A copy of the frames alone still scans, under the campaign's settings. + state = ( + json.loads(state_path.read_text(encoding="utf-8")) if state_path.is_file() else None + ) + base = ( + campaign.reaction_run + if state is None + else ReactionRunConfig.from_dict(state["config"]) + ) settings = _scan_settings(parser, arguments, base) result = scan_frames( run_frames(root, timestep_fs=trajectory.timestep_fs), occurrences=recorded_occurrences(root), + last_frame=0 if state is None else state["global_frame"], **settings, ) out = root / "scan" if arguments.out is None else arguments.out / trajectory.id _print_scan(result.write(out), result, trajectory_id=trajectory.id) + if "missing_frames" in result.summary: + print( + f"warning: {trajectory.id} lacks saved frames " + f"{frame_ranges(result.summary['missing_frames'])}, so these tables leave them " + "out, and the timing of any change across a gap is not known", + file=sys.stderr, + ) elif arguments.command == "validate": payload = { "campaign": str(campaign.source), diff --git a/src/reactionflow/scan.py b/src/reactionflow/scan.py index 3710b6d..81d0d59 100644 --- a/src/reactionflow/scan.py +++ b/src/reactionflow/scan.py @@ -80,8 +80,9 @@ def write(self, directory: str | Path) -> Path: def run_frames(root: str | Path, *, timestep_fs: float | None = None) -> Iterator[Frame]: """Yield each saved observation of a ``ReactionRun`` once, in order. - A new trajectory generation starts with the frame it resumed from, which is skipped. Reading - stops at the first frame that cannot be read, such as the one a running job is writing. + A new trajectory generation starts with the frame it resumed from, which is skipped. A segment + is read up to its first frame that cannot be read, such as the one a running job is writing, and + a segment that is missing, empty or cut is passed over, so that the later ones still count. """ segments = sorted( @@ -92,13 +93,17 @@ def run_frames(root: str | Path, *, timestep_fs: float | None = None) -> Iterato for segment in segments: path = segment / "trajectory.traj" if not path.is_file() or path.stat().st_size == 0: - return - with Trajectory(path) as trajectory: + continue + try: + trajectory = Trajectory(path) + except Exception: # a cut that took out the offset table + continue + with trajectory: for index in range(len(trajectory)): try: atoms = trajectory[index] except Exception: - return + break frame = int(atoms.info["reactionflow_global_frame"]) if frame <= last: continue @@ -158,6 +163,7 @@ def scan_frames( region_radius: int | None = 2, species_every: int = 10, occurrences: Mapping[OccurrenceKey, str] | None = None, + last_frame: int | None = None, ) -> ScanResult: """Replay ordered frames through a new detector and tracker and tabulate what they find. @@ -165,6 +171,8 @@ def scan_frames( numbers, as ``run_frames`` and ``file_frames`` do. Molecules are the bonded components of the first frame, and candidates fall into the local classes ``reactionflow status`` reports. ``occurrences``, from ``recorded_occurrences``, names the candidates a run recorded itself. + ``last_frame``, a run's last checkpointed observation, makes the summary list the frames from 0 + to it that were not read, as ``missing_frames``. Timing across such a gap is not known. """ if species_every < 1: @@ -181,6 +189,8 @@ def scan_frames( species: list[dict[str, Any]] = [] waiting: frozenset[tuple[Bond, str]] = frozenset() reverted = net_changes = 0 + missing: list[list[int]] = [] + read_to = -1 def add(candidate: ReactionCandidate, *, at_end: bool) -> None: nonlocal net_changes @@ -211,6 +221,9 @@ def add(candidate: ReactionCandidate, *, at_end: bool) -> None: ) for frame, time_ps, atoms in frames: + if last_frame is not None and frame > read_to + 1: + missing.append([read_to + 1, frame - 1]) + read_to = frame times[frame] = time_ps events = detector.process(atoms, frame=frame) if not symbols: @@ -263,6 +276,8 @@ def add(candidate: ReactionCandidate, *, at_end: bool) -> None: ) for candidate in tracker.finish(): add(candidate, at_end=True) + if last_frame and read_to < last_frame: + missing.append([read_to + 1, last_frame]) resolved = sum(row["resolved"] for row in candidates) summary: dict[str, Any] = { @@ -287,6 +302,8 @@ def add(candidate: ReactionCandidate, *, at_end: bool) -> None: if occurrences is not None: summary["recorded_by_run"] = len(occurrences) summary["matching_recorded"] = sum(bool(row["occurrence_id"]) for row in candidates) + if missing: + summary["missing_frames"] = missing return ScanResult(bond_changes, candidates, species, summary) From e90d4ee6f3c5a5c7f5e68d72ab7098a8c7952cfc Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:40:52 -0500 Subject: [PATCH 03/25] Keep a trajectory's reactions in status when one of its files cannot be read One unreadable bundle or result.json dropped all of a trajectory's reactions and showed events 0, so a copy of a run's results alone looked empty, and a missing database looked like no reactions at all. Each file now costs only what depends on it. Another occurrence's bundle stands for its class, a class with none keeps its counts and barriers in its own row, and the error column and a line under the table say what was left out. Lost saved frames show there too. --- docs/campaigns.md | 7 ++- src/reactionflow/status.py | 111 +++++++++++++++++++++++++++++-------- 2 files changed, 94 insertions(+), 24 deletions(-) diff --git a/docs/campaigns.md b/docs/campaigns.md index 4687201..08d76b8 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -249,8 +249,11 @@ imports the MLIP. The first table has one row per trajectory: model profile, pressure, phase, step, detected events, converged pathways out of those refined, and the latest error. That error is the permanent failure for a `failed` trajectory, and otherwise the contents of `last-error.json`. Trajectories -without output are listed as `not started`. A trajectory whose files cannot be read is reported -with that error while the rest of the report still prints. +without output are listed as `not started`. A file that cannot be read costs only what depends on it: +the error column counts what was left out, `--json` lists it as `unreadable`, and a line under the +table says the report is incomplete. A class whose candidate bundles cannot all be read, as in a copy +of a run's results alone, keeps its counts and barriers in a row of its own, because classes merge +only by their graphs. Saved frames that are gone appear in the error column and as `missing_frames`. Pathway counts include every refined occurrence. An event refined in steps counts each step as a refined pathway, and only its steps can converge. diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index a144f27..a685ec9 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -17,6 +17,7 @@ from .candidates import ReactionCandidate, ReactionClasses, _reaction_graph from .detection import _symbols_by_id from .pathway import CONVERGED_STATUSES +from .run import frame_ranges, missing_frames from .store import _read_bundle # Converged bands, including those refined before the saddle check, which have no verdict. @@ -64,10 +65,13 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: """Fill one trajectory's status row and return its classes. Each class is (candidate, events, results, occurrence ID of the candidate). Each step of an - event refined in steps is one more, an event of its own with the step's bonds. + event refined in steps is one more, an event of its own with the step's bonds. The candidate is + None when no bundle of the class can be read, as in a copy of a run's results alone, so that its + counts and barriers still show. The row is filled as reading progresses, so a damaged file later on still leaves the phase and - step from state.json in place. + step from state.json in place. Files that cannot be read are named in ``row["unreadable"]``, + and each costs only what depends on it. """ state = _read_json(root / "state.json") @@ -76,35 +80,77 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: if state["phase"] == "failed" and failure: row["error"] = f"{failure.get('type')}: {failure.get('message')}" records = _occurrences(root) + unreadable: Counter[str] = Counter() + lost = missing_frames(root, state["global_frame"]) + if lost: + row["missing_frames"] = lost members: dict[str, list[sqlite3.Row]] = {} for record in records: members.setdefault(record["class_id"], []).append(record) classes = [] steps = [] pathways: Counter[str] = Counter() + + read: dict[str, ReactionCandidate | None] = {} + + def bundle(occurrence: sqlite3.Row) -> ReactionCandidate | None: + name = occurrence["bundle"] + if name not in read: + try: + read[name] = _read_bundle(root / name)[1] + except Exception: + read[name] = None + unreadable["candidate bundle"] += 1 + return read[name] + + def result_file(path: Path) -> dict[str, Any] | None: + try: + result = _read_json(path) + pathways[result["status"]] += 1 + return result + except Exception: + unreadable["result file"] += 1 + return None + for occurrences in members.values(): - _, candidate = _read_bundle(root / occurrences[0]["bundle"]) + # Any occurrence can stand for its class. + candidate, example = None, occurrences[0]["occurrence_id"] + for occurrence in occurrences: + candidate = bundle(occurrence) + if candidate is not None: + example = occurrence["occurrence_id"] + break results = [] for occurrence in occurrences: path = root / "pathways" / occurrence["occurrence_id"] / "result.json" - if path.is_file(): - result = _read_json(path) - results.append(result) - pathways[result["status"]] += 1 - if result.get("steps"): - _, event = _read_bundle(root / occurrence["bundle"]) - for step in result["steps"]: - step_path = path.parent / f"step-{step['step']}" / "result.json" - step_result = _read_json(step_path) - bonds = { + result = result_file(path) if path.is_file() else None + if result is None: + continue + results.append(result) + if result.get("steps"): + event = bundle(occurrence) + for step in result["steps"]: + step_path = path.parent / f"step-{step['step']}" / "result.json" + step_result = result_file(step_path) + if step_result is None: + continue + part = event and replace( + event, + **{ key: frozenset(map(tuple, step_result[key])) for key in ("reactant_bonds", "product_bonds") - } - example = occurrence["occurrence_id"] - steps.append((replace(event, **bonds), 1, [step_result], example)) - pathways[step_result["status"]] += 1 - classes.append((candidate, len(occurrences), results, occurrences[0]["occurrence_id"])) + }, + ) + steps.append((part, 1, [step_result], occurrence["occurrence_id"])) + classes.append((candidate, len(occurrences), results, example)) row.update(events=len(records), pathways=dict(pathways)) + notes = [f"{count} {what}{'' if count == 1 else 's'}" for what, count in unreadable.items()] + if not (root / "reactions.sqlite3").is_file(): + bundles = root.glob("candidates/[!.]*/candidate.json") + if count := sum(1 for _ in bundles): + notes.append(f"no reactions.sqlite3 to count {count} candidates") + if notes: + row["unreadable"] = notes return classes + steps @@ -149,14 +195,27 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A row["error"] = _read_json(root / "last-error.json").get("error") except Exception as error: row["error"] = f"unreadable last-error.json: {type(error).__name__}" + notes = [row["error"]] if row["error"] else [] + if "missing_frames" in row: + lost = frame_ranges(row["missing_frames"]) + notes.insert(0, f"saved frames {lost} missing or unreadable") + if "unreadable" in row: + notes.insert(0, "unreadable: " + ", ".join(row["unreadable"])) + row["error"] = "; ".join(notes) or None merged = known.setdefault((model, spec.pressure_GPa), ReactionClasses(radius)) for candidate, events, results, occurrence_id in classes: - name = merged.find(candidate) + # A class with no readable candidate has no graph to merge by, so it keeps its own row. + name = None if candidate is None else merged.find(candidate) if name is None: name = str(len(entries)) - merged.add(name, candidate) + if candidate is not None: + merged.add(name, candidate) entries[name] = { - "reaction": _describe(candidate, radius), + "reaction": ( + f"class of {occurrence_id} (no readable candidate)" + if candidate is None + else _describe(candidate, radius) + ), "example": {"trajectory": spec.id, "occurrence_id": occurrence_id}, "model": model, "pressure_GPa": spec.pressure_GPa, @@ -251,8 +310,16 @@ def format_status(report: dict[str, Any]) -> str: + ", ".join(f"{count} {phase}" for phase, count in sorted(phases.items())), "", ] + incomplete = sum("unreadable: " in (row["error"] or "") for row in report["trajectories"]) + if incomplete: + lines += [ + f"Incomplete: files of {incomplete} " + f"{'trajectory' if incomplete == 1 else 'trajectories'} cannot be read, so what " + "follows leaves them out; the error column and --json say which.", + "", + ] if not report["reactions"]: - lines.append("No reactions recorded yet.") + lines.append("No reactions could be read." if incomplete else "No reactions recorded yet.") return "\n".join(lines) reaction_rows = [] for entry in report["reactions"]: From ceb87339fb91a8639dd506a8982cb85d85a23f48 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:43:15 -0500 Subject: [PATCH 04/25] Say why EON searches fail, in status and as a run goes A run whose searches all failed printed nothing, and status showed only a failed count at confidence 0; the reason was only in attempts/*.json. A run that ended in an error left no trace in status, and status said not started when records existed without search-state.json. Now status counts failed searches by reason, run prints each new reason on stderr up to ten, and the error that ended a run is kept in last-error.json and shown by status until the next run starts. Whether a run should stop by itself is a separate decision; this only reports. --- docs/eon-search.md | 4 +- src/reactionflow/search/cli.py | 154 ++++++++++++++++++++---------- src/reactionflow/search/runner.py | 13 +++ 3 files changed, 122 insertions(+), 49 deletions(-) diff --git a/docs/eon-search.md b/docs/eon-search.md index 889d12b..04c450d 100644 --- a/docs/eon-search.md +++ b/docs/eon-search.md @@ -245,7 +245,9 @@ output directory. show up in the step list. `basins: "bonds"` merges states only when no bond changes between them; a low-barrier reaction is still a separate step. - A state is searched until it reaches the confidence target, with no cap on searches; `status` - shows the searches and confidence of the current state. + shows the searches and confidence of the current state, and counts the failed searches by reason + (`failure_reasons` in `--json`). `run` prints each new reason on standard error as it first + appears, up to ten. A confidence that stays at 0 while searches fail will not recover by itself. - Searches use a fixed, nonsingular cell, fully periodic or fully nonperiodic, with `FixAtoms` as the only constraint. - By default EON converges saddles and minima on the norm of the whole force vector diff --git a/src/reactionflow/search/cli.py b/src/reactionflow/search/cli.py index 9cf4b20..e2ede30 100644 --- a/src/reactionflow/search/cli.py +++ b/src/reactionflow/search/cli.py @@ -9,7 +9,7 @@ import time from collections import Counter from collections.abc import Iterator -from contextlib import contextmanager +from contextlib import contextmanager, suppress from dataclasses import asdict, dataclass from pathlib import Path from typing import Any @@ -131,13 +131,34 @@ def starting_atoms(self) -> Atoms: ) +def _last_error(root: Path) -> dict[str, str]: + """The error that ended the last run, as a report entry, or nothing.""" + + try: + return {"last_error": json.loads((root / "last-error.json").read_text())["error"]} + except FileNotFoundError: + return {} + except Exception as error: + return {"last_error": f"unreadable last-error.json: {type(error).__name__}"} + + def search_status(campaign: SearchCampaign) -> dict[str, Any]: """Read progress and records without creating directories or loading resources.""" root = campaign.output_root - report: dict[str, Any] = {"mode": "eon", "output_root": str(root)} + report: dict[str, Any] = {"mode": "eon", "output_root": str(root), **_last_error(root)} if not (root / "search-state.json").exists(): - return report | {"status": "not_started", "steps": 0, "time_s": 0.0, "attempts": 0} + records = sum( + 1 + for kind in ("states", "processes", "attempts", "steps") + for _ in (root / kind).glob("*") + ) + status = ( + {"status": "state_missing", "records": records} + if records + else {"status": "not_started"} + ) + return report | status | {"steps": 0, "time_s": 0.0, "attempts": 0} state = read_search_state(root) config = AKMCConfig(**json.loads((root / "search-contract.json").read_text())["akmc"]) store = SearchStore(root) @@ -149,6 +170,11 @@ def search_status(campaign: SearchCampaign) -> dict[str, Any]: table = [{**processes[pid], "id": pid} for pid in counters["processes"]] kept = {process["id"] for process in kept_processes(table, config)} saddles = [item["minima"] for item in attempts if "minima" in item] + reasons = Counter( + str(item.get("message") or "no message") for item in attempts if item["status"] == "failed" + ) + if reasons: + report["failure_reasons"] = dict(reasons.most_common()) return report | { "status": "stopped" if state["stop_reason"] else "incomplete", "stop_reason": state["stop_reason"], @@ -197,8 +223,12 @@ def search_status(campaign: SearchCampaign) -> dict[str, Any]: def _format_status(report: dict[str, Any]) -> str: + error = [f"Last error: {report['last_error']}"] if "last_error" in report else [] if report["status"] == "not_started": - return "EON AKMC: not started" + return "\n".join(["EON AKMC: not started", *error]) + if report["status"] == "state_missing": + line = f"EON AKMC: search-state.json is missing, with {report['records']} records present" + return "\n".join([line, *error]) chemical = report["chemical_steps"] steps = report["steps"] if chemical is None else f"{chemical} chemical of {report['steps']}" lines = [ @@ -211,6 +241,12 @@ def _format_status(report: dict[str, Any]) -> str: f"States: {report['states']}; processes: {report['processes']}; searches: " f"{report['attempts']} {json.dumps(report['attempt_statuses'], sort_keys=True)}", ] + lines += error + reasons = report.get("failure_reasons", {}) + if reasons: + shown = [f"{count} {reason}" for reason, count in list(reasons.items())[:3]] + more = f"; and {len(reasons) - 3} other reasons" if len(reasons) > 3 else "" + lines.append(f"Failed searches: {'; '.join(shown)}{more}") if report["failed_with_saddle"]: lines.append( f"Failed searches that kept a saddle: {report['failed_with_saddle']}, of which " @@ -252,49 +288,49 @@ def _eon_output(path: Path | None) -> Iterator[None]: os.close(saved) -def handle_command(arguments: argparse.Namespace) -> int: - campaign = SearchCampaign.load(arguments.campaign) - command = arguments.command - if command == "prepare": - if arguments.install: - extra = (f"pyeonclient[ase]=={EON_VERSION}",) - return install_and_prepare(campaign, index=arguments.index, extra_packages=extra) - models = prepare_campaign(campaign, index=arguments.index) - print(json.dumps({"campaign": str(campaign.source), "models": models}, sort_keys=True)) - return 0 - if command == "status": - report = search_status(campaign) - if arguments.json: - print(json.dumps(report, indent=2, sort_keys=True)) - else: - print(_format_status(report)) - return 0 - if command == "run": - from ..cli import resolve_task_index, visible_gpu - - resolve_task_index(1, arguments.index) - if campaign.require_gpu: - visible_gpu() - if arguments.download: - prepare_adapter(campaign.adapter) +@contextmanager +def _keeping_last_error(root: Path) -> Iterator[None]: + """Clear the last run's error, and keep this run's for `status` if it ends in one.""" + + path = root / "last-error.json" + with suppress(OSError): + path.unlink() + try: + yield + except Exception as error: + with suppress(OSError): + root.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps({"error": f"{type(error).__name__}: {error}"}) + "\n") + raise + + +def _validate(campaign: SearchCampaign) -> None: + atoms = campaign.starting_atoms() + report = { + "campaign": str(campaign.source), + "mode": "eon", + "atoms": len(atoms), + "output_root": str(campaign.output_root), + "tasks": 1, + "kmc_steps": campaign.akmc.steps, + "temperature_K": campaign.akmc.temperature_K, + "adapter_factory": campaign.adapter.factory, + "require_gpu": campaign.require_gpu, + "rigid_rotation": campaign.eon.rotation_for(atoms), + "warnings": [PINNED] if free_rotation(atoms) == "pinned" else [], + } + print(json.dumps(report, sort_keys=True)) + + +def _run(campaign: SearchCampaign, arguments: argparse.Namespace) -> None: + from ..cli import resolve_task_index, visible_gpu + + resolve_task_index(1, arguments.index) + if campaign.require_gpu: + visible_gpu() + if arguments.download: + prepare_adapter(campaign.adapter) atoms = campaign.starting_atoms() - warnings = [PINNED] if free_rotation(atoms) == "pinned" else [] - if command == "validate": - report = { - "campaign": str(campaign.source), - "mode": "eon", - "atoms": len(atoms), - "output_root": str(campaign.output_root), - "tasks": 1, - "kmc_steps": campaign.akmc.steps, - "temperature_K": campaign.akmc.temperature_K, - "adapter_factory": campaign.adapter.factory, - "require_gpu": campaign.require_gpu, - "rigid_rotation": campaign.eon.rotation_for(atoms), - "warnings": warnings, - } - print(json.dumps(report, sort_keys=True)) - return 0 adapter = load_mlip_adapter(campaign.adapter, campaign.trajectories[0]) preflight = getattr(adapter, "preflight", None) if callable(preflight): @@ -305,8 +341,8 @@ def handle_command(arguments: argparse.Namespace) -> int: lease = getattr(adapter, "_lease", None) if lease is None: raise TypeError("EON needs a built-in model adapter or reactionflow.adapters.ase") - for warning in warnings: - print(f"warning: {warning}", file=sys.stderr) + if free_rotation(atoms) == "pinned": + print(f"warning: {PINNED}", file=sys.stderr) log = campaign.output_root / "eon.log" if campaign.eon.log == "file" else None # One model load gives the calculator and the same identity an MD restart records. with lease() as (calculator, identity), _eon_output(log): @@ -320,4 +356,26 @@ def handle_command(arguments: argparse.Namespace) -> int: campaign.output_root, atoms, backend=backend, contract=contract, config=campaign.akmc ) print(json.dumps(search_status(campaign), sort_keys=True)) + + +def handle_command(arguments: argparse.Namespace) -> int: + campaign = SearchCampaign.load(arguments.campaign) + command = arguments.command + if command == "prepare": + if arguments.install: + extra = (f"pyeonclient[ase]=={EON_VERSION}",) + return install_and_prepare(campaign, index=arguments.index, extra_packages=extra) + models = prepare_campaign(campaign, index=arguments.index) + print(json.dumps({"campaign": str(campaign.source), "models": models}, sort_keys=True)) + elif command == "status": + report = search_status(campaign) + if arguments.json: + print(json.dumps(report, indent=2, sort_keys=True)) + else: + print(_format_status(report)) + elif command == "validate": + _validate(campaign) + else: + with _keeping_last_error(campaign.output_root): + _run(campaign, arguments) return 0 diff --git a/src/reactionflow/search/runner.py b/src/reactionflow/search/runner.py index f729bac..3999caf 100644 --- a/src/reactionflow/search/runner.py +++ b/src/reactionflow/search/runner.py @@ -9,6 +9,7 @@ import fcntl import hashlib import json +import sys from dataclasses import asdict, replace from pathlib import Path from typing import Any @@ -24,6 +25,10 @@ from .labels import bond_changes, reversed_changes from .records import SearchStore +# Distinct reasons for a failed search that `run` prints, so that a model failing in a new way each +# time cannot flood a job's log. `status` counts them all. +_PRINTED_REASONS = 10 + def _save_state(root: Path, state: dict[str, Any]) -> None: temporary = root / ".search-state.json.tmp" @@ -411,6 +416,7 @@ def run_exploration( _save_state(root, state) counted = "chemical_steps" if config.basins == "bonds" else "steps" stop_reason = "steps_completed" + printed: set[str] = set() while state[counted] < config.steps: if state["active"] is None: counters = state["states"][state["current"]] @@ -439,6 +445,13 @@ def run_exploration( else: result, structures = _attempt(store, state, active, backend, config) store.write("attempts", active["id"], result, structures=structures) + if ( + result["status"] == "failed" + and result["message"] not in printed + and len(printed) < _PRINTED_REASONS + ): + printed.add(result["message"]) + print(f"{active['id']} failed: {result['message']}", file=sys.stderr) _commit_attempt(state, result) if result["status"] == "no_push": # Every later search from this state would select nothing as well. From 0ded79f13f992288b952c2c60e9ec319de3676c8 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:45:53 -0500 Subject: [PATCH 05/25] Refuse a band whose images hold different atoms in check_saddle A last image with 3 atoms or a middle image with O2 validated with no complaint, and a saddle image in another atom order displaced the wrong atoms: a real saddle came back no_saddle. A NaN energy on the first image gave saddle_validated with a NaN barrier. Images now must hold the same atoms, elements and periodicity, images with IDs are put in the first image's order, a non-finite energy ends as saddle_unchecked, and a pressure on a band with no cell says so. --- docs/pathway-refinement.md | 6 ++++-- src/reactionflow/pathway.py | 32 ++++++++++++++++++++++++++++++-- 2 files changed, 34 insertions(+), 4 deletions(-) diff --git a/docs/pathway-refinement.md b/docs/pathway-refinement.md index aad5c7d..4436be1 100644 --- a/docs/pathway-refinement.md +++ b/docs/pathway-refinement.md @@ -144,7 +144,7 @@ The check ends a converged band as one of these: | `saddle_validated` | one mode, and the two sides relax to the band's two ends | | `no_saddle` | an endpoint is the highest point, or the saddle image has no mode | | `saddle_rejected` | more than one mode, or the sides do not reach both ends | -| `saddle_unchecked` | the Hessian failed, or the connectivity check failed or was inconclusive | +| `saddle_unchecked` | an image's energy was not finite, the Hessian failed, or the connectivity check failed or was inconclusive | The outcome's message says which case applies. A saddle lies above both endpoints, so when an endpoint is the highest point of the band, there is no saddle to test, and the outcome has no @@ -200,7 +200,9 @@ pathway's `images.traj`, a band from another code, or an EON saddle written as t (minimum, saddle, minimum). The campaign supplies the model, the bond thresholds, and the pathway settings. `--index` picks the trajectory whose model and pressure are used; a campaign with several trajectories needs it, even inside a Slurm job. `--pressure` replaces that pressure. The verdict is written beside the band as `.check.json`, in the form -of a pathway's `result.json`, and its status, barrier, and message are printed. +of a pathway's `result.json`, and its status, barrier, and message are printed. The images must hold +the same atoms, elements, and periodicity, or the check refuses the band; images with atom IDs in +another order are put in the first image's order, and a file without IDs is taken in its own order. Campaign and state files written before the saddle check may name `active_radius` and `imaginary_frequency_cutoff_cm1`. They still open, and those two settings are ignored. A trajectory diff --git a/src/reactionflow/pathway.py b/src/reactionflow/pathway.py index cd01eca..e17d9ba 100644 --- a/src/reactionflow/pathway.py +++ b/src/reactionflow/pathway.py @@ -576,6 +576,26 @@ def _verdict( return "saddle_unchecked", f"the connectivity check was inconclusive{detail}" +def _one_atom_list(images: Sequence[Atoms]) -> list[Atoms]: + """The images in the first one's atom order, or a refusal if they do not hold the same atoms.""" + + first = images[0] + ids = atom_ids(first) + band = [first] + for number, image in enumerate(images[1:], 1): + order = {atom_id: index for index, atom_id in enumerate(atom_ids(image))} + if order.keys() != set(ids): + raise ValueError(f"band image {number} holds different atoms from image 0") + if list(order) != list(ids): + image = image[[order[atom_id] for atom_id in ids]] + if image.get_chemical_symbols() != first.get_chemical_symbols(): + raise ValueError(f"band image {number} holds different elements from image 0") + if image.pbc.tolist() != first.pbc.tolist(): + raise ValueError(f"band image {number} has different periodicity from image 0") + band.append(image) + return band + + def check_saddle( images: Sequence[Atoms], calculator: Calculator, @@ -593,9 +613,14 @@ def check_saddle( """ options = config or PathwayConfig() - band = [canonical_copy(image) for image in images] - if len(band) < 3: + if len(images) < 3: raise ValueError("a band needs its two ends and at least one image between them") + carrying = {"atom_id" in image.arrays or "atom_ids" in image.info for image in images} + if len(carrying) > 1: + raise ValueError("some band images carry atom IDs and others do not, so atoms cannot match") + band = _one_atom_list([canonical_copy(image) for image in images]) + if pressure_GPa is not None and band[0].cell.rank < 3: + raise ValueError("a band checked at a pressure needs a full cell, for its volume") energies = [] for image in band: image.calc = calculator @@ -619,6 +644,9 @@ def check_saddle( "images": tuple(band), "pressure_GPa": pressure_GPa, } + if not np.isfinite(profile).all(): + message = "the model returned a non-finite energy for a band image" + return PathwayOutcome("saddle_unchecked", message=message, **common) # A saddle between two minima lies above both. When an endpoint is the highest point, the band # has not found one, so there is no barrier or saddle to check. if max(profile[1:-1]) <= max(profile[0], profile[-1]): From 0e0364fe95284e1bb0759a502d4f84488cd0a2dc Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:45:53 -0500 Subject: [PATCH 06/25] Report the residual forces of a checked band MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit check said saddle_validated for a highest image with a force of 0.545 eV/Å and for unrelaxed ends that made a barrier 6.7 times too low. The first force sat deep in the JSON and the second nowhere. check now records the largest atomic force on the first, highest and last images, and a validated band's message says when they pass neb_fmax or relax_fmax. The verdict and barrier do not change: whether check should withhold them is a scientific acceptance threshold, left open. --- docs/pathway-refinement.md | 7 +++++++ src/reactionflow/cli.py | 23 +++++++++++++++++++++-- src/reactionflow/pathway.py | 25 +++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/docs/pathway-refinement.md b/docs/pathway-refinement.md index 4436be1..fd4d444 100644 --- a/docs/pathway-refinement.md +++ b/docs/pathway-refinement.md @@ -204,6 +204,13 @@ of a pathway's `result.json`, and its status, barrier, and message are printed. the same atoms, elements, and periodicity, or the check refuses the band; images with atom IDs in another order are put in the first image's order, and a file without IDs is taken in its own order. +The verdict rests on curvature and connectivity, and the barrier is the highest image minus the first +image as given, so it means what it says only for relaxed ends and a stationary highest image, which +a band from another code may not have. `largest_force_eV_A` in the verdict gives the largest atomic +force on the first, highest, and last images, and the message of a validated band says when these +pass `neb_fmax` or `relax_fmax`. A converged band of this package can pass `neb_fmax` slightly, +because that tolerance applies to the band's projected forces. Neither changes the verdict. + Campaign and state files written before the saddle check may name `active_radius` and `imaginary_frequency_cutoff_cm1`. They still open, and those two settings are ignored. A trajectory started before the check cannot resume under it, because its pathway settings, and so its campaign diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index ac11d5a..a85789d 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -25,7 +25,7 @@ from .detection import _is_count from .mlip import load_mlip_adapter from .model_setup import install_and_prepare, model_catalog, prepare_adapter, prepare_campaign -from .pathway import check_saddle +from .pathway import check_saddle, largest_forces from .run import ReactionRun, ReactionRunConfig, RunSummary, frame_ranges, missing_frames from .scan import ScanResult, file_frames, recorded_occurrences, run_frames, scan_frames from .status import campaign_status, format_status @@ -198,13 +198,31 @@ def check_band( config=campaign.reaction_run.pathway, pressure_GPa=pressure_GPa, ) + forces = largest_forces(outcome, calculator) result: dict[str, object] = { "band": str(band.resolve()), "trajectory_id": trajectory.id, "reactionflow_version": __version__, **outcome.to_dict(), + "largest_force_eV_A": forces, } del result["method"] # how the band was made is not known here + # This reports and does not decide: the verdict rests on curvature and connectivity alone. + if outcome.barrier is not None and None not in forces.values(): + options = campaign.reaction_run.pathway + concerns = [] + if forces["highest"] > options.neb_fmax: + concerns.append( + f"the highest image carries a force of {forces['highest']:.3g} eV/Å, where a " + f"stationary point would carry about neb_fmax, {options.neb_fmax:g}" + ) + if max(forces["first"], forces["last"]) > options.relax_fmax: + concerns.append( + f"the ends carry forces of {forces['first']:.3g} and {forces['last']:.3g} eV/Å, " + f"where relaxed ends would carry about relax_fmax, {options.relax_fmax:g}, and " + "the barrier is measured from the first image as given" + ) + result["message"] = "; ".join(concerns) path = band.with_name(f"{band.stem}.check.json") temporary = path.with_name(f".{path.name}-{uuid4().hex}.tmp") try: @@ -482,7 +500,8 @@ def main(argv: Sequence[str] | None = None) -> int: elif not math.isfinite(pressure): parser.error("--pressure must be a finite number") path, result = check_band(campaign, arguments.band, index=index, pressure_GPa=pressure) - summary = {key: result[key] for key in ("status", "barrier_eV", "message")} + keys = ("status", "barrier_eV", "message", "largest_force_eV_A") + summary = {key: result[key] for key in keys} print(json.dumps({"check": str(path), **summary}, sort_keys=True)) else: _run_and_record(campaign, arguments) diff --git a/src/reactionflow/pathway.py b/src/reactionflow/pathway.py index e17d9ba..6802eaa 100644 --- a/src/reactionflow/pathway.py +++ b/src/reactionflow/pathway.py @@ -681,6 +681,31 @@ def check_saddle( ) +def largest_forces(outcome: PathwayOutcome, calculator: Calculator) -> dict[str, float | None]: + """The largest atomic force in eV/Å on a checked band's first, highest and last images, or None + for an image where the model fails or gives no finite one. + + ``check_saddle`` asks only about curvature and connectivity, and takes the barrier from the + first image as given. These say whether the images it read are stationary points at all. + """ + + profile = outcome.enthalpies or outcome.energies + forces: dict[str, float | None] = {} + for name, index in ( + ("first", 0), + ("highest", 1 + int(np.argmax(profile[1:-1]))), + ("last", len(profile) - 1), + ): + atoms = outcome.images[index].copy() + atoms.calc = calculator + try: + largest = float(np.linalg.norm(atoms.get_forces(apply_constraint=False), axis=1).max()) + except Exception: # the verdict is already made, and this must not cost it + largest = float("nan") + forces[name] = largest if np.isfinite(largest) else None + return forces + + def refine_pathway( candidate: ReactionCandidate, *, From 5abadec8edb627057541376ab956dc617a3c7729 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:46:44 -0500 Subject: [PATCH 07/25] Refuse region radius 0, which finds no reactions Radius 0 passed validate and a whole run then found nothing, with exit 0: a region of single atoms has no bond in it to change. Real frames gave 0 candidates against 36. Making 0 work would only move the emptiness to refinement, whose endpoint check compares bonds a region of changed atoms leaves out. status --radius 0 and same_reaction(radius=0) are a different radius, for classes, and stay. --- docs/candidate-tracking.md | 5 +++-- docs/offline-scan.md | 3 ++- src/reactionflow/candidates.py | 7 +++++-- src/reactionflow/cli.py | 4 ++-- src/reactionflow/run.py | 7 +++++-- 5 files changed, 17 insertions(+), 9 deletions(-) diff --git a/docs/candidate-tracking.md b/docs/candidate-tracking.md index 030108b..fa07a4a 100644 --- a/docs/candidate-tracking.md +++ b/docs/candidate-tracking.md @@ -31,8 +31,9 @@ Start the detector and tracker together so the tracker sees a baseline before an - The first stable topology establishes the accepted reactant basin. - Each changed region has its own stability window. A region is the changed atoms and the atoms - within `region_radius` bonds of them (2 by default), over the union of accepted, stable, and - provisional bonds. Regions that share an atom merge, so interacting changes stay together. + within `region_radius` bonds of them (2 by default, at least 1), over the union of accepted, + stable, and provisional bonds. Regions that share an atom merge, so interacting changes stay + together. - `region_radius=None` makes a region the whole bonded component instead. In a polymer, network, or slab that component holds every change in it, so unrelated changes merge and a crossing anywhere holds back a finished one. A campaign sets the radius as diff --git a/docs/offline-scan.md b/docs/offline-scan.md index 834c503..dcc9659 100644 --- a/docs/offline-scan.md +++ b/docs/offline-scan.md @@ -24,7 +24,8 @@ the run recorded. `candidates.csv` gives each one's occurrence ID, and in `summa `matching_recorded` equals `recorded_by_run` for a finished run. Options change the replay and nothing else: -- `--radius R` sets the region radius in bonds, and `--radius none` takes whole bonded components. +- `--radius R` sets the region radius in bonds, from 1, and `--radius none` takes whole bonded + components. - `--persistence N` sets how many observations a crossing must last to become a bond change. - `--stability N` sets how many observations a region's new bonds must hold. - `--pair A-B FORM BREAK` sets the formation and break distances in Å for one element pair. It can diff --git a/src/reactionflow/candidates.py b/src/reactionflow/candidates.py index ad806d2..64e50b9 100644 --- a/src/reactionflow/candidates.py +++ b/src/reactionflow/candidates.py @@ -217,8 +217,11 @@ class ReactionTracker: def __init__(self, *, stability_frames: int = 3, region_radius: int | None = 2) -> None: if not _is_count(stability_frames, 1): raise ValueError("stability_frames must be a positive integer") - if region_radius is not None and not _is_count(region_radius, 0): - raise ValueError("region_radius must be a non-negative integer or None") + if region_radius is not None and not _is_count(region_radius, 1): + raise ValueError( + "region_radius must be a positive integer or None: at 0 a region holds single " + "atoms, so no bond can change in it" + ) self.stability_frames = int(stability_frames) self.region_radius = None if region_radius is None else int(region_radius) self._symbols: dict[int, str] | None = None diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index a85789d..d59e3af 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -241,8 +241,8 @@ def _scan_settings( if arguments.radius == "none": radius = None elif arguments.radius is not None: - if not arguments.radius.isdigit(): - parser.error("--radius takes a non-negative integer or none") + if not arguments.radius.isdigit() or int(arguments.radius) < 1: + parser.error("--radius takes a positive integer or none") radius = int(arguments.radius) persistence = arguments.persistence if persistence is None: diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index e53f26a..6eb0527 100644 --- a/src/reactionflow/run.py +++ b/src/reactionflow/run.py @@ -57,9 +57,12 @@ def __post_init__(self) -> None: raise ValueError(f"{name} must be a positive integer") radius = self.candidate_region_radius if radius is not None and ( - isinstance(radius, bool) or not isinstance(radius, Integral) or radius < 0 + isinstance(radius, bool) or not isinstance(radius, Integral) or radius < 1 ): - raise ValueError("candidate_region_radius must be a non-negative integer or null") + raise ValueError( + "candidate_region_radius must be a positive integer or null: at 0 a region holds " + "single atoms, so no bond can change in it and the run finds no reactions" + ) def to_dict(self) -> dict[str, object]: return { From d74778940a328ef451ce31d32763becf14636177 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:48:59 -0500 Subject: [PATCH 08/25] Split a status row's barriers by direction A class holds both directions of a reaction, and each barrier is measured from its own event's reactant, so min, median and max ran over both under one label: 0.072 / 0.350 / 0.627 at 20 GPa, where the truth is 0.628 one way and 0.073 the other, and the median is the barrier of nothing. Now the barrier is for the direction the row names, and a further column, reverse_barriers_eV in --json, for the other. A symmetric reaction has the same graph both ways and stays one list, so its row does not change. A barrier whose event cannot be read goes to unknown_direction_barriers_eV. --- docs/campaigns.md | 6 ++- src/reactionflow/status.py | 87 +++++++++++++++++++++++++++++--------- 2 files changed, 71 insertions(+), 22 deletions(-) diff --git a/docs/campaigns.md b/docs/campaigns.md index 08d76b8..20f8d4f 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -262,8 +262,10 @@ pressure, because only those barriers are comparable. Classes merge by the local changed bonds and the atoms within `--radius` bonds of them (default 2), including forward/reverse equivalence ([local classes](reaction-identity.md#local-classes)). Each row lists the bond changes with the formula of those atoms, the number of trajectories and events, converged bands out of -those refined, validated saddles, and the minimum, median, and maximum barrier. Each step of an -event refined in steps is an event of its own local class +those refined, validated saddles, and the minimum, median, and maximum barrier. A barrier is measured +from each event's own reactant, so a class that runs both ways has two barriers: the row's barrier +is for the direction it names, and a further column, `reverse_barriers_eV` in `--json`, is for the +other. Each step of an event refined in steps is an event of its own local class ([steps](pathway-refinement.md#reactant-frame-and-steps)). Barriers come only from validated saddles ([saddle check](pathway-refinement.md#saddle-check)). A line below the table totals the saddle check's verdicts across the campaign (validated, no saddle, rejected, unchecked) diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index a685ec9..964e3ee 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -14,7 +14,7 @@ from ase.formula import Formula from .campaign import CampaignConfig -from .candidates import ReactionCandidate, ReactionClasses, _reaction_graph +from .candidates import ReactionCandidate, ReactionClasses, _isomorphic, _reaction_graph from .detection import _symbols_by_id from .pathway import CONVERGED_STATUSES from .run import frame_ranges, missing_frames @@ -64,10 +64,10 @@ def _describe(candidate: ReactionCandidate, radius: int) -> str: def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: """Fill one trajectory's status row and return its classes. - Each class is (candidate, events, results, occurrence ID of the candidate). Each step of an - event refined in steps is one more, an event of its own with the step's bonds. The candidate is - None when no bundle of the class can be read, as in a copy of a run's results alone, so that its - counts and barriers still show. + Each class is (candidate, events, results, occurrence ID of the candidate, and for each result + the event it refined, for its direction). Each step of an event refined in steps is one more, an + event of its own with the step's bonds. The candidate is None when no bundle of the class can be + read, as in a copy of a run's results alone, so that its counts and barriers still show. The row is filled as reading progresses, so a damaged file later on still leaves the phase and step from state.json in place. Files that cannot be read are named in ``row["unreadable"]``, @@ -103,6 +103,26 @@ def bundle(occurrence: sqlite3.Row) -> ReactionCandidate | None: unreadable["candidate bundle"] += 1 return read[name] + def as_refined( + occurrence: sqlite3.Row, candidate: ReactionCandidate | None + ) -> ReactionCandidate | None: + """The class's candidate with this occurrence's bonds, read without its structures, since + each atom's element is the same throughout a trajectory.""" + + if candidate is None: + return None + try: + data = _read_json(root / occurrence["bundle"] / "candidate.json")["candidate"] + return replace( + candidate, + atom_ids=tuple(data["atom_ids"]), + reactant_bonds=frozenset(map(tuple, data["reactant_bonds"])), + product_bonds=frozenset(map(tuple, data["product_bonds"])), + ) + except Exception: + unreadable["candidate bundle"] += 1 + return None + def result_file(path: Path) -> dict[str, Any] | None: try: result = _read_json(path) @@ -120,13 +140,15 @@ def result_file(path: Path) -> dict[str, Any] | None: if candidate is not None: example = occurrence["occurrence_id"] break - results = [] + results, refined = [], [] for occurrence in occurrences: path = root / "pathways" / occurrence["occurrence_id"] / "result.json" result = result_file(path) if path.is_file() else None if result is None: continue results.append(result) + validated = result["status"] == "saddle_validated" + refined.append(as_refined(occurrence, candidate) if validated else None) if result.get("steps"): event = bundle(occurrence) for step in result["steps"]: @@ -141,14 +163,14 @@ def result_file(path: Path) -> dict[str, Any] | None: for key in ("reactant_bonds", "product_bonds") }, ) - steps.append((part, 1, [step_result], occurrence["occurrence_id"])) - classes.append((candidate, len(occurrences), results, example)) + steps.append((part, 1, [step_result], occurrence["occurrence_id"], [part])) + classes.append((candidate, len(occurrences), results, example, refined)) row.update(events=len(records), pathways=dict(pathways)) notes = [f"{count} {what}{'' if count == 1 else 's'}" for what, count in unreadable.items()] if not (root / "reactions.sqlite3").is_file(): bundles = root.glob("candidates/[!.]*/candidate.json") if count := sum(1 for _ in bundles): - notes.append(f"no reactions.sqlite3 to count {count} candidates") + notes.append(f"reactions.sqlite3 is missing, with {count} on disk") if notes: row["unreadable"] = notes return classes + steps @@ -166,6 +188,7 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A # Only barriers from the same model and pressure are comparable, so classes merge within those. known: dict[tuple[Any, ...], ReactionClasses] = {} entries: dict[str, dict[str, Any]] = {} + labelled: dict[str, Any] = {} # the graph of the direction each class is named for for index, spec in enumerate(campaign.trajectories): model = campaign.adapter_profile_for(index) row: dict[str, Any] = { @@ -203,16 +226,17 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A notes.insert(0, "unreadable: " + ", ".join(row["unreadable"])) row["error"] = "; ".join(notes) or None merged = known.setdefault((model, spec.pressure_GPa), ReactionClasses(radius)) - for candidate, events, results, occurrence_id in classes: + for candidate, events, results, occurrence_id, refined in classes: # A class with no readable candidate has no graph to merge by, so it keeps its own row. name = None if candidate is None else merged.find(candidate) if name is None: name = str(len(entries)) if candidate is not None: merged.add(name, candidate) + labelled[name] = _reaction_graph(candidate, radius=radius) entries[name] = { "reaction": ( - f"class of {occurrence_id} (no readable candidate)" + f"class of {occurrence_id} (no readable candidate, direction unknown)" if candidate is None else _describe(candidate, radius) ), @@ -233,7 +257,7 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A if spec.id not in entry["trajectories"]: entry["trajectories"].append(spec.id) entry["events"] += events - for result in results: + for result, event in zip(results, refined, strict=True): entry["pathways"][result["status"]] += 1 validation = result.get("frequency_validation") if validation: @@ -242,7 +266,15 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A if connectivity: entry["connectivity"][connectivity["status"]] += 1 if result["status"] == "saddle_validated": - entry["barriers_eV"].append(result["barrier_eV"]) + key = "barriers_eV" + # A barrier is measured from each event's own reactant, so a class that runs + # both ways has two barriers. A row with no candidate names no direction. + if candidate is not None: + if event is None: + key = "unknown_direction_barriers_eV" + elif not _isomorphic(labelled[name], _reaction_graph(event, radius=radius)): + key = "reverse_barriers_eV" + entry.setdefault(key, []).append(result["barrier_eV"]) for entry in entries.values(): for key in ("pathways", "frequency", "connectivity"): entry[key] = dict(entry[key]) @@ -276,6 +308,13 @@ def _pressure(value: float | None) -> str: return "NVT" if value is None else f"{value:g}" +def _spread(barriers: list[float]) -> str: + if not barriers: + return "-" + values = (min(barriers), statistics.median(barriers), max(barriers)) + return " / ".join(f"{value:.3f}" for value in values) + + def _converged(pathways: dict[str, int]) -> str: # An event refined in steps has no band of its own; its steps are counted. total = sum(count for status, count in pathways.items() if status != "stepwise") @@ -321,9 +360,9 @@ def format_status(report: dict[str, Any]) -> str: if not report["reactions"]: lines.append("No reactions could be read." if incomplete else "No reactions recorded yet.") return "\n".join(lines) + both_ways = any("reverse_barriers_eV" in entry for entry in report["reactions"]) reaction_rows = [] for entry in report["reactions"]: - barriers = entry["barriers_eV"] reaction_rows.append( [ entry["reaction"], @@ -333,12 +372,8 @@ def format_status(report: dict[str, Any]) -> str: str(entry["events"]), _converged(entry["pathways"]), str(entry["pathways"].get("saddle_validated", 0)), - "-" - if not barriers - else " / ".join( - f"{value:.3f}" - for value in (min(barriers), statistics.median(barriers), max(barriers)) - ), + _spread(entry["barriers_eV"]), + *([_spread(entry.get("reverse_barriers_eV", []))] if both_ways else []), ] ) lines.append( @@ -352,11 +387,23 @@ def format_status(report: dict[str, Any]) -> str: "converged", "validated", "barrier eV (min / median / max)", + *(["reverse barrier eV (min / median / max)"] if both_ways else []), ], reaction_rows, ) ) lines.append("") + if both_ways: + lines.append( + "A class holds both directions of a reaction. Barrier is for the direction the row " + "names, and reverse barrier for the other." + ) + unknown = sum(len(e.get("unknown_direction_barriers_eV", [])) for e in report["reactions"]) + if unknown: + lines.append( + f"{unknown} barrier{'' if unknown == 1 else 's'} of unknown direction, whose candidate " + "files cannot be read, left out of the table and listed in --json." + ) verdicts = sum((Counter(entry["pathways"]) for entry in report["reactions"]), Counter()) if any(verdicts[status] for status in _CONVERGED): line = ( From 9b0d7808067a149c88f8157a5d872019e4e732e7 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:50:12 -0500 Subject: [PATCH 09/25] Refuse to extend a completed run instead of returning its old summary run() on a completed run with a larger total_steps returned the old summary with no message and ran no MD, so a Python caller could read it as the longer run. The command line never reached this, because the campaign contract refuses a changed total_steps first. --- src/reactionflow/run.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index 6eb0527..9bf416b 100644 --- a/src/reactionflow/run.py +++ b/src/reactionflow/run.py @@ -709,6 +709,12 @@ def run( self._start(atoms) elif atoms is not None: raise ValueError("initial atoms may only be supplied to a new run") + if self._phase == "completed" and total_steps > self._global_step: + raise ValueError( + f"this run completed at step {self._global_step} and cannot be extended to " + f"{total_steps}: its last candidates are closed. Start a new trajectory from a " + "saved frame, which carries its momenta" + ) self._refine_after_md = refine_after_md lost = missing_frames(self.root, self._global_frame) if lost: From edd6c7318083ff0b6965b153cce4eed591fc4870 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:50:18 -0500 Subject: [PATCH 10/25] Refuse a runtime that restores other atoms than its checkpoint holds MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit docs/exact-restart.md says such a provider is refused, but only the step counter was compared at a restore: a provider shifting atoms by 0.123 Å ran to completion and looked like an exact continuation. The refusal leaves the checkpoint in place, so the run resumes once the provider is fixed. --- src/reactionflow/run.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index 9bf416b..284fd12 100644 --- a/src/reactionflow/run.py +++ b/src/reactionflow/run.py @@ -438,6 +438,8 @@ def _write_frame(self, trajectory: Trajectory, atoms: Atoms, written: int) -> No def _run_segment(self, total_steps: int, runtime_provider: ExactRuntimeProvider) -> None: """Run MD from the current checkpoint until refinement is due or the run ends.""" + # A restore may hand back the snapshot's own atoms, so compare with a copy from before it. + saved_atoms = None if self._snapshot is None else self._snapshot.atoms.copy() manager = ( runtime_provider.start(self._initial) if self._snapshot is None @@ -449,6 +451,8 @@ def _run_segment(self, total_steps: int, runtime_provider: ExactRuntimeProvider) self._checkpoint(runtime) elif int(runtime.nsteps) != self._global_step: raise ValueError("exact runtime step counter does not match the durable run") + elif not _same_atomic_state(runtime.atoms, saved_atoms): + raise ValueError("exact runtime restored other atoms than its checkpoint holds") else: self._record_environment_change(runtime.snapshot()) From 1849e2d8161ad00984f1ca7dc53f18ae4331c061 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:50:24 -0500 Subject: [PATCH 11/25] Stop for a NaN energy in the model check, as the docs say The check after a refinement with no barrier tested forces only, so a model returning NaN energy with finite forces completed with exit 0 and every event neb_failed, with no mention of NaN. The campaign docs say the job stops if the model fails or returns NaN. --- src/reactionflow/run.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index 284fd12..391dbc6 100644 --- a/src/reactionflow/run.py +++ b/src/reactionflow/run.py @@ -653,8 +653,13 @@ def _check_model(self, runtime_provider: ExactRuntimeProvider) -> None: assert self._active_checkpoint is not None snapshot = ExactRestartSnapshot.read(self._active_checkpoint / "exact-restart") with runtime_provider.restore(snapshot) as runtime: - if runtime.atoms.calc is not None and not np.isfinite(runtime.atoms.get_forces()).all(): + atoms = runtime.atoms + if atoms.calc is None: + return + if not np.isfinite(atoms.get_forces()).all(): raise ValueError("the model returned non-finite forces at the last checkpoint") + if not np.isfinite(atoms.get_potential_energy()): + raise ValueError("the model returned a non-finite energy at the last checkpoint") def _refine_due(self) -> bool: """Whether to refine now: the queue waits only for running MD under ``refine_after_md``.""" From f500cd480d344fd17f74d18ba9325c9d5f39b0ff Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:50:30 -0500 Subject: [PATCH 12/25] Keep a bad saved frame from failing the whole run, and store the right reactant The search for the reactant frame ran outside the refinement's try, so a saved frame with a singular cell or no IDs marked the whole trajectory failed for good, not that one event. And after falling back from an earlier frame to the tracker's, a later failure stored the earlier frame as the reactant. --- src/reactionflow/pathway.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/reactionflow/pathway.py b/src/reactionflow/pathway.py index 6802eaa..151d525 100644 --- a/src/reactionflow/pathway.py +++ b/src/reactionflow/pathway.py @@ -753,20 +753,20 @@ def _refine_pathway( changed = candidate.reactant_bonds ^ candidate.product_bonds first, last = candidate.reactant_frame, candidate.product_frame saved = {**frames, first: candidate.reactant, last: candidate.product} - # By the tracker's reactant frame a changed bond can already be inside its hysteresis gap, and - # relaxing that frame can carry the change on, so the reactant comes from the latest frame - # before the bonds began to change. - start = next( - ( - frame - for frame in range(first, first - options.reactant_lookback - 1, -1) - if frame in saved and _on_reactant_side(saved[frame], candidate, detector_config) - ), - first, - ) - window = {frame: saved[frame] for frame in sorted(saved) if start <= frame <= last} images: list[Atoms] = [] try: + # By the tracker's reactant frame a changed bond can already be inside its hysteresis gap, + # and relaxing that frame can carry the change on, so the reactant comes from the latest + # frame before the bonds began to change. + start = next( + ( + frame + for frame in range(first, first - options.reactant_lookback - 1, -1) + if frame in saved and _on_reactant_side(saved[frame], candidate, detector_config) + ), + first, + ) + window = {frame: saved[frame] for frame in sorted(saved) if start <= frame <= last} prepared = _prepare_frames(window, options, variable_cell=setup.pressure_GPa is not None) if not set(candidate.atom_ids) <= set(atom_ids(prepared[first])): raise ValueError("candidate atom IDs are missing from its endpoints") @@ -783,6 +783,7 @@ def relaxed(frame: int, stage: str) -> tuple[Atoms, bool]: ): # Bonds away from the event differ between the two frames, so keep the tracker's. start = first + images[0] = prepared[start] reactant, converged = relaxed(start, "relax_reactant") state = _made(reactant, candidate, detector_config) if converged else None # The intermediates in the order the trajectory reached them, each with the first frame From d246284f00ff3b942d0da83a33f160270a637984 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:53:33 -0500 Subject: [PATCH 13/25] Warn when a trajectory's settings quietly shape its result These all ran with exit 0 and no message: a multi-frame structure started from its last frame; momenta in the structure gave 34845 K at a 10 K target; FixAtoms held in MD and was dropped in refinement; a left-handed or partly periodic NPT cell ran all MD and then refined nothing; and a run too short to confirm any reaction reported none. validate and a trajectory's first start now say so. Nothing here stops a run. --- docs/campaigns.md | 7 ++++++ src/reactionflow/cli.py | 49 +++++++++++++++++++++++++++++++++++++++++ tests/test_campaign.py | 4 +++- 3 files changed, 59 insertions(+), 1 deletion(-) diff --git a/docs/campaigns.md b/docs/campaigns.md index 20f8d4f..8af804b 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -240,6 +240,13 @@ with internal state, such as ASE's EMT and its neighbor list, can differ at roun does after any restart. The setting is not part of a trajectory's contract, so you can change it when you resubmit. A trajectory already paused for refinement finishes that refinement first. +`validate`, and `run` when it starts a trajectory, warn about settings that quietly shape the result, +on standard error for `run` and under `warnings` for `validate`: a structure with several frames, of +which the run takes the last; momenta in the structure, which replace `temperature_K` for the +starting velocities; constraints, which hold in MD and are removed in refinement; a cell at a pressure +that refinement cannot use; and a run with too few observations to confirm any reaction. None of +these stops the run. + ## Campaign status `reactionflow status campaign.json` summarizes a campaign without changing it, so it is safe to diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index d59e3af..9fa6076 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -138,6 +138,49 @@ def visible_gpu(environment: Mapping[str, str] | None = None) -> str: return devices[0] +def _start_warnings(campaign: CampaignConfig, index: int) -> list[str]: + """What a trajectory's settings will silently do to its result. None of this stops it.""" + + trajectory = campaign.trajectory(index) + path = campaign.structure_for(index) + name = os.path.relpath(path, campaign.source.parent) + atoms = read(path) + notes = [] + try: + several = len(read(path, index=slice(0, 2))) > 1 + except Exception: # a format that reads one image and takes no slice + several = False + if several: + notes.append(f"{name} holds several frames; a run starts from the last one") + if atoms.has("momenta"): + notes.append( + f"{name} carries momenta, which set the starting velocities in place of temperature_K; " + "the thermostat then moves them toward it" + ) + if atoms.constraints: + notes.append( + f"{name} has constraints, which hold in MD and are removed in refinement, so its " + "barriers are for a system with every atom free" + ) + if trajectory.pressure_GPa is not None and ( + not atoms.pbc.all() or atoms.cell.rank < 3 or atoms.cell.handedness <= 0 + ): + notes.append( + f"{name} has no fully periodic right-handed cell, which refinement at a pressure " + "needs, so MD can run but every event would end unresolved, with no barrier" + ) + run = campaign.reaction_run + observations = trajectory.total_steps // run.observation_interval + needed = run.detector.persistence_frames + run.candidate_stability_frames - 1 + if observations < needed: + notes.append( + f"{trajectory.total_steps} steps at one observation per {run.observation_interval} " + f"make {observations}, and confirming a reaction takes at least {needed}, so such a " + "run cannot report one" + ) + return notes + + def run_selected_trajectory( campaign: CampaignConfig, *, @@ -160,6 +203,8 @@ def run_selected_trajectory( run = ReactionRun.open(root) atoms = read(campaign.structure_for(index)) if run.phase == "new" else None else: + for note in _start_warnings(campaign, index): + print(f"warning: {note}", file=sys.stderr) atoms = read(campaign.structure_for(index)) adapter = load_mlip_adapter(adapter_spec, trajectory) preflight = getattr(adapter, "preflight", None) @@ -478,6 +523,7 @@ def main(argv: Sequence[str] | None = None) -> int: file=sys.stderr, ) elif arguments.command == "validate": + count = len(campaign.trajectories) payload = { "campaign": str(campaign.source), "trajectories": len(campaign.trajectories), @@ -488,6 +534,9 @@ def main(argv: Sequence[str] | None = None) -> int: }, "adapter_profile_counts": dict(Counter(campaign.trajectory_profiles)), } + notes = (note for index in range(count) for note in _start_warnings(campaign, index)) + if warnings := list(dict.fromkeys(notes)): + payload["warnings"] = warnings print(json.dumps(payload, sort_keys=True)) elif arguments.command == "check": if arguments.index is None and len(campaign.trajectories) > 1: diff --git a/tests/test_campaign.py b/tests/test_campaign.py index 8a78344..bf06724 100644 --- a/tests/test_campaign.py +++ b/tests/test_campaign.py @@ -313,7 +313,9 @@ def test_cli_validate_has_no_campaign_size_ceiling(tmp_path, capsys) -> None: path = _campaign(tmp_path, count=130) assert main(["validate", str(path)]) == 0 - assert json.loads(capsys.readouterr().out) == { + report = json.loads(capsys.readouterr().out) + assert len(report.pop("warnings")) == 2 # this toy has no cell at its pressure, and few steps + assert report == { "adapter_profile_counts": {"test": 130}, "structure_atoms": {"structure.extxyz": 1}, "campaign": str(path), From ae1aeb525c506f9d593cb8f32b935528041b00e0 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:55:36 -0500 Subject: [PATCH 14/25] Warn when a scanned file names atoms only by LAMMPS type, or has no cell A LAMMPS dump with id, type and coordinates reads carbon as hydrogen, and reported 0 bond changes with exit 0 where the file holds a C-C break: a false zero on the commonest dump layout. A file with no cell gave a false formed bond when an atom wrapped at a boundary. The counts do not change. --- docs/offline-scan.md | 5 +++- src/reactionflow/cli.py | 52 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 54 insertions(+), 3 deletions(-) diff --git a/docs/offline-scan.md b/docs/offline-scan.md index dcc9659..2ee9aea 100644 --- a/docs/offline-scan.md +++ b/docs/offline-scan.md @@ -40,7 +40,10 @@ Any source that is not a `.json` campaign file is read as a trajectory file, in reads, for example an XDATCAR, an extxyz file, or a LAMMPS dump. Atoms keep the IDs the file carries, such as an `atom_id` column in extxyz, and otherwise get IDs 0 to N−1 in file order, so the file must keep its atom order. ASE sorts a LAMMPS dump by its `id` column, and it names the -elements from an `element` or `mass` column, because types alone do not. The settings start from +elements from an `element` or `mass` column, because types alone do not: a dump with only types is +read as H, He, and so on, and the scan warns, on standard error and under `warnings` in the summary. +It also warns on standard error about a file with no periodic cell, where a wrapped atom reads as a +bond change. The settings start from the defaults of a campaign (thresholds from covalent radii, persistence 3, stability 3, radius 2) and take the same options. `--frame-ps` gives the time between frames for the time columns, which are empty without it. The tables go to `-scan/` beside the file unless `--out` names another diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 9fa6076..2f3ae4b 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -10,14 +10,16 @@ import socket import sys from collections import Counter -from collections.abc import Mapping, Sequence +from collections.abc import Iterator, Mapping, Sequence from contextlib import suppress from dataclasses import asdict, replace from numbers import Integral from pathlib import Path from uuid import uuid4 +from ase import Atoms from ase.io import read +from ase.io.formats import filetype from ._durable import ensure_directory, file_digest, publish, sync_directory from ._version import __version__ @@ -318,6 +320,46 @@ def _scan_settings( } +def _file_warnings(path: Path, first: Atoms) -> list[str]: + """What a trajectory file's form does to a scan of it, read off its first frame, to keep with + the scan's summary. A file with no cell is often a gas-phase run, so that note is only printed. + """ + + if not first.pbc.any(): + print( + f"warning: {path.name} has no periodic cell, so an atom wrapped at a box boundary " + "reads as a bond change. A gas-phase trajectory is fine; otherwise use a format " + "that keeps it", + file=sys.stderr, + ) + notes = [] + try: + dump = filetype(str(path)) == "lammps-dump-text" + except Exception: # the frames were read, so this only costs the note below + dump = False + if dump: + with path.open(encoding="utf-8", errors="replace") as stream: + columns = next( + (line.split()[2:] for line in stream if line.startswith("ITEM: ATOMS")), [] + ) + if "element" not in columns and "mass" not in columns: + notes.append( + f"{path.name} names atoms only by LAMMPS type, which ASE reads as atomic numbers " + f"(type 1 as H, type 2 as He), so these results are for " + f"{first.get_chemical_formula()}. Add an element or mass column to the dump" + ) + return notes + + +def _noting_first_frame( + path: Path, frames: Iterator[tuple[int, float | None, Atoms]], notes: list[str] +) -> Iterator[tuple[int, float | None, Atoms]]: + for number, item in enumerate(frames): + if number == 0: + notes.extend(_file_warnings(path, item[2])) + yield item + + def _print_scan(path: Path, result: ScanResult, **identity: object) -> None: counts = {key: value for key, value in result.summary.items() if key != "settings"} print(json.dumps({"scan": str(path), **identity, **counts}, sort_keys=True)) @@ -451,9 +493,15 @@ def main(argv: Sequence[str] | None = None) -> int: if frame_ps is not None and not (math.isfinite(frame_ps) and frame_ps > 0): parser.error("--frame-ps must be a positive number") settings = _scan_settings(parser, arguments, ReactionRunConfig()) - result = scan_frames(file_frames(arguments.campaign, frame_ps=frame_ps), **settings) source = arguments.campaign + notes: list[str] = [] + frames = _noting_first_frame(source, file_frames(source, frame_ps=frame_ps), notes) + result = scan_frames(frames, **settings) + if notes: + result.summary["warnings"] = notes _print_scan(result.write(arguments.out or source.with_name(f"{source.name}-scan")), result) + for note in notes: + print(f"warning: {note}", file=sys.stderr) return 0 raw = json.loads(arguments.campaign.read_text(encoding="utf-8")) if not isinstance(raw, dict): From 133425ab61c3e3d45889159be5ab04a1280063f1 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:56:02 -0500 Subject: [PATCH 15/25] Refuse eon.displace_centers on a push that ignores it displace_centers applies only to a local push, and displace defaults to all, so a campaign setting only the centers pushed every atom without a word, and its exploration was not the restricted one it asked for. --- src/reactionflow/search/eon.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/reactionflow/search/eon.py b/src/reactionflow/search/eon.py index ae9d84f..f2a07ff 100644 --- a/src/reactionflow/search/eon.py +++ b/src/reactionflow/search/eon.py @@ -148,6 +148,11 @@ def __post_init__(self) -> None: raise ValueError("eon.displace_centers holds element symbols or atom indices") if isinstance(center, Integral) and center < 0: raise ValueError("eon.displace_centers atom indices must be non-negative") + if centers and self.displace != "local" and self.guide is None: + raise ValueError( + "eon.displace_centers applies only to a local push; " + "set eon.displace to 'local', or leave the centers out" + ) object.__setattr__(self, "displace_centers", tuple(centers)) if type(self.equivalent_atoms) is not bool: raise ValueError("eon.equivalent_atoms must be true or false") From 76846a009c67d47f92775740c4bfb0cc728d884b Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:57:19 -0500 Subject: [PATCH 16/25] Refuse a model file replaced during a run on the generic adapter The generic adapter hashed its model files once per process, so a file replaced between MD and refinement mixed two potentials under one recorded identity with no message. The Torch adapters already rehash before every calculator and refuse this. --- src/reactionflow/adapters/ase.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/reactionflow/adapters/ase.py b/src/reactionflow/adapters/ase.py index bee7051..ed69a57 100644 --- a/src/reactionflow/adapters/ase.py +++ b/src/reactionflow/adapters/ase.py @@ -244,14 +244,17 @@ def __init__(self, *, trajectory: TrajectorySpec, options: Mapping[str, Any]) -> self._contract: ComponentState | None = None def _model_file_hashes(self) -> dict[str, str]: - if self._model_hashes is None: - hashes: dict[str, str] = {} - for path in self.model_files: - if not path.is_file(): - raise FileNotFoundError(f"MLIP model file is missing: {path}") - hashes[str(path)] = file_digest(path) - self._model_hashes = hashes - return self._model_hashes + # Rehash before every new calculator: replacing a file mid-run must not let MD and + # refinement use different potentials under one recorded identity. + hashes: dict[str, str] = {} + for path in self.model_files: + if not path.is_file(): + raise FileNotFoundError(f"MLIP model file is missing: {path}") + hashes[str(path)] = file_digest(path) + if self._model_hashes is not None and hashes != self._model_hashes: + raise ValueError("model files changed during the run") + self._model_hashes = hashes + return hashes def _package_versions(self, module_name: str) -> dict[str, str]: names = set(self.packages) From 6b350cf98993c60a98b0a9a7621b56f366abdeb8 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:58:37 -0500 Subject: [PATCH 17/25] Say in status why events ended with no band status read converged 2/9 and its verdict line was silent about 6 collapsed and 1 unresolved, or about a transient model failure, whose text was only in that event's result.json. A left-handed cell's 'SSNEB requires fully periodic right-handed cells' was found the same way. Now they are counted by status and message, up to five in the text and all under no_band in --json. --- docs/campaigns.md | 2 ++ src/reactionflow/status.py | 19 ++++++++++++++++++- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/docs/campaigns.md b/docs/campaigns.md index 8af804b..4b0e0e1 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -277,6 +277,8 @@ other. Each step of an event refined in steps is an event of its own local class saddles ([saddle check](pathway-refinement.md#saddle-check)). A line below the table totals the saddle check's verdicts across the campaign (validated, no saddle, rejected, unchecked) and counts separately the bands that converged before the check existed, which carry no verdict. +Events that ended with no converged band are counted by status and message, up to five in the text +and all of them as `no_band` in `--json`. NVT barriers are potential energies; NPT barriers are enthalpies. `--json` prints the same data, including every barrier value, the counts of each pathway, frequency and connectivity status, trajectory ID, and one example occurrence per class, for scripted analysis. diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index 964e3ee..63b9d8e 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -189,6 +189,7 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A known: dict[tuple[Any, ...], ReactionClasses] = {} entries: dict[str, dict[str, Any]] = {} labelled: dict[str, Any] = {} # the graph of the direction each class is named for + no_band: Counter[tuple[str, str]] = Counter() for index, spec in enumerate(campaign.trajectories): model = campaign.adapter_profile_for(index) row: dict[str, Any] = { @@ -259,6 +260,8 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A entry["events"] += events for result, event in zip(results, refined, strict=True): entry["pathways"][result["status"]] += 1 + if result["status"] not in _CONVERGED | {"stepwise", "saddle_validated"}: + no_band[result["status"], result.get("message") or "no message"] += 1 validation = result.get("frequency_validation") if validation: entry["frequency"][validation["status"]] += 1 @@ -288,12 +291,18 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A item["reaction"], ), ) - return { + report = { "campaign": str(campaign.source), "radius": radius, "trajectories": trajectories, "reactions": reactions, } + if no_band: + report["no_band"] = [ + {"status": status, "message": message, "events": count} + for (status, message), count in no_band.most_common() + ] + return report def _table(header: list[str], rows: list[list[str]]) -> str: @@ -425,6 +434,14 @@ def format_status(report: dict[str, Any]) -> str: f"{verdicts['time_limit']} refinement{'' if verdicts['time_limit'] == 1 else 's'}" " stopped at max_hours_per_barrier, with no barrier." ) + ended = [item for item in report.get("no_band", []) if item["status"] != "time_limit"] + if ended: + lines.append( + f"Ended with no converged band: {sum(item['events'] for item in ended)} events." + ) + lines += [f" {i['events']} {i['status']}: {i['message'][:120]}" for i in ended[:5]] + if len(ended) > 5: + lines.append(f" and {len(ended) - 5} other reasons, listed in --json as no_band.") radius = report["radius"] lines.append( f"A class is a bond change with the atoms within {radius} bond{'' if radius == 1 else 's'}" From db6fc54ed061946b4611e09c97de609f698c78b9 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:58:45 -0500 Subject: [PATCH 18/25] Count events waiting to be refined in status status showed converged 2/2 while four events sat in the queue, because it read only the phase, step and failure from state.json. With refine_after_md a whole run's events wait there, so converged x/y could read as everything done. --- docs/campaigns.md | 3 ++- src/reactionflow/status.py | 7 +++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/docs/campaigns.md b/docs/campaigns.md index 4b0e0e1..930ee4c 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -262,7 +262,8 @@ table says the report is incomplete. A class whose candidate bundles cannot all of a run's results alone, keeps its counts and barriers in a row of its own, because classes merge only by their graphs. Saved frames that are gone appear in the error column and as `missing_frames`. Pathway counts include every refined occurrence. An event refined in steps counts each step as a -refined pathway, and only its steps can converge. +refined pathway, and only its steps can converge. Events still waiting to be refined are not among +them: a line under the table counts them, and `--json` gives `queued`. The second table merges reaction classes across trajectories that share a model profile and pressure, because only those barriers are comparable. Classes merge by the local change: the diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index 63b9d8e..a00c5c5 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -84,6 +84,8 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: lost = missing_frames(root, state["global_frame"]) if lost: row["missing_frames"] = lost + if state.get("pending_pathway_ids"): + row["queued"] = len(state["pending_pathway_ids"]) members: dict[str, list[sqlite3.Row]] = {} for record in records: members.setdefault(record["class_id"], []).append(record) @@ -358,6 +360,11 @@ def format_status(report: dict[str, Any]) -> str: + ", ".join(f"{count} {phase}" for phase, count in sorted(phases.items())), "", ] + queued = [row for row in report["trajectories"] if row.get("queued")] + if queued: + where = ", ".join(f"{row['queued']} in {row['trajectory']}" for row in queued[:5]) + more = f", and {len(queued) - 5} more trajectories" if len(queued) > 5 else "" + lines += [f"Waiting to be refined, not yet in the converged column: {where}{more}.", ""] incomplete = sum("unreadable: " in (row["error"] or "") for row in report["trajectories"]) if incomplete: lines += [ From 6fd9152356b86627f9211436cc3de9467859a124 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:59:24 -0500 Subject: [PATCH 19/25] Show a resume under another environment in status A resume with a changed package version or calculator source warned once in a job log and added a line to restart-environment.jsonl, and status stayed clean. A run whose calculator source was edited to double its energy and forces at step 340 read the same as one that was never touched. --- src/reactionflow/status.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index a00c5c5..ddb8977 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -86,6 +86,14 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: row["missing_frames"] = lost if state.get("pending_pathway_ids"): row["queued"] = len(state["pending_pathway_ids"]) + changes = root / "restart-environment.jsonl" + if changes.is_file(): + try: + lines = [json.loads(line) for line in changes.read_text().splitlines()] + fields = sorted({field for line in lines for field in line["changed"]}) + row["environment_changes"] = {"resumes": len(lines), "fields": fields} + except Exception: + unreadable["restart-environment.jsonl"] += 1 members: dict[str, list[sqlite3.Row]] = {} for record in records: members.setdefault(record["class_id"], []).append(record) @@ -168,7 +176,10 @@ def result_file(path: Path) -> dict[str, Any] | None: steps.append((part, 1, [step_result], occurrence["occurrence_id"], [part])) classes.append((candidate, len(occurrences), results, example, refined)) row.update(events=len(records), pathways=dict(pathways)) - notes = [f"{count} {what}{'' if count == 1 else 's'}" for what, count in unreadable.items()] + notes = [ + what if what.endswith("jsonl") else f"{count} {what}{'' if count == 1 else 's'}" + for what, count in unreadable.items() + ] if not (root / "reactions.sqlite3").is_file(): bundles = root.glob("candidates/[!.]*/candidate.json") if count := sum(1 for _ in bundles): @@ -365,6 +376,15 @@ def format_status(report: dict[str, Any]) -> str: where = ", ".join(f"{row['queued']} in {row['trajectory']}" for row in queued[:5]) more = f", and {len(queued) - 5} more trajectories" if len(queued) > 5 else "" lines += [f"Waiting to be refined, not yet in the converged column: {where}{more}.", ""] + for row in report["trajectories"]: + if "environment_changes" in row: + changes = row["environment_changes"] + lines += [ + f"{row['trajectory']} resumed {changes['resumes']} time(s) under another " + f"environment, so it is no longer bitwise exact: {', '.join(changes['fields'])}. " + "A changed calculator source can mean a changed potential.", + "", + ] incomplete = sum("unreadable: " in (row["error"] or "") for row in report["trajectories"]) if incomplete: lines += [ From 867b1a3b4547634d65074f23a4b4a25d5c1264d7 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 10:59:39 -0500 Subject: [PATCH 20/25] Label each scan row by its own direction candidates.csv gave a class's first-seen direction as every row's label, so broken C0-O1 could sit under formed C-O. The class number still merges both directions. --- docs/offline-scan.md | 5 +++-- src/reactionflow/scan.py | 11 ++++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/docs/offline-scan.md b/docs/offline-scan.md index 2ee9aea..dfb0880 100644 --- a/docs/offline-scan.md +++ b/docs/offline-scan.md @@ -57,8 +57,9 @@ directory. first frame. - `candidates.csv` has one row per candidate: the reactant, product and observed frames and the times of the first two, resolved or not, whether it was still changing at the end, its number of - atoms, its bond changes, its local class and label as `reactionflow status` names them (radius 2), - and the run's occurrence ID when the run recorded the same candidate. + atoms, its bond changes, its local class as `reactionflow status` merges them (radius 2), the + label of that class in this row's direction, and the run's occurrence ID when the run recorded the + same candidate. - `species.csv` counts each formula among the bonded components every `--every` frames, with the number of starting molecules in each component, which is the chain length in a polymer. - `summary.json` holds the settings and the counts, which are also printed: diff --git a/src/reactionflow/scan.py b/src/reactionflow/scan.py index 81d0d59..9872bd0 100644 --- a/src/reactionflow/scan.py +++ b/src/reactionflow/scan.py @@ -180,7 +180,7 @@ def scan_frames( detector = BondChangeDetector(detector_config) tracker = ReactionTracker(stability_frames=stability_frames, region_radius=region_radius) classes = ReactionClasses(radius=2) - labels: list[str] = [] + class_count = 0 times: dict[int, float | None] = {} symbols: dict[int, str] = {} molecule = nx.utils.UnionFind() @@ -193,11 +193,11 @@ def scan_frames( read_to = -1 def add(candidate: ReactionCandidate, *, at_end: bool) -> None: - nonlocal net_changes + nonlocal net_changes, class_count name = classes.find(candidate) if name is None: - labels.append(_describe(candidate, 2)) - name = str(len(labels)) + class_count += 1 + name = str(class_count) classes.add(name, candidate) changed = candidate.reactant_bonds ^ candidate.product_bonds if candidate.resolved: @@ -215,7 +215,8 @@ def add(candidate: ReactionCandidate, *, at_end: bool) -> None: "atoms": len(candidate.atom_ids), "changes": _changes(candidate, symbols), "class": int(name), - "class_label": labels[int(name) - 1], + # A class holds both directions, so the label is this row's own. + "class_label": _describe(candidate, 2), "occurrence_id": (occurrences or {}).get((window, frozenset(changed)), ""), } ) From 1f8bc72651b55a16a95dea3cfcbc5dfa6960ec11 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 11:00:28 -0500 Subject: [PATCH 21/25] Correct what the docs promise about exact restart and foreign bands docs/exact-restart.md named a strict mode that does not exist and did not say that nothing tests a calculator's exactness, while the README's own xTB example leaves it after each pause or kill. The campaign docs called that a rounding-level difference; docs/models.md said a changed device continues when the campaign's device option is refused; and the README did not say check ignores residual forces. --- README.md | 7 +++++-- docs/campaigns.md | 5 +++-- docs/exact-restart.md | 11 ++++++++--- docs/models.md | 7 ++++--- examples/xtb/README.md | 4 +++- 5 files changed, 23 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 2019df1..b4ab62f 100644 --- a/README.md +++ b/README.md @@ -66,12 +66,15 @@ wait in a queue until then ([when barriers are refined](docs/campaigns.md#when-b Two more commands work on finished output. `reactionflow scan` replays the saved frames through bond detection, with the run's settings or new ones, and tabulates every bond change and species ([offline scan](docs/offline-scan.md)). `reactionflow check campaign.json band.traj` tests whether -the highest image of any band, from ReactionFlow or another code, is a saddle between its ends. +the highest image of any band, from ReactionFlow or another code, is a saddle between its ends. It +reports the band's residual forces and does not require them to be small, so a barrier from another +code's band means what it says only for relaxed ends and a stationary highest image. A number for `pressure_GPa` runs NPT, and the pathway is then refined with variable-cell SSNEB and gives an enthalpy barrier. `null` runs NVT with ordinary NEB. Each trajectory writes its frames, every detected occurrence, every refinement with its band and saddle check, and a -`reactions.sqlite3` index under `output_root//`. Nothing is overwritten, and a +`reactions.sqlite3` index under `output_root//`. Published events and refinements are +never overwritten, though a rerun of `scan` or `check` replaces its own tables or verdict, and a resubmitted job continues each unfinished trajectory from its last checkpoint. Refining one event stops after `max_hours_per_barrier` hours, 4 by default, and records diff --git a/docs/campaigns.md b/docs/campaigns.md index 930ee4c..e06f96f 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -236,8 +236,9 @@ marked failed refines its queue when resubmitted, then reports the failure again MD restarts exactly after a pause, and refinement reads only saved frames, so with a model that keeps no state between evaluations both orders give the same frames, events, and barriers. A model -with internal state, such as ASE's EMT and its neighbor list, can differ at rounding level, as it -does after any restart. The setting is not part of a trajectory's contract, so you can change it +with internal state can differ, as it does after any restart: ASE's EMT and its neighbor list at +rounding level, and tblite's xTB, which reuses its wavefunction, by a whole different path +([exact restart](exact-restart.md)). The setting is not part of a trajectory's contract, so you can change it when you resubmit. A trajectory already paused for refinement finishes that refinement first. `validate`, and `run` when it starts a trajectory, warn about settings that quietly shape the result, diff --git a/docs/exact-restart.md b/docs/exact-restart.md index 10cd1e3..88762a1 100644 --- a/docs/exact-restart.md +++ b/docs/exact-restart.md @@ -40,13 +40,18 @@ Bitwise continuation requires the same ASE and NumPy versions, so the codec reco do not influence the propagated state. Calculator state is not guessed. An MLIP adapter must either serialize all mutable state or -declare that deterministic reconstruction from pinned model/configuration data is exact. An -adapter that cannot satisfy that contract cannot be used in strict exact-restart mode. +declare that deterministic reconstruction from pinned model/configuration data is exact. The generic ASE calculator adapter is for deterministic, stateless inference. It binds the calculator factory and kwargs, model-file SHA-256 values, calculator source, and selected installed package versions to the checkpoint. Calculators with mutable inference state or their own RNG must -use a custom adapter that snapshots that state explicitly. +use a custom adapter that snapshots that state explicitly. Nothing tests this at run time: the +adapter labels every calculator exact and gives no warning when one is not. tblite reuses its +previous wavefunction, so each pause or kill of the xTB example changes its path: 6e-8 Å by step 100, +growing to 6.3 Å within 20 ps, with 7 events uninterrupted, 8 with `refine_after_md`, and 10 after +three kills, on one machine. Each is still a valid sample, but not the same trajectory. ASE's EMT and +LJ differ at rounding level, at most 7e-15 Å over 3 ps. ANI-1xnr, AIMNet2-RXN, MACE and UMA on a CPU +stayed bitwise exact in the cases tried; no GPU has been tested. The built-in Torch adapters, ANI-1xnr included, use deterministic PyTorch execution, disable TF32, and fix the cuBLAS workspace configuration. They record the model-file SHA-256 values, options, diff --git a/docs/models.md b/docs/models.md index 99b2ff9..3989405 100644 --- a/docs/models.md +++ b/docs/models.md @@ -303,9 +303,10 @@ every pathway stage. For CPU runs, set `device: "cpu"` and campaign `require_gpu worker must see exactly one usable GPU. Each checkpoint records the model files, options, package versions, adapter/runtime source, and the -Torch/platform/device environment. A changed model file or option refuses the resume. A changed -version, source hash, or device continues with a warning and a line in -`restart-environment.jsonl`, because the run is then no longer bitwise exact +Torch/platform/device environment. A changed model file or option refuses the resume, and that +includes the campaign's `device` option, since it is part of the trajectory's contract. A changed +version or source hash, or other hardware behind the same option, continues with a warning and a line +in `restart-environment.jsonl`, because the run is then no longer bitwise exact ([exact restart](exact-restart.md)). Keep the environment of an active campaign when you can. For Torch backends, deterministic algorithms are enabled and TF32 is disabled; unsupported nondeterministic operations fail. UMA uses fixed float32 inference without compilation or MoLE merging, so its diff --git a/examples/xtb/README.md b/examples/xtb/README.md index f563563..3cc264c 100644 --- a/examples/xtb/README.md +++ b/examples/xtb/README.md @@ -31,7 +31,9 @@ Saddle check: 9 validated, 0 no saddle, 0 rejected, 0 unchecked. Barriers come o Each proton transfer is one occurrence of the same reaction class, the O-H bond that breaks and the O-H bond that forms in malonaldehyde. Every saddle has one imaginary mode near 1150i cm⁻¹ and connects the two enols. The number of events and their barriers change with the machine, because -floating-point differences grow along a chaotic trajectory. Eight 20 ps runs from 400 to 800 K +floating-point differences grow along a chaotic trajectory, and with every pause or kill, because +tblite reuses its wavefunction ([exact restart](../../docs/exact-restart.md)). On one machine this +example gave 7 events uninterrupted and 10 after three kills. Eight 20 ps runs from 400 to 800 K gave 55 transfers, all validated, at 0.04 to 0.19 eV. Each refinement is in `pathways//`: `images.traj` holds the band, and `result.json` From 0c551bca5aad69bee7807b805972b0ea287b3477 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 13:42:04 -0500 Subject: [PATCH 22/25] Stop an EON run when the model itself has stopped working Every exception in a search became a failed search with no limit, so a model that broke after the start relaxed wrote failed records until the job's time limit: 880 failed searches and 3.7 MB in 8 s with a toy calculator on real EON, and no reason anywhere but the attempt records. If the model raises or returns a non-finite value during a search, it is now evaluated once on that search's starting structure, which it evaluated before. If that fails too, the run exits 1 without recording the search, and a rerun repeats it with its saved seed. If that succeeds, the failure belonged to the geometry: the search is recorded as failed, with what the model raised, and the run goes on. There is no count or threshold. The call into the model records the failure itself, because EON's own code can report it as an ordinary ending. For the same reason a Ctrl-C inside a model call is raised again once EON returns, rather than written down as a failed search that changes every later record. --- docs/eon-search.md | 7 ++- src/reactionflow/search/eon.py | 81 +++++++++++++++++++++++++++++----- 2 files changed, 75 insertions(+), 13 deletions(-) diff --git a/docs/eon-search.md b/docs/eon-search.md index 04c450d..5d44f63 100644 --- a/docs/eon-search.md +++ b/docs/eon-search.md @@ -222,7 +222,8 @@ exactly once, and an interrupted search is repeated with its saved seed. Keeping does not promise bitwise agreement across EON builds or numerical environments. The output directory holds the model/settings contract in `search-contract.json`, progress in -`search-state.json`, and immutable checksummed records under `states/`, `processes/`, +`search-state.json`, the error that ended the last run, if any, in `last-error.json`, and immutable +checksummed records under `states/`, `processes/`, `attempts/`, and `steps/`. A process record keeps its saddle, product, barrier, prefactors, and bond changes; the reverse of each process is recorded for the product state. An attempt that EON rejected after its saddle converged keeps the saddle, both minima with their energies, and the @@ -248,6 +249,10 @@ output directory. shows the searches and confidence of the current state, and counts the failed searches by reason (`failure_reasons` in `--json`). `run` prints each new reason on standard error as it first appears, up to ten. A confidence that stays at 0 while searches fail will not recover by itself. + The one automatic stop: if the model raises or returns a non-finite value during a search, it is + evaluated once on that search's starting structure. If that fails too, the run exits with status + 1 without recording the search, and a rerun repeats it with its saved seed. If that succeeds, the + failure belonged to the geometry, so the search is recorded as failed and the run goes on. - Searches use a fixed, nonsingular cell, fully periodic or fully nonperiodic, with `FixAtoms` as the only constraint. - By default EON converges saddles and minima on the norm of the whole force vector diff --git a/src/reactionflow/search/eon.py b/src/reactionflow/search/eon.py index f2a07ff..f07ce59 100644 --- a/src/reactionflow/search/eon.py +++ b/src/reactionflow/search/eon.py @@ -15,7 +15,7 @@ import math import tempfile from collections.abc import Mapping -from dataclasses import asdict, dataclass, field +from dataclasses import asdict, dataclass, field, replace from importlib import import_module from importlib.metadata import version from numbers import Integral, Real @@ -352,14 +352,18 @@ class _PreservingCalculator(Calculator): implemented_properties = ("energy", "forces") - def __init__(self, calculator: Calculator, template: Atoms): + def __init__(self, calculator: Calculator, template: Atoms, failures: list[BaseException]): super().__init__() self.calculator = calculator self.template = template.copy() self.template.info = copy.deepcopy(template.info) + self.failures = failures def calculate(self, atoms=None, properties=None, system_changes=all_changes): super().calculate(atoms, properties, system_changes) + for failure in self.failures: + if not isinstance(failure, Exception): # Ctrl-C: end this search without the model + raise failure if not np.array_equal(atoms.numbers, self.template.numbers): raise ValueError("EON changed atom order or species") # EON passes a zero force-evaluation cell for isolated molecules while @@ -370,13 +374,19 @@ def calculate(self, atoms=None, properties=None, system_changes=all_changes): raise ValueError("EON changed the fixed cell") physical = _moved(self.template, atoms.positions) physical.calc = self.calculator - energy = physical.get_potential_energy() - # EON applies the fixed-coordinate mask itself. - forces = physical.get_forces(apply_constraint=False) - if not np.isfinite(energy) or forces.shape != physical.positions.shape: - raise ValueError("calculator returned invalid energy or force shape") - if not np.isfinite(forces).all(): - raise ValueError("calculator returned nonfinite forces") + try: + energy = physical.get_potential_energy() + # EON applies the fixed-coordinate mask itself. + forces = physical.get_forces(apply_constraint=False) + if not np.isfinite(energy) or forces.shape != physical.positions.shape: + raise ValueError("calculator returned invalid energy or force shape") + if not np.isfinite(forces).all(): + raise ValueError("calculator returned nonfinite forces") + except BaseException as error: + # EON's own code can take this for an ordinary failed stage and carry on, so it is + # noted here for the search to look at afterwards. See EONBackend.search. + self.failures.append(error) + raise self.results = {"energy": energy, "forces": forces} @@ -415,6 +425,7 @@ def __init__(self, calculator: Calculator, settings: EONSettings, temperature_K: self.calculator = calculator self.settings = settings self.temperature_K = float(temperature_K) + self._failures: list[BaseException] = [] self._eon = _load_eon() # Comparisons never evaluate the potential; they need only the comparison settings. self._comparison = { @@ -474,7 +485,7 @@ def _parameters( return parameters def _matter(self, atoms: Atoms, fixed: np.ndarray, parameters): - bridge = _PreservingCalculator(self.calculator, atoms) + bridge = _PreservingCalculator(self.calculator, atoms, self._failures) potential = self._eon.potential_from_ase(bridge) geometry = atoms.copy() geometry.set_constraint() @@ -518,9 +529,54 @@ def relax(self, atoms: Atoms) -> tuple[Atoms, float]: raise ValueError("starting structure did not relax to a minimum") return _moved(source, relaxed.positions), float(relaxed.potential_energy) - def search(self, atoms: Atoms, *, seed: int) -> ProcessResult: - """Run one process search from a relaxed state; failures are results, not errors.""" + def _raise_interrupt(self) -> None: + """Ctrl-C inside a model call, which EON's own code can take for a failed stage.""" + + for failure in self._failures: + if not isinstance(failure, Exception): + raise failure + + def _check_model(self, atoms: Atoms) -> None: + """After the model failed in a search, stop the run if it fails on the search's start too. + + The start is a structure it has evaluated before, so this tells a model that has stopped + working from a geometry it cannot handle. A failing search would otherwise be recorded, and + the next one started, for as long as the job lasts. + """ + physical = canonical_copy(atoms) + physical.calc = self.calculator + try: + energy = physical.get_potential_energy() + forces = physical.get_forces(apply_constraint=False) + if not (np.isfinite(energy) and np.isfinite(forces).all()): + raise ValueError("calculator returned nonfinite energy or forces") + except Exception as error: + raise RuntimeError( + "the model fails on this search's starting structure, which it evaluated before, " + "so the search was not recorded and a rerun repeats it: " + f"{type(error).__name__}: {error}" + ) from error + + def search(self, atoms: Atoms, *, seed: int) -> ProcessResult: + """Run one process search from a relaxed state. + + Failures are results, not errors, unless the model itself has stopped working or the search + was interrupted. Then this raises, and the search is repeated with its seed on a rerun. + """ + + self._failures.clear() + result = self._search(atoms, seed) + self._raise_interrupt() + if self._failures: + self._check_model(atoms) + # EON may have reported this as an ordinary ending, such as too many iterations. + raised = f"{type(self._failures[0]).__name__}: {self._failures[0]}" + if raised not in result.message: + result = replace(result, message=f"{result.message} (the model raised {raised})") + return result + + def _search(self, atoms: Atoms, seed: int) -> ProcessResult: metadata: dict[str, Any] = {"seed": int(seed)} try: source = canonical_copy(atoms) @@ -658,6 +714,7 @@ def prefactors(self, first: Atoms, saddle: Atoms, second: Atoms) -> tuple[float, with tempfile.TemporaryDirectory() as work, contextlib.chdir(work): values = tuple(float(v) for v in self._eon.get_prefactors(parameters, *matters)) except Exception: + self._raise_interrupt() return None low, high = parameters.prefactor_min_value, parameters.prefactor_max_value return values if all(low <= value <= high for value in values) else None From d525749fab04f2e5c0c18e14be0307cbafacd2a6 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 13:53:38 -0500 Subject: [PATCH 23/25] Report what a damaged copy of a run lacks more exactly in status A step's result.json without its bonds, or a state.json from before global_frame, still dropped the whole trajectory's reactions. One lost bundle was counted twice. A copy with results but no database or candidates read as empty, and a trajectory with records but no state.json as not started. A last-error.json holding something other than a string crashed the whole report. --- src/reactionflow/status.py | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index ddb8977..520e18f 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -81,7 +81,7 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: row["error"] = f"{failure.get('type')}: {failure.get('message')}" records = _occurrences(root) unreadable: Counter[str] = Counter() - lost = missing_frames(root, state["global_frame"]) + lost = missing_frames(root, state.get("global_frame", 0)) if lost: row["missing_frames"] = lost if state.get("pending_pathway_ids"): @@ -102,6 +102,7 @@ def _trajectory(root: Path, row: dict[str, Any]) -> list[tuple[Any, ...]]: pathways: Counter[str] = Counter() read: dict[str, ReactionCandidate | None] = {} + lost_bundles: set[str] = set() def bundle(occurrence: sqlite3.Row) -> ReactionCandidate | None: name = occurrence["bundle"] @@ -110,7 +111,7 @@ def bundle(occurrence: sqlite3.Row) -> ReactionCandidate | None: read[name] = _read_bundle(root / name)[1] except Exception: read[name] = None - unreadable["candidate bundle"] += 1 + lost_bundles.add(name) return read[name] def as_refined( @@ -130,7 +131,7 @@ def as_refined( product_bonds=frozenset(map(tuple, data["product_bonds"])), ) except Exception: - unreadable["candidate bundle"] += 1 + lost_bundles.add(occurrence["bundle"]) return None def result_file(path: Path) -> dict[str, Any] | None: @@ -166,24 +167,32 @@ def result_file(path: Path) -> dict[str, Any] | None: step_result = result_file(step_path) if step_result is None: continue - part = event and replace( - event, - **{ + try: + bonds = { key: frozenset(map(tuple, step_result[key])) for key in ("reactant_bonds", "product_bonds") - }, - ) + } + except Exception: + unreadable["result file"] += 1 + continue + part = event and replace(event, **bonds) steps.append((part, 1, [step_result], occurrence["occurrence_id"], [part])) classes.append((candidate, len(occurrences), results, example, refined)) row.update(events=len(records), pathways=dict(pathways)) + if lost_bundles: + unreadable["candidate bundle"] = len(lost_bundles) notes = [ what if what.endswith("jsonl") else f"{count} {what}{'' if count == 1 else 's'}" for what, count in unreadable.items() ] if not (root / "reactions.sqlite3").is_file(): - bundles = root.glob("candidates/[!.]*/candidate.json") - if count := sum(1 for _ in bundles): - notes.append(f"reactions.sqlite3 is missing, with {count} on disk") + count = sum(1 for _ in root.glob("candidates/[!.]*/candidate.json")) + results = sum(1 for _ in root.glob("pathways/[!.]*/result.json")) + if count or results: + notes.append( + f"reactions.sqlite3 is missing, with {count} candidates and {results} results " + "on disk" + ) if notes: row["unreadable"] = notes return classes + steps @@ -225,6 +234,11 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A classes = _trajectory(root, row) except Exception as error: # a damaged trajectory must not hide the rest row["error"] = f"unreadable: {type(error).__name__}: {error}" + elif any( + any(root.glob(f"{kind}/[!.]*")) for kind in ("candidates", "pathways", "segments") + ): + row["phase"] = "unknown" + row["unreadable"] = ["state.json is missing, with records on disk"] # A worker can fail before it writes state.json, so its last error is read for every # trajectory, and read last so that a damaged file cannot hide the reactions above. if row["error"] is None and (root / "last-error.json").is_file(): @@ -232,7 +246,7 @@ def campaign_status(campaign: CampaignConfig, *, radius: int = 2) -> dict[str, A row["error"] = _read_json(root / "last-error.json").get("error") except Exception as error: row["error"] = f"unreadable last-error.json: {type(error).__name__}" - notes = [row["error"]] if row["error"] else [] + notes = [str(row["error"])] if row["error"] else [] if "missing_frames" in row: lost = frame_ranges(row["missing_frames"]) notes.insert(0, f"saved frames {lost} missing or unreadable") From 96223b551ccf1bd0fba03f3343a76c9e4df48fa3 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 14:17:52 -0500 Subject: [PATCH 24/25] Make the start warnings exact, and read each structure once The observation count rounded down, though MD also observes the end of a last, shorter interval, so 25 steps at one per 10 were told they could not report what they then reported. A .con structure was told it has constraints, because ASE's reader attaches a FixAtoms that fixes nothing. A cell-less structure at pressure was told MD can run, and it fails at once. And validate read a structure twice per trajectory, 6.8 s for 2000 trajectories of 288 atoms against 0.36 s. --- src/reactionflow/cli.py | 97 ++++++++++++++++++++++------------------- 1 file changed, 53 insertions(+), 44 deletions(-) diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 2f3ae4b..0e10e4c 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -10,7 +10,7 @@ import socket import sys from collections import Counter -from collections.abc import Iterator, Mapping, Sequence +from collections.abc import Iterable, Iterator, Mapping, Sequence from contextlib import suppress from dataclasses import asdict, replace from numbers import Integral @@ -18,6 +18,7 @@ from uuid import uuid4 from ase import Atoms +from ase.constraints import FixAtoms from ase.io import read from ase.io.formats import filetype @@ -140,47 +141,56 @@ def visible_gpu(environment: Mapping[str, str] | None = None) -> str: return devices[0] -def _start_warnings(campaign: CampaignConfig, index: int) -> list[str]: - """What a trajectory's settings will silently do to its result. None of this stops it.""" +def _start_warnings(campaign: CampaignConfig, indices: Iterable[int]) -> list[str]: + """What these trajectories' settings will silently do to their results. None of this stops them. - trajectory = campaign.trajectory(index) - path = campaign.structure_for(index) - name = os.path.relpath(path, campaign.source.parent) - atoms = read(path) - notes = [] - try: - several = len(read(path, index=slice(0, 2))) > 1 - except Exception: # a format that reads one image and takes no slice - several = False - if several: - notes.append(f"{name} holds several frames; a run starts from the last one") - if atoms.has("momenta"): - notes.append( - f"{name} carries momenta, which set the starting velocities in place of temperature_K; " - "the thermostat then moves them toward it" - ) - if atoms.constraints: - notes.append( - f"{name} has constraints, which hold in MD and are removed in refinement, so its " - "barriers are for a system with every atom free" - ) - if trajectory.pressure_GPa is not None and ( - not atoms.pbc.all() or atoms.cell.rank < 3 or atoms.cell.handedness <= 0 - ): - notes.append( - f"{name} has no fully periodic right-handed cell, which refinement at a pressure " - "needs, so MD can run but every event would end unresolved, with no barrier" - ) - run = campaign.reaction_run - observations = trajectory.total_steps // run.observation_interval - needed = run.detector.persistence_frames + run.candidate_stability_frames - 1 - if observations < needed: - notes.append( - f"{trajectory.total_steps} steps at one observation per {run.observation_interval} " - f"make {observations}, and confirming a reaction takes at least {needed}, so such a " - "run cannot report one" - ) - return notes + Each structure is read once, since a large campaign can name one for every trajectory. + """ + + notes: dict[str, None] = {} + structures: dict[Path, tuple[Atoms, bool]] = {} + for index in indices: + trajectory = campaign.trajectory(index) + path = campaign.structure_for(index) + if path not in structures: + try: + several = len(read(path, index=slice(0, 2))) > 1 + except Exception: # a format that reads one image and takes no slice + several = False + structures[path] = read(path), several + atoms, several = structures[path] + name = os.path.relpath(path, campaign.source.parent) + if several: + notes[f"{name} holds several frames; a run starts from the last one"] = None + if atoms.has("momenta"): + notes[ + f"{name} carries momenta, which set the starting velocities in place of " + "temperature_K; the thermostat then moves them toward it" + ] = None + # ASE's .con reader attaches a FixAtoms that fixes nothing. + if any(not isinstance(c, FixAtoms) or len(c.index) for c in atoms.constraints): + notes[ + f"{name} has constraints, which hold in MD and are removed in refinement, so its " + "barriers are for a system with every atom free" + ] = None + if trajectory.pressure_GPa is not None and ( + not atoms.pbc.all() or atoms.cell.rank < 3 or atoms.cell.handedness <= 0 + ): + notes[ + f"{name} has no fully periodic right-handed cell, which refinement at a pressure " + "needs, so every event would end unresolved, with no barrier" + ] = None + run = campaign.reaction_run + # MD also observes the end of a last, shorter interval. + observations = -(-trajectory.total_steps // run.observation_interval) + needed = run.detector.persistence_frames + run.candidate_stability_frames - 1 + if observations < needed: + notes[ + f"{trajectory.total_steps} steps at one observation per " + f"{run.observation_interval} make {observations}, and confirming a reaction takes " + f"at least {needed}, so such a run cannot report one" + ] = None + return list(notes) def run_selected_trajectory( @@ -205,7 +215,7 @@ def run_selected_trajectory( run = ReactionRun.open(root) atoms = read(campaign.structure_for(index)) if run.phase == "new" else None else: - for note in _start_warnings(campaign, index): + for note in _start_warnings(campaign, [index]): print(f"warning: {note}", file=sys.stderr) atoms = read(campaign.structure_for(index)) adapter = load_mlip_adapter(adapter_spec, trajectory) @@ -582,8 +592,7 @@ def main(argv: Sequence[str] | None = None) -> int: }, "adapter_profile_counts": dict(Counter(campaign.trajectory_profiles)), } - notes = (note for index in range(count) for note in _start_warnings(campaign, index)) - if warnings := list(dict.fromkeys(notes)): + if warnings := _start_warnings(campaign, range(count)): payload["warnings"] = warnings print(json.dumps(payload, sort_keys=True)) elif arguments.command == "check": From 964c4d12a2e94170e1fba35033ce913a47b1b9e6 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Fri, 6 Nov 2026 14:17:52 -0500 Subject: [PATCH 25/25] Read a compressed dump for the LAMMPS type check ASE sees through .gz, .bz2 and .xz, while the check read the compressed bytes as text, found no columns and told the user to add an element column the file already had. --- src/reactionflow/cli.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 0e10e4c..802d4c2 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -20,7 +20,7 @@ from ase import Atoms from ase.constraints import FixAtoms from ase.io import read -from ase.io.formats import filetype +from ase.io.formats import filetype, open_with_compression from ._durable import ensure_directory, file_digest, publish, sync_directory from ._version import __version__ @@ -348,7 +348,7 @@ def _file_warnings(path: Path, first: Atoms) -> list[str]: except Exception: # the frames were read, so this only costs the note below dump = False if dump: - with path.open(encoding="utf-8", errors="replace") as stream: + with open_with_compression(str(path)) as stream: columns = next( (line.split()[2:] for line in stream if line.startswith("ITEM: ATOMS")), [] )