"""Score `mepd channels` output against the g-xTB reference saddle. python -m benchmark.score_channels --runs RUNS --dataset SET.csv --inputs mepd_profile.toml --out scores.json The question: does `channels` find a *lower* saddle connecting the same two endpoints than the dataset's reference saddle? Both are measured on the same surface (g-xTB): the reference saddle is the dataset's TS re-converged at g-xTB and IRC-verified to connect these endpoints. Both are also measured from the same floor: the lowest reactant-side energy found anywhere in the run. Referencing each saddle to whatever conformer its own search started from would be gameable -- start from a strained conformer and the climb looks short -- and wrong, since reaching that conformer from the lowest one costs energy that precedes the barrier. See `lowest_reactant_energy`. For each run, every TS in `channels/` and `alternate-channels/` is a candidate: those are the ones whose IRC actually reconnects the requested pair. Off-target exits are excluded -- they go somewhere else, so a low barrier among them means nothing. A multistep alternate channel's barrier is its *highest* step, since that is what rate-limits the route. """ from __future__ import annotations import argparse import collections import csv import json import statistics from pathlib import Path from typing import Optional from benchmark import reaction_qm as rqm _REPO_ROOT = Path(__file__).resolve().parents[2] def _nodes_of_xyz(path: Path, charge: int, multiplicity: int) -> list: from mepd.chain import Chain from mepd.inputs import ChainInputs try: return list(Chain.from_xyz(path, ChainInputs(), charge=charge, spinmult=multiplicity)) except Exception: return [] def _energy(node, engine) -> Optional[float]: try: return float(engine.compute_energies([node])[0]) except Exception: return None def _energy_of_xyz(path: Path, engine, charge: int, multiplicity: int) -> Optional[float]: nodes = _nodes_of_xyz(path, charge, multiplicity) return _energy(nodes[0], engine) if nodes else None def lowest_reactant_energy(run_dirs, engine, charge: int, multiplicity: int, e_separated: Optional[float], ref_reactant_node) -> tuple: """Energy of the lowest reactant-side structure anywhere in this run. A barrier measured from whichever conformer a search happened to start from is gameable: begin from a strained conformer and the climb to the saddle looks small. It is not physically the barrier either, because getting from the lowest conformer to that strained one costs energy that must be paid first. So every barrier here -- the reference saddle's and every discovered saddle's -- is measured from one common floor: the lowest reactant-side energy the run turned up. Candidates are the minimized start conformers `channels` sampled, any IRC endpoint that relaxed back to the reactant, and the separately optimized fragments. Only structures whose connectivity matches the reactant count, so a product-side well can never become the floor. `run_dirs` is one run directory or several runs of the same reaction (e.g. one per conformer backend); the floor is then the lowest over all of them, so runs compared against each other share one zero. Returns (energy, provenance).""" import mepd.cli as cli best, source = None, None if e_separated is not None: best, source = e_separated, "separated_fragments" if isinstance(run_dirs, (str, Path)): run_dirs = [run_dirs] run_dirs = [Path(d) for d in run_dirs] tag = (lambda d, prov: f"{prov}@{d.parent.name}") if len(run_dirs) > 1 else (lambda d, prov: prov) candidates = [] for run_dir in run_dirs: conf = run_dir / "conformers" / "start.xyz" if conf.is_file(): candidates += [(n, tag(run_dir, "start_conformer")) for n in _nodes_of_xyz(conf, charge, multiplicity)] ts_dir = run_dir / "ts" if ts_dir.is_dir(): for irc_fp in sorted(ts_dir.glob("*_irc.xyz")): nodes = _nodes_of_xyz(irc_fp, charge, multiplicity) if len(nodes) >= 2: candidates += [(nodes[0], tag(run_dir, "irc_endpoint")), (nodes[-1], tag(run_dir, "irc_endpoint"))] for node, prov in candidates: try: if not cli._connectivity_matches(node, ref_reactant_node): continue except Exception: continue e = _energy(node, engine) if e is not None and (best is None or e < best): best, source = e, prov return best, source def score(runs: Path, dataset: Path, inputs: Path, floor_roots=()) -> list: """`floor_roots`: other run roots holding runs of the same reactions (matched by directory name). Their reactant-side structures join this run's in setting the floor -- see `lowest_reactant_energy`.""" from mepd.inputs import RunInputs engine = RunInputs.open(inputs).engine bench = {r["reaction_id"]: r for r in csv.DictReader(open(dataset))} out = [] for run_dir in sorted(Path(runs).glob("rxn_*")): status_fp = run_dir / "_driver_status.json" if not status_fp.is_file(): continue st = json.loads(status_fp.read_text()) rid = st.get("reaction_id") row = bench.get(rid) if row is None: continue charge, mult = int(row["charge"]), int(row["multiplicity"]) from benchmark.screen_gxtb import _node, _separated_energy if rid.startswith("RXN_"): import h5py with h5py.File(rqm.H5_PATH, "r") as fh: rxn = rqm.load_reaction(fh, rid, rqm.build_index()) e_sep = _separated_energy(engine, rxn.reactants, {}) ref_reactant_node = _node( *zip(*rqm._separate_fragments(rxn.reactants)), charge, mult ) else: # Not a Reaction-QM entry (e.g. the Geodesic MLP set): the dataset's # reactant geometry defines the reactant, and there is no # separated-fragments reference. e_sep = None ref_reactant_node = _nodes_of_xyz(_REPO_ROOT / "data" / "reaction_qm" / row["start_xyz"], charge, mult)[0] floor_dirs = [run_dir] + [ Path(r) / run_dir.name for r in floor_roots if (Path(r) / run_dir.name).is_dir() ] e_ref, e_ref_source = lowest_reactant_energy( floor_dirs, engine, charge, mult, e_sep, ref_reactant_node ) if e_ref is None: continue # The reference saddle is re-scored against the same floor rather than # taken from the dataset table, which is referenced to separated # fragments. Mixing the two would make any "improvement" partly an # artifact of the two barriers having different zeroes. ref_ts_fp = _REPO_ROOT / "data" / "reaction_qm" / row["ts_reference_xyz"] e_ref_ts = _energy_of_xyz(ref_ts_fp, engine, charge, mult) if e_ref_ts is None: continue ref_barrier = (e_ref_ts - e_ref) * rqm.HARTREE_TO_KCAL found = [] for sub, kind in (("channels", "channel"), ("alternate-channels", "alternate")): root = run_dir / sub if not root.is_dir(): continue for group in sorted(root.iterdir()): if not group.is_dir(): continue if kind == "channel": ts_files = [group / "ts.xyz"] else: # A route's barrier is its rate-limiting step. ts_files = sorted(group.glob("step_*/ts.xyz")) energies = [ _energy_of_xyz(fp, engine, charge, mult) for fp in ts_files if fp.is_file() ] energies = [e for e in energies if e is not None] if not energies: continue barrier = (max(energies) - e_ref) * rqm.HARTREE_TO_KCAL found.append({"kind": kind, "group": group.name, "barrier_kcal": round(barrier, 3)}) rec = { "reaction_id": rid, "n_atoms": st.get("n_atoms"), "n_reactant_frags": st.get("n_reactant_frags"), "ref_barrier_kcal": round(ref_barrier, 3), "ref_barrier_sep_kcal": float(row["barrier_kcal"]), "ref_b3lyp_barrier_kcal": float(row["ref_b3lyp_barrier_kcal"]), "reactant_floor_source": e_ref_source, "reactant_floor_hartree": e_ref, "returncode": st.get("returncode"), "wall_seconds": st.get("wall_seconds"), "floor_below_separated_kcal": ( round((e_sep - e_ref) * rqm.HARTREE_TO_KCAL, 3) if e_sep is not None else None ), "n_found": len(found), "found": found, } if found: best = min(found, key=lambda f: f["barrier_kcal"]) rec["best_barrier_kcal"] = best["barrier_kcal"] rec["best_kind"] = best["kind"] rec["improvement_kcal"] = round(ref_barrier - best["barrier_kcal"], 3) out.append(rec) return out def report(recs: list) -> None: n = len(recs) if not n: print("no scored runs yet") return withch = [r for r in recs if r.get("n_found")] print(f"{n} run(s) scored; {len(withch)} produced at least one endpoint-connecting channel") if not withch: return imp = [r["improvement_kcal"] for r in withch] lower = [r for r in withch if r["improvement_kcal"] > 1.0] same = [r for r in withch if abs(r["improvement_kcal"]) <= 1.0] higher = [r for r in withch if r["improvement_kcal"] < -1.0] print(f"\nvs the g-xTB reference saddle (>1 kcal/mol counts as a real difference):") print(f" found LOWER : {len(lower):>4} ({100*len(lower)/len(withch):.0f}%)") print(f" matched : {len(same):>4} ({100*len(same)/len(withch):.0f}%)") print(f" only higher : {len(higher):>4} ({100*len(higher)/len(withch):.0f}%)") if lower: li = sorted((r["improvement_kcal"] for r in lower), reverse=True) print(f"\n improvement among the {len(lower)} lower: median {statistics.median(li):.1f} " f"max {li[0]:.1f} kcal/mol") print(" largest improvements:") for r in sorted(lower, key=lambda r: -r["improvement_kcal"])[:8]: print(f" {r['reaction_id']} ref {r['ref_barrier_kcal']:>8.1f} -> " f"{r['best_barrier_kcal']:>8.1f} ({r['improvement_kcal']:+.1f}, {r['best_kind']})") src = collections.Counter(r.get("reactant_floor_source") for r in recs) print(f"\n reactant floor came from: {dict(src)}") drops = [r["floor_below_separated_kcal"] for r in recs if r.get("floor_below_separated_kcal")] if drops: print(f" floor sat {statistics.median(drops):.1f} kcal/mol below separated fragments " f"(median over {len(drops)}); max {max(drops):.1f}") kinds = [f["kind"] for r in withch for f in r["found"]] print(f"\n channel objects found: {kinds.count('channel')} direct, {kinds.count('alternate')} multistep") multi = [r for r in withch if r["n_found"] > 1] print(f" runs with >1 distinct connecting channel: {len(multi)}") def main(argv=None) -> int: p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) p.add_argument("--runs", type=Path, required=True) p.add_argument("--dataset", type=Path, required=True) p.add_argument("--inputs", type=Path, required=True) p.add_argument("--out", type=Path, help="write per-reaction scores as JSON") p.add_argument("--floor-roots", type=Path, nargs="*", default=(), help="other run roots of the same reactions; their reactant-side structures share the floor") args = p.parse_args(argv) recs = score(args.runs, args.dataset, args.inputs, floor_roots=args.floor_roots) report(recs) if args.out: args.out.write_text(json.dumps(recs, indent=2)) print(f"\nwrote {args.out}") return 0 if __name__ == "__main__": raise SystemExit(main())