/neb-irc-activation-energy

name: "neb-irc-activation-energy"
description: "NEB-IRC activation energy pipeline for reaction barriers using GFN2-xTB and pysisyphus. Optimize reactant and product geometries, run CI-NEB path search, optimize the transition state with a Hessian, verify with IRC (one imaginary mode, endpoints matching reactant/product, single NEB maximum), and report the electronic and Gibbs barriers. Use when you need a transition state, reaction barrier, activation energy, minimum energy path, or intrinsic reaction coordinate. Covers reactant/product atom-ordering pitfalls, feasibility sizing for single-core runs, and thermochemistry corrections. Renders an IRC energy-profile plot and an animated TS imaginary-mode HTML viewer. For 2D reaction scheme drawing use rdkit-chemdraw-cdxml."
license: "CC-BY-4.0"

NEB-IRC activation energy pipeline

Overview

Computes a reaction activation energy end to end — optimize reactant and product, find the minimum energy path with climbing-image NEB, refine the transition state with a Hessian, and confirm it with IRC — using GFN2-xTB through pysisyphus. Outputs a verified TS geometry, the barrier (ΔE‡; ΔG‡ after thermal corrections), an IRC energy profile (`irc_energy_profile.png`), and an animated TS imaginary-mode viewer (`ts_imaginary_mode.html`). Verification is required, not optional: a converged TS is meaningless until its single imaginary mode and IRC endpoints are checked.

When to Use

• Finding the transition state for an elementary reaction step and its activation energy
• Computing a reaction barrier (ΔE‡ or ΔG‡) to rank a series of related reactions
• Running a climbing-image NEB / minimum energy path between a reactant and product
• Verifying a candidate TS with IRC — does it connect the intended reactant and product?
• Screening barriers at a cheap semi-empirical level before committing DFT time

Reach for DFT on a multi-core node instead when you need quantitative agreement with experiment; GFN2-xTB barriers are semi-quantitative (see `references/energetics.md`). For a 2D scheme figure of the reaction, use the `rdkit-chemdraw-cdxml` skill instead.

Prerequisites

• **Tools**: `xtb` (GFN2-xTB engine), `pysisyphus` (`pysis` CLI for the path pipeline)
• **Python**: `matplotlib` for the IRC plot (the TS animation HTML needs no packages)
• **Input**: `reactant.xyz` and `product.xyz` with **identical atom ordering** (see Step 1)
• **Environment**: single core suffices; set `OMP_NUM_THREADS` to match physical cores

**Work in a local scratch dir** (e.g. `/tmp/rxn/`), not a mounted/networked workspace: pysisyphus creates and deletes symlinks and throws `PermissionError` mid-run on s3fs/FUSE. Copy results out at the end.

**Materialize the bundled scripts into the scratch dir first.** They can't be run in place from the skill directory, so use your file tools to read each one and save it into your working dir before running it. The scripts live in this skill's `scripts/` folder (next to this SKILL.md):

• `scripts/setup_env.sh`
• `scripts/pipeline.yaml`
• `scripts/check_result.py`
• `scripts/plot_irc.py`

The TS imaginary-mode animation is **not** produced here — read the **molecular-visualization-3dmol** skill and use its `mol_viewer.py` (Step 6).

Check for the tools; install only if missing (inside pixi/conda, invoke via `pixi run xtb`):

cd /tmp/rxn
command -v xtb && command -v pysis || bash setup_env.sh   # xtb binary + pysisyphus, ~2-3 min
source "${ROOT:-${HOME:-/tmp}/xtbenv}/env.sh"              # re-source in every new shell

Workflow

Step 1: Prepare reactant and product geometries

Most failures originate here, not in the NEB. Build the product by editing a **copy** of the reactant so atom ordering is identical — a permuted order gives a path that is geometrically valid and chemically meaningless. Align non-reacting groups so a spectator conformational change does not fold into the barrier. For bimolecular reactions use a pre-reaction complex as the reactant, not separated fragments (NEB converges poorly from infinite separation, and the reference state changes the reported barrier — record it).

# Build product from a copy of the reactant, moving only the reacting atoms.
from pathlib import Path

lines = Path("reactant.xyz").read_text().splitlines()
natoms = int(lines[0])
atoms = [ln.split() for ln in lines[2:2 + natoms]]   # [symbol, x, y, z] per atom
# ... edit ONLY the coordinates of atoms that move; keep order + symbols ...
Path("product.xyz").write_text("\n".join([str(natoms), "product"] + [" ".join(a) for a in atoms]) + "\n")

Step 2: Check the endpoints sit in different basins

If the reacting groups start too close, preoptimization carries the reactant over the barrier and both endpoints relax to the same structure; the NEB then returns a flat profile and the TS search aborts. This looks like success in the log until it fails minutes later, so compare the two pre-optimized endpoints on the key reacting bond.

import numpy as np

def load_xyz(fn):
    lines = open(fn).read().splitlines()
    return np.array([[float(v) for v in ln.split()[1:4]] for ln in lines[2:2 + int(lines[0])]])

r, p = load_xyz("first_pre_opt.xyz"), load_xyz("last_pre_opt.xyz")   # written by preopt
i, j = 0, 5                                     # indices of the atoms whose bond changes
dr, dp = np.linalg.norm(r[i] - r[j]), np.linalg.norm(p[i] - p[j])
assert abs(dr - dp) > 0.3, "endpoints nearly identical: move the reacting fragment further out"

Step 3: Run the pipeline

The template chains preopt → IDPP interpolation → CI-NEB → RS-I-RFO TS optimization with Hessian → IRC both directions → endpoint reoptimization. **Set `charge` and `mult` in `pipeline.yaml` before running** — the default `charge: 0` is wrong for any ion and converges silently to a meaningless TS. Add `alpb: <solvent>` for solution reactions. Whatever you set here must match every standalone `xtb` call in Step 5. Raise `max_cycles