From 56852a31c594d1218250a9275770f0f8e540cf48 Mon Sep 17 00:00:00 2001 From: Austin Ellis Date: Sat, 7 Nov 2026 12:00:27 -0500 Subject: [PATCH] Export a refined band as VASP frame folders for ElecTrace reactionflow export-electrace --out writes one numbered folder per image of a band, each with only a POSCAR, plus an electrace.toml and a reactionflow.json. The configuration reads the folders as one NEB path with index atom identity, sets pressure_gpa for an NPT band so that ElecTrace compares the same enthalpy as barrier_eV, and lists the atoms of the bonds that differ between the band's ends as focus atoms. ReactionFlow keeps bonds by stable atom ID, so those are mapped through the written atom order; reactionflow.json records that order with the pathway's result.json. Every image is sorted by element in the same stable order, so each species is one POSCAR block and an atom keeps its index. A band whose images do not share their atoms and order, or that has no cell, is refused, as is an output folder that is not empty. Constraints and the model's energies are not written, so ElecTrace reads only the energies of the calculations run on the frames. A band that did not converge or failed its saddle check is exported with its status printed. An event refined in steps points to its step-N folders, which export with their own bonds. --- README.md | 6 +- docs/pathway-refinement.md | 25 +++++ src/reactionflow/cli.py | 14 +++ src/reactionflow/electrace_export.py | 140 +++++++++++++++++++++++++++ 4 files changed, 184 insertions(+), 1 deletion(-) create mode 100644 src/reactionflow/electrace_export.py diff --git a/README.md b/README.md index 2019df1..6ac6a0a 100644 --- a/README.md +++ b/README.md @@ -63,10 +63,14 @@ By default a trajectory pauses for each reaction and refines it before MD goes o `"refine_after_md": true` in the campaign, MD runs to its last step without pausing, and the events wait in a queue until then ([when barriers are refined](docs/campaigns.md#when-barriers-are-refined)). -Two more commands work on finished output. `reactionflow scan` replays the saved frames through +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. +`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 +([electronic analysis](docs/pathway-refinement.md#electronic-analysis-with-electrace)). 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, diff --git a/docs/pathway-refinement.md b/docs/pathway-refinement.md index aad5c7d..d7f7d36 100644 --- a/docs/pathway-refinement.md +++ b/docs/pathway-refinement.md @@ -208,6 +208,31 @@ started before the check cannot resume under it, because its pathway settings, a contract, have changed. Finish it with the version that started it, and test its bands with `reactionflow check`. Its results keep their `ci_neb_converged` status. +## Electronic analysis with ElecTrace + +```bash +reactionflow export-electrace runs//pathways/ --out band-dft +``` + +`reactionflow export-electrace` lays out one band for +[ElecTrace](https://github.com/Austin243/ElecTrace), which follows charges, critical points and +bonds from image to image. `band-dft` must be new or empty. It gets one numbered folder per image +(`00`, `01`, ...) holding only a `POSCAR`, plus two files: + +- `electrace.toml` analyzes the folders as one NEB path with `atom_identity = "index"`. An NPT band + sets `pressure_gpa`, so that ElecTrace compares H = E + PV at the band's pressure, the same + enthalpy as `barrier_eV`. `focus_atoms` lists the atoms of the bonds that differ between the + band's two ends. +- `reactionflow.json` keeps the pathway's `result.json` and the changed bonds, plus `atom_ids`, + the stable ID of the atom at each POSCAR index. + +Each species is one POSCAR block, and an atom has the same index in every frame. Constraints are +not written, and neither are the model's energies, which stay in `reactionflow.json`: run the +electronic-structure calculations and analyses in the frame folders, then `electrace analyze +band-dft/electrace.toml`. A band whose saddle check failed, or that did not converge, is still +exported, and its status is printed. An event refined in steps has no band of its own; export one +of its `step-N` folders, whose focus atoms come from that step's bonds. + ## Constant-pressure SSNEB Campaigns pass `pressure_GPa` automatically. Direct calls to refinement or run methods must supply diff --git a/src/reactionflow/cli.py b/src/reactionflow/cli.py index 35f47e7..f69a43a 100644 --- a/src/reactionflow/cli.py +++ b/src/reactionflow/cli.py @@ -22,6 +22,7 @@ from ._version import __version__ from .campaign import CampaignConfig from .detection import _is_count +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 @@ -323,6 +324,16 @@ def _parser() -> argparse.ArgumentParser: scan.add_argument("--frame-ps", type=float, help="ps between the frames of a trajectory file") scan.add_argument("--out", type=Path, help="output directory") + export = commands.add_parser( + "export-electrace", help="lay out one refined band as VASP frame folders for ElecTrace" + ) + export.add_argument( + "pathway", + type=Path, + help="a band's folder: /pathways/, or one of its step-N folders", + ) + export.add_argument("--out", type=Path, required=True, help="a new or empty folder") + status = commands.add_parser( "status", help="summarize trajectories and reactions without changing any run" ) @@ -377,6 +388,9 @@ def main(argv: Sequence[str] | None = None) -> int: parser.error(str(error)) print(json.dumps({"backends": catalog}, indent=2, sort_keys=True)) return 0 + if arguments.command == "export-electrace": + print(json.dumps(export_band(arguments.pathway, arguments.out), sort_keys=True)) + return 0 if arguments.command == "scan" and arguments.campaign.suffix != ".json": if arguments.index is not None: parser.error("--index picks a campaign trajectory") diff --git a/src/reactionflow/electrace_export.py b/src/reactionflow/electrace_export.py new file mode 100644 index 0000000..e00657a --- /dev/null +++ b/src/reactionflow/electrace_export.py @@ -0,0 +1,140 @@ +"""Lay out one refined band as numbered VASP frame folders for ElecTrace.""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any + +import numpy as np +from ase import Atoms +from ase.io import read, write + +from .detection import atom_ids + + +def _bonds(values: Any) -> frozenset[tuple[int, int]]: + return frozenset((min(map(int, pair)), max(map(int, pair))) for pair in values) + + +def _changed_bonds(pathway: Path, result: dict[str, Any]) -> frozenset[tuple[int, int]]: + """Bonds, as atom-ID pairs, that differ between the band's two ends.""" + + if "step" in result: + # A step of a stepwise event records the bonds of its own two ends. + return _bonds(result["reactant_bonds"]) ^ _bonds(result["product_bonds"]) + if pathway.parent.name != "pathways": + raise ValueError(f"{pathway} is not a folder under a trajectory's pathways/") + bundle = pathway.parent.parent / "candidates" / result["occurrence_id"] / "candidate.json" + candidate = json.loads(bundle.read_text(encoding="utf-8"))["candidate"] + return _bonds(candidate["reactant_bonds"]) ^ _bonds(candidate["product_bonds"]) + + +def _toml_string(value: str) -> str: + return json.dumps(value) # a JSON string with ASCII escapes is also a TOML basic string + + +def export_band(pathway: str | Path, out: str | Path) -> dict[str, Any]: + """Write the band in ``pathway`` to ``out`` as ``00``, ``01``, ... folders of POSCARs. + + ``out`` also gets an ``electrace.toml`` for that path, with index atom identity, the band's + pressure, and the atoms of its changed bonds as focus atoms, and a ``reactionflow.json`` + with the refinement record and the atom order. Energies stay in ``reactionflow.json``: they + are the model's, and ElecTrace should read the energies of the calculations run on the frames. + """ + + pathway, out = Path(pathway), Path(out) + result = json.loads((pathway / "result.json").read_text(encoding="utf-8")) + if result.get("schema_version") != 1: + raise ValueError(f"unsupported pathway result in {pathway}") + images: list[Atoms] = read(pathway / "images.traj", index=":") # type: ignore[assignment] + if len(images) < 2: + steps = sorted(path.name for path in pathway.glob("step-*") if path.is_dir()) + if steps: + raise ValueError( + f"{pathway} was refined in steps and has no band of its own; export one of " + + ", ".join(str(pathway / step) for step in steps) + ) + raise ValueError(f"{pathway} holds no band (status {result.get('status')!r})") + + ids = atom_ids(images[0]) + symbols = images[0].get_chemical_symbols() + for number, image in enumerate(images): + if atom_ids(image) != ids or image.get_chemical_symbols() != symbols: + raise ValueError(f"image {number} of {pathway} does not keep the first image's atoms") + if abs(np.linalg.det(image.cell.array)) < 1e-6: + raise ValueError(f"image {number} of {pathway} has no cell, which a POSCAR needs") + + if out.exists() and (not out.is_dir() or any(out.iterdir())): + raise FileExistsError(f"{out} already exists and is not an empty folder") + out.mkdir(parents=True, exist_ok=True) + + # One stable element sort for every image, so each species is one POSCAR block and an atom + # has the same index in every frame. + order = sorted(range(len(symbols)), key=symbols.__getitem__) + index_of = {ids[original]: new for new, original in enumerate(order)} + width = max(2, len(str(len(images) - 1))) + labels = [f"{number:0{width}d}" for number in range(len(images))] + for label, image in zip(labels, images, strict=True): + frame = image[order] + frame.set_constraint() + frame.calc = None + (out / label).mkdir() + write(out / label / "POSCAR", frame, format="vasp", direct=True, vasp5=True) + + changed = sorted(_changed_bonds(pathway, result)) + focus = sorted({index_of[atom_id] for bond in changed for atom_id in bond}) + pressure = result.get("pressure_GPa") + name = f"ReactionFlow {result['occurrence_id']}" + if "step" in result: + name += f" step {result['step']}" + config = [ + "# Written by reactionflow export-electrace. Frames hold only POSCARs: run the", + "# calculations and descriptor analyses in them, then run electrace analyze.", + "[project]", + f"name = {_toml_string(name)}", + 'series_kind = "neb"', + *(() if pressure is None else (f"pressure_gpa = {float(pressure)!r}",)), + "", + "[inputs]", + 'root = "."', + "", + "[matching]", + 'atom_identity = "index"', + "", + "[outputs]", + "# Atoms of the bonds that differ between the band's two ends, as POSCAR indices.", + f"focus_atoms = {focus}", + "", + ] + (out / "electrace.toml").write_text("\n".join(config), encoding="utf-8") + + record = { + "schema_version": 1, + "source": str(pathway.resolve()), + "frames": labels, + "periodic": [bool(value) for value in images[0].pbc], + "atom_ids": [ids[original] for original in order], + "changed_bonds": { + "atom_ids": [list(bond) for bond in changed], + "indices": [[index_of[first], index_of[second]] for first, second in changed], + }, + "result": result, + } + (out / "reactionflow.json").write_text( + json.dumps(record, indent=2, sort_keys=True) + "\n", encoding="utf-8" + ) + validation = result.get("frequency_validation") or {} + return { + "export": str(out), + "frames": len(labels), + "focus_atoms": focus, + "status": result.get("status"), + "saddle_check": validation.get("status"), + "barrier_eV": result.get("barrier_eV"), + "barrier_quantity": result.get("barrier_quantity"), + "pressure_GPa": pressure, + } + + +__all__ = ["export_band"]