--- name: orca-irc description: Generate ORCA input files for IRC (Intrinsic Reaction Coordinate) calculations and post-process the results. Use this skill whenever the user asks about IRC calculations, reaction path following, confirming transition state connectivity, or tracing a minimum energy path from a TS in ORCA. Also trigger when the user mentions IRC endpoints, forward/backward reaction paths, or needs to verify that a TS connects to expected reactants and products. --- # ORCA IRC Calculation Skill This skill handles the full IRC workflow: generating an ORCA input file from a confirmed transition state geometry, and post-processing the output to extract the reaction-path energy profile and endpoint geometries. **Scope:** Input generation, local post-processing, and (optionally) HPC submission via the CatGo workflow engine. The "Submitting to HPC" section below covers the proven Expanse flow. If the user is running on their own non-CatGo infrastructure, just generate the input file from the template in Stage 1 and skip the submission section. **Target version:** ORCA 6.x. The parser and defaults below are written against the ORCA 6 output layout used by CatGo's own parser (`server/catgo/utils/orca_output.py::OrcaIrcOutput`). ## Prerequisites IRC requires a transition state geometry that has been confirmed by a frequency calculation showing exactly one imaginary frequency corresponding to the reaction coordinate. If the user has not run a frequency calculation on the candidate TS, suggest the `orca-freq` skill before generating an IRC input. Running IRC on a geometry that isn't a true first-order saddle point produces meaningless paths. This is a suggestion, not a hard stop — if the user confirms they know what they're doing and still wants the input, generate it. ## Stage 1: Input Generation Before generating the input, gather from the user: - **TS geometry** — file path or pasted XYZ coordinates - **Method and basis set** — must match whatever was used for the TS optimization and frequency calculation. Method/basis consistency across Opt → Freq → IRC is critical; mixing levels means the TS may not be a stationary point on the new surface. - **Charge and multiplicity** - **Number of cores** (`%pal nprocs`) - **Memory per core** (`%maxcore`, in MB) Optional (apply only if user specifies): - Solvent model (CPCM). Gas-phase IRC is standard practice for confirming connectivity; solvation mainly affects barrier heights, not topology. - Dispersion correction — again, match whatever was used for the TS search. ### Hessian Handling (important — read all of this) IRC moves downhill from the TS along the imaginary-frequency mode. To identify that mode, ORCA needs a Hessian **at the TS geometry**. Getting this step right is the single biggest factor in whether the IRC gives a sensible path. **What the Hessian is used for, concretely:** 1. Diagonalize to find the mode with a negative eigenvalue (the imaginary frequency). That eigenvector is the initial descent direction. 2. Build a mass-weighted quadratic model around the TS. The very first displacement along the mode is scaled by `Scale_Displ_SD`, chosen (via `Init_Displ_DE`) to correspond to a target drop in energy. 3. Subsequent steps use an adaptive integrator in mass-weighted coordinates; the Hessian is not re-diagonalized every step (that's what `SD` vs. `LQA` flavors are about), but the first-step quality dominates the final path. **Consequences:** - A Hessian with more than one negative eigenvalue means your geometry is not a clean TS. IRC will pick *one* of the negative modes and the result is arbitrary. Fix the TS first. - A numerically noisy Hessian (bad SCF convergence, too-loose grid) can flip the sign of a small imaginary mode or rotate the eigenvector into a nearby low-frequency mode. The IRC then follows the wrong coordinate — often tumbling, rotation, or a methyl rock instead of the bond-making/breaking mode. `calc_numfreq` or `InitHess read` from a converged `AnFreq` job is the fix. - Reusing a `.hess` file from a calculation at a *different* method/basis/grid/solvation is a silent footgun. The eigenmode numbering can differ, so "the imaginary mode" may not be the same mode on the IRC surface. **Only use `InitHess read` if the `.hess` was produced at the same level as this IRC job.** **The three options:** **Option A — `InitHess read` (preferred when applicable):** Read a pre-computed `.hess` file from the frequency calculation at the same level of theory. Fastest and most accurate, because the Hessian is already well-converged and you skip recomputing it. Ask the user for the filename. ``` %irc InitHess read Hess_Filename "previous_freq.hess" end ``` **Option B — `InitHess calc_anfreq`:** Compute analytical frequencies at the start of the IRC. Use when no `.hess` file is available and the method supports analytical Hessians (common DFT functionals, HF). Cleanest and still fairly fast for small-to-medium systems. ``` %irc InitHess calc_anfreq end ``` **Option C — `InitHess calc_numfreq`:** Compute numerical frequencies. Use when analytical Hessians aren't available or aren't reliable for the chosen method (some double hybrids, multireference, custom-parameter methods). Much slower (scales with 6N SCFs) but works everywhere. ``` %irc InitHess calc_numfreq end ``` **Decision flow for the skill:** ask "Do you have a `.hess` file from your frequency calculation at this level of theory?" If yes → Option A with filename. If no and the method supports analytical Hessians → Option B. If no and it doesn't → Option C. ### IRC Parameters and ORCA Defaults ORCA 6's documented defaults: - `MaxIter` = **20** iterations per arm - `Step` = **adaptive**, governed by `Scale_Displ_SD` (initial value derived from `Init_Displ_DE`) and dynamically rescaled — bounded to 1/16× – 4× of the initial value. There is no single fixed step-size default to quote. The default 20 iterations is often too few to reach a clear minimum, so the recommended override is: ``` %irc MaxIter 70 PrintLevel 1 Direction both end ``` Do **not** set `Step` unless the user specifically asks for a fixed step — the adaptive algorithm is usually better than any constant. If the path behaves badly, the right lever is usually `MaxIter`, `InitHess`, or the TS quality, not a fixed step size. Rationale for the overrides: - `MaxIter 70` — ORCA's default of 20 is almost always too few for a path that reaches a clear minimum. Bump to 100–150 for flat PES regions. - `PrintLevel 1` — prints energy and gradient at each step so the run can be monitored and the path summary is written in full. - `Direction both` — follows the path in both directions from the TS. Almost always what you want when confirming connectivity. ### Composite Methods For composite methods like `PBEh-3c` or `r²SCAN-3c`, do NOT specify a separate basis set — the basis is built into the method keyword. ``` ! PBEh-3c IRC ``` This automatically satisfies method/basis consistency since the composite method has a fixed internal basis. ### Template ``` # ORCA IRC Calculation # Method: {method} {basis} # TS confirmed by frequency calculation: {num_imaginary} imaginary mode(s) %pal nprocs {nprocs} end %maxcore {maxcore} ! {method} {basis} {dispersion} IRC %irc InitHess {inithess_option} {hess_filename_line} MaxIter 70 PrintLevel 1 Direction both end %output jsongbwfile True jsonpropfile True end * xyzfile {charge} {multiplicity} {ts_geometry_file} ``` When generating the input: - If the user provides inline coordinates, use `* xyz {charge} {mult}` followed by the coordinates and closing `*`. - Omit `{dispersion}` if already included in a composite method. - Omit `{basis}` for composite methods. - Omit the `Hess_Filename` line unless using `InitHess read`. - Keep the `%output jsongbwfile True jsonpropfile True end` line — it makes ORCA emit the JSON files OPI's `Output.parse()` consumes during post-processing (Stage 2). ### About OPI input builders OPI (`pip install orca-pi`) ships typed `BlockIrc` builders that validate keys at construction. **However, the catgo backend's `irc` node already emits its own `%irc` block from node params (`irc_max_iter`, `irc_direction`, `init_hess`, etc.).** Pasting an OPI-built `%irc` block via `extra_blocks` would produce **two `%irc` blocks** in the same `.inp`, which is undefined behavior. For this skill, **stick with node params** for `%irc` content and use `extra_blocks` only for `%output`. The OPI parsing wins (Stage 2) still apply. If you need a knob `BlockIrc` exposes that the node params don't (`scale_init_displ`, `tolmaxg`, `monitor_internals`, etc.), open that as a node-def gap rather than dual-emitting blocks. ## Submitting to HPC (Expanse) — proven flow Use this when the user wants the CatGo workflow engine to run the IRC on Expanse. Skip if they only want the input file. > **Use `catgo_workflow` (graph-based), NOT `catgo_workflow_engine` (task-based).** > The graph-based tool auto-captures the viewer structure on `create` and lets > you wire opt → freq → IRC explicitly. Param keys differ from the task-based > API: graph-based uses `method`/`basis`, task-based uses `orca_method`/`orca_basis`. ### 1. Confirm TS structure is loaded and find session_id ```json catgo_view(action: "get_state") ``` ```bash curl -s http://localhost:8000/api/hpc/connections ``` Copy the `session_id` for `host: login.expanse.sdsc.edu`. ### 2. Create the workflow (auto-captures the TS structure) ```json catgo_workflow(action: "create", name: "IRC from TS") ``` ### 3. Add the IRC node ```json catgo_workflow(action: "batch", workflow_id: "", operations: [ {"op": "add_node", "node_type": "irc", "label": "irc1", "params": { "software": "orca", "method": "B3LYP", "basis": "def2-SVP", "charge": 0, "multiplicity": 1, "dispersion": "D3BJ", "max_iterations": 70 }}, {"op": "connect", "from_id": "", "to_id": "irc1", "from_handle": "structure", "to_handle": "structure"} ]) ``` ### Canonical IRC node params (what the engine actually reads) | Parameter | Default | Description | |---|---|---| | `method` | r2SCAN-3c | DFT functional | | `basis` | 6-31G | Basis set (omit for composite methods) | | `charge` / `multiplicity` | 0 / 1 | | | `dispersion` | (none) | `D4` \| `D3BJ` \| `D3` — **use this field, NOT `extra_keywords`** | | `grid` | DefGrid2 | `DefGrid1/2/3` | | `max_iterations` | 30 | per-arm iteration cap (legacy alias: `max_irc_iterations`) | | `num_cores` / `max_core_mb` | 4 / 4000 | `%pal nprocs` / `%maxcore` | > ⚠️ **Hessian-source / step-control params (`init_hess`, `irc_direction`, `Init_Displ_DE`, etc.) are not currently read by the IRC engine code path** ([server/workflow/engines/orca.py:218](server/workflow/engines/orca.py#L218)). The engine emits a minimal `%irc MaxIter ... end` block from `max_iterations` only. To control Hessian source, direction, or initial displacement, either (a) extend the node-def, or (b) skip the workflow engine and run ORCA directly on a hand-edited input file using the Stage 1 template above. ### 4. Run with the full HPC run_config IRC chains many SCFs along the path — bump `walltime` and prefer `shared`/`compute` over `debug`. Read `server/templates/orca_generic.sh` and pass its contents as `default_template`. ```json catgo_workflow(action: "run", workflow_id: "", run_config: { "execution_mode": "hpc", "default_session_id": "", "base_work_dir": "/expanse/lustre/projects/sdp126/jyang25/ORCA/catgo", "default_job_params": { "nodes": 1, "ntasks": 8, "cpus_per_task": 1, "walltime": "04:00:00", "partition": "shared" }, "cluster_configs": { "": { "account": "sdp126", "partition": "shared", "module_loads": "module load cpu/0.17.3b\nmodule load gcc/10.2.0/npcyll4\nexport PATH=$HOME/openmpi-4.1.8/bin:$PATH\nexport LD_LIBRARY_PATH=$HOME/openmpi-4.1.8/lib:$LD_LIBRARY_PATH", "orca_dir": "/home/jyang25/orca_6_1_1_RRP8", "default_template": "", "default_job_params": { "nodes": 1, "ntasks": 8, "cpus_per_task": 1, "walltime": "04:00:00", "partition": "shared" } } } }) ``` The local-scratch template stages I/O to `$TMPDIR/orca_$SLURM_JOB_ID` and copies results back. Required on Expanse — Lustre kills ORCA's many-small-file I/O during per-step SCFs. ### 5. Monitor ```json catgo_workflow(action: "status", workflow_id: "") ``` ### 6. Pull files for post-processing When status is COMPLETED, pull the output and trajectory files referenced in Stage 2 below: ```bash mkdir -p ./local_run for f in ORCA.out ORCA_IRC.dat ORCA_IRC.xyz ORCA_IRC_F.xyz ORCA_IRC_B.xyz \ ORCA.property.json ORCA.json; do curl -s -X POST http://localhost:8000/api/hpc/files/read-content \ -H 'Content-Type: application/json' \ -d "{\"session_id\":\"\",\"file_path\":\"/$f\"}" \ > ./local_run/$f done ``` Then run the parser/plotter from Stage 2 against the local copy. ### Submission gotchas - `catgo_workflow_engine.add_task` doesn't auto-attach the viewer structure → "No input structure provided". - `partition=workq` (Shaheen default) is invalid on Expanse → use `debug`/`shared`/`compute`. - `partition=debug` capped at 30 min — IRC almost always needs more. - Missing `account=sdp126` → "Invalid account or account/partition combination". - Missing `module_loads` + `orca_dir` → `orca` not on PATH; per-step SCFs silently produce nothing. - After re-connecting to Expanse, the session_id changes — re-discover via `/api/hpc/connections` and update both `default_session_id` and the `cluster_configs` key. - A `.hess` from a different method/basis/grid is a silent footgun for `init_hess: "read"` — only reuse a `.hess` produced at exactly this IRC's level of theory. ## Stage 2: Post-Processing After the IRC completes, ORCA produces: | File | Contents | |------|----------| | `{basename}.out` | Main output with per-step energies/gradients and the `IRC PATH SUMMARY` table | | `{basename}_IRC.dat` | Clean tabular file: step, energy (Eh), path length — easiest programmatic source when present | | `{basename}_IRC.xyz` | Full trajectory: concatenated geometries from backward end → TS → forward end | | `{basename}_IRC_F_trj.xyz` | Forward-arm trajectory | | `{basename}_IRC_B_trj.xyz` | Backward-arm trajectory | | `{basename}_IRC_F.xyz` | Final endpoint geometry (forward direction) | | `{basename}_IRC_B.xyz` | Final endpoint geometry (backward direction) | The forward/backward labeling is arbitrary — it depends on the sign of the initial displacement. The user identifies which endpoint is reactant vs. product by looking at the geometries. ### Parsing the Energy Profile Parse the `IRC PATH SUMMARY` table in the main output file. OPI does not model this block as a typed object, but its grepper recipe replaces the bespoke "find marker, slice 10 KB, walk lines" boundary handling with a single call. Requires `pip install orca-pi`. Each row has the shape: ``` [<= TS] ``` Reference parser using OPI's grepper: ```python import re import sys sys.path.insert(0, ".claude/skills") # for the _shared helper from _shared.orca_opi import grep_block ROW_PATTERN = re.compile( r"(\d+)\s+([-\d.]+)\s+([-\d.]+)\s+([\d.]+)\s+([\d.]+)\s*(.*)" ) def parse_irc_path(output_file): # offset=2 skips the header + ruler; count is a safe upper bound. lines = grep_block(output_file, "IRC PATH SUMMARY", offset=2, count=400) steps = [] ts_raw_step = None for line in lines: line = line.strip() if not line: continue m = ROW_PATTERN.match(line) if not m: continue raw_step = int(m.group(1)) is_ts = "<= TS" in m.group(6) if is_ts: ts_raw_step = raw_step steps.append({ "step": raw_step, "energy_eh": float(m.group(2)), "dE_kcal_mol": float(m.group(3)), "max_gradient": float(m.group(4)), "rms_gradient": float(m.group(5)), "is_ts": is_ts, }) # Re-index so backward steps are negative, TS = 0, forward positive if ts_raw_step is not None: for s in steps: s["step"] -= ts_raw_step return steps # Sanity check before parsing — replaces hand-grep of "ORCA TERMINATED NORMALLY" from opi.output.grepper.recipes import has_terminated_normally assert has_terminated_normally(output_file) ``` ### Plotting the Energy Profile Plotting matches the style used by CatGo's `IrcPathPlot.svelte`: ΔE in kcal/mol on the y-axis, IRC step number (re-indexed relative to the TS) on the x-axis, backward arm in purple, forward arm in green, TS marked with a red point and an amber dashed vertical line. ```python import matplotlib.pyplot as plt BACKWARD_COLOR = "#8b5cf6" FORWARD_COLOR = "#10b981" TS_COLOR = "#ef4444" TS_GUIDE = "#f59e0b" def plot_irc_profile(steps, output_png="irc_profile.png"): xs = [s["step"] for s in steps] ys = [s["dE_kcal_mol"] for s in steps] ts_idx = next((i for i, s in enumerate(steps) if s["is_ts"]), None) fig, ax = plt.subplots(figsize=(8, 5)) if ts_idx is not None: ax.plot(xs[: ts_idx + 1], ys[: ts_idx + 1], color=BACKWARD_COLOR, linewidth=1.5, marker="o", markersize=3, label="Backward") ax.plot(xs[ts_idx:], ys[ts_idx:], color=FORWARD_COLOR, linewidth=1.5, marker="o", markersize=3, label="Forward") ax.scatter([xs[ts_idx]], [ys[ts_idx]], color=TS_COLOR, s=60, zorder=5, label="TS") ax.axvline(xs[ts_idx], color=TS_GUIDE, linestyle="--", linewidth=1, alpha=0.7) else: ax.plot(xs, ys, color=FORWARD_COLOR, linewidth=1.5, marker="o", markersize=3) ax.set_xlabel("IRC Step (relative to TS)") ax.set_ylabel("ΔE (kcal/mol)") ax.legend(loc="best", frameon=False) fig.tight_layout() fig.savefig(output_png, dpi=150) print(f"Profile saved to {output_png}") if __name__ == "__main__": import sys steps = parse_irc_path(sys.argv[1]) out = sys.argv[2] if len(sys.argv) > 2 else "irc_profile.png" plot_irc_profile(steps, out) ``` ### Viewing the IRC profile in the IDE After `plot_irc_profile(...)` writes the PNG, surface it inline with the shared helper: ```python sys.path.insert(0, ".claude/skills") from _shared.orca_opi import show_png show_png("irc_profile.png", "IRC energy profile") # prints `![IRC energy profile](irc_profile.png)` ``` Then **reply to the user with that markdown link** so Claude Code renders the figure inline in chat. Notes: - `dE_kcal_mol` in the PATH SUMMARY is already relative to the TS, so TS sits at ΔE = 0 and both arms descend to negative values. No extra conversion needed. - If the user wants path length on the x-axis instead of step number, use `{basename}_IRC.dat` as the source — it's cleaner than regex-parsing the main output. ### Extracting Endpoint Geometries Read `{basename}_IRC_F.xyz` and `{basename}_IRC_B.xyz` to show the user the endpoint structures. Remind the user that these endpoint geometries are **not** fully optimized — they should run a follow-up geometry optimization on each endpoint to get proper reactant and product structures. ## Common Issues and Troubleshooting **IRC terminates early / doesn't reach a minimum:** Increase `MaxIter` — try 100 or 150 for flat PES regions. Recompute with `InitHess calc_numfreq` if the initial direction looked wrong. **IRC follows the wrong mode:** The initial Hessian is bad or the TS has extra low-lying imaginary frequencies. Use `InitHess read` with a `.hess` file from a converged frequency job at the same level, or switch from `calc_anfreq` to `calc_numfreq`. If the TS has more than one imaginary frequency, the geometry is not a true TS — fix it first. **IRC connects to unexpected structures:** The TS probably corresponds to a different reaction than intended. Common with NEB-TS when the algorithm finds a stepwise mechanism instead of a concerted one. Re-run the TS search with a different guess geometry or a tighter NEB. **Energy goes up instead of down:** Sign of a bad Hessian or the geometry not being a true TS. Recompute the Hessian with `NumFreq` and verify exactly one imaginary mode before rerunning IRC. **Solvation:** Gas-phase IRC is standard for connectivity confirmation. If the user needs a solvated energy profile, single-point CPCM calculations on the gas-phase IRC geometries are more practical than running the full IRC with CPCM.