Golden paths: ferric in applications

End-to-end workflows for real questions, written so an agent can execute them without reading the rest of the repo. Each one states what it answers, what it cannot, and how to tell it went wrong.

Read For agents first if you are automated. It lists the failure modes that cost the most time. Individual steps live in Recipes.

Golden path A — Reaction energetics of a cation cascade

Answers: which intermediate is the bottleneck, and which branch point controls product selectivity.

Worked example: terpene synthase carbocation cascades (geranyl → linalyl → α-terpinyl → product). The same shape applies to any closed-shell cation mechanism.

Step 0 — gate before you spend hours

pip install ferric
git clone https://github.com/mgoldey/ferric && cd ferric   # examples/, testdata/, scripts/
ferric examples/water-rhf.toml    # expect converged = true, energy = -74.9631468000

Don't skip this, and don't build from source first. The wheel installs in about a minute, and a source build takes about 10 minutes. If the install is broken, every later failure will look like a chemistry problem. The clone is needed because the wheel doesn't ship examples/, testdata/ or scripts/.

Run the study itself under scripts/ferric-limited --max=8G --high=7G --. These can be multi-hour jobs, and an OOM kill leaves a truncated log with no error message.

Step 1 — build the species, and check parity

Generate starting geometries however you like; RDKit ETKDG + MMFF is fine. Use force fields for geometry seeds only, never for a reported energy.

Then, before any QM, verify electron parity for every species:

n_elec = sum(ATOMIC_NUMBER[sym] for sym in symbols) - charge
assert (n_elec % 2 == 0) == (multiplicity % 2 == 1), (
    f"{name}: {n_elec} electrons cannot have multiplicity {multiplicity}"
)

A cation is charge=1, multiplicity=1. Getting this wrong is the most common first failure. ferric catches it with a clear message, but checking in your own loop fails faster and names the species.

Check geometry as well as parity. A structure can have a correct formula and still be physically impossible. Measured: a hand-built C10H17+ passed parity (76 electrons, singlet) while carrying a 0.902 A C-H contact, shorter than a real bond, and a hydrogen 1.030 A from a second carbon. The SCF stalled at 37 iterations, which is the correct behaviour for overlapping nuclei.

# formula check and geometry check fail on different mistakes; keep both
assert (n_elec % 2 == 0) == (multiplicity % 2 == 1)   # bad stoichiometry
assert min_interatomic_distance(xyz) > 0.9            # gross overlap only (your own helper, Angstrom)
# The failure above PASSES that screen (closest contact 0.902 A): its H sat
# ~1 A from two carbons at once. Also reject any H bonded to two heavy atoms:
assert max_heavy_neighbours_of_hydrogen(xyz, cutoff=1.3) <= 1

Step 1b: relax with xtb before any DFT

xtb is the middle tier of the funnel. GFN2-xTB fixes a bad structure in seconds:

LD_LIBRARY_PATH=~/.local/lib/x86_64-linux-gnu \
  xtb mol.xyz --opt --chrg 1 --uhf 0 --gfn 2
# -> xtbopt.xyz

A local install needs LD_LIBRARY_PATH. Without it, xtb fails with libxtb.so.6: cannot open shared object file; the install itself is fine.

Measured on the same molecule, method and basis:

hand-built geometry : closest contact 0.902 A -> SCF Stalled @37, E = -389.5110690
xtb-relaxed (54 its): closest contact 1.084 A -> SCF converged,   E = -390.3794234

The relaxed structure is 0.868 Ha (545 kcal/mol) lower and converges; the hand-built one does not.

Step 2 — optimize each intermediate

[molecule]
xyz = "alpha_terpinyl.xyz"
charge = 1
multiplicity = 1

[basis]
name = "6-31g"           # def2-SVP did not fit on a workstation here; see below

[method]
kind = "ksdft"
task = "optimize"

[dft]
functional = "PBE"

[scf]
max_iter = 400           # 173 iterations was needed; the default is 200
energy_conv = 1e-7       # a loose, reachable bound; see below
# density_conv defaults to 1e-6 and is the criterion that does the real work.
# Tighten it if you need more than the energy; see below.

[optimize]
max_steps = 60

These settings follow from the measurements below. Don't raise the basis or tighten energy_conv without measuring one species first.

Delocalised cations need iterations, and a reachable energy_conv. Measured on the same 27-atom cation, changing only these two knobs:

max_iter 100, energy_conv 1e-8  ->  MaxIter @100,  29 min, 3.30 GB, no result
max_iter 400, energy_conv 1e-7  ->  CONVERGED @173, 4.9 min, 1.81 GB

173 iterations is normal for a delocalised carbocation; the default cap of 200 is not sized for this. energy_conv = 1e-8 was unreachable: under density fitting the energy change floors with naux, so a dE bound tighter than the RI noise floor is never met and the run spends its whole iteration budget on it. density_conv is the criterion that converges the density. Leave the dE bound loose enough to be satisfiable.

The exit reason is diagnostic. Stalled means the SCF could not make progress: suspect the geometry. MaxIter means it was progressing and ran out of iterations: raise the cap.

Set density_conv as well as energy_conv. E_HF is variational and tolerates a loose density; anything depending linearly on the MO coefficients does not. A correlation energy computed from an unconverged density can look reproducible and still shift when anything upstream changes, because the density was never converged.

Checkpoint every species as it finishes. A crash at species 6 of 8 should cost one species, not the whole run.

The method and basis decide whether the study is feasible. Measured on the same 27-atom cation, single point (both the basis and the functional differ between the two rows, so this does not isolate the cost of the basis):

def2-SVP / B3LYP : >22 min, 6.04 GB, SIGKILLed before converging
6-31G   / PBE    :  4.9 min, 1.81 GB, CONVERGED

6-31G/PBE used a third of the memory and finished. Don't plan a def2-SVP cascade on a workstation without measuring one species first.

Treat the optimization time as an estimate. An optimization takes tens of gradient steps. Estimate: at ~5 min per 6-31G single point, 30–60 steps would take roughly 2.5–5 hours per species if each gradient step costs about as much as a single point. Neither a gradient step nor a full optimization of this system has been timed. Later steps start from the previous density and may converge faster than a cold start. Time one optimization before budgeting a cascade.

Use a release build. Measured: the same cation at B3LYP/def2-SVP ran for more than 34 minutes on a debug binary, sharing the machine with one other job, without converging a single point. The wheel and cargo build --release are both release builds.

A 27-atom cation at def2-SVP needs ~6 GB. Budget the study with that in mind. Running eight species of that size at once exhausts a workstation. Always wrap the run in scripts/ferric-limited --max=8G --high=7G --. For agents explains why the in-process budget alone is not enough.

Size it before you start. Measured: the tert-butyl cation (C4H9+, 101 basis functions at def2-SVP) is a single point that takes seconds. A C10 cascade cation (C10H17+) has 225 basis functions. Estimated at N^3.5 scaling, that is about 16x the single-point cost. That is cheap for an energy, but an optimization multiplies it by the step count, so the optimizations of the C10 species dominate the cost of the study. Run one species end to end and time it before you queue eight.

Step 3 — confirm the ordering survives the functional

Run single points on the optimized geometries with at least two more functionals, for example B3LYP and wB97X-V (both available; wB97X-V resolves to libxc's HYB_GGA_XC_WB97X_V), at a basis that fitted in Step 2.

If the functionals disagree on the ordering of intermediates, report the disagreement. Do not average them or pick the one that matches your hypothesis. A functional-dependent ordering is a finding about the system.

Step 4 — properties at the key intermediate

Hirshfeld and Löwdin charges, the ESP at the nuclei, and the static polarizability of the π-stabilized cation. From Python these are ferric.hirshfeld_charges, ferric.lowdin_charges and ferric.esp_at_atoms.

Cost warning: the polarizability comes from ferric-rpa and is an RPA-level calculation. Budget it separately from the DFT run. It can cost more than the optimization before it.

Step 5 (optional) — barriers

Intermediates alone give you a thermodynamic profile. If you need barriers, ferric has a Python-only transition-state toolchain. None of it is wired to the CLI's method.task:

One imaginary frequency is necessary but not sufficient for a transition state: a methyl rotor also gives one. Use run_irc to confirm that the saddle connects the two intermediates you meant.

What this cannot tell you

State these in any write-up. They define the scope of the claim.

• Gas-phase cluster models. There's no enzyme environment unless you add QM/MM (QM/MM).
• Analytic Hessians for RHF, UHF and ROHF only. Their frequencies use the analytic Hessian (RHF and UHF also with an ECP or fixed external charges or field; ROHF with neither); KS-DFT (ROKS included), implicit solvent and polarizable embedding frequencies, and every TS search, use finite differences of analytic gradients (Step 5). That costs 6N gradients per Hessian, so a C10 cation Hessian is hundreds of gradient evaluations. Without Step 5 you have intermediate energies but no barriers, and kcat depends on barriers.
• No thermochemistry. There's no entropy, enthalpy or free-energy correction anywhere, from the CLI or Python. You can compute a zero-point energy yourself from the frequencies (½Σhν over the real modes). Python's FrequencyResult doesn't provide one.
• Relative energies only. Absolute totals carry basis-set and functional errors far larger than the differences you are interpreting.

Reading the output

Before publishing any number

Check the capability's grade in Capabilities and validation. A method being available from the CLI doesn't mean its numbers are production-grade. The grades are Proven, Proven (narrow), Smoke and Spike; a few kinds are not graded.