diff --git a/README.md b/README.md index 6ac6a0a..21430c1 100644 --- a/README.md +++ b/README.md @@ -66,7 +66,9 @@ wait in a queue until then ([when barriers are refined](docs/campaigns.md#when-b Three 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. `reactionflow export-electrace --out ` lays out a band as VASP frame folders and an [ElecTrace](https://github.com/Austin243/ElecTrace) configuration, to follow charges, critical points and bonds along it @@ -75,7 +77,8 @@ points and bonds along it 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 4687201..e06f96f 100644 --- a/docs/campaigns.md +++ b/docs/campaigns.md @@ -236,10 +236,18 @@ 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, +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 @@ -249,22 +257,30 @@ 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. +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 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) 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/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/eon-search.md b/docs/eon-search.md index 889d12b..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 @@ -245,7 +246,13 @@ 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. + 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/docs/exact-restart.md b/docs/exact-restart.md index 1346c80..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, @@ -87,7 +92,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/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/docs/offline-scan.md b/docs/offline-scan.md index c4bdcb6..dfb0880 100644 --- a/docs/offline-scan.md +++ b/docs/offline-scan.md @@ -13,16 +13,19 @@ 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` `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 @@ -37,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 @@ -51,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: @@ -62,7 +69,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/docs/pathway-refinement.md b/docs/pathway-refinement.md index d7f7d36..985ce03 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,16 @@ 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. + +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 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` 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) 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 f69a43a..b0a6cad 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -8,15 +8,19 @@ import math import os import socket +import sys from collections import Counter -from collections.abc import Mapping, Sequence +from collections.abc import Iterable, 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.constraints import FixAtoms from ase.io import read +from ase.io.formats import filetype, open_with_compression from ._durable import ensure_directory, file_digest, publish, sync_directory from ._version import __version__ @@ -25,8 +29,8 @@ from .electrace_export import export_band 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 .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 from .wizard import create_campaign @@ -138,6 +142,58 @@ def visible_gpu(environment: Mapping[str, str] | None = None) -> str: return devices[0] +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. + + 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( campaign: CampaignConfig, *, @@ -160,6 +216,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) @@ -198,13 +256,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: @@ -223,8 +299,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: @@ -255,6 +331,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 open_with_compression(str(path)) 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)) @@ -373,6 +489,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)) @@ -398,9 +517,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): @@ -438,21 +563,39 @@ 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": + count = len(campaign.trajectories) payload = { "campaign": str(campaign.source), "trajectories": len(campaign.trajectories), @@ -463,6 +606,8 @@ def main(argv: Sequence[str] | None = None) -> int: }, "adapter_profile_counts": dict(Counter(campaign.trajectory_profiles)), } + if warnings := _start_warnings(campaign, range(count)): + payload["warnings"] = warnings print(json.dumps(payload, sort_keys=True)) elif arguments.command == "check": if arguments.index is None and len(campaign.trajectories) > 1: @@ -475,7 +620,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 cd01eca..151d525 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]): @@ -653,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, *, @@ -700,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") @@ -730,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 diff --git a/src/reactionflow/run.py b/src/reactionflow/run.py index 1a36d02..391dbc6 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 { @@ -108,6 +111,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) @@ -391,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 @@ -402,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()) @@ -409,11 +460,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 +528,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)) @@ -580,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``.""" @@ -640,7 +718,21 @@ 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: + 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. diff --git a/src/reactionflow/scan.py b/src/reactionflow/scan.py index 3710b6d..9872bd0 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: @@ -172,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() @@ -181,13 +189,15 @@ 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 + 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: @@ -205,12 +215,16 @@ 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)), ""), } ) 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 +277,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 +303,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) 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/eon.py b/src/reactionflow/search/eon.py index ae9d84f..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 @@ -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") @@ -347,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 @@ -365,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} @@ -410,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 = { @@ -469,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() @@ -513,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) @@ -653,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 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. diff --git a/src/reactionflow/status.py b/src/reactionflow/status.py index a144f27..520e18f 100644 --- a/src/reactionflow/status.py +++ b/src/reactionflow/status.py @@ -14,9 +14,10 @@ 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 from .store import _read_bundle # Converged bands, including those refined before the saddle check, which have no verdict. @@ -63,11 +64,14 @@ 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. + 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. + 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,121 @@ 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.get("global_frame", 0)) + if lost: + 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) classes = [] steps = [] pathways: Counter[str] = Counter() + + read: dict[str, ReactionCandidate | None] = {} + lost_bundles: set[str] = set() + + 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 + lost_bundles.add(name) + 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: + lost_bundles.add(occurrence["bundle"]) + return None + + 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"]) - results = [] + # 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, refined = [], [] 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) + 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"]: + step_path = path.parent / f"step-{step['step']}" / "result.json" + step_result = result_file(step_path) + if step_result is None: + continue + try: bonds = { 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"])) + 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(): + 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 @@ -120,6 +210,8 @@ 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 + no_band: Counter[tuple[str, str]] = Counter() for index, spec in enumerate(campaign.trajectories): model = campaign.adapter_profile_for(index) row: dict[str, Any] = { @@ -142,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(): @@ -149,14 +246,28 @@ 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 = [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") + 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) + 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)) - merged.add(name, candidate) + if candidate is not None: + merged.add(name, candidate) + labelled[name] = _reaction_graph(candidate, radius=radius) entries[name] = { - "reaction": _describe(candidate, radius), + "reaction": ( + f"class of {occurrence_id} (no readable candidate, direction unknown)" + if candidate is None + else _describe(candidate, radius) + ), "example": {"trajectory": spec.id, "occurrence_id": occurrence_id}, "model": model, "pressure_GPa": spec.pressure_GPa, @@ -174,8 +285,10 @@ 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 + 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 @@ -183,7 +296,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]) @@ -197,12 +318,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: @@ -217,6 +344,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") @@ -251,12 +385,34 @@ 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}.", ""] + 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 += [ + 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) + 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"], @@ -266,12 +422,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( @@ -285,11 +437,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 = ( @@ -311,6 +475,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'}" 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),