neb-irc-activation-energy git:20260815.88171adA

Immutable. This exact content never changes and is served at /api/v1/blob/c55e2dd890c61dc4.

---
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 there first.** They can be read from the skill path but are
not on the sandbox's path, and are not on GitHub. Read each with the **leading-slash** skill
path (without the leading `/` the read hits the empty workdir) and write it locally:

```python
import os
os.makedirs("/tmp/rxn", exist_ok=True); os.chdir("/tmp/rxn")
_SKILL = "/SciAgent-Skills/skills/scientific-computing/neb-irc-activation-energy/scripts"
for name in ("setup_env.sh", "pipeline.yaml", "check_result.py", "make_visuals.py"):
    open(name, "w").write(read_file(f"{_SKILL}/{name}"))   # read_file = your file tool
```

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

```bash
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).

```python
# 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.

```python
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` to 150–200 above
~50 atoms.

```bash
pysis pipeline.yaml > pipeline.log 2>&1
tail -40 pipeline.log            # confirm it reached the endopt/IRC stage
```

ΔE‡ and the reaction energy come from the pipeline's `| BARRIERS |` block, referenced to the
reactant endpoint (`check_result.py` in Step 4 reports the same numbers):

```
  Left:     0.00 kJ mol-1      # reactant endpoint
    TS:   107.84 kJ mol-1      # dE‡ = TS - Left = 107.84 kJ/mol
 Right:    10.19 kJ mol-1      # dE_rxn = Right - Left = 10.19 kJ/mol
```

### Step 4: Verify the transition state

Never report a barrier before this passes. The checker applies all three gates (see Key
Concepts) and exits 0 only if every one passes.

```bash
python3 check_result.py pipeline.log
# [PASS] imaginary frequencies  exactly 1 at -621.8 cm-1
# [PASS] IRC endpoints          forward and backward matched distinct inputs
# [PASS] NEB profile            elementary, single barrier, span 107.5 kJ/mol
#   electronic barrier dE‡ (GFN2-xTB): 107.8 kJ/mol = 25.77 kcal/mol
```

### Step 5: Compute thermochemistry and report the barrier

Report ΔE‡ and the level of theory; for comparison against experimental rates, report ΔG‡,
which needs Hessians on the TS and reactant. **Every standalone `xtb` call must use the same
`--chrg`, `--uhf`, and solvent as `pipeline.yaml`, on the pipeline's optimized geometries — not
the raw input.** A gas-phase Hessian against a solvated barrier, or a missing `--chrg`, yields
a nonsensical (often negative) ΔG‡.

```bash
# flags must match pipeline.yaml (here: charge -1, aqueous)
xtb forward_end_final_geometry.xyz --hess --gfn 2 --alpb water --chrg -1 --uhf 0 > r_hess.log 2>&1
xtb ts_final_geometry.xyz          --hess --gfn 2 --alpb water --chrg -1 --uhf 0 > ts_hess.log 2>&1
# dG‡ = G(TS) - G(reactant), from the "TOTAL FREE ENERGY" lines
grep -i "TOTAL FREE ENERGY" r_hess.log ts_hess.log
```

### Step 6: Generate the two visual deliverables

Deliver both, alongside the numbers. `make_visuals.py` builds them in one call and reads
gfn/charge/mult/solvent from `pipeline.yaml`, so the IRC energies are recomputed at the
pipeline's exact level (override with `--charge/--mult/--alpb/--gfn` only if you edited the
pipeline after running).

```bash
python3 make_visuals.py            # -> irc_energy_profile.png + ts_imaginary_mode.html
```

Per-frame single points are recomputed because `forward_irc.trj` / `backward_irc.trj` carry no
energy in their comment lines. The animation reuses `ts_imaginary_mode_000.trj` (written by
`tsopt: do_hess: True`); open the HTML in a browser (it loads 3Dmol.js from a CDN) — it plays
the mode back and forth with a play/pause button and a speed slider.

## Key Parameters

Set in `pipeline.yaml` unless noted.

| Parameter | Section | Default | Range / Options | Effect |
|-----------|---------|---------|-----------------|--------|
| `gfn` | `calc` | `2` | `0`, `1`, `2` | GFN parametrization; `2` is the standard choice |
| `charge` / `mult` | `calc` | `0` / `1` | integers | Set explicitly; wrong values silently give a wrong TS |
| `alpb` | `calc` | off | solvent name (`water`, ...) | Implicit solvation; changes the barrier |
| `pal` | `calc` | `1` | physical cores | Threads; xTB parallel efficiency is modest |
| `between` | `interpol` | `8` | `6`–`12` | Intermediate NEB images (`between`+2 total) |
| `climb` | `cos` | `True` | `True`/`False` | Climbing image; gives a usable TS guess |
| `max_cycles` | `opt` | `80` | `80`–`200` | Raise above ~50 atoms |
| `do_hess` | `tsopt` | `True` | `True`/`False` | Required — the imaginary-frequency gate depends on it |

## Key Concepts

**The three verification gates** and what a failure means:

| Gate | Failure | Meaning and fix |
|---|---|---|
| Imaginary frequencies | 0 found | Optimizer fell into a minimum; perturb along the NEB tangent and rerun tsopt |
| Imaginary frequencies | ≥2 found | Higher-order saddle; displace along the lowest non-reactive mode and reoptimize |
| Imaginary frequencies | 1, wrong mode | Confirm it is the bonds breaking/forming — a methyl rotation also shows exactly one |
| IRC endpoints | mismatch | TS connects other species; may be a real TS for a different reaction — do not report |
| NEB profile | multiple maxima | Not an elementary step; split at the intervening minimum and run each segment |

**Barrier definitions** differ by reference and corrections: ΔE‡ (electronic), ΔE‡+ZPE, ΔH‡,
ΔG‡ (for kinetics). For bimolecular reactions the reference state (separated reactants vs
pre-reaction complex) shifts the number — record which was used. Full table in
`references/energetics.md`.

**Submerged barriers.** For an ion + neutral in the gas phase (e.g. anionic SN2) the TS can
sit *below* the separated reactants — a negative barrier against separated reactants is real,
not a bug. Reference to the pre-reaction complex (the pipeline's endpoint) and/or add `alpb`
solvation, which lifts it to a positive, experiment-comparable value.

## Common Recipes

### Recipe: Split a multi-step reaction into elementary steps

When to use: `check_result.py` reports multiple NEB maxima — it names the intermediate image
index(es). Extract that image from the NEB path and run the pipeline on each half.

```python
from pathlib import Path

N = 5   # interior-minimum image index reported by check_result.py's NEB profile gate
lines = Path("final_geometries.trj").read_text().splitlines()   # one frame per NEB image
n = int(lines[0].split()[0])
Path("intermediate.xyz").write_text("\n".join(lines[N * (n + 2):(N + 1) * (n + 2)]) + "\n")
# then: pipeline on reactant.xyz + intermediate.xyz, and intermediate.xyz + product.xyz
```

### Recipe: Refine the electronic energy with a DFT single point

When to use: the xTB geometry is fine but you need a better barrier. Keep the xTB geometry and
thermal corrections; replace only the electronic energy on 3 structures (≤30 atoms).

```python
from pyscf import gto, dft

def electronic_energy(xyz_fn, basis="def2-svp", xc="b3lyp"):
    mol = gto.M(atom=xyz_fn, basis=basis)     # xyz file path accepted directly
    mf = dft.RKS(mol); mf.xc = xc
    return mf.kernel()                          # Hartree

dE = (electronic_energy("ts_final_geometry.xyz") - electronic_energy("reactant.xyz")) * 2625.4996
print(f"DFT dE‡ = {dE:.1f} kJ/mol (add xTB G_corr for dG‡)")
```

## Expected Outputs

- `pipeline.log` — full run log; parsed by `check_result.py`
- `ts_final_geometry.xyz` — the optimized transition state
- `final_geometries.trj` — the converged NEB path (per-image energies in comment lines)
- `ts_imaginary_mode_000.trj` — TS displaced along the imaginary mode (input to the animation)
- `irc_energy_profile.png`, `ts_imaginary_mode.html` — the two deliverables (Step 6)
- A barrier: ΔE‡ from the BARRIERS block, ΔG‡ after Hessian thermal corrections

## Troubleshooting

| Problem | Cause | Solution |
|---------|-------|----------|
| `xtb: command not found` | Env not sourced | `source "${ROOT:-${HOME:-/tmp}/xtbenv}/env.sh"` in every new shell |
| `scripts/…: No such file` / GitHub 404 | Bundled scripts not on sandbox path | Materialize via the leading-slash skill path + `read_file` (Prerequisites); not on GitHub |
| `setup_env.sh: HOME: unbound variable` | `HOME` unset under `set -u` | Fixed in the shipped script; if patching, `export HOME="${HOME:-/tmp}"` first |
| `PermissionError` on a symlink mid-run | pysisyphus symlinks on a mounted/s3fs dir | Run in a local dir (`/tmp/rxn/`), copy results back |
| ΔG‡ negative or absurd | Hessian charge/solvent ≠ pipeline, or raw input geometry used | Match `--chrg`/`--uhf`/`alpb` to `pipeline.yaml`; use the optimized endpoint geometry (Step 5) |
| Barrier below separated reactants | Submerged barrier for ion + neutral | Expected in gas phase; reference the pre-reaction complex and/or add `alpb` |
| NEB profile nearly flat, TS aborts | Endpoints in the same basin | Move the reacting fragment further out (Step 2) |
| Path is chemically nonsensical | Permuted atom order between endpoints | Rebuild product from a copy of the reactant (Step 1) |
| Job killed / never finishes | System too large for the budget | Check `references/feasibility.md`; split stages or shrink the model |
| `qm` optimizer crashes on cycle 1 | QuickMin instability | Use `type: lbfgs` under `opt` (the template default) |
| `No imaginary-mode trajectory found` | `ts_imaginary_mode_000.trj` absent | Ensure `tsopt: do_hess: True` ran; re-run tsopt |
| Animation HTML blank | 3Dmol.js blocked (offline / strict CSP) | Open with network access; the viewer loads 3Dmol from a CDN |

## Bundled Resources

- `scripts/setup_env.sh` — installs xtb (GitHub release) + pysisyphus (PyPI), writes `env.sh`
- `scripts/pipeline.yaml` — full preopt→NEB→TSopt→IRC→endopt template with inline comments
- `scripts/check_result.py` — verification of the three gates (exit 0 = all pass); prints ΔE‡
- `scripts/make_visuals.py` — builds `irc_energy_profile.png` + `ts_imaginary_mode.html`
- `references/feasibility.md` — measured timings, atom-count sizing, what DFT can/can't do here
- `references/energetics.md` — ΔE‡/ΔH‡/ΔG‡ definitions, thermochemistry, reporting conventions

## References

- [xtb documentation](https://xtb-docs.readthedocs.io/) — GFN2-xTB methods and CLI
- [pysisyphus documentation](https://pysisyphus.readthedocs.io/) — COS/NEB, TS optimizers, IRC
- Bannwarth, Ehlert, Grimme, *J. Chem. Theory Comput.* 2019, 15, 1652 — GFN2-xTB method paper
- Steinmetzer, Kupfer, Gräfe, *Int. J. Quantum Chem.* 2021, 121, e26550 — pysisyphus paper
- [ASE](https://wiki.fysik.dtu.dk/ase/) / [PySCF](https://pyscf.org/) — geometry I/O and DFT single points