ferric

A quantum chemistry engine written in Rust, with Python bindings and a TOML-driven command line. It computes Hartree–Fock and DFT energies and gradients, the MP2 family, coupled cluster, RPA and GW, and constrained DFT, with libint2 for the integrals.

pip install ferric      # Linux x86_64, Python 3.10–3.13

Start here

Implemented ≠ validated

Working code is not a checked number. These pages describe what exists; how far each capability's numbers have been checked against an independent reference differs a lot between methods. Each one is graded individually in Capabilities and validation. Read it before relying on a result.

The idea

ferric is organized around electronic response: how the density reacts to a perturbation. That object appears as the polarizability \( \alpha \), the dielectric function \( \varepsilon \) and the susceptibility \( \chi \), and three of the method families here are different ways of getting it right where standard methods get it wrong:

• Attenuated MP2: MP2 builds dispersion from an uncoupled polarizability, which overbinds some systems (π-stacking is the classic case). Attenuating the correlation operator removes the long-range part where that error lives.
• PDEP-RPA / GW: the dielectric matrix is the density–density response. PDEP keeps only its dominant eigenmodes, a low-rank representation of the screening used by RPA correlation and the GW interaction.
• Constrained DFT: a constraint couples to the density and reads its response, building charge-localized diabatic states and their electron-transfer couplings.

Electronic response develops the argument and says which parts of it are demonstrated and which remain a design premise.

Source, license, citation

github.com/mgoldey/ferric, dual-licensed MIT / Apache-2.0. To cite ferric and the methods you used, see References and citing.

Installation

There are two ways in. Most people want the first.

Fastest: the prebuilt wheel

Wheels are published to PyPI for Linux x86_64 (manylinux_2_28), CPython 3.10–3.13. libint2 and libxc are statically linked into the extension, and OpenBLAS ships inside the wheel as a bundled shared library, so nothing needs compiling. The wheel's libint2 carries second derivatives, so analytic Hessians work from a plain pip install; see What the libint2 build carries.

pip install ferric        # or: uv pip install ferric

The only releases so far are pre-releases (0.1.0rc*). pip and uv both select them when no final release exists, so the plain command above works. Pin a version (ferric==0.1.0rc5) if you need a reproducible environment.

The wheel gives you two things:

• the Python module: import ferric
• a ferric command on PATH that runs TOML input files (the same CLI as cargo run --bin ferric)

Check it works

python -c "import ferric; print(ferric.__file__)"

Then run your first calculation. It takes under a second.

What the wheel does not contain

The wheel holds the compiled library and nothing else. The examples/ input files, the testdata/ molecules and the tools/ pipeline scripts live in the git repository. Pages that use them say so. To have them locally:

git clone https://github.com/mgoldey/ferric
cd ferric
ferric examples/water-rhf.toml

The CLI resolves [molecule] xyz = "..." relative to the current directory, so run example files from the repository root.

Building from source

Build from source if you are changing ferric, need the MPI build, or are on a platform without a wheel. ferric links libint2, a C++ integral library; scripts/install-libint.sh installs a prebuilt copy in seconds.

Prerequisites

• Rust 1.75+: install via rustup
• libint2 2.13.1: the conda-forge Linux x86-64 build, which includes second-derivative integrals; scripts/install-libint.sh downloads it (checksum pinned) and needs curl, python3, zstd and patchelf
• libxc, OpenBLAS, LAPACK, Eigen3 and Boost headers
• Python 3.10+ and maturin: optional, for the Python bindings

Steps (Ubuntu 22.04+)

The apt list matches what CI installs (.github/workflows/ci.yml).

# 1. System dependencies
sudo apt-get install -y build-essential cmake g++ gfortran wget git \
    libeigen3-dev libopenblas-dev liblapack-dev pkg-config \
    libxc-dev libboost-dev patchelf zstd \
    python3-dev python3-pip python3-venv

# 2. Get ferric and install libint2 (seconds)
git clone https://github.com/mgoldey/ferric
cd ferric
scripts/install-libint.sh ~/.local/libint2-2.13.1
export LIBINT2_PREFIX=~/.local/libint2-2.13.1

# 3. Build ferric
cargo build --release

# 4. Check it: seconds, not minutes
OPENBLAS_NUM_THREADS=1 ./target/release/ferric examples/water-rhf.toml

The last line should end with converged = true and energy = -74.9631468000 Hartree, as shown in the first calculation.

LIBINT2_PREFIX is read by crates/ferric-integrals/build.rs; unset, it defaults to the installer's ~/.local/libint2-2.13.1 when that exists, and to ~/.local otherwise, so the export above is optional for the default path. The install script records the library's absolute path in every binary that links it, so moving or deleting that directory breaks the build's executables until they are rebuilt. A prefix holding a static libint2.a (a from-source libint2 build) also works, but a library generated without second derivatives has no analytic Hessians; frequencies then use finite differences of the analytic gradient.

The workspace has three binaries (ferric, ferric-cli and ferric-batch), so cargo run needs to be told which one: cargo run --release --bin ferric -- input.toml.

Running the test suite

OPENBLAS_NUM_THREADS=1 cargo test --workspace

.cargo/config.toml sets OPENBLAS_NUM_THREADS=1 for cargo-invoked processes only when the variable is unset (its [env] entry has no force = true, so an exported value wins). The prefix above makes this command use 1 regardless of your shell. Do not raise it above 1; see Threading. The Python binding tests are pytest, not cargo; see CONTRIBUTING.md.

Debug vs release

• Debug (cargo build): fast to compile, slow to run. Use it while iterating on Rust code. It catches debug_assert! violations and integer overflow that release builds silently allow.
• Release (cargo build --release): slow to compile, fast to run. Use it for anything you will actually wait on: real molecules, benchmarks, and any RPA/GW/CC job.

Python bindings from source

uv sync
uv run maturin develop --release
uv run python -c "import ferric; print(ferric.__file__)"

Use uv run maturin develop, not a bare maturin develop. A bare one can install into a different interpreter from the one uv run python loads, and the stale build keeps getting imported.

What the libint2 build carries

Integral classes, derivative orders and angular-momentum limits are fixed when libint2's source is generated, so no build flag changes them. Highest angular momentum for energy / 1st / 2nd derivatives (– = not generated):

With either build, analytic Hessians cover orbital bases up to f functions; a basis with g or higher functions uses finite differences of the analytic gradient. The wheel has no G12 class, so F12 methods need a source build. scripts/generate-libint-small.sh regenerates the wheel's export. The upstream mpqc4 tarball (libint2 2.7.2) has no second derivatives and no G12 class: built against it, ferric uses finite-difference Hessians and the G12 tests skip with an explicit message.

MPI

MPI is a source build only. There is no MPI wheel on PyPI: the wheels workflow deliberately does not publish one.

sudo apt-get install -y libopenmpi-dev openmpi-bin libclang-dev
cargo build --release --workspace --features mpi
mpirun -np 4 -x OPENBLAS_NUM_THREADS=1 ./target/release/ferric examples/water-rhf.toml

libclang-dev is needed by mpi-sys's bindgen step. MPICH and Intel MPI are not supported: their libmpi.so.12 ABI differs from OpenMPI 4.x's libmpi.so.40.

The CLI is the MPI entry point; the Python API is not

The CLI is SPMD by design: mpirun -np N ferric input.toml runs one rank per process.

mpirun -np N python script.py is not supported. The bindings expose no rank or world-size accessor, so if rank == 0 cannot be written. Every rank runs the whole script, prints N times, and races on the same output files.

Threading

The ferric binary and import ferric pin OpenBLAS to one thread when OPENBLAS_NUM_THREADS is unset, and honour the variable when it is set (cargo sets it to 1 through .cargo/config.toml, but only when it is unset, so an exported value reaches cargo-launched processes). Leave it unset or at 1. ferric uses rayon for its parallelism; a threaded OpenBLAS on top of that oversubscribes the machine, and its LU routines can crash when called from rayon workers. For throughput across many independent jobs, prefer many single-threaded processes over one multi-threaded job.

Your first calculation

This page takes you from an installed ferric to one checked number, then shows where to go next. It assumes you have run pip install ferric (Installation). No git clone is needed until the last section.

We compute the Hartree–Fock energy of water in the minimal STO-3G basis. It takes well under a second and has a known answer, so you can tell immediately whether your installation is right.

1. From Python

Save this as water.py:

import ferric

# XYZ format: atom count, a comment line, then symbol x y z in Ångström.
water = ferric.Molecule.from_xyz_string("""3
water
O   0.000000   0.000000   0.117790
H   0.000000   0.755453  -0.471161
H   0.000000  -0.755453  -0.471161
""", 0, 1)                       # charge 0, spin multiplicity 1 (singlet)

basis = ferric.BasisSet.bundled("sto-3g")
rhf = ferric.run_rhf(water, basis)

print(rhf.converged, f"{rhf.energy:.10f}")

Run it:

OPENBLAS_NUM_THREADS=1 python water.py

You should see:

True -74.9631468000

Two things to notice:

• Always read converged. When run_rhf runs out of iterations it still returns an energy; it sets the flag and does not raise. A number with converged == False is not a result. See Sharp bits.
• Energies are in Hartree. Geometries go in as Ångström; see Sharp bits for which accessors return Bohr.

2. The same thing from the command line

The wheel also installs a ferric command that reads a TOML input file. Write the geometry to water.xyz:

printf '3\nwater\nO 0.000000 0.000000 0.117790\nH 0.000000 0.755453 -0.471161\nH 0.000000 -0.755453 -0.471161\n' > water.xyz

and write water-rhf.toml next to it:

[molecule]
xyz = "water.xyz"      # resolved relative to the directory you run ferric from

[basis]
name = "sto-3g"

[method]
kind = "rhf"

OPENBLAS_NUM_THREADS=1 ferric water-rhf.toml

The output ends with:

RHF/sto-3g on water.xyz
  nbasis     = 7
  iterations = 8
  converged  = true
  energy     = -74.9631468000 Hartree

Same molecule, same number. The CLI rejects any key it does not recognise, so a typo is an error rather than a silently ignored setting. Every key is listed in the input reference.

3. Add correlation

To go beyond Hartree–Fock, change the method and the basis. RI-MP2 needs an auxiliary (fitting) basis alongside the orbital basis, and STO-3G is too small for a meaningful correlation energy, so this example moves to cc-pVDZ:

bs  = ferric.BasisSet.bundled("cc-pvdz")
aux = ferric.BasisSet.bundled("cc-pvdz-ri")
mp2 = ferric.run_rimp2(water, bs, aux)
print(f"RI-MP2 total energy: {mp2.total_energy:.10f} Ha")

RI-MP2 total energy: -76.2308014541 Ha

In TOML the same calculation changes the [basis] name and the method kind, and adds an [mp2] section naming the auxiliary basis:

[basis]
name = "cc-pvdz"

[method]
kind = "rimp2"

[mp2]
auxbasis = "cc-pvdz-ri"

If auxbasis is omitted the CLI uses cc-pvdz-ri, but name it explicitly: the auxiliary basis should match the orbital basis. This is examples/water-rimp2.toml.

4. Where to next

Running the bundled examples

The repository has one input file per method under examples/, indexed in Examples. They refer to molecules under testdata/, which the wheel does not include, so run them from a clone:

git clone https://github.com/mgoldey/ferric && cd ferric
OPENBLAS_NUM_THREADS=1 ferric examples/water-rimp2.toml
OPENBLAS_NUM_THREADS=1 ferric examples/water-attmp2.toml   # attenuated MP2, ω = 0.420 Å⁻¹

If you built from source rather than installing the wheel, replace ferric with cargo run --release --bin ferric --.

Sharp bits

Things that behave differently from what you might expect, each of which has cost someone a run. This page follows the pattern of JAX's "Sharp Bits": each entry says what surprises people, why ferric does it that way, and what to do.

Leave OPENBLAS_NUM_THREADS unset

What happens: the ferric CLI and import ferric pin OpenBLAS to one thread when OPENBLAS_NUM_THREADS is unset, and honour the variable when it is set. (.cargo/config.toml sets it to 1 for cargo-launched processes, also only when it is unset.) A value above 1 is therefore a deliberate choice, and its effect depends on the workload: it can speed up one calculation and slow another down several-fold.

Why: ferric's parallelism comes from rayon. OpenBLAS threads run inside rayon's worker threads, so the thread count is roughly BLAS threads × rayon threads. Measured on benzene with the ferric CLI, single runs on a heavily loaded machine with 6 physical cores (12 hardware threads), so read these as order-of-magnitude:

"unset" and 1 run the same configuration, so the gap between those two columns is the machine's load, not a setting. B3LYP timings were noisy, with no clear trend. None of these runs crashed.

Three effects matter beyond speed:

• Oversubscription can stall a run. Once BLAS threads × rayon threads exceed the physical cores, a job can slow by orders of magnitude, as the 12-thread RI-JK run shows.
• Results change slightly. The BLAS thread count changes the order of floating-point summation. The RI-JK energy at 6 threads differs from the 1-thread energy by 6e-6 Ha; expect differences at the 1e-6 to 1e-5 Ha level.
• LU routines can crash. Threaded OpenBLAS LU factorization (dgetrf) can overflow the 2 MB stacks of rayon's worker threads. ferric calls LU solves in BSE, PCM, COSMO, MBD and orbital-rotation steps, among others.

Do: leave OPENBLAS_NUM_THREADS unset, and ferric runs with one BLAS thread. If you set it higher for a particular job, never above the number of physical cores, and check the timing and the energy against a one-thread run of the same input. For throughput across many molecules, run many single-threaded processes.

converged is a flag, and it does not mean "ground state"

What happens: run_rhf, run_uhf and run_rohf return a result object even when they hit max_iter; they do not raise. (run_dft, run_ksdft and the correlated drivers behave differently: they raise if their reference SCF does not converge.)

Why: a non-converged state is still useful for diagnosis, and raising would discard it.

Do: check result.converged in Python, or converged = true in CLI output, before using any number, and treat False as "no result". True means the SCF reached a stationary point, not necessarily the lowest one; open-shell molecules especially can land on an excited solution. Comparing UHF with ROHF, or a different starting guess, is a cheap check. See Python bindings.

energy_conv is a sanity bound, not a target

What happens: setting [scf] energy_conv very tight (say 1e-10) can make a KS-DFT calculation never converge.

Why: convergence is driven by the density criterion (density_conv). The energy change between iterations has a noise floor set by the integration grid and density fitting, and a tight energy bound can sit below it. The CLI default is deliberately loose (1e-3).

Do: tighten density_conv, not energy_conv.

Density fitting is on in run_dft but off in run_rhf

What happens: comparing a ferric DFT energy with a code that uses exact Coulomb gives a discrepancy that grows with molecule size.

Why: run_dft, run_ksdft and the CLI's ksdft density-fit Coulomb (RI-J) with def2-universal-jkfit by default, and, for functionals with exact exchange (hybrids and range-separated hybrids), exchange (RI-K) with the same basis. run_rhf builds exact four-centre J and K unless you pass df_j_aux/df_k_aux. The difference is the fitting error that density fitting always carries, not a bug. Against exact J at PBE/STO-3G it is 0.28 kcal/mol for water, 1.16 for benzene and 9.5 for a 71-atom drug molecule.

Do:

• To compare with an exact-Coulomb code (ORCA with NORI, PySCF without density_fit()), pass df_j_aux="exact" to run_dft or run_ksdft, and df_k_aux="exact" too for a hybrid. "", "none", "off" and "conventional" mean the same. In the CLI, set [scf] df_j_aux = "" and df_k_aux = "". The CLI's SCF J/K log line then reads RI-JK via with a blank name; the run uses exact J and K.
• Or fit on both sides: density_fit(auxbasis="def2-universal-jkfit") in PySCF, or run_rhf(..., df_j_aux="def2-universal-jkfit", df_k_aux="def2-universal-jkfit") in ferric, so J and K are both fitted on both sides. run_rhf takes a basis name or "" here, not "exact".
• run_qmmm (KS methods), run_gw and run_u_gw with xc, run_tddft with a functional, run_tdhf_static_polarizability, run_double_hybrid and run_rs_mp2_rpa always density-fit their reference SCF with def2-universal-jkfit and have no opt-out.

Units differ by interface and by accessor

Do: check the docstring of any coordinate accessor before converting it. Converting both QM/MM accessors "to be safe" makes one of them wrong by a factor of 1.89.

Energies take a BasisSet; geometry changes take a basis name

What happens: run_rhf(mol, ferric.BasisSet.bundled("sto-3g")) works, but run_optimize wants the string "sto-3g".

Why: entry points that move atoms (optimization, frequencies, saddle search) build the basis themselves for each geometry, so they take its name.

Do: check the signature in Python bindings.

Not everything is on the CLI

The CLI runs the method.kinds in the capability matrix, with task = energy, optimize or frequencies. It also has [qmmm], [cosmo], [pcm], [external_potential], [scf] stability_descent, [mp2] att_operator and [dft] grid_prune. These capabilities are Python-only:

Examples need the repository

examples/*.toml point at molecules under testdata/, and the pipeline tools live under tools/. Neither ships in the wheel. Clone the repository and run from its root; see Installation.

[memory] budget_gb does not cap the whole process

What happens: a job's resident memory can exceed budget_gb, and a job whose budget is too small usually still runs, slower, rather than failing. Benzene PBE/def2-SVP with budget_gb = 0.002 finishes with the same SCF iterations, a peak RSS of 157 MB, and 11.2 s of wall time instead of 3.5 s.

Why: the CLI turns budget_gb (or, when unset, 0.8 × available RAM) into one process-wide ledger. The large, size-dependent allocations reserve their bytes from it before allocating: three-index RI tensors, the DFT grid's AO cache, MO-transformed RI blocks, and the large tensors of the MP2, RPA, GW and CC methods. Two such allocations alive at the same time therefore cannot each claim the whole budget. Nothing else is charged: basis-sized matrices (Fock, density, DIIS history), integral engines, BLAS and per-thread scratch, allocator overhead. The spill path's two scratch blocks are reported but not refused, and together they can reach about twice the budget.

When a charged allocation does not fit, the result depends on the allocation:

• The SCF's three-index tensor is spilled to a file in $TMPDIR (/tmp by default) and reread on every iteration.
• The DFT grid AO cache is recomputed at every Fock build instead of stored. The energy is bit-identical.
• An allocation with no fallback stops the job with an error naming it. For some methods this check comes after the SCF. Benzene RI-MP2/def2-SVP with budget_gb = 0.001 completes the SCF, then exits with RI-MP2 MO-side blocks ... requires 0.01 GB; budget is 0.00 GB.

Some stages print a warning when resident memory passes 110% of the budget. The warning never stops the run.

From Python, memory_budget_gb= sets the same per-allocation limits, but no shared ledger is installed. Each check compares its own allocation with the whole budget, not with what the other allocations have left. Two checks instead subtract the process's current RSS and allow 90% of the remainder: the KS-DFT grid AO cache (store or recompute) and the UKS Newton/TRAH f_xc kernel's second grid cache. The f_xc check subtracts RSS in the CLI too.

Do: leave room below the machine's real limit for the uncharged part. On a shared machine, also cap the process externally so a runaway job dies in its own cgroup, for example with scripts/ferric-limited -- ferric input.toml (defaults MemoryMax=12G, MemoryHigh=10G, no swap) or systemd-run --user --scope -p MemoryMax=12G -- ferric input.toml.

Implemented is not validated

A method listed on these pages exists and runs. How closely its numbers have been checked against an independent reference differs a lot from method to method. Look up the grade in Capabilities and validation before relying on a number.

ferric recipes

Copy-paste workflows that end in a number. Each gives its expected output and rough runtime, so you can tell "still running" from "silently wrong". That confusion is the most common failure when you drive a QC code you haven't used before.

Recipes 0, 1, 3 and 4 were executed as written for this page, on 2026-09-23. Recipe 2's output is quoted from an earlier run. Recipe 5 is a sketch, marked as such where it appears. Numbers labelled MEASURED came out of the program.

Prerequisite: a working ferric. The fast path is the prebuilt wheel (about a minute). A source build takes ~30 minutes, most of it libint2:

pip install ferric

See Installation. Build from source only if you are changing ferric itself or need MPI.

Which recipes need a git clone. The wheel contains the compiled library and the ferric command, nothing else. examples/, testdata/, scripts/ and the tools/ Python package live in the repository (what the wheel does not contain). Run those recipes from the repository root:

On a shared machine, run anything real under a memory cap. ferric's memory budget charges only the large, size-dependent tensors; basis-sized matrices, integral engines, BLAS scratch and allocator overhead are not charged, so a job's resident memory can exceed the budget. Only a cgroup puts a hard ceiling on the process, and without one an overshoot can trigger the system-wide OOM killer, which may kill unrelated processes. The details are in Sharp bits and For agents. scripts/ferric-limited (repository, Linux with systemd) wraps a command in a systemd-run --user scope:

scripts/ferric-limited --max=8G --high=7G -- ferric input.toml

0. From a SMILES to an energy

import ferric
from tools.structure import from_smiles        # needs a clone + RDKit

mol = from_smiles("CCO")                        # ethanol
res = ferric.run_dft(mol, ferric.BasisSet.bundled("def2-svp"),
                     functional="PBE", dispersion="d3bj")
assert res.converged                            # always check this first
print(res.total_energy, res.e_scf, res.e_dispersion)

MEASURED (2026-09-23, 4 threads on a loaded machine, about 4 s):

-154.72507297  -154.72071374  -0.00435923

• from_smiles returns an RDKit ETKDG embedding plus an MMFF cleanup, seeded so it is reproducible. That is a starting structure, not a minimum. For anything you will report, relax it (xtb, then a ferric optimization) first. See recipe 3 and Golden paths Step 1b.
• total_energy is e_scf + e_dispersion. dispersion=None (the default) leaves e_dispersion as None, meaning "not evaluated", never 0.0.
• run_dft density-fits the Coulomb term by default. That matters when you compare against a code using exact Coulomb (recipe 6).
• from_smiles reads the charge from the SMILES. Spin is never inferred: pass multiplicity=2 for a radical.

For a file instead of a SMILES, tools.structure.read("x.sdf") returns the same kind of Molecule (PDB, mmCIF, PQR, SDF, mol2, .gro and XYZ are supported). A PDB must already carry explicit hydrogens.

1. Single-point energy from a TOML file

The first success that confirms your install works. From the repository root:

ferric examples/water-rhf.toml

MEASURED (about 2 s including start-up):

RHF/sto-3g on testdata/molecules/water.xyz
  nbasis     = 7
  iterations = 9
  converged  = true
  energy     = -74.9631468000 Hartree

From a source checkout, the same run is cargo run --release --bin ferric -- examples/water-rhf.toml.

For your own molecule you need an .xyz (Å) and a TOML file:

[molecule]
xyz = "mymol.xyz"        # relative to the directory you run ferric from

[basis]
name = "def2-svp"

[method]
kind = "rhf"

For every method.kind, see Capabilities and validation. For every TOML key, see the input reference.

2. Charged and open-shell species

Read this before you run any ion, radical or metal center. charge and multiplicity are [molecule] keys. Most files in examples/ leave them at their defaults (0 and 1). The open-shell exceptions are examples/h_uhf.toml and examples/oh-ugw.toml.

[molecule]
xyz = "tBu.xyz"
charge = 1          # cation
multiplicity = 1    # closed-shell singlet

[basis]
name = "def2-svp"

[method]
kind = "ksdft"

[dft]
functional = "B3LYP"

Relax the geometry with xtb first. A correct formula doesn't guarantee a physical structure, and the parity check below can't tell the difference. MEASURED: a hand-built C10H17+ with a 0.902 Å C–H contact stalled the SCF. After xtb mol.xyz --opt --chrg 1 --uhf 0 --gfn 2 it converged, 545 kcal/mol lower. See Golden paths Step 1b.

MEASURED on the tert-butyl cation C4H9+ (13 atoms), B3LYP/def2-SVP:

converged  = true
energy     = -157.4301362681 Hartree

Expect seconds to about a minute at this size.

The error you will hit first

If the electron count and the multiplicity disagree, you get this error. It's worth learning to recognise:

error: inconsistent charge/multiplicity: 35 electrons with multiplicity 1
implies n_alpha = (35 + 1 - 1) / 2 = 35/2, which is not a non-negative
integer... An odd electron count needs an even multiplicity (2, 4, ...)
and vice versa

It means your geometry, your charge or your multiplicity is wrong. ferric can handle ions and open shells. The message prints the arithmetic and the rule, so check the atom count in your .xyz against the first line of the file, then the charge and multiplicity you set.

A radical (open-shell doublet) uses the same keys:

[molecule]
xyz = "radical.xyz"
charge = 0
multiplicity = 2    # one unpaired electron -> UHF/UKS

From Python, multiplicity belongs to the molecule, not to the run_* call: ferric.Molecule.from_xyz("x.xyz", charge=0, multiplicity=2).

3. Optimize a geometry

ferric examples/h2-lda-opt.toml

MEASURED (about 2 s): converged = true, steps = 2, final E = -1.1212649781 Hartree.

The pattern is task = "optimize" next to any kind that has analytic gradients:

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

[dft]
functional = "B3LYP"

[optimize]
max_steps = 30

Capabilities and validation lists which methods have analytic gradients. Harmonic frequencies use task = "frequencies" (finite differences of the analytic gradient; see examples/water-frequencies.toml).

Runtime depends on the system, so measure before you plan. A whole optimization at drug-like size has not been timed on this page. For scale, two single-point measurements:

• 32 atoms, PBE/def2-SVP: 96 s (the figure recorded in the tools.pipeline.tiers.tier4_dft docstring).
• A 27-atom delocalised cation, PBE/6-31G: 4.9 min, because it needed 173 SCF iterations. At B3LYP/def2-SVP the same molecule was killed for memory after 22 min (Golden paths, Step 2).

An optimization multiplies the single-point cost by the number of steps. Measure one species before you queue many.

4. Ligand screening: force field → xtb → DFT

tools.pipeline.run_funnel runs a tiered funnel: cheap scoring on many candidates, and expensive QM only on the few that survive. Use it when you have N molecules and want DFT numbers on the good ones. Don't write your own loop.

This was run exactly as shown, from the repository root:

from tools.campaign.hierarchy import Tier
from tools.isomers.model import Isomer
from tools.pipeline import Stage, run_funnel
from tools.pipeline.tiers import tier2_forcefield, tier3_gfn2, tier4_dft

# Three isomers of C3H6O2, so comparing absolute energies is meaningful.
smiles = ["CCC(=O)O", "COC(C)=O", "CCOC=O"]
candidates = [
    Isomer(smiles=s, kind="structural", transform=s, parent_smiles=smiles[0])
    for s in smiles
]
stages = [
    Stage(Tier.FORCE_FIELD,   tier2_forcefield, keep=3, name="ff"),
    Stage(Tier.SEMIEMPIRICAL, tier3_gfn2,       keep=2, name="xtb"),
    Stage(Tier.QUANTUM,       tier4_dft,        keep=1, name="dft"),
]
rep = run_funnel(candidates, stages, {"seed": 0xF00D, "basis": "sto-3g"})
print(rep.table())
for iso in rep.survivors:
    print(iso.canonical, rep.value("dft", iso.canonical))

MEASURED (2026-09-23, one process):

tier  stage           in   out  failed     secs   s/cand  note
   2  ff               3     3       0      0.4     0.13  ff: kept 3 of 3 scored
   3  xtb              3     2       0      0.1     0.04  xtb: kept 2 of 3 scored
   4  dft              2     1       0      4.8     2.38  dft: kept 1 of 2 scored
      TOTAL                                 5.3
      dominant tier 4 (dft) = 90% of wall
COC(C)=O -264.5628932123858

STO-3G is a smoke-test basis, so don't read chemistry into which isomer survived. The run shows the plumbing works. Change "basis" in the context for real work.

What the funnel does that a hand-written loop usually doesn't:

• Ranks ascending at every tier. Lower is better, because every tier reports an energy or an energy-like score. For the same reason, rank only candidates that share a molecular formula. An absolute energy of a larger molecule is lower because it has more electrons, not because it is better. For substituent series, rank against the parent (parent_smiles). The pipeline notes §0b cover that.
• Drops failures instead of ranking them. A candidate that a tier failed on is counted as failed. It is never treated as having scored well. That's easy to get wrong by hand, and it silently corrupts a screen.
• Times each tier. The tier that actually costs the run is often not the one the cost table predicts.
• Stops early when the population is empty, instead of running an expensive tier on nothing.

tier4_dft adds D3(BJ) dispersion by default. Pass context["dispersion"] = None for the bare SCF energy. Docking (tier1_dock) needs the ferric[docking] extra, a receptor (context["receptor_pdbqt"]) and a box centre (context["box_center"]). The pipeline notes §0b show the substituent version, with liability flags and parent-relative gating.

Before you rank anything by the DFT tier, read the noise measurement. MEASURED on a real campaign: the best available ΔΔE noise over a pose ensemble was 4.07 kcal/mol, against substituent effects of 1–2 kcal/mol. See the pharma coverage notes.

Related modules in tools/active_site/: ligand_embedding, pose_relaxation, binding_energy, prescreen, pocket_charges, pocket_field.

5. Rank residues for mutation (electrostatic pre-screen)

ferric does not design mutations. It can rank which active-site residues most influence a reactive center, which gives you a short list worth testing instead of a guess.

Sketch, not executed for this page. The two library calls are real. The per-residue split is ordinary Python written for this page, and derive_pocket_charges needs pdb2pqr and a pocket PDB.

from collections import defaultdict
from tools.active_site.pocket_charges import PocketCharges, derive_pocket_charges
from tools.active_site.pocket_field import pocket_field_at_atoms

pocket = derive_pocket_charges("pocket.pdb")   # runs pdb2pqr; fills residue_ids
site_xyz = [(x, y, z)]                         # reactive-center atom(s), Angstrom

# pocket_field_at_atoms returns the TOTAL field. To rank residues, split the
# charges by residue and evaluate each group on its own.
by_res = defaultdict(list)
for q, rid in zip(pocket.charges, pocket.residue_ids):
    by_res[rid].append(q)
contrib = {
    rid: pocket_field_at_atoms(PocketCharges(qs, pocket.source_pdb, pocket.ff), site_xyz)
    for rid, qs in by_res.items()
}   # each value: (N_sites, 4) array of [phi, Ex, Ey, Ez], atomic units

The method:

1. Compute the field each residue produces at the reactive center.
2. Rank residues by their contribution. That's your candidate list.
3. Check the top few with QM/MM (QM/MM; the embedding matches pyscf.qmmm.mm_charge to <1e-8 Ha). This step is what turns the ranking into physics rather than electrostatic hand-waving.

Limits. Include these in any write-up that uses this method. It's a classical point-charge pre-screen, so it has no polarization response of the protein, no sterics, no conformational change on mutation, and no ΔΔG. It ranks hypotheses. A residue it flags is a candidate for QM/MM, not a designed mutation.

6. Comparing against another code

Three causes account for most spurious "ferric disagrees" reports:

Density fitting. ferric's KS-DFT (kind = "ksdft", run_dft) fits the Coulomb term by default, and the fitting error grows with system size. MEASURED at PBE/STO-3G against exact Coulomb: water 0.28, benzene 1.16, and a 71-atom drug molecule 9.5 kcal/mol. Compare like with like. run_dft(..., df_j_aux="exact") turns fitting off, or you can turn it on in the other code.

Grid. ferric's KS-DFT default grid is (75, 110) (radial, angular), flat, with no pruning. PySCF's default grid.level=3 is roughly (75, 302). The difference is worth ~1e-5 Ha on water and grows with atom count, because grid error scales with the number of atoms, not the basis size. A small basis on a big molecule can be worse than a big basis on a small one. Match grids before you conclude anything.

Convergence criteria. Setting energy_conv alone can leave the density loosely converged. Variational quantities (E_HF) don't mind. Anything that depends linearly on the MO coefficients (correlation energies, properties) does. density_conv is the criterion that actually converges the density, so tighten it when you care about those quantities. Keep energy_conv loose: with density fitting, a very tight energy_conv may be unreachable (Golden paths, Step 2).

[scf]
density_conv = 1e-9

Choosing a method

Start from what you want to compute. Each row names a method that ferric implements for that task, where to run it (CLI method.kind or Python only), how far it is validated, and a shipped example to copy. Grades, and the anchors behind them, are on Capabilities and validation. Cost is given as formal scaling with system size N, not as timings.

This page only recommends what the code and its tests support. Where ferric has no validated option for a task, the row says so.

By task

Notes on the table:

• "Proven" for an attenuated MP2 method means the implementation reproduces its reference, not that the method is accurate on your system. Its parameters are fitted to interaction-energy benchmarks in one basis (the aug-cc-pVTZ sets without counterpoise, with frozen core); in another basis, with counterpoise correction, or with frozen_core = 0 (the default), you are extrapolating. See The MP2 family.
• KS-DFT uses density fitting by default (def2-universal-jkfit), and the fitting error grows with system size. Turn it off or match it before comparing with an exact-Coulomb code; see SCF and DFT.
• CC results are sensitive to the aux basis: at CC accuracy the RI error depends on the auxiliary set. See Coupled cluster.

Pairing an orbital basis with its auxiliary basis

Every correlated method needs an RI auxiliary basis ([mp2] auxbasis or [rpa] auxbasis; the auxbasis argument in Python). Use the set built for your orbital basis. The bundled pairs are:

def2-tzvpp-rifit and def2-qzvpp-rifit are also bundled as auxiliary sets. The ECP orbital sets aug-cc-pvdz-pp and aug-cc-pvtz-pp have no bundled -pp-rifit partner. STO-3G and 6-31G have no matching RI set; the examples pair them with cc-pvdz-ri or a larger RI set. Only the names above (plus cc-pvdz-f12 and cc-pvdz-f12-optri) load by name; any other set goes in through [basis] path as a Gaussian-94 file (orbital basis only). def2-universal-jkfit is the default JK set for KS-DFT.

When nothing here fits

• Analytic Hessians cover RHF and UHF (exact J/K, no ECP, basis up to f functions); ROHF, KS-DFT and correlated methods' frequencies are by finite differences of the analytic gradient.
• Open-shell coupled cluster and open-shell TS search are not available (the one open-shell CC-type method, LinLCCD(hh), is library-only). Open-shell KS-DFT is available from the CLI: ksdft with multiplicity > 1 runs UKS, and uhf/rohf with [dft] functional run UKS/ROKS. See Capabilities and validation: open shells in the CLI.
• TDDFT with the XC kernel exists only as a library-only spike; the CLI and Python TDDFT omit it. See RPA, GW and excited states.

For PySCF users

ferric's Python API will look familiar if you use PySCF, but it is built differently. This page maps the PySCF calls you already know onto ferric's, then lists the differences that change numbers or cause errors.

The table only lists pairs where both sides exist. If a PySCF feature has no row, assume ferric does not have it; check the full function reference to be sure. The PySCF names were checked against PySCF 2.12.1.

Rosetta table

Assume import ferric and, on the PySCF side, the matching from pyscf import ... (gto, scf, dft, mp, cc, tdscf, gw, lo, df, qmmm, solvent, geomopt, hessian).

Molecule and basis

SCF and DFT

Geometry and vibrations

Correlation, response and properties

How ferric differs

Stateless calls, not mutable method objects. PySCF builds an mf object, lets you set attributes on it, then runs it. ferric has one function per calculation. Settings are keyword arguments, and the call returns a result object. There is nothing to reconfigure and rerun; call the function again. Correlated drivers run their own reference SCF, so you pass them the molecule and basis, not a converged mf.

Spin is multiplicity, and it lives on the molecule. PySCF's spin is 2S; ferric's multiplicity is 2S+1. It is set when the molecule is built (from_xyz_string(..., multiplicity=3)). run_uhf and run_rohf have no spin keyword.

The auxiliary basis is always explicit. PySCF picks an auxiliary basis for you when you call .density_fit() without one. ferric's MP2, CC, RPA, GW and TDDFT drivers take the RI basis as a required argument. Only SCF-level fitting has a default: run_dft and run_ksdft fit J (and K for hybrids) with def2-universal-jkfit unless you pass df_j_aux/df_k_aux. These fit their reference SCF with def2-universal-jkfit and take no override: run_qmmm with a KS method, run_gw and run_u_gw with xc, run_tddft with a functional, run_tdhf_static_polarizability, run_double_hybrid, and run_rs_mp2_rpa (its HF reference).

Density fitting is on in run_dft and off in run_rhf. run_dft uses RI-J by default, and RI-K for hybrids. dft.RKS in PySCF uses exact Coulomb unless you call .density_fit(). run_rhf builds exact four-centre J and K unless you pass df_j_aux/df_k_aux. The difference between a fitted and an exact energy is the fitting error, not a bug; the ferric source records it at PBE/STO-3G against conventional J as 0.28 (water), 1.16 (benzene) and 9.5 kcal/mol (a 71-atom drug molecule). Pass df_j_aux="exact" to run_dft when you compare against exact-Coulomb PySCF.

The DFT grid is smaller and not pruned. ferric's default is 75 radial × 110 Lebedev angular points on every atom, with no pruning. PySCF's default (grids.level = 3) uses 302 angular points for H through Ne and 434 from Na on, with 50 radial points for H and 75 for C–Ne, and prunes with nwchem_prune. Expect grid-level differences in DFT energies.

Correlated methods use RI integrals throughout. ferric's CCSD and CCSD(T) build their integrals from the auxiliary basis you pass. The fair PySCF comparison is cc.CCSD(mf).density_fit(...), not plain cc.CCSD(mf), which uses exact four-index integrals.

Range-separation ω is in Å⁻¹ in the run_* drivers. omega on run_attenuated_rimp2 and run_rs_mp2_rpa is in Å⁻¹ (default 0.420). PySCF, and ferric's low-level compute_eri3_mo and compute_metric_2c, take Bohr⁻¹.

Point charges are in Bohr; QmmmSystem coordinates are in Å. point_charges= on run_rhf and the other SCF drivers takes Bohr. PySCF's qmmm.mm_charge defaults to mol.unit, which is Ångström unless you changed it. ferric.QmmmSystem takes Ångström coordinates, and its point_charges() accessor returns Bohr, ready to pass to run_rhf.

An unconverged SCF can raise. PySCF's kernel() returns the energy either way. ferric's run_rhf, run_uhf, run_rohf and run_qmmm do the same, so check .converged. run_dft and every correlated driver raise instead when their SCF does not converge.

TDDFT refuses some functionals. ferric's run_tddft includes the f_xc response and matches PySCF TDA/TDDFT to 1e-3 eV for HF, LDA, PBE and B3LYP on the tested systems. It refuses meta-GGA, VV10 and range-separated functionals, which PySCF accepts, and it runs closed-shell references only.

GW defaults to a Hartree–Fock reference. PySCF's gw.GW(mf) runs on whatever mf you pass, usually a DFT one. ferric.run_gw runs its own HF reference unless you pass xc="pbe" or another functional.

The PCM defaults differ. ferric's solvent= is IEF-PCM with 110 tesserae per atomic sphere (pcm_lebedev_order=110). PySCF's PCM defaults to C-PCM with Lebedev order 29 (302 points per sphere) and SWIG discretization. The named solvent "water" is ε = 78.4 in ferric; PySCF's default ε is 78.3553.

AO matrices do not match element by element. ferric's AO basis-function conventions come from libint2 and are not the same as PySCF's (libcint). MEASURED on CO/cc-pVDZ: RHF total energies agree to 3e-12 Ha and orbital energies to 8e-9 Ha, while the two AO density matrices differ element by element by up to 1.5. Compare invariant quantities (energies, orbital energies, charges), not raw AO matrices. Axis order differs too: compute_eri3 returns (naux, n_bf, n_bf), with the auxiliary index first.

Basis sets come from a fixed bundled list. BasisSet.bundled(name) knows 25 names (see Bundled basis sets). There is no Python loader for a basis file or for a per-element basis dictionary.

Results come back as numpy arrays or lists, not live objects. Matrices are numpy.ndarray. Per-atom quantities such as charges are Python lists.

QM/MM

Run a quantum region inside a classical environment: a ligand in a protein pocket, a reacting site in an enzyme, a solute in explicit solvent. The QM atoms get a real wavefunction; everything else becomes point charges that polarize it.

What you get: QM region selection (by index, by radius, or by whole residue), link atoms across covalent cuts with four boundary-charge schemes, Gaussian-smeared charges, Thole polarizable embedding, an AMBER-form MM crate, analytic QM and MM forces, a full-structure gradient across the cut, geometry optimization, and a TIP3P solvation droplet. Structures come from PDB, mmCIF, PQR, SDF, mol2, GROMACS .gro, XYZ or SMILES.

What is missing: no AMBER prmtop reader (go through OpenMM), and no periodic boundary conditions — the solvation droplet is finite, with a vacuum boundary.

Before you start: relax the geometry with xtb

The single highest-value habit in this whole workflow. DFT on an unrelaxed structure wastes hours and often simply fails.

Measured on this project, a 27-atom carbocation built by hand:

hand-built    C7-H24 = 0.902 Å      SCF: Stalled@37, never converged
xtb GFN2      C7-H24 = 1.084 Å      SCF: converged
                                    545 kcal/mol lower

A 0.9 Å C–H bond is not a slightly-off geometry; it is a structure whose SCF has no reason to converge. xtb finds that in seconds where DFT spends hours failing.

xtb input.xyz --opt --gfn 2 --chrg 1        # seconds
ferric single_point.toml                     # then the expensive step

Use it for the initial geometry of anything you did not get from a crystal structure or a previous optimization — hand-built ligands, docked poses, edited residues, anything from a SMILES string. See Applications for the full docking → xtb → DFT funnel.

Build note: xtb built with -O3 miscompiles gradients on this platform. Build with -Doptimization=2.

The smallest working example

import ferric

# The FULL structure: every atom, QM and MM alike.
symbols = ["O", "H", "H",   "O", "H", "H"]     # two waters
coords  = [[0.0, 0.0, 0.0], [0.76, 0.59, 0.0], [-0.76, 0.59, 0.0],
           [0.0, 0.0, 3.0], [0.76, 0.59, 3.0], [-0.76, 0.59, 3.0]]
charges = [-0.834, 0.417, 0.417,  -0.834, 0.417, 0.417]   # MM partial charges (e)

# QM = the first water; the second becomes point charges.
sys = ferric.QmmmSystem(symbols, coords, charges, qm_indices=[0, 1, 2])

res = ferric.run_qmmm(sys, "cc-pvdz", method="rhf")
print(res.energy)

charges are ignored for atoms that end up in the QM region, so you can pass a full force-field charge set and let the selection decide.

Selecting the QM region

Three ways, and the choice matters more than it looks:

# 1. Explicit indices — full control.
ferric.QmmmSystem(symbols, coords, charges, qm_indices=[0, 1, 2])

# 2. Everything within a radius of some seed atoms.
ferric.QmmmSystem(symbols, coords, charges,
                  qm_seeds=[0], qm_radius_angstrom=5.0)

# 3. Whole residues — a residue joins if ANY of its atoms is in range.
#    Prefer this for proteins: a radius cut alone will slice through a
#    residue and leave you with a chemically meaningless fragment.
ferric.QmmmSystem(symbols, coords, charges,
                  qm_seeds=[0], qm_radius_angstrom=5.0,
                  residue_ids=residue_ids)

residue_ids together with qm_indices is an error, not a silent no-op — an explicit index list and a residue expansion are contradictory instructions.

Cutting a covalent bond

If the QM/MM boundary crosses a bond, you need a link atom and a scheme for what to do with the host's charge:

sys = (ferric.QmmmSystem(symbols, coords, charges, qm_indices=qm)
       .with_link_atoms(bonds)                        # bonds crossing the cut
       .with_boundary_charges(bonds, "rcd"))          # keep | delete-host | rc | rcd

The parser is strict: an unknown string is an error, never a silent default.

Energies, forces, and optimization

res = ferric.run_qmmm(sys, "cc-pvdz", method="uhf")
# KS-DFT: method="rks" or "uks" AND xc, e.g. method="rks", xc="PBE".
# xc without rks/uks (or rks/uks without xc) raises ValueError.
res.energy
res.qm_gradient()      # dE/dR on the QM atoms
res.mm_forces()        # forces on the MM sites
res.full_gradient()    # across the cut, link rows folded back onto the frontier

opt = ferric.run_optimize_qmmm(sys, "cc-pvdz", method="rhf", move_mm="none")

move_mm chooses which MM atoms relax: "none" (the default, where only QM atoms move), ("within", radius_angstrom), ("residues", [ids]) (the system must have been built with residue_ids=), or "all". Any value other than "none" requires mm_topology=, because moving MM atoms without a force field to hold their own geometry together is meaningless.

Polarizable embedding

Thole polarizable sites, if fixed charges are not enough:

sys = ferric.QmmmSystem(symbols, coords, charges,
                        qm_indices=qm,
                        polarizabilities_angstrom3=alphas)
res = ferric.run_qmmm(sys, "cc-pvdz", method="rhf", thole_a=2.1304)
res.e_pol                  # the polarization energy
res.induced_dipoles()

What is validated, and how far

Checked against pyscf.qmmm.mm_charge: energy shift agrees to <1e-8 Ha, MM forces to 3e-10. An empty or all-zero MM region is bit-identical to a gas phase calculation — the trivial limit is a genuine no-op, not an approximation that happens to be small.

Not validated: periodic boundary conditions (absent), and the solvation droplet, which is a hard-sphere packing at roughly bulk density rather than an equilibrated box.

Reading a structure

Every format lands in the same place, so a PDB and an XYZ of one molecule give a bit-identical Molecule:

from tools.structure import read, from_smiles

mol = read("ligand.pdb")        # or .cif .pqr .sdf .mol2 .gro .xyz
mol = from_smiles("CCO")        # ETKDG geometry -- tier-2 grade, NOT optimized

A PQR carries MM charges as well as coordinates, which is why the CLI's [qmmm] section reads one. Note that a docked pose from Vina is united-atom — nonpolar hydrogens are merged into their carbons — so it is not a QM geometry until those are restored; tier1_dock does that for you.

From the CLI

A [qmmm] section runs QM/MM from a TOML file, with no Python needed. The QM region becomes the molecule that is solved, and the MM region becomes the external potential it is solved in. A complete, runnable input is examples/water-qmmm.toml (repository, not the wheel):

[molecule]
xyz = "testdata/molecules/water.xyz"   # required by the parser, NOT read with [qmmm]
charge = 0                             # apply to the QM region
multiplicity = 1

[basis]
name = "sto-3g"

[method]
kind = "rhf"
task = "energy"

[qmmm]
pqr = "testdata/molecules/water_na.pqr"
qm_indices = [0, 1, 2]            # zero-based; or: qm_seeds = [0], qm_radius_angstrom = 1.5
# link_bonds = [[0, 3]]           # required when the cut crosses a covalent bond
# boundary_scheme = "delete-host" # default; also "keep", "rc", "rcd"

MEASURED on that file (water QM, one Na⁺ 4 Å away as MM, STO-3G): vacuum −74.9629466809, embedded −74.9653197421 Ha (−1.489 kcal/mol). That matches ferric.run_rhf(point_charges=...) to all ten printed digits.

The geometry and the MM charges both come from the PQR. An xyz has no partial charges, and an MM region without charges is just a set of ignored coordinates. [molecule] xyz must still be present, because the parser requires it, but it is not read when [qmmm] is present. That matters when you compute a vacuum reference: delete [qmmm] and ferric falls back to the xyz. If the xyz holds a different geometry, you are comparing two different molecules. In the example above, the xyz is an optimized water that gives −74.9631468000, which is 0.13 kcal/mol of error in the difference. Compare at one geometry.

boundary_scheme defaults to "delete-host", not "keep". Keeping the host charge across a covalent cut puts a bare point charge inside the link atom's bond length (MEASURED 0.443 Å on an ethane C–C cut). A geometry optimization in that field diverges instead of failing with an error. "keep" is still available and is the right choice when the cut isn't covalent. An unknown scheme name is an error.

[qmmm] and [external_potential] together are refused. The MM region is an external potential, so combining them would silently count a contribution twice.

Scope of the CLI section. It reads only a PQR, because it needs charges and geometry together. Use Python (tools.structure, below) for the other formats. tools.active_site.solvate can write a solvated system straight to a PQR for this section. The CLI section does electrostatic embedding only. Polarizable sites, smeared charges, the MM force field and MM relaxation are Python-only.

Solvating a solute

from tools.active_site.solvate import solvate, write_pqr

drop = solvate(symbols, coords_angstrom, radius_angstrom=12.0)
write_pqr("solvated.pqr", symbols, coords_angstrom, charges, drop)

The solute is written first, so its indices are 0 .. n-1 and can go straight into [qmmm] qm_indices. Waters are TIP3P at roughly bulk density; this is a starting structure, not an equilibrated one, and a droplet has a vacuum boundary. Repeat over several seed= values before trusting a difference — dE_statistics does that and refuses fewer than two seeds.

Known limits, stated plainly

• No AMBER prmtop reader. Go through OpenMM (active_site.mm_topology.topology_from_openmm), or build the arrays.
• No periodic boundary conditions. solvate() gives a finite droplet with a vacuum boundary — adequate for a local environment, not for bulk.
• The pocket field is fixed unless you ask otherwise. move_mm="none" is the default in run_optimize_qmmm; the MM sites do not relax with the QM region until you widen it.
• No QM/MM dispersion. D3/D4/XDM/VV10 are all QM-atom-pairwise, so dispersion between the QM region and the MM charges is absent. The MM crate supplies Lennard-Jones terms for the MM-MM part only.
• The MM crate is AMBER-form (harmonic bonds and angles, periodic torsions, Lennard-Jones, Coulomb), validated against OpenMM.

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 has 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/, tools/, 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 ~30 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/, scripts/ or tools/.

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 — force fields are 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 single most common first failure; ferric catches it with a clear message, but checking in your own loop fails faster and names the species.

The parity gate is NOT a geometry check, and you need both. A structure can have a perfectly 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 sitting 1.030 A from a second carbon. The SCF stalled at 37 iterations. That is the correct behaviour for overlapping nuclei, and it cost a run to discover.

# 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

This is the middle tier of the funnel and skipping it is expensive. 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

LD_LIBRARY_PATH is not optional on a local install: bare xtb fails with libxtb.so.6: cannot open shared object file, which reads like a broken install and is not one.

MEASURED on the same molecule, same method, same 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

0.868 Ha = 545 kcal/mol lower, and the difference between a result and a failure. Seconds of xtb bought that.

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 100
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

Two things to take from that. 173 iterations is normal for a delocalised carbocation -- the default cap is not sized for this. And energy_conv: 1e-8 was unreachable: under density fitting the energy change floors with naux, so a tighter dE than the RI noise floor can never be met and the run burns its whole iteration budget chasing it. density_conv is the criterion that does the real work in this codebase; leave the dE bound loose enough to be satisfiable.

Note the failure MODE is diagnostic. Stalled means the SCF could not make progress -- suspect the geometry. MaxIter means it was progressing and ran out of room -- raise the cap. They are different problems and the exit reason names which one you have.

Set density_conv, not just 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 possible at all. 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. ESTIMATED: 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 is an outage, not a plan. 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's still cheap for an energy, but an optimization multiplies it by the step count. So the optimizations of the C10 species are the dominant cost of the whole 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 and do not pick the one that matches your hypothesis. A cascade whose ordering is functional-dependent is a finding about the system, not a number to be cleaned up.

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, not a cheap add-on to the DFT run. Budget it separately. 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 aren't hedging. They define the scope of the claim.

• Gas-phase cluster models. There's no enzyme environment unless you add QM/MM (golden path B, and QM/MM).
• Analytic Hessians for RHF and UHF only. RHF and UHF frequencies use the analytic Hessian; ROHF, KS-DFT and embedded 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, not 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.

Golden path B — Which residue should I mutate?

Answers: which active-site residues most influence a reactive center, so you have a ranked short list.

Does not answer: what mutation to make. ferric ranks hypotheses. It does not design, and it does not predict ΔΔG.

Step 1 — classical pre-screen (cheap, all residues)

derive_pocket_charges (in tools/active_site/pocket_charges.py) runs pdb2pqr on the pocket and records each charge's residue (residue_ids, res_names). pocket_field_at_atoms returns the potential and field [phi, Ex, Ey, Ez] (atomic units) at the sites you give it. It returns the total. To rank residues, group the charges by residue and evaluate each group separately. Recipes §5 has a sketch of the code.

Step 2 — QM/MM the top few (expensive, short list only)

Use QM/MM (ferric.QmmmSystem + ferric.run_qmmm). The embedding energy shift matches pyscf.qmmm.mm_charge to <1e-8 Ha. This step is what turns the answer into physics rather than electrostatics.

The pre-screen exists to make this step affordable. Running QM/MM on every residue is the thing the funnel pattern is designed to avoid.

Limits

Classical point charges: no protein polarization response, no sterics, no conformational change on mutation, no ΔΔG. A flagged residue is a candidate for QM/MM, not a designed mutation.

Golden path C — Screen many ligands down to DFT

Answers: of N candidates, which few deserve expensive QM.

Use tools.pipeline.run_funnel (repository tools/, not the wheel). Don't write your own loop. Recipes §4 has a funnel you can run (force field → xtb → DFT on three isomers, with its measured output). The pipeline notes §0b show the substituent version: parent-relative gating, liability flags and the optional docking tier.

Why use it instead of your own loop:

• Failed candidates are dropped and counted, never ranked. A hand-written screen that sorts ascending on a sentinel value silently promotes its failures to the top. The funnel exists to prevent that bug.
• Per-tier wall times. The tier that actually dominates is often not the one the cost table predicts. MEASURED in both recorded runs: the DFT tier took 90–96% of the wall time.
• Early stop on an empty population.
• Ascending rank at every tier, since each tier reports an energy-like score. Rank only candidates with the same formula, or rank relative to a parent.

The funnel will produce an ordering that the noise doesn't support. MEASURED on a real campaign: the best available ΔΔE noise over a pose ensemble was 4.07 kcal/mol, against substituent effects of 1–2. Read the pharma coverage notes before you rank anything.

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.

Toxicity screening

Structural-alert and predicted-liability readouts for a molecule, from the command line. This is a repository tool (tools/tox), not part of the wheel. Run it from a git clone of ferric with RDKit installed.

python -m tools.tox --offline "CC(=O)Oc1ccccc1C(=O)O"

The screen combines a local RDKit pass with two optional web providers. The local pass checks several hundred compiled SMARTS patterns from six published alert catalogs (Brenk, PAINS, NIH, and the Glaxo, Dundee and BMS sets via ChEMBL), plus Lipinski and Veber rules, and needs no network.

What the web providers return (probed 2026-09-23):

• ADMETlab 3.0's documented POST /api/admet returns HTTP 404. The client uses the live but undocumented POST /api/single/admet instead, one molecule per request. The host rate-limits hard, so expect HTTP 429, reported as an unavailable provider (exit 3), on larger batches.
• ProTox-3.0 has no documented JSON API, and by design the provider never scrapes the HTML results page. It never contributes endpoints: every online run lists it as unsupported, and it never affects the exit status.

For a batch, --offline is still the safer default. The local screen needs no network, and it is the part whose output does not depend on a third-party service being up.

Usage

python -m tools.tox [--offline | --require-online] [--fail-on-alerts]
                    [--timeout SECONDS] [--json] SMILES|FILE [SMILES|FILE ...]

An input is read as a file when it exists and ends in .smi, .smiles or .txt; otherwise it is treated as a literal SMILES. A file holds one <smiles> [label] per line, and # starts a comment.

# candidates.smi
CC(=O)Oc1ccccc1C(=O)O    aspirin
CN1C=NC2=C1C(=O)N(C)C(=O)N2C   caffeine

python -m tools.tox --offline candidates.smi

Exit status

When several apply, the precedence is 1 > 2 > 4 > 3 > 0.

These are distinct on purpose. "No alerts found" and "the alert screen did not run" produce similar-looking output, and they mean opposite things — so a provider failure is never folded into success. A web-service outage is neither a usage error nor a verdict, so it has its own status: 3 means the local results are usable and the online predictions are missing. If a pipeline needs the online predictions, pass --require-online. Without --fail-on-alerts, alerts are reported but do not change the exit status.

Each provider's outcome is in the JSON as provider_status: ok, unavailable (outage, HTTP error, rate limit, timeout, unparseable response), unsupported (ProTox), no_result, error or contract_violation. The human output shows the same tag in brackets on its !! lines.

Reading the output

Abridged, for aspirin with --offline (each line is followed by a one-line explanation, omitted here):

  alert_brenk                            0.3333 probability  [higher=worse] rdkit-alerts
  alert_pains                                 0 probability  [higher=worse] rdkit-alerts
  alert_total_count                           2 count        [higher=worse] rdkit-alerts
  desc_clogp                               1.31 log10        [higher=worse] rdkit-alerts
  desc_mw                                 180.2 Da           [higher=worse] rdkit-alerts
  lipinski_violation_fraction                 0 probability  [higher=worse] rdkit-alerts

Every line states its polarity, and the JSON carries the same flag as higher_is_worse. An aggregator that guesses the direction will invert a safety ranking, so read the flag rather than assuming one. ADMETlab's F20%/F30% columns are an example of how easy the guess is to get wrong: they are the probability of low oral bioavailability (below 20% or 30%), so they are reported as low_bioavailability_20pct and low_bioavailability_30pct, higher is worse. Reading them as bioavailability, higher is better, inverts them.

A value of None means unknown, never zero. For a probability-valued endpoint, 0.0 means "confidently predicted negative", which is the opposite of "no information".

What the alert scores are not

The alert_* endpoints are scaled hit counts (n/3, capped at 1.0). The output labels them, like lipinski_violation_fraction and veber_violation_fraction, with the unit probability. That label marks a 0–1, higher-is-worse scale, and it is what puts these endpoints into the mean that ToxAssessment.liability_score (the rank-only aggregate on the object assess_smiles returns) takes over every probability-unit endpoint; alert_total_count (unit count) stays out of it. The number itself is a rank-only liability density, not a probability of toxicity, and each alert line's explanation says so. A molecule with zero alerts is not thereby safe: danuglipron screens clean across all six catalogs and was discontinued for a liver signal. Structural alerts catch known problem substructures; they say nothing about dose, exposure or on-target pharmacology.

From Python

from tools.tox.assess import assess_smiles, assess_many

a = assess_smiles("CCO", include_web=False)
for e in a.endpoints:
    if e.known:
        print(e.name, e.value, e.units, "higher_is_worse" if e.higher_is_worse else "")

# A batch reuses ONE provider list: compiling the SMARTS catalogs dominates
# the runtime of a per-molecule loop.
results = assess_many({"aspirin": "CC(=O)Oc1ccccc1C(=O)O"}, include_web=False)

assess_many returns results in input order and does not rank them — ranking depends on what else is being traded off, so it is the caller's job.

ferric for agents

Read this first if you are an automated agent driving ferric. It's a routing page: what exists, where each thing is documented, and which mistakes cost the most time. The linked pages hold the details.

Start here

Don't build from source unless you need to. The prebuilt wheel installs in about a minute. A source build takes ~30 minutes, most of it libint2.

pip install ferric
git clone https://github.com/mgoldey/ferric && cd ferric   # for examples/, testdata/, tools/, scripts/
ferric examples/water-rhf.toml

Expect converged = true and energy = -74.9631468000 Hartree. If you don't get that, stop and fix the install (Installation). Nothing else will work, and the failures will be confusing.

pip install puts a ferric command on PATH that runs the same CLI as the source build. The wheel holds only the compiled library and that command. The examples/, testdata/, scripts/ and tools/ directories need the clone.

Which wheel

Use the nightly artifact when you need a fix that has landed on main but hasn't been tagged. Don't build from source just to get an unreleased fix.

Build from source when you are changing ferric or need MPI. MPI builds are source-only (see Installation).

Capability → entry point

Before you trust a number

Check the grade. Not every capability is equally validated. A method being available from the CLI doesn't mean its numbers are production-grade. Read its row in Capabilities and validation before quoting a value.

Check convergence, not the exit code. converged = true (CLI) or .converged (Python) is the signal. A run that hit max_iter still prints an energy. In the run log, read run_end.converged, and for correlated methods read the result record (Run logs).

When comparing to another code, check density fitting and the grid. Recipes §6 covers both.

Failure modes that cost the most time

1. Charge/multiplicity parity. An odd electron count needs an even multiplicity. The error message prints the arithmetic and the rule, and it means your .xyz or your charge is wrong. ferric can handle ions.
2. Debug instead of release. A debug binary is far slower. MEASURED: a 27-atom cation didn't converge a single point in 34 minutes on one. If a job seems hung, check which binary is running. cargo run --release --bin ferric and the wheel are both release builds. ./target/debug/ferric is not.
3. MPI is a workspace feature, not a CLI one. ferric-cli has no mpi feature of its own, so -p ferric-cli --features mpi fails. The working build is cargo build --release --workspace --features mpi.
4. Threading. Don't set OPENBLAS_NUM_THREADS above 1. The CLI and import ferric pin OpenBLAS to one thread when the variable is unset and honour an explicit value, and cargo sets it to 1 only when it is unset (an exported value wins). Multithreaded BLAS under ferric's rayon parallelism can crash (LU routines) or oversubscribe the machine. It's a correctness setting, not a performance tip.
5. Under-converged density. energy_conv alone doesn't converge the density. Correlation energies and properties inherit the full first-order density error. E_HF doesn't, because it's variational. Set density_conv.
6. Serialize QM jobs. Two concurrent ferric runs don't just halve throughput. MEASURED: a pair of them pushed a single point that normally takes seconds past a 15-minute timeout. Each run wants all the cores and several GB. Run one at a time, or use the funnel, which times each tier.

Memory: the budget predicts, the cgroup enforces

ferric works out a memory budget from [memory] budget_gb, FERRIC_MEM_BUDGET_GB, or, if neither is set, 0.8 × available RAM. The CLI prints that budget at startup and installs it as one shared, debited pool for the whole process: the three-index RI tensors, the DFT grid's AO cache and the large tensors of the MP2, RPA, GW and CC methods reserve their bytes from it, so two allocations alive at the same time cannot each claim the whole budget. From Python, memory_budget_gb= sets the same per-allocation limits but installs no shared pool, so each check compares its own allocation with the whole budget. In both, an allocation that does not fit is spilled to disk, recomputed on demand, or refused with an error naming it, depending on the allocation (see Sharp bits).

The budget isn't a cap on process memory (RSS). It covers the dominant tensors, not every allocation. Basis-sized matrices, integral engines (for example libint2's C++-side pools), BLAS and per-thread scratch and allocator overhead aren't charged, so treat the printed number as a floor on what the process needs. Lowering FERRIC_MEM_BUDGET_GB changes how the charged work is split into blocks. It doesn't shrink the uncharged part.

Size the budget from a gradient, not a single point. MEASURED on benzene/cc-pVDZ/PBE/RI-J: the SCF alone peaked at 0.616 GB, and the SCF plus gradient at 3.033 GB. The KS gradient holds the AO second derivatives on the grid. A budget sized from an energy run under-provisions a geometry optimization by roughly 4x, and the first optimization step is where it fails.

The cgroup is the only hard ceiling. On Linux with systemd:

scripts/ferric-limited --max=8G --high=7G -- ferric input.toml

An unbounded overshoot triggers the system-wide OOM killer, which can kill unrelated processes. Inside a cgroup, only the job dies. dmesg | grep -i oom confirms a kill. A truncated log doesn't prove one. As a rule of thumb, plan on ~6 GB for 27 atoms at def2-SVP.

Telling a running job from a dead one

• Look for evidence of work, not a process name. pgrep -f also matches your own grep, a bash -c wrapper, or a defunct entry at 0% CPU. Sort by CPU instead: ps -eo pid,stat,pcpu,rss,etimes,comm --sort=-pcpu | head.
• A log that has stopped growing doesn't mean the job is dead. stdout is buffered, so a long QM run can legitimately print nothing for tens of minutes. The run log (Run logs) is flushed record by record, so watch it instead.
• Terminated usually means timeout(1), not the OOM killer. Check whether a timeout wrapped the command before you look in dmesg.
• A grep filter that matches nothing looks exactly like "the job produced nothing". Test the filter against known-good output first.

Honest scope

ferric is a quantum chemistry engine: energies, gradients, finite-difference Hessians and properties. It isn't a protein-engineering tool, a docking program or an MD code. The tools/ layer connects it to RDKit, xtb and docking for screening, but the QM is the product.

Its analytic Hessian covers RHF and UHF only (other frequencies are finite differences of the gradient), and it has no thermochemistry (entropy, enthalpy, free energy). Transition-state search and IRC are Python-only and closed-shell only (Golden paths Step 5). If a task asks for something ferric can't do, such as designing mutations or predicting ΔΔG, say so instead of approximating it with something adjacent.

Methods overview

What exists, grouped by family. Maturity varies a great deal between these. Check the grade and its anchors on Capabilities and validation before trusting a number. To go from a task to a method, start with Choosing a method.

QM/MM embedding is on its own page: QM/MM.

Exact and local correlation

A correlated calculation is two choices: the method (what is computed: RI-MP2, dRPA, LinLCCD) and the approximation (whether, and how, its amplitudes are truncated). ferric keeps them apart. method.kind (and the Python function) names the method, and the method is computed exactly. A local approximation is asked for explicitly, with the CLI's [local] section or the Python local= and eps= keywords, and only rimp2, drpa and linlccd have one.

The one local scheme is the amplitude threshold: in the Boys-localized basis, a pair amplitude whose localized integral is at or below eps is dropped (single threshold, Wang et al. 2023). Three facts govern its use:

• eps is part of the model. It has no default and must be written, and every printout and run-log record of a local run carries it, with the fraction of amplitudes kept ("local": null marks an exact run). Report it with any number you quote.
• eps = 0 is the exact method. It keeps every amplitude and reproduces the exact energy (RI-MP2 to ≤ 1e-9 Ha; dRPA to the canonical plasmon formula; LinLCCD of the same variant), which is how each local path is anchored.
• The error is one-sided and grows faster than linearly in eps: dropped amplitudes give less correlation energy, measured with the same sign at every point, and each decade of eps raises the error by more than a decade (the fraction of amplitudes dropped grows with it). MP2 has Hylleraas stationarity; dRPA does not, so its error is first order in the dropped amplitudes and the one-sidedness is measured, not guaranteed.

Riccati, plasmon and PDEP are not approximations of each other: they are algorithms for the same exact dRPA. The drCCD Riccati solve at eps = 0, the plasmon formula and full-rank PDEP agree to ≤ 2.6e-14 Ha (H2/STO-3G and water/6-31G, PDEP at 64 quadrature points). Which one to use is a cost question: see RPA, GW and excited states.

Formal scaling with system size N, for orientation (these are the textbook exponents, not measured timings): SCF and DFT N⁴ for exact four-centre integrals, lower with density fitting and screening; RI-MP2 and its variants N⁵; CCSD N⁶; (T) N⁷.

Infrastructure

Shared machinery underneath all of the above:

• Screening: Schwarz bounds on every four-centre path, CSB/CSAM (Thompson & Ochsenfeld 2017) as options, LinK exchange (Ochsenfeld, White & Head-Gordon 1998). QQR (Maurer, Lambrecht & Ochsenfeld 2012) is implemented and validated as a bound but not used in production.
• Spherical and Cartesian basis support, with BSE-JSON and Gaussian-94 parsers and a set of bundled orbital, RI, JK and ECP bases.
• einsum!: a tensor-contraction macro that routes contractions through BLAS3 GEMMs, used throughout the CC and MP3 code.
• Memory budgets: allocation ceilings that refuse an oversized job with a named breakdown instead of letting it be OOM-killed.
• Python bindings (pyo3) and a TOML-driven CLI. Not every method has both; Capabilities and validation says which.

Properties

ESP at nuclei and on the molecular surface, electric field, static and atom-partitioned polarizabilities, Mulliken, Löwdin, Hirshfeld, CHELPG and RESP charges, density matrices, and NPZ export of ML-ready features.

Negative results

Known limits and measured negatives are kept in one place, Capabilities and validation: known limits, so a correction has one row to change. Two are worth knowing before you pick a method:

• TDHF/RPAx \( C_6 \) stays about 63% low regardless of gap. Use that kernel for static polarizabilities only.
• Local MP2: the integral-direct local MP2 (rimp2 with [local] integral_direct = true) is measured sub-quadratic (about N^1.24 erfc, N^1.4 Coulomb) on alkanes C20–C48; it is at about parity with exact RI-MP2 at C20 and about 6× faster at C32. Local MP2 without integral_direct makes no scaling claim. See The MP2 family.

The AO-sparse Laplace SOS-MP2 truncation is not a negative: the radius it needs does not track the molecular diameter. See the MP2 page.

SCF and DFT

Ground-state self-consistent field methods, their nuclear gradients, and the things built on them: geometry optimization, harmonic frequencies, transition-state search, implicit solvation and dispersion correction. The [scf] and [dft] keys are in Input file; grades are on Capabilities and validation.

Hartree–Fock

What it is. RHF (closed shell), UHF and ROHF (open shell), with DIIS, Schwarz screening, and a choice of direct, LinK, density-fitted (RI-J / RI-K) or seminumerical (COSX) Fock builds (see Choosing how exchange is built below). The open-shell solvers add per-spin DIIS, a virtual-block level shift, augmented-Hessian Newton, and Maximum-Overlap-Method (MOM) orbital tracking for near-degenerate cases (Gilbert, Besley & Gill 2008). These exist because specific systems failed without them.

Run it. method.kind = "rhf" / "uhf" / "rohf" (examples/water-rhf.toml, examples/h_uhf.toml); Python ferric.run_rhf, run_uhf, run_rohf. Spin and charge are set on the molecule (Molecule.from_xyz(path, charge, multiplicity)), not on the solver call.

Defaults. Exact four-index J and K (no density fitting) unless [scf] df_j_aux / df_k_aux are set. Convergence gates on the density: density_conv = 1e-6 by default; energy_conv (default 1e-3) is a sanity bound, not the convergence target.

Accuracy. Proven (rhf, uhf, rohf).

Kohn–Sham DFT

What it is. Kohn–Sham DFT through libxc: LDA, GGA, hybrid and range-separated hybrid functionals (for example PBE, B3LYP, ωB97X-V), VV10 nonlocal correlation, and the SCAN / r2SCAN meta-GGAs (τ-dependent; no density Laplacian). The solvers cover RKS, UKS and ROKS.

Run it. method.kind = "ksdft" with [dft] functional = "PBE" (examples/water-wb97xv.toml, examples/benzene-dfb3lyp.toml); Python ferric.run_dft or run_ksdft. The two Python entry points are closed shell (RKS) only. In the CLI, ksdft on a molecule with multiplicity > 1 runs UKS, and kind = "uhf"/"rohf" with [dft] functional run UKS/ROKS. Open-shell KS is also the reference inside pdep-rpa and gw ([rpa] xc with multiplicity > 1), and in Python through run_frequencies(reference="uhf", xc=...) and QM/MM (run_qmmm(method="uks")); see Capabilities and validation.

Defaults.

• Functional: LDA if [dft] functional is omitted.
• Grid: Becke partitioning of atom-centred Treutler–Ahlrichs radial × Lebedev angular grids, 75 × 110 points per atom, unpruned.
• Pruning is opt-in: [dft] grid_prune = "nwchem" (examples/water-pbe-pruned-grid.toml) removes about 23% of points at 75 × 110. It is accepted for task = "energy" only; optimization and frequency runs refuse a pruned grid because the gradient's grid-response term is built on the unpruned grid.
• Density fitting: RI-J and RI-K are on by default for ksdft and run_dft (def2-universal-jkfit), while rhf is exact by default. In Python, df_j_aux="" selects exact Coulomb. The fitting error grows with size: measured at PBE/STO-3G against exact J it was 0.28 kcal/mol for water, 1.16 for benzene and 9.5 for a 71-atom drug molecule. Compare with an exact-Coulomb code only after turning fitting off, or at matched fitting.
• Meta-GGAs get an automatic 0.5 Ha virtual-block level shift (ramped to zero at convergence) when you set none, because τ amplifies grid noise and plain DIIS limit-cycles.

Accuracy. ksdft is Proven. SCAN/r2SCAN energies are pinned against PySCF (ferric-scf/tests/dft_scan.rs); their closed-shell gradients against PySCF with grid response (dft_gradient_mgga.rs). Other meta-GGA variants (deorbitalized, +VV10, hybrid meta-GGA) are refused, not approximated.

Dispersion correction: D3(BJ)

What it is. Grimme's D3 with Becke–Johnson damping, a post-SCF pairwise correction, implemented natively in Rust. The two-body term only; the three-body (ATM) term is not implemented.

Run it. [dft] dispersion = "d3bj" on method.kind = "ksdft" (examples/water-pbe-d3bj.toml); "d3bj(b3lyp)" borrows another functional's fitted damping parameters. Python: run_dft(..., dispersion="d3bj"), or ferric.d3bj_energy(mol, functional) for the correction alone.

Scope. The damping parameters are fitted per functional, so the key is refused on any method other than ksdft. task = "optimize" works (the D3 gradient is implemented; ferric-d3/tests/gradient_vs_fd.rs), and so does closed-shell task = "frequencies": see Frequencies with dispersion.

Accuracy. Agrees with simple-dftd3 1.6.0 to below 1e-12 Ha on water, methane, benzene and an argon dimer with the three-body term off on both sides (ferric-d3/tests/vs_reference_dftd3.rs). The omitted three-body term is 0.1% of the two-body energy at benzene and grows with size.

Dispersion correction: MBD@rsSCS

What it is. The many-body dispersion energy of Ambrosetti et al. 2014 (range-separated self-consistent screening), added post-SCF. Its per-atom inputs are Tkatchenko–Scheffler free-atom α, C6 and R_vdW scaled by Hirshfeld volume ratios v_A / v_A^free of the converged SCF density; the free-atom volumes come from live free-atom SCFs in the same basis and SCF settings, solved for the isolated atom: point charges, external fields, implicit solvent, polarizable sites and cDFT constraints of the molecular run are not applied to it. The model itself is described under Polarizabilities and dispersion coefficients.

Run it. [dft] dispersion = "mbd" on a Kohn–Sham SCF (examples/water-pbe-mbd.toml) uses the β published for [dft] functional (PBE 0.83, PBE0 0.85, HSE06 0.85); "mbd(pbe0)" uses another functional's β. Any other functional is an error, not a default β. The printout gives E(KS-DFT), E(MBD@rsSCS) with β and the functional, and the corrected total; the JSON run log's run_end record carries the corrected total as energy, with scf_energy and a dispersion object (model, params, energy, beta, volume_ratios) in extra (D3(BJ) runs log the same object without beta and volume_ratios). UKS/ROKS energy runs, which write no run_end record, write a dispersion record with the same fields. Python: run_dft(..., dispersion="mbd"), which also reports DftResult.volume_ratios.

Gradient. task = "optimize" on an RKS, UKS or ROKS reference and run_dft(with_gradient=True) add the exact analytic MBD@rsSCS gradient: the explicit dependence on the nuclear positions, and the dependence through the Hirshfeld volumes, including how the SCF density itself responds to the displacement. That response is the orbital relaxation, obtained from one coupled-perturbed Kohn–Sham (Z-vector) solve per gradient; it costs about as much as a few SCF iterations. On a UKS reference the solve is coupled across the α and β orbital rotations (Coulomb couples the spins; exchange and the XC kernel are spin-resolved), and the orthonormality of each spin's occupied orbitals enters separately. On a ROKS reference the solve runs over the three rotation blocks of the shared orbitals (closed→virtual, open→virtual, closed→open) with the exact restricted-open-shell orbital Hessian, and the closed and open orbitals are kept orthonormal as one set, which adds a closed–open cross term to the orthonormality contribution. The volumes are integrated on a lattice anchored to the molecular centroid, and its motion is part of the gradient. Against finite differences of the full SCF + MBD calculation at 6-31G the gradient agrees to 6e-9 Hartree/Bohr for H2O with PBE, PBE0, HSE06 and PBE with RI-J, and to 9e-10 for NH3; the orbital relaxation alone is 1.0e-5 for H2O (11.5% of the largest MBD component). For UKS doublets and a triplet (NH2, OH, O2 at 6-31G with PBE, PBE0 and HSE06) it agrees to ≤ 1.9e-9, and to 6e-9 for OH with PBE + RI-J and PBE0 + RI-JK, against an orbital relaxation of 3e-6–1e-5. For ROKS doublets and triplets (HCO, NH2, CH2, O2 at 6-31G) it agrees to ≤ 2.1e-11 with PBE and PBE0, ≤ 4.5e-10 with HSE06 and 7.5e-9 for HCO with PBE + RI-J, against an orbital relaxation of 2.7e-6–9.8e-6. The relaxation term needs an LDA, GGA or hybrid-GGA functional without VV10, no implicit solvation, polarizable embedding or cDFT constraints, and integer aufbau occupation (no MOM); UKS and ROKS additionally need exchange that is not COSX. Other setups are refused rather than given an approximate gradient. Open-shell frequencies with dispersion are refused. Python's run_dft is closed-shell only.

Frequencies with dispersion

What it is. Harmonic frequencies on the dispersion-corrected surface E(KS) + E(disp), closed-shell Kohn–Sham only. The Hessian is the central difference of the corrected analytic gradient: at each of the 6N displaced geometries the SCF is converged, the KS gradient and the dispersion gradient are evaluated there, and their sum is differenced. For MBD@rsSCS the dispersion gradient is the exact one above, so the response of the Hirshfeld volumes to the displacement, through the SCF density, is in the Hessian. There is no analytic Hessian on this path: [frequencies] hessian = "analytic" is an error and "auto" runs finite differences.

Run it. [dft] dispersion with method.task = "frequencies" on a closed-shell KS SCF. The printout gives the corrected energy, E(KS-DFT) and the dispersion energy at the input geometry, and the JSON run log gets a dispersion record with the same fields as an energy run. Python: run_frequencies(mol, basis, xc="PBE", dispersion="d3bj"), which reports .e_dispersion and the corrected .energy. Open-shell (UKS/ROKS) frequencies with dispersion are refused.

Accuracy. For PBE/STO-3G water, the D3(BJ) part of the Hessian agrees with 4-point second differences of the D3(BJ) energy to 5.3e-9 Hartree/Bohr² (largest element 1.9e-5), and the resulting frequencies with those of the KS Hessian plus that independent D3 Hessian to 1e-5 cm⁻¹. The MBD@rsSCS part agrees along three fixed directions with second differences of the full SCF + MBD energy to 1.5e-6 Hartree/Bohr² (0.8% of the largest curvature, 2.0e-4). That figure is set by noise in the energy differences, not by the Hessian construction. A gradient without the orbital-relaxation term misses by 3.2e-5. Dispersion leaves the Hessian's asymmetry (7.6e-6 Hartree/Bohr² for PBE/6-31G water) and the projected translation/rotation modes (below 1e-4 cm⁻¹) where they are without it.

The shifts are small for a single molecule. For PBE/6-31G water at its uncorrected PBE minimum (1594.35, 3500.27 and 3669.29 cm⁻¹), D3(BJ) shifts the three modes by +0.006, −0.077 and −0.095 cm⁻¹, and MBD@rsSCS by +0.107, −0.406 and −0.516 cm⁻¹.

Implicit solvation

Two independent implementations. Both are threaded uniformly through RHF, UHF, ROHF and the KS variants, and both leave the energy bit-identical to vacuum when unset.

Neither has cavitation, dispersion or repulsion terms, and neither has an analytic gradient. An unknown solvent name is an error, not a silent vacuum.

Gradients, optimization, frequencies, transition states

Gradients. Analytic nuclear gradients for RHF, UHF, ROHF and KS-DFT (including meta-GGA, with grid response), checked against finite differences and, for DFT, against PySCF. RI-MP2 also has an analytic gradient, closed shell only (there is no unrestricted MP2 nuclear gradient).

Optimization. method.task = "optimize" for rhf, uhf, rohf, ksdft, rimp2 and pdep-rpa (examples/h2_opt.toml, examples/h2-lda-opt.toml); on an open-shell molecule only uhf, rohf and ksdft. Python ferric.run_optimize is RHF.

Harmonic frequencies. method.task = "frequencies" for rhf, uhf, rohf and ksdft (examples/water-frequencies.toml); Python ferric.run_frequencies(mol, basis_name, reference="rhf", xc=...). Mass-weighted, translations and rotations projected out. The Hessian is:

• Analytic for closed-shell RHF, and for UHF of any multiplicity, with exact four-centre J/K, no ECP, no external potential or solvent, and a basis up to f functions: one SCF plus a coupled-perturbed HF solve per nuclear coordinate (for UHF, one solve coupling the α and β orbital rotations). This needs a libint2 with second derivatives (the build scripts/install-libint.sh installs).
• Central finite differences of the analytic gradient everywhere else (ROHF, KS-DFT, RI J/K, ECPs, embedding, g functions): 6N gradient evaluations. The displacement [frequencies] delta (default 5e-3 Bohr) is a real accuracy knob; the printed Hessian asymmetry, zero in exact arithmetic, is the check that it and the SCF thresholds suit the system.

[frequencies] hessian (Python hessian=) selects it: "auto" (default, analytic where it applies), "analytic" (an error where it does not) or "fd". The output names the one that ran (Hessian = analytic; Python .hessian_source).

Transition states and IRC (Python only, closed shell). ferric.run_saddle(mol, basis_name, xc=...) searches for a first-order saddle by P-RFO, using two finite-difference Hessians with Bofill updates between them; it costs 2(6N + 1) + (n_steps + 1) gradient evaluations and refuses to start from a geometry with no negative mode. result.is_transition_state() requires both convergence and exactly one imaginary frequency. ferric.run_irc(mol, basis_name, result.imaginary_mode) then follows the reaction path both ways to show which minima the saddle connects (measured about 71 gradients per branch on NH3 inversion). Both refuse multiplicity ≠ 1.

Choosing how exchange is built

Four ways to build K, and the choice is a real one, measured on this code. Butane, one thread; TZ is def2-TZVP (184 functions), QZ is def2-QZVP (528).

Time for one exchange build (seconds, butane, one thread):

Time for a full SCF (seconds, butane, one thread):

Energy error against exact exchange:

Start with density fitting whenever its three-index tensor fits in memory (n_aux × n_bf² × 8 bytes: 1 GB for butane/QZVP, 7 GB for octane/QZVP). It spills to disk when it does not. It is two to three orders of magnitude faster per build than anything else here. Its error is a fitting error, not zero, and it grows with system size.

Need exact exchange? Direct or LinK; both are exact to the screening threshold as K builders (butane/def2-SVP: link == direct to 9e-12 Ha). LinK's cost is being re-measured, so no LinK timing is quoted here. k_builder is honoured by UHF and ROHF as well as RHF: the open-shell solvers build K_α and K_β from one builder instance, refreshing its density-dependent state per spin. Whether LinK is the faster choice for a given system is a separate question from whether it is honoured — see the cost note above, which is being re-measured. link is skipped with a warning, never silently, whenever density-fitted J/K is active; both link and cosx are skipped with a warning when the functional uses no exact exchange or is range-separated (exchange then comes from the SR/LR fitters).

RIJCOSX. k_builder = "cosx" combines with density-fitted Coulomb: when RI-J is active (df_j_aux named, or the Kohn–Sham default), J comes from RI-J and K from COSX, as in ORCA's RIJCOSX and Psi4's DFDIRJ+COSX. COSX then replaces RI-K: the Kohn–Sham RI-K default is not applied, and an explicitly named df_k_aux next to k_builder = "cosx" is an error. df_j_aux = "" keeps exact J with COSX K. The two approximations add: on water/6-31G (HF) the COSX error is 7.140e-6 Ha with exact J and 7.141e-6 Ha with RI-J, and the cross term is second order (2.8–5.8 × the product of the two errors; 7.8e-10 Ha at the (50,110) grid, 2.3e-11 Ha at (75,302)). Gradients take J from the RI-J derivative and K from the COSX derivative (FD agreement 1.7e-9 Ha/Bohr RHF, 3.6e-9 UHF).

COSX is for large basis sets on systems too big for RI-JK. Its cost per grid point barely moves with angular momentum while analytic exchange grows roughly tenfold from SVP to QZVP, so it wins at high angular momentum, not at large system size. Measured on one thread at the flat (50,110) grid:

• On butane, against exact direct exchange: at def2-TZVP the full COSX SCF is 3.7× slower (358 s vs 98 s); at def2-QZVP the COSX K build takes 90 s, against about 400 s for the direct build's single J+K sweep, at a relative K error of 5.4e-5.
• On n-alkanes at def2-SVP, against LinK: slower at every size measured, 1.59× (C20), 1.09× (C32) and 1.22× (C48), with no trend toward parity. At def2-TZVP on C20 it is faster (0.67×).

Below quadruple zeta it is the wrong tool unless RI-JK's three-index tensor does not fit in memory.

A density-driven pair screen (on the product of the integral bound and the local half-transformed density) keeps the K error below 2e-6 Ha at the default threshold, and the half-transforms D·X and X·Gᵀ run over per-batch sparse AO lists (cosx_half_transform = "sparse", the default). At def2-SVP the K build grows as N^1.29–N^1.32 between C20 and C48; that is one family of molecules in one basis.

On the flat (50,110) grid COSX's energy error is about 5e-6 Ha on water/cc-pVDZ (4.9e-6 in cosx_k_anchors.rs, 5.1e-6 in the ORCA comparison run), 1.7e-4 Ha on butane/def2-SVP and 1.2e-4 Ha on butane/def2-TZVP (against exact exchange). ORCA 6.1.1's COSX at its own default grid gives 4.9e-6, 3.8e-5 and 1.3e-5 Ha on the same systems and bases with about half as many points: ORCA's grids are pruned around a 194-point valence shell and it evaluates its final energy once on a finer grid. ferric's default does the same (cosx_grid pruning and cosx_final_pass, below). Refining the flat grid converges water to 3.4e-8 Ha, but butane/def2-SVP stays at 3.4e-5 Ha at Lebedev-302 for every radial grid. That residual is angular: an independent COSX (PySCF SGX, 75 radial shells) on the same system goes from 5.5e-5 Ha at 302 points per shell to -8.5e-6 at 434 and 2.7e-6 at 590. ferric's COSX follows the same path at 75 radial shells: -3.4e-5 at 302, -5.6e-6 at 434 and +1.8e-6 at 590, for 1.8x and 2.2x the 302-point wall time. Reaction energies cancel most of the error (0.02 kcal/mol on an isodesmic alkane reaction at the flat (50,110) grid); absolute energies do not. Four knobs, all optional:

• cosx_grid = { radial = 35, angular = 194, prune = "sgx" } is the default; { radial = 50, angular = 110 } is the flat grid. The angular order matters most: on butane/def2-TZVP the error falls from 1.2e-4 Ha at 110 points per shell to 4.0e-6 Ha at 302, while going from 50 to 100 radial shells changes it by 1e-6 or less. Cost grows with the number of points. angular must be one of 6, 14, 26, 50, 110, 194, 302, 434 or 590. prune = "sgx" prunes the grid the way ORCA's GridX and PySCF's SGX do: five radial regions per atom (NWChem boundaries on Bragg radii) with Lebedev orders from one row, picked by the peak angular — at 194 the regions get 26/50/110/194/110 points. ferric's pruned COSX energy equals PySCF SGX on the same grid to 1.8e-12 Ha. A table without prune is flat.
• cosx_final_pass = true (the default; false turns it off) re-evaluates exchange once on a larger grid at the converged density and reports that energy; the SCF-grid energy is printed and logged next to it. cosx_final_grid picks the grid (default { radial = 50, angular = 302, prune = "sgx" }). The pass is not self-consistent, but it lands within 1.5e-7 Ha of an SCF converged on the final grid (table below). It is energy-only: gradients and geometry tasks run without it.
• cosx_overlap_fit = true (default) applies the Izsák–Neese overlap correction. At the default grid it helps; on coarser grids it makes things worse, and its benefit is strongly molecule-dependent — large on water, nil to negative on ethane — so do not expect the factor quoted in the literature.
• cosx_backend = "md3c1e" (default) is the batched McMurchie–Davidson integral kernel. "cosx-a" is the per-point libint2 path, about three times slower and kept only as the cross-check the kernel is anchored against.
• cosx_screen_thresh = 1e-7 (default) is the density-driven pair-screening threshold. 0.0 disables screening bit-identically; 1e-6 already fails a 1e-6 Ha K-error bar on butane. Not available with the "cosx-a" backend.

Setting any cosx_* key without k_builder = "cosx", or k_builder = "cosx" together with a named df_k_aux, is an error.

Which COSX grid. Measured against exact exchange (exact J on both sides, RHF, overlap fit on unless noted), one run each, six threads; times are whole SCFs, indicative only. Errors, E − E_exact in Ha:

Whole-SCF wall time relative to exact-K RHF on the same system (COSX is slower than exact exchange at all of these sizes and bases; see above for where it wins):

What the table shows:

• The pruned sgx (35,194) grid uses 0.65× the points of flat (50,110) and is more accurate on seven of the eight molecules, by 1.7× (benzene) to 29× (ethane). On methane it is 1.15× worse (9.2e-6 against 8.0e-6 Ha). This is the grid geometry tasks (optimize, frequencies) use, since they run without the final pass.
• The final pass on sgx (50,302) reproduces an SCF converged on that grid to 4e-8 Ha on seven molecules and 1.5e-7 Ha on sulfamethoxazole, at 0.45–0.59× of that SCF's cost. Together with the pruned SCF grid it is more accurate than flat (50,110) on all eight molecules and on the reaction (0.003 against 0.074 kcal/mol), and faster on all eight (0.80–0.91× flat (50,110)'s wall time). This pair is the default for energies; ROHF/ROKS has no final pass and skips it with a note.
• More radial shells (35 → 50) at a 194 peak buy little; removing the pruning (flat 194) costs 2.7× the points for errors of the same size. The overlap fit helps on every molecule but water (1.8× to 21×; on water it is 1.2× worse).

COSX gradients.

task = "optimize" and task = "frequencies" with k_builder = "cosx" differentiate the COSX energy itself (grid-function, ESP-integral and Becke-weight derivatives). With the default overlap fit the fitted exchange is not variational in the orbitals, so the gradient adds an orbital-response (Z-vector) term. The gradient is exact for RHF, RKS and UHF with cosx_overlap_fit = false, and for RHF and UHF with the default fit; measured against finite differences of the COSX energy it agrees to 2e-9–4e-9 Ha/Bohr, on flat and pruned (sgx) grids alike, and with RI-J (RIJCOSX). The gradient differentiates the SCF-grid energy: with cosx_final_pass the reported energy is the final-grid one, so geometry tasks run without the pass and a gradient of a final-pass result is refused (an ORCA-style gradient evaluated on the final grid misses finite differences of the final-grid energy by 4.6e-7–7.3e-7 Ha/Bohr on water/6-31G, so it is not used). Fitted COSX with a Kohn–Sham functional, UKS and ROHF/ROKS are refused for gradient tasks before the SCF runs.

Convergence

A few things worth knowing before debugging a stubborn SCF:

Check converged. These routines return a result whether or not they converged; a non-converged SCF is a result with converged = false, not an error. Downstream code that ignores the flag will happily consume a half-converged density.

Multiple solutions are real. For systems like alkane chains, different initial guesses converge to genuinely different SCF solutions — not a convergence failure but a different basin. The guess picks the basin.

Density-fitting has a noise floor. DF-JK introduces an error floor that makes energy-based convergence criteria below roughly 1e-9 meaningless; SCF gates on the density RMS change instead.

Near-linear-dependence. Diffuse (aug-) basis sets on close-packed systems can drive the overlap matrix near-singular; the canonical-orthogonalization threshold is tunable via FERRIC_LINDEP_THRESH.

Screening

• Schwarz bounds on every 4-centre path, built so they can never underestimate (zero-valued table entries are floored; left unfloored, one such entry was measured to cost 1.5e-4 Ha)
• LinK (Ochsenfeld, White & Head-Gordon 1998) — exchange via significant-pair and density-pair lists; its scaling is being re-measured
• QQR (Maurer, Lambrecht & Ochsenfeld 2012) is implemented and validated as a bound but is not used in production: on LinK it screened only 0.009% more quartets than Schwarz at alkane_16 (measured against a LinK pair-list implementation that skipped quartets; not yet repeated on the current lists)
• COSX shell-pair screening uses a primitive-level Hölder bound that provably never underestimates; an overlap-based bound can underestimate, which silently corrupts K

QM/MM embedding

QM/MM (electrostatic, polarizable and Gaussian-smeared embedding, link atoms, boundary-charge schemes, and a [qmmm] CLI section that reads a PQR) has its own page: QM/MM, including its CLI section.

Determinism

The Fock build's reduction folds partial matrices in a strict ascending group order, independent of thread count and of the memory band width. Results are bit-identical across RAYON_NUM_THREADS — a property pinned by tests, not just intended.

This matters more than it might seem: a tree-fold reduction would be equally deterministic but would produce different bits, since floating-point addition is not associative. The ascending order is load-bearing.

Cite

DIIS: Pulay 1980. MOM: Gilbert, Besley & Gill 2008. LinK: Ochsenfeld, White & Head-Gordon 1998. COSX: Neese et al. 2009; overlap fit: Izsák & Neese 2011. CSB/CSAM screening: Thompson & Ochsenfeld 2017. DFT grids: Becke 1988, Treutler & Ahlrichs 1995, Lebedev & Laikov 1999. Functionals through libxc (Lehtola et al. 2018): PBE, B3LYP, ωB97X-V, SCAN, r2SCAN, VV10. D3(BJ): Grimme et al. 2010, Grimme, Ehrlich & Goerigk 2011. IEF-PCM: Cancès, Mennucci & Tomasi 1997; COSMO: Klamt & Schüürmann 1993. P-RFO: Banerjee et al. 1985; Bofill 1994. Full entries in References.

The MP2 family

The largest method family here. Every variant is density-fitted (RI): each needs an orbital basis and a matching RI auxiliary basis ([mp2] auxbasis, or the auxbasis argument in Python). Formal cost is O(N⁵) for the RI-MP2 transformation. Grades per method.kind are on Capabilities and validation. The [mp2] keys are in Input file.

RI-MP2

What it is. Second-order Møller–Plesset theory on an RHF reference (UHF for an open shell, below) with 3-centre/2-centre density fitting. Canonical (non-RI) MP2 is also implemented, for cross-validation, not production.

Open shells. rimp2 and oo-rimp2 accept an open-shell molecule (multiplicity > 1) from the CLI, for task = "energy" only: they solve the same plain UHF that kind = "uhf" runs, then take unrestricted RI-MP2 (UMP2, as PySCF mp.MP2(uhf)) or unrestricted OO-RI-MP2. There is no unrestricted MP2 nuclear gradient, and [mp2] kappa is refused on an open shell. From Python, run_rimp2 does the same (UHF + UMP2 when multiplicity > 1; result.reference says which, and kappa raises); run_oo_rimp2 is closed shell only. Both unrestricted methods are validated on the CH3 doublet at cc-pVDZ: U-RI-MP2 against ORCA RI-MP2 NoRI, U-OO-RI-MP2 against an independent numpy OO-RI-MP2 (see the anchors).

The other MP2-family kinds (mp3, att-rimp2, scs-mp2, scs-mp2-2terfc, laplace-mp2, laplace-sos-mp2, rs-mp2-rpa) and their Python drivers are closed shell only: set multiplicity = 1. mp2-v switches to a UHF reference when multiplicity > 1 (CLI, energy only); Python run_mp2_v is closed shell only.

Run it. method.kind = "rimp2" (examples/water-rimp2.toml, examples/water-rimp2-frozen-core.toml); Python ferric.run_rimp2(mol, bs, aux).

Accuracy. Proven. RI-MP2 is size-extensive to 2e-12 Ha for a well-separated dimer against twice the monomer (ferric-mp2/tests/rimp2_size_extensivity.rs, which asserts 1e-7).

Knobs. frozen_core defaults to 0 (all electrons correlated). Several published parameterizations below assume frozen core, so set it when you use them.

Attenuated MP2

What it is. MP2 with the Coulomb operator in the correlation energy replaced by a short-range one. Two forms are implemented:

• erfc: \( \mathrm{erfc}(\omega r)/r \), which is \( 1/r \) at short range and decays to zero beyond roughly \( 1/\omega \).
• terfc: a short-range operator that is Coulombic inside a cutoff radius \( r_0 \) and attenuated outside it, with a controlled curvature (Dutoi & Head-Gordon 2008). It needs precomputed interpolation tables (FERRIC_TERF_TABLE_DIR, generated under terf-tables/); without them the run errors rather than substituting another operator.

Why. MP2's error for non-covalent interactions has two parts that partly cancel in small basis sets: basis-set superposition error, and an error in the long-range part of the correlation energy. Attenuation removes the long-range part. That also means attenuated MP2 has no long-range dispersion at all: its asymptotic \( C_6 \) is zero. The method works in the basis and on the systems it was fitted for; MP2-V (below) adds long-range dispersion back through VV10.

Parameterization. The attenuation parameters are fitted, not derived, and each fit belongs to one basis. The aug-cc-pVTZ sets that ferric ships as defaults were fitted on S66 without counterpoise and with frozen core (stated in crates/ferric-mp2/src/scs.rs):

Running a fitted variant in another basis, with counterpoise, or with frozen_core = 0 (the default) is extrapolation outside the fit.

Run it.

ω is in Å⁻¹ in the CLI and Python and in Bohr⁻¹ in the Rust API. The SCS-MP2 defaults are Grimme's c_OS = 1.2, c_SS = 0.333.

Accuracy. att-rimp2, scs-mp2 and scs-mp2-2terfc are Proven (the erfc energy is pinned against PySCF/libcint in testdata/reference/h2o_cc-pvdz_attenuated-rimp2-erfc0p420.json). mp2-v is Smoke: its VV10 half is bit-identical to the ωB97X-V code path, but there is no published MP2-V total energy to compare with, and ferric evaluates the VV10 term on the converged HF density (the post-HF variant) on its own 50×50 grid rather than SG-1. Read the header of examples/water-mp2v.toml before quoting a number.

Robust fitting. When the RI metric differs from the operator in the integrals, as it does for attenuated operators, robust (Dunlap) density fitting is required, not optional.

RS-MP2 + LR-RPA

What it is. Short-range MP2 plus long-range direct RPA: the attenuated-MP2 idea with the missing long-range correlation supplied by response rather than by a fitted dispersion term. Two formulations:

• delta-lr (default): \( E_{MP2}[\text{Coulomb}] + (E_{dRPA}[\text{erf}] - 2E_{OS}[\text{erf}]) \), one dRPA call.
• coupled-rings: \( E_{MP2}[\text{Coulomb}] + \Delta dRPA[\text{Coulomb}] - \Delta dRPA[\text{erfc}] \), two dRPA calls; includes the mixed short/long-range ring diagrams.

Run it. method.kind = "rs-mp2-rpa" (examples/water-rs-mp2-rpa.toml, [mp2] formulation); Python run_rs_mp2_rpa(..., omega=0.420, formulation="delta-lr").

Accuracy. Smoke. The ω→0 and ω→∞ limits reduce exactly to MP2 and to MP2 + dRPA; numbers at production ω are not established on new systems.

Orbital-optimized MP2 and MP3

OO-RI-MP2 (oo-rimp2, examples/water-oo-rimp2.toml, run_oo_rimp2): orbitals optimized for the MP2 Lagrangian, with a level-shifted Newton step, orbital DIIS, Cayley rotations and backtracking. Proven (narrow): the energy matches an independent numpy OO-RI-MP2 to 7.5e-13 Ha (H2O, NH3, UHF CH3 at cc-pVDZ), and the closed-shell analytic gradient matches a finite difference of its own energy to 8e-9 Ha/Bohr. ORCA 6.1.1's OO-RI-MP2 stops 3.7e-8 to 7.0e-8 Ha above the same minimum, so it is only a loose cross-check.

MP3 (mp3, examples/water-mp3.toml, run_mp3): spin-orbital third-order Møller–Plesset through the einsum! framework. Proven.

Laplace formulations

RI-Laplace MP2 (laplace-mp2, examples/water-laplace-rimp2.toml, run_laplace_mp2): MP2 through a Laplace transform of the energy denominator, using pseudo-density matrices in the AO basis. The implementation is dense. It is the correctness reference for that formulation, not a reduced-scaling path, and no reduced scaling has been measured.

Laplace SOS-MP2 (laplace-sos-mp2, examples/water-laplace-sos-mp2.toml, run_laplace_sos_mp2): opposite-spin-only MP2 (Jung et al. 2004, c_OS = 1.3 by default) with a minimax Laplace quadrature (n_quad must be 3, 5 or 7; anything else is an error). There is deliberately no c_SS: dropping the same-spin term is what lets the denominator factorize. With c_OS = 1.0 it reproduces the opposite-spin component of RI-MP2 to quadrature error, which is the test anchor. Three formulations (sos_formulation):

• mo (default) and ao are both exact and agree to round-off.
• ao-sparse restricts each Boys-localized orbital's pseudo-density to an AO domain of radius domain_cutoff_bohr (required, no default). It is the one approximate variant.

What is measured for ao-sparse:

• Against the exact AO path on n-alkanes (cc-pVDZ, c_OS = 1, n_quad = 7, 2026-07-28), chemical accuracy (absolute error below 1.6 mHa) needs a domain radius of 3, 3, 3, 4, 5 and 5 Bohr for C2, C4, C6, C8, C10 and C12. The diameter grows fivefold over that series, so radius/diameter falls from 0.52 to 0.17.
• In the STO-3G tests, a 12 Bohr domain that is exact for butane (10.5 Bohr across) is also exact for octane (19.9 Bohr across), and a 4 Bohr domain on butane is already within 0.1%.
• A 71-atom drug molecule (danuglipron, 31.3 Bohr across, STO-3G) is within 0.05% at 4 Bohr.

The test sos_ao_sparse_truncation_radius_is_transferable_across_sizes in crates/ferric-mp2/src/laplace.rs pins the STO-3G butane/octane comparison: at 12 Bohr butane is exact (relative error below 1e-9) and octane is within 1e-6, and octane at 3 Bohr is worse than at 12 Bohr. The 4 Bohr butane figure and the danuglipron run (recorded in the test's doc comment) are measurements, not assertions. The alkane sweep is kept in the project's working notes.

How to read it (provisional): the radius needed grows, but far more slowly than the molecule. That points to a finite decay length rather than strict saturation: no single radius is shown to suffice at every size.

What is not claimed: any speedup. The domains discard contributions but the tensor algebra is still dense, so there are no timings to report.

Exact and local MP2

method.kind = "rimp2" is the method; it is computed exactly unless a [local] section asks for a local approximation (see Exact and local correlation). The local approximation is closed shell and energy only.

Amplitude-threshold local MP2 (examples/water-rimp2-local.toml; [local] scheme = "amplitude-threshold", eps = 1e-4; Python run_rimp2(..., local="amplitude-threshold", eps=1e-4)): the single-threshold local MP2 of Wang, Aldossary, Shi, Liu, Li & Head-Gordon (2023), with localized virtuals and per-pair domain-local RI fits. eps has no default: it is part of the model. eps = 0 reproduces RI-MP2 exactly; a finite eps carries a one-sided truncation error that grows faster than linearly in eps. Every printout and run-log record of a local run states eps and the fraction of amplitudes kept.

The exact RI-MP2 reference that measures that error is opt-in. It is a full canonical RI-MP2 and forms the global (naux, nocc·nvir) tensor, so a run with it switched on is not reduced-cost. In the CLI, [local] reference = true computes it and prints the local error against it; without it the output reads E_corr(canonical RI) = not computed (opt-in: set [local] reference = true) and the run log's e_corr_canonical_ri is null. In Python, compute_reference=True computes it; by default result.local["e_corr_canonical_ri"] is None.

Integral-direct local MP2 ([local] integral_direct = true, examples/alkane8-rimp2-local-direct.toml; Python integral_direct=True) is the reduced-cost path. Its correlation assembly never forms the global 3-index tensor; only the opt-in reference does. Locality comes from an integral-free pair gate (gate_cal), per-occupied auxiliary-fit and virtual domains (aux_radius, virt_radius), and truncation of each orbital's AO support (ao_tail). With every map at its trivial setting it reproduces the global 3-index path and canonical RI-MP2 (tests/lmp2_direct.rs).

What is measured for the integral-direct local MP2 (n-alkanes C20 → C48, 6-31G with cc-pVDZ-RI, a quiet machine, 2026-09-07; benchmark bench_direct_alkane_series in crates/ferric-mp2/tests/lmp2_direct.rs). The run froze the carbon cores (frozen_core = number of carbons) and calibrated the pair gate (0.7 Coulomb, 0.02 erfc with ω = 1.0) at eps = 1e-4; the defaults are all-electron with no pair gate. Timings are the correlation stage (assembly + solve) and exclude the canonical reference:

• Correlation-stage cost grows as about N^1.24 with the erfc kernel and N^1.4 with Coulomb, fitted to the last three sizes.
• It is at about parity with canonical RI-MP2 at C20 (about 1.1× slower) and about 6× faster at C32 (5.7–6.3×).

How to read it (provisional): sub-quadratic on this series, but the fit has three points on one family of molecules in one basis. It is not shown to be linear, and it is not measured on 3-D or diffuse systems.

Local MP2 without integral_direct makes no scaling claim. Its assembly is pair-local (no dense J is formed), but it still builds the global 3-index tensor. Use it as the reference implementation and for small systems.

Cite

RI-MP2 auxiliary sets: Weigend et al. 1998. SCS-MP2: Grimme 2003. SOS-MP2: Jung et al. 2004. OO-MP2: Lochan & Head-Gordon 2007; Bozkaya et al. 2011. Laplace MP2: Häser & Almlöf 1992. Attenuated MP2: Goldey & Head-Gordon 2012; terfc: Dutoi & Head-Gordon 2008, Goldey, Dutoi & Head-Gordon 2013; SCS-MP2(2terfc): Goldey & Head-Gordon 2014; MP2-V: Goldey, Belzunces & Head-Gordon 2015. LMP2: Wang et al. 2023. Robust fitting: Dunlap 2000. Full entries in References.

Coupled cluster

"RI-CC" here means the two-electron integrals in the amplitude equations come from density fitting (three-centre B tensors from an RI auxiliary basis), not from exact four-centre integrals. Everything below needs an orbital basis and an RI auxiliary basis. Formal cost is O(N⁶) for CCSD and O(N⁷) for (T).

CCSD and CCSD(T)

What it is. RI-CCSD and the perturbative triples correction (T), in two implementations:

• Spin-adapted, closed shell (ccsd_closed_shell, ccsd_t_closed_shell): amplitudes over spatial orbitals, the algorithm PySCF's cc.CCSD uses. This is what both the CLI and Python run for an RHF reference. Its VVVV block is 16× smaller than the spin-orbital one; measured about 8–10× faster than the spin-orbital CCSD at cc-pVDZ, and the (T) step 9.6–42× faster.
• Spin-orbital (ccsd, ccsd_t): built from the same RHF spatial orbitals and kept as the cross-check of the spin-adapted solvers, so it is closed shell only as well. Its (T) streams one occupied triple at a time, so memory is O(n_o·n_v³)-class rather than the dense six-index tensor.

No open-shell CCD, CCSD or CCSD(T) exists in ferric: every CC solver reads restricted orbitals and refuses an unrestricted reference. The only open-shell coupled-cluster-type method is LinLCCD(hh): the CLI linlccd kind is closed shell, and its open-shell (UHF) version is library-only (ferric_cc::linlccd_u::u_linlccd).

Run it.

• CCSD: method.kind = "ccsd" (examples/water-ccsd.toml, water/cc-pVDZ); Python ferric.run_ccsd(mol, bs, aux).
• CCSD(T): method.kind = "ccsd(t)" (examples/water-ccsd-t.toml); Python ferric.run_ccsd_t(mol, bs, aux).
• CCD: method.kind = "ccd" (examples/water-ccd.toml); Python ferric.run_ccd(mol, bs, aux).
• All three kinds are closed shell only and refuse multiplicity > 1.

Aux basis. The RI error is not negligible at CC accuracy: on water / cc-pVDZ with cc-pvdz-ri, RI-CCSD differs from exact-integral CCSD by about 1.3e-4 Ha. Pick the aux for the accuracy you need, not by habit.

Accuracy. ccsd is Proven.

Limits. Closed-shell entry points only in the CLI and Python. The spin-adapted (T) rejects spin-orbital amplitudes with a typed error rather than mixing conventions.

LinLCCD and ωB97X-L-V

LinLCCD(hh) is linearized coupled-cluster doubles with the hole–hole ladder kept to all orders, closed shell only. The ladder keeps the correlation energy finite as the HOMO–LUMO gap closes, where MP2 diverges. method.kind = "linlccd" (examples/water-linlccd.toml). Proven (narrow): no other quantum chemistry code implements LinLCCD(hh), so its energy is checked against an independent numpy solve of the same equations on PySCF density-fitted integrals (H2O and NH3 through the closed-shell path, and UHF OH through the library-only open-shell path; agreement ≤1.2e-12 Ha). With the ladder off it reduces exactly to RI-MP2, and with exact integrals its driver terms reproduce canonical MP2.

[mp2] linlccd_variant (Python run_linlccd(variant=...)) selects the method's ladder terms: hh (default), drivers-only (no ladder, equal to RI-MP2) or full (hole–hole plus particle–particle, with CCD-like VVVV memory). Every variant is computed exactly by default. All three variants match an independent numpy solve on the same density-fitted integrals (H2O and NH3, 6-31G and cc-pVDZ, ≤ 4.4e-13 Ha); drivers-only also matches PySCF DFMP2.

Local LinLCCD. With [local] scheme = "amplitude-threshold" and eps (examples/water-linlccd-local.toml; Python run_linlccd(..., local="amplitude-threshold", eps=1e-4)) the same variant is solved in the Boys-localized basis with pair amplitudes at or below eps dropped (see Exact and local correlation). eps has no default and is printed and logged with the kept fraction; eps = 0 reproduces the exact LinLCCD of the same variant and the numpy reference (≤ 4.4e-13 Ha, every variant), which is how the local path is anchored. [local] reference = true also runs the exact LinLCCD and prints the local error. Closed shell and energy only.

ωB97X-L-V is a double-hybrid functional that uses short-range LinLCCD(hh) instead of MP2 for its correlation term. It converges its own ωB97X-L Kohn–Sham reference (a non-converged reference is an error), then adds the LinLCCD(hh) correction on those orbitals. method.kind = "wb97x-l-v" (examples/water-wb97xlv.toml). [dft] lambda and omega override the published 0.6 and 0.1 Bohr⁻¹; omitting them gives the published values. λ enters the amplitude equations as well as the energy (the paper's eqn 22), so the correlation term is quadratic in λ at leading order. Validated against PySCF with the published parameters plus a numpy LinLCCD(hh) on water and OH, and against the paper's Be₂ bond energy (2.3 kcal/mol; ferric 2.299).

MP2-based double hybrids

B2PLYP and DSD-PBEP86: a KS reference with weighted exchange and correlation components, plus scaled (SCS-)RI-MP2 correlation. method.kind = "b2plyp" / "dsd-pbep86" (examples/water-b2plyp.toml); Python ferric.run_double_hybrid(mol, bs, aux, kind="b2plyp"). Spike: no comparison to a reference code yet.

Implementation

All contractions go through einsum!, a macro that maps tensor contractions onto BLAS3 GEMMs. The permutation copies that feed those GEMMs are parallelized, because for a strided permutation the copy can dominate the contraction it feeds: measured at 47% at n_v = 40 and 70% at n_v = 80. The copies are bit-identical regardless of thread count, since a permutation writes each output element exactly once; a test pins this and has been checked to fail when deliberately broken.

Memory

The amplitude tensors dominate and grow as \( n_o^2 n_v^2 \), or \( (2n_o)^2 (2n_v)^2 \) in the spin-orbital drivers. Memory budgets are enforced: an oversized job is refused with a breakdown naming the dominant term instead of being OOM-killed partway through.

Cite

CCSD: Scuseria, Janssen & Schaefer 1988; spin-adapted closed-shell CCSD equations: Hirata et al. 2004. (T): Raghavachari et al. 1989; closed-shell (T) algorithm: Rendell, Lee & Komornicki 1991. LinLCCD(hh): Carter-Fenk 2025. ωB97X-L-V: Ransford & Carter-Fenk 2026. B2PLYP: Grimme 2006. DSD-PBEP86: Kozuch & Martin 2011. Review: Bartlett & Musiał 2007. Full entries in References.

RPA, GW and excited states

Methods built on the density–density response function: RPA correlation energies, GW quasiparticle energies, BSE and TDDFT excitations, and polarizabilities and \( C_6 \) coefficients. All need an RI auxiliary basis ([rpa] auxbasis). The [rpa], [gw] and [tddft] keys are in Input file; grades are on Capabilities and validation.

What PDEP does in ferric

The independent-particle response \( \chi_0(i\omega) \) is built in the RI auxiliary basis from the three-centre B tensors, as an explicit (Adler–Wiser) sum over every occupied–virtual pair \( ia \), weighted by \( 4\varepsilon_{ia}/(\omega^2 + \varepsilon_{ia}^2) \) (sternheimer::dielectric_matrix). The sum over empty states is still there.

PDEP (projective dielectric eigenpotentials) then works in the eigenbasis of the static dielectric matrix in that RI space. Eigenpotentials whose eigenvalue is within trunc_thresh of 1 (default 1e-4) carry almost no screening and are dropped, so the frequency-dependent work runs in a smaller basis. That compression, not the removal of the empty-state sum, is what PDEP contributes here. How much it saves depends on the threshold; runs that need the full-rank answer set trunc_thresh = 0.0, as the GW, BSE and \( C_6 \) examples do.

RPA correlation energy

What it is. Direct RPA (dRPA) correlation from the dielectric eigenvalues on an imaginary-frequency quadrature.

• PDEP-RPA, closed shell: method.kind = "pdep-rpa" (examples/water-pdep-rpa.toml); Python ferric.run_pdep_rpa. Proven.
• U-PDEP-RPA, open shell over a spin-summed dielectric. From the CLI, set method.kind = "pdep-rpa" with multiplicity > 1 and task = "energy": the CLI solves UHF (UKS with [rpa] xc) with MOM after 5 iterations and runs U-PDEP-RPA on it. It is CLI-only: Python run_pdep_rpa is closed shell only. The library (ferric_rpa::run_u_pdep_rpa) also accepts a ROHF (or ROKS) reference, which it semi-canonicalizes first: each spin uses the orbitals and energies of its own Fock matrix, diagonalized in its occupied and virtual blocks (see the anchors).
• Attenuated RPA: short-range correlation with an erfc operator.
• RS-MP2 + LR-RPA: short-range MP2 plus long-range dRPA, on the MP2 page.

The static eigensolve defaults to Lanczos, with a dense path for small problems. Geometry optimization with pdep-rpa is supported (task = "optimize") on a closed-shell RHF reference.

Exact and local dRPA

What it is. method.kind = "drpa" (Python ferric.run_drpa) is dRPA@HF computed by the drCCD Riccati equations in the Boys-localized basis, closed shell, energy only. It is exact by default: no amplitude is truncated (examples/water-drpa.toml), and the energy equals the canonical plasmon formula to ≤ 1e-12 Ha. Riccati, plasmon and full-rank PDEP-RPA are algorithms for the same exact dRPA; measured agreement is ≤ 2.6e-14 Ha (H2/STO-3G, water/6-31G with and without frozen core, PDEP with trunc_thresh = 0 at 64 Gauss–Legendre points). Proven (narrow).

Which exact algorithm. The Riccati solve holds a ring-product plan of no³·nv² numbers, no times the size of the amplitudes, so it is the small-system path: C12 thrashed and was then killed for memory. A run that cannot fit the memory budget is refused before the SCF and pointed at pdep-rpa with [rpa] trunc_thresh = 0, which gives the same energy (to its frequency-quadrature error) at far lower memory and is faster at every size measured (n-alkanes C4–C16).

The local approximation ([local] scheme = "amplitude-threshold" with eps, examples/water-drpa-local.toml; Python run_drpa(..., local="amplitude-threshold", eps=1e-4)) drops pair amplitudes whose localized |2(ia|jb)| is at or below eps. eps has no default and is printed and logged with the kept fraction; eps = 0 is the exact method. dRPA is not variational, so the error is first order in what is dropped. It is positive (less correlation) at every point measured, which is a measurement, not a guarantee, and grows faster than linearly in eps: on n-octane / 6-31G it is 6.7e-7, 3.9e-5, 5.5e-4 and 1.0e-2 Ha at eps = 1e-6, 1e-5, 1e-4 and 1e-3, keeping 82%, 50%, 18% and 3% of the amplitudes. [local] reference = true (Python compute_reference=True) also computes the canonical plasmon dRPA and prints the error against it. [local] eps_sweep (Python run_drpa_scan) evaluates several eps on one SCF and one localized assembly. Proven (narrow) through its exact limit.

What is measured (n-alkanes, 6-31G / cc-pVDZ-RI, frozen carbon cores, Coulomb, single thread; local error against the plasmon formula, PDEP error against full-rank PDEP at 32 quadrature points). The fraction of the truncated object retained at 1 kcal/mol error:

The amplitude threshold compresses more as the molecule grows, and PDEP's fraction stays flat. The two fractions are of different objects (no²nv² amplitudes against naux modes), so this is not a cost comparison. In wall time at a matched error of about 1 kcal/mol (reference off), truncated PDEP stays about 8× faster from C4 to C16 (C16: 165.7 s local at eps = 1e-4 against about 21 s PDEP), and both grow at the same rate in that range. No speedup over PDEP is claimed for local dRPA.

GW

What it is. Quasiparticle energies from the GW self-energy: G0W0, COHSEX, evGW0 and evGW, closed shell, plus unrestricted U-GW. The starting point is HF by default or a KS functional ([rpa] xc). From a KS starting point the static term Σx − v_xc enters the G0W0, evGW₀ and evGW quasiparticle equation (for U-GW, each spin's own equation with that spin's v_xc), so Σc is evaluated at the shifted root; COHSEX and U-COHSEX, which are static, add it to the quasiparticle energy.

Run it. method.kind = "gw" with [gw] method = "g0w0" (examples/water-g0w0-pbe.toml, open shell examples/oh-ugw.toml); Python ferric.run_gw, run_u_gw. The open-shell reference is UHF by default; [gw] reference = "rohf" (Python run_u_gw(reference="rohf")) uses ROHF instead, or ROKS with [rpa] xc (examples/oh-ugw-rohf.toml), semi-canonicalized per spin as for U-PDEP-RPA.

Accuracy. Smoke; treat results as about ±0.3 eV.

Limits. The quasiparticle equation is solved by a Newton root search on the self-energy, which is fragile near \( \Sigma_c \) poles. Runs report whether each root and each eigenvalue-self-consistency loop converged, and warn when one did not; check those flags.

BSE-TDA

What it is. Bethe–Salpeter excitation energies in the Tamm–Dancoff approximation on top of G0W0@HF quasiparticle energies, closed shell.

Run it. method.kind = "bse-tda" (examples/water-bse-tda.toml, and a set of *-bse-tda-augdz.toml examples for small organics); Python ferric.run_bse_tda.

Accuracy. Smoke. Given the same quasiparticle energies, the lowest five singlet excitation energies match an independent numpy BSE-TDA (PySCF density-fitted integrals, static RPA W) to 1.8e-10 Ha and their oscillator strengths to 2.6e-9, for H2O at cc-pVDZ and aug-cc-pVDZ and NH3 and CH2O at cc-pVDZ. The quasiparticle energies come from the internal G0W0@HF, which matches PySCF only at [rpa] n_quad = 100 and trunc_thresh = 0; the defaults are coarser. The core and high-virtual quasiparticle energies are ill-conditioned, which moves the lowest five excitations by at most 2.2e-7 Ha.

TDDFT and TDA

What it is. Linear-response excitations in the Tamm–Dancoff approximation (TDA, which is CIS for an HF reference) and the full Casida equations, closed shell.

Run it. method.kind = "tda" or "tddft" with [tddft] n_roots and optionally xc (examples/water-tda.toml, examples/water-tddft-pbe.toml); Python ferric.run_tddft(mol, bs, aux, functional=..., method="tda") or method="casida".

Scope. Closed-shell references, singlet excitations. With a DFT reference the \( (ia|f_{xc}|jb) \) XC-kernel term is included; with no functional the result is CIS/TDHF. Meta-GGA, VV10 and range-separated functionals are refused. Grade: Proven (narrow, closed shell): water, formaldehyde and NH3 at 6-31G and aug-cc-pVDZ with HF, LDA, PBE and B3LYP match PySCF TDA/TDDFT to at most 6.5e-4 eV, test bar 1e-3 eV (ferric-tddft/tests/validation_tddft.rs).

A separate, library-only TDA-DFT in ferric-gw/src/tddft.rs uses the same kernel (ferric_dft::lr_kernel) and is pinned against PySCF (ferric-gw/tests/tda_dft_vs_pyscf.rs). The user-facing TDA reproduces it to 1e-8 Ha for HF, PBE and B3LYP.

Polarizabilities and dispersion coefficients

What exists. Static molecular and atom-partitioned polarizabilities, Casimir–Polder \( C_6 \) coefficients from per-atom dynamic polarizabilities \( \alpha^A(i\omega) \), and many-body dispersion (MBD). Three sources feed the \( C_6 \) contraction ([rpa] c6_source):

• ts: the Tkatchenko–Scheffler single-pole model (default).
• mbd: self-consistent dipole screening of the TS polarizabilities. This is the full-range screening of Tkatchenko et al. 2012 (Gaussian-damped dipole tensor, no Fermi range separation), not the range-separated screening of MBD@rsSCS.
• pdep: dynamic PDEP-RPA polarizabilities (examples/water-c6-pdep.toml, examples/argon-c6-rpa-pbe.toml).

MBD@rsSCS dispersion energy. ferric.mbd_rsscs_energy(mol, volume_ratios, functional=... | beta=...) (Rust: ferric_rpa::dispersion::mbd_rsscs_energy) computes the MBD@rsSCS energy of Ambrosetti et al. 2014 for a finite molecule, standalone (no SCF). Inputs are per-atom Hirshfeld volume ratios, which scale the free-atom α, C6 and R_vdW (Z = 1–54). The polarizabilities are screened with the short-range part (1 − f) of the Gaussian-damped dipole tensor on libMBD's 15-point imaginary-frequency grid; the energy couples the screened oscillators through the long-range part f of the bare dipole tensor, with the Fermi function f(R) = 1/(1 + exp(−6(R/(β(R_A + R_B)) − 1))). β is functional dependent: PBE 0.83, PBE0 and HSE06 0.85; any other functional needs an explicit beta. It returns the energy and the screened α₀, C6, R_vdW and ω per atom, and raises on a polarization catastrophe (a non-positive coupled-oscillator eigenvalue). No periodic systems. As a dispersion correction on a Kohn–Sham SCF it is [dft] dispersion = "mbd" / run_dft(dispersion="mbd"), which take the volume ratios from the converged density (reported as DftResult.volume_ratios) and add the energy and its analytic gradient (see SCF and DFT).

The argon-c6-rpa-pbe.toml header records C6(Ar–Ar) = 56.4 a.u. at RPA@PBE/aug-cc-pVTZ against the DOSD value 64.3 (−12%). Which of TS and PDEP-RPA gives better molecular \( C_6 \) is not established.

Use an augmented basis for any polarizability or \( C_6 \): without diffuse functions the dipole response is badly underestimated.

TDHF/RPAx \( C_6 \) is a measured negative. \( C_6 \) built on the RPAx@PBE kernel stays about 63% low regardless of the gap. Its static polarizability is not established either: at a physical scissor (0.36 Ha) water/cc-pVDZ gives an isotropic α of 5.20 a.u. against the DOSD 9.64 a.u. (−46%). The 9.24 a.u. quoted in the example's header comes from scissor = 0, where the tensor has a negative diagonal component, and that setting is refused. method.kind = "tdhf-static-polarizability" computes static α only. It needs a KS reference ([rpa] xc), and at the default [gw] scissor = 0 it can hit an excitonic instability, which is reported as an error rather than a negative α; examples/water-tdhf-static-alpha.toml hits it as shipped, so set scissor to about 0.3–0.4 Ha.

Cite

PDEP: Wilson, Gygi & Galli 2008. RI-RPA quadrature: Eshuis, Yarkony & Furche 2010; minimax grids: Kaltak, Klimeš & Kresse 2014. GW: Hedin 1965; GW100: van Setten et al. 2015. TDDFT review: Dreuw & Head-Gordon 2005. TS: Tkatchenko & Scheffler 2009; MBD: Tkatchenko et al. 2012; MBD@rsSCS: Ambrosetti et al. 2014. Full entries in References.

Constrained DFT

Charge- and spin-constrained DFT, and the electron-transfer couplings that follow from it.

Run it

Python and the Rust library; there is no CLI section.

• run_cdft(mol, basis_set, constraints, functional=None, ...) returns a CdftResult: a constrained UHF solve, or UKS when functional names a libxc functional other than "HF" (None and "HF", any case, give UHF).
• CdftConstraint(atoms, target, kind="charge") defines one fragment constraint. atoms are 0-based atom indices. target is the electron population on the fragment, the Becke-weighted trace \( \mathrm{Tr}[W D] \), not a net charge: \( N_\alpha + N_\beta \) for kind="charge" and \( N_\alpha - N_\beta \) for kind="spin". A neutral He atom has a charge population of 2.0; He⁺ has 1.0.
• cdft_coupling(state_a, state_b) returns a CdftCouplingResult with the Wu–Van Voorhis coupling h_ab, the determinant overlap s_ab and the two diabat energies e_a, e_b. The raw element uses each state's free energy F = E + λN, which makes the coupling independent of a constant shift of the constraint operator.

In Rust the entry points are ferric_scf::cdft_driver::solve_cdft_uhf and ferric_scf::cdft_coupling::coupling_hab; see the Rust API.

This example reproduces the HeNe⁺ constrained solution that crates/ferric-scf/tests/cdft_outer_loop.rs pins (E = −130.4021906 Ha, λ = −2.754 Ha per electron), with the He fragment held at 2.0 electrons:

import ferric

# HeNe+ doublet at 2.0 Å; hold the He atom (index 0) at 2.0 electrons.
mol = ferric.Molecule.from_xyz_string(
    "2\nHeNe+\nHe 0.0 0.0 0.0\nNe 0.0 0.0 2.0\n", 1, 2
)
r = ferric.run_cdft(
    mol,
    ferric.BasisSet.bundled("def2-svp"),
    [ferric.CdftConstraint([0], 2.0, kind="charge")],
    guess="hcore",
    level_shift=0.5,
    max_iter=400,
    lambda_tol=1e-5,
    max_outer=40,
    stability_descent=False,
    grid_radial=99,
    grid_angular=302,
)
print(r.converged)
print(f"{r.energy:.6f}")
print(f"{r.populations[0]:.5f}")

True
-130.402190
2.00000

r.energy is the ordinary UHF/UKS energy at the constrained density, without the constraint term. r.lambdas holds the multipliers (Ha per electron) and r.weight_matrix(i) the AO weight operator of constraint i.

What to know before using it:

• A returned result has a converged λ loop. When the outer loop exceeds max_outer, run_cdft raises RuntimeError. The inner SCF at the final λ can still be unconverged, so check r.converged: it is true only when the inner SCF converged and every constraint is met to lambda_tol.
• stability_descent defaults to True here (in run_uhf it defaults to False). After the λ loop converges, a constrained saddle point is followed downhill to a lower state that still meets the constraint. The example turns it off to reproduce the pinned solution; with it on, this HeNe⁺ case ends about 0.0245 Ha lower. The descent is skipped for a KS reference.
• The weight grid defaults to 99 × 302. It must resolve populations below lambda_tol; 75 × 110 resolves them only to about 1e-4. Setting grid_radial or grid_angular uses that grid for the XC quadrature too.
• The UKS path is the one compared against another code. Constrained UKS/PBE matches NWChem 7.2.2 (see Accuracy). The UHF path (functional=None or "HF") has only internal checks: NWChem's standalone SCF (Hartree–Fock) module has no cDFT, because the Becke weight operator is built on the XC grid of its DFT module. A Hartree–Fock comparison through that DFT module (xc HFexch) has not been made.
• One constraint is the tested case. For one constraint the outer loop is a Newton step kept inside a sign-change bracket. It backs off from a λ whose inner SCF does not converge, and it discards a finite-difference derivative whose two inner solves landed in different SCF states. With several constraints the outer loop is a plain k × k Newton step without these safeguards.
• cdft_coupling has strict preconditions and raises ValueError when one fails. Each state carries exactly one kind="charge" constraint, both are converged, and both come from the same molecule, geometry, basis, charge and multiplicity and the same Hamiltonian: functional, df_j_aux/df_k_aux, k_builder, XC grid, point charges and external field. Two states that are the same determinant (\( |S_{ab}| \to 1 \)) also raise. The sign of h_ab is a determinant-phase convention; compare \( |H_{ab}| \).

Accuracy

Constrained UKS/PBE energies and multipliers are compared against NWChem 7.2.2 cdft ... pop becke on LiH, HF and H2O⁺ (charge and spin constraints) at 6-31G and def2-SVP (ferric-scf/tests/validation_cdft.rs): E(N) − E_unconstrained to 2.0e-7 Ha and λ to 1.1e-6 against NWChem's grid limit, and dE/dN = −λ to 2.3e-12 Ha (this identity is checked at def2-SVP, on the first target of each constraint kind); numbers are on Capabilities and validation. No external reference value is stated for UHF-cDFT energies. The He₂⁺ coupling ingredients (determinant overlap, one- and two-electron transition elements, |V|) match NWChem's et module to ≤ 6e-11 Ha on NWChem's own determinants, and the KS diabats and couplings match end to end when started from NWChem's λ, except def2-SVP at 3.50 Å, where inner SCF solves within 1e-9 in λ of the root do not converge in 100 iterations (see Capabilities and validation). These diabats are validated only from NWChem's λ, set through the Rust RhfConfig::cdft_lambda_init; run_cdft starts from λ = 0, where the symmetric pair is delocalized and the inner SCF is bistable. The tests also check the coupling kernel on synthetic matrices and He₂⁺ identities (ferric-scf/tests/cdft_coupling.rs), probe HeNe⁺ over a distance series (cdft_coupling_hene.rs), and check exact identities on LiH/def2-SVP (cdft_uhf.rs): the constraint is satisfied, λ = 0 reproduces plain UHF, and the constraint composes with an external point charge. The Python tests (crates/ferric-python/tests/test_cdft.py) rerun those configurations through the bindings and check cdft_coupling against an independent transition-density construction. Treat UHF-cDFT energies as unvalidated against other codes; see Capabilities and validation.

The response connection

A cDFT constraint couples a Lagrange multiplier \( \lambda \) to a fragment-weighted density operator. The derivative

\[ \frac{\partial N}{\partial \lambda} \]

— how much charge moves per unit constraint potential — is a susceptibility. So cDFT probes the same object as RPA and GW and attenuated MP2, through a different coupling.

Implementation

• Fragment charge and spin constraints via a grid-Becke weight operator
• A nested Lagrange-multiplier solve (Wu–Van Voorhis): an inner SCF at fixed \( \lambda \), an outer Newton iteration on \( \lambda \) itself

The nesting is what makes cDFT more expensive than a plain SCF — each outer step is a full converged inner solve.

Electron-transfer coupling

Once you have two charge-localized diabatic states, the coupling \( H_{ab} \) between them follows from a non-orthogonal determinant overlap, computed via Löwdin biorthogonalization.

That gives the matrix element governing electron-transfer rates in Marcus theory, from states that are constructed rather than guessed.

A caveat

The λ convergence tolerance (lambda_tol in Python, cdft_lambda_tol in Rust, default 1e-5 electrons) interacts with the coupling calculation in a way worth checking: a loosely converged \( \lambda \) produces diabatic states that are not quite the ones you asked for, and \( H_{ab} \) inherits that error. Tighten it before trusting a coupling. The exception is a constraint whose population barely responds to \( \lambda \), such as He₂⁺ at the localized-hole plateau: there the tolerance has to be loosened to match that flatness (the tests use 1e-2), or the Newton step drives \( \lambda \) off a cliff.

Cite

Wu & Van Voorhis 2006 (electron-transfer coupling from cDFT); Becke 1988 (the fragment weight partition). Full entries in References.

Electronic response

ferric is organized around electronic response: how the electron density reacts to a perturbation. Standard quantum-chemistry codes are usually organized around a hierarchy of wavefunction ansätze (HF → MP2 → CCSD → CCSD(T)). That is a perfectly good organizing principle. It is not the one used here.

The central object appears under several names depending on which coupling you look at:

These are the same physics viewed through different couplings. A code that computes one well should be able to compute the others, and errors in one should be diagnosable as errors in the others.

The claim

The premise behind the architecture is that response is:

1. Local in real space: a density fluctuation here does not much affect the density far away, so the response should be sparse in a localized basis.
2. Low-rank in its eigenspectrum: the dielectric matrix has a small number of dominant eigenmodes, so it can be compressed without losing the physics.

If both hold, organizing the computation around response should make it cheaper: attenuate the operator, keep the dominant dielectric modes.

What is actually demonstrated

Low rank: demonstrated, and narrower than it sounds. PDEP compresses the dielectric matrix to its dominant eigenpotentials in the RI auxiliary space, and that works in production paths; see RPA and GW. It does not remove the sum over empty states: ferric builds \( \chi_0 \) by summing over every occupied–virtual pair, and PDEP compresses what comes out of that sum.

Locality: one positive result, still without a speedup.

• AO-sparse Laplace SOS-MP2. Restricting each localized orbital's pseudo-density to an AO domain works: the radius needed grows far more slowly than the molecule (chemical accuracy at 3 to 5 Bohr from ethane to dodecane, while radius/diameter falls from 0.52 to 0.17; within 0.05% at 4 Bohr on a 71-atom drug molecule). The tensor algebra is still dense, so no timing gain is claimed. This result is specific to that formulation and does not carry over to the other locality lanes below. The test sos_ao_sparse_truncation_radius_is_transferable_across_sizes pins the STO-3G butane/octane comparison (12 Bohr exact on both; octane worse at 3 Bohr). The C2-C12 sweep and the drug-molecule figure are measurements, not regression tests.
• Local MP2 (amplitude threshold) has localized virtuals and per-pair domain-local RI fits. The integral-direct variant (rimp2 with [local] integral_direct = true) is measured at about N^1.24 (erfc) to N^1.4 (Coulomb) on alkanes C20–C48, three points in one basis, so the reading is provisional. Local MP2 without integral_direct still builds the global 3-index tensor and makes no scaling claim.
• RI-Laplace MP2 is dense; it is the correctness reference for the AO formulation, not a reduced-scaling path.

The locality premise is therefore partly supported and not yet cashed in as cost. Measured limits are kept on Capabilities and validation.

Why this framing is useful anyway

Even where the scaling payoff has not arrived, the response framing makes the error in one method diagnosable through another.

MP2's dispersion is the clearest case. Its dispersion energy is built from an uncoupled (uncoupled Hartree–Fock) response, in which the density fluctuation does not feel the field it creates. For systems with low-lying, highly polarizable excitations, such as π-stacked and other π systems, that uncoupled response over-polarizes, and MP2 overbinds. This is documented in the literature that replaces MP2's uncoupled dispersion with a coupled one (Cybulski & Lytle 2007; Heßelmann 2008; Pitoňák & Heßelmann 2010). The size of that coupling correction varies from system to system, and ferric has no coupled-dispersion (MP2C-style) implementation of its own, so read this as the literature's diagnosis, not a ferric measurement.

Attenuated MP2 takes a blunter route: it removes the long-range correlation altogether and relies on a fitted short-range operator plus the basis-set error it is fitted in. That is why it has zero asymptotic \( C_6 \), why its parameters belong to one basis, and why MP2-V adds long-range dispersion back through VV10. RS-MP2 + LR-RPA puts the long-range part back through response instead. See The MP2 family.

Where the methods come from

The method families in ferric are not an arbitrary selection. Each one works on the response function from a different direction. The definitions and parameters live on the method pages; this page is the map.

Attenuated MP2: removing the long-range part

MP2's dispersion comes from an uncoupled response, which over-polarizes for systems with low-lying, highly polarizable excitations (π-stacked aromatics are the classic case); it is not a general property of polarizable molecules. In small basis sets the resulting overbinding is partly cancelled by basis-set superposition error, which disguises it.

Attenuated MP2 replaces \( 1/r \) in the correlation energy with a short-range operator (erfc or terfc) and fits its range to interaction energies in a chosen basis. It removes the long-range correlation entirely rather than correcting it, so it has no asymptotic \( C_6 \), and its parameters are specific to the basis and protocol they were fitted in. MP2-V restores long-range dispersion with VV10; RS-MP2 + LR-RPA restores it with long-range RPA.

Goldey & Head-Gordon (JPCL 2012) introduced it in aug-cc-pVDZ; Goldey, Dutoi & Head-Gordon (PCCP 2013) introduced terfc in aug-cc-pVTZ; the dual-attenuated SCS variant is Goldey & Head-Gordon (JPCB 2014). Operators, parameters and fitting protocol: The MP2 family.

PDEP-RPA and GW: compressing the response

The dielectric matrix is built from the density–density response function. ferric forms the independent-particle response in the RI auxiliary basis by summing over occupied–virtual pairs, then works in the eigenbasis of the static dielectric matrix, dropping eigenpotentials that carry almost no screening. PDEP (projective dielectric eigenpotentials) is that compression. It is the demonstrated part of the low-rank premise; the empty-state sum that feeds it is still there. See RPA, GW and excited states.

Constrained DFT: reading the response

A cDFT constraint couples a Lagrange multiplier \( \lambda \) to a fragment-weighted density operator. The derivative \( \partial N / \partial \lambda \), how much charge moves per unit constraint potential, is a susceptibility.

That makes cDFT a direct probe of the same object. It also yields charge-localized diabatic states whose electron-transfer couplings \( H_{ab} \) follow from non-orthogonal determinant overlaps. See Constrained DFT.

What this buys

Three families, one object. An error in the polarizability shows up as an error in dispersion, in screening, and in charge-transfer coupling, so a fix validated in one place has predictable consequences in the others.

That is the design bet. Whether it pays off in cost is still open; see Electronic response for what has and has not been measured.

Architecture

A Cargo workspace of focused crates, layered so that method crates depend on integrals and core, not on each other.

                          +------------------+
                          |   ferric-cli     |   TOML config -> most methods
                          +--------+---------+   (+ ferric-python: pyo3 bindings)
                                   |
   +-----------+-----------+-------+------+-----------+------------+
   |           |           |              |           |            |
+--v----+ +----v----+ +----v----+   +-----v----+ +----v-----+ +---v------+
|ferric | |ferric   | |ferric   |   |ferric    | |ferric    | |ferric    |
|-scf   | |-mp2     | |-dft     |   |-rpa      | |-gw       | |-cc       |
|RHF/UHF| |RI-MP2,  | |RKS/UKS/ |   |PDEP-RPA, | |G0W0,     | |CCD/CCSD/ |
|/ROHF, | |OO,att,  | |ROKS,    |   |U-PDEP,   | |COHSEX,   | |(T)       |
|KS-DFT,| |SCS,     | |libxc,   |   |response  | |evGW,     | +----------+
|DIIS,  | |2terfc,  | |Becke    |   |props,    | |U-GW      |
|MOM,AH,| |Laplace  | |grids,   |   |ESP/Hirsh/| +-----+----+
|cDFT,  | +----+----+ |VV10     |   |NPZ export|       |
|grads  |      |      +----+----+   +-----+----+       |
+---+---+      |           |              |            |
    |          +-----+-----+------+-------+------------+
    |                |     |      |
    |   +------------v--+ +v------v-----+   ferric-tensors (sparse),
    |   |ferric-export | |ferric-      |   ferric-quadrature (Laplace/grid roots)
    |   |cube,NPZ,GTO  | |integrals    |   support crates
    |   +--------------+ |libint2 FFI  |
    |                    |shim/shim.cc |   Coulomb/erf/erfc, 1e/2e/3c/2c, derivs
    +--------+-----------+------+------+
             |                  |
        +----v------------------v----+
        |        ferric-core         |   Molecule, BasisSet, Shell, elements,
        |                            |   BSE-JSON / G94 parsers, bundled bases
        +----------------------------+

Layers

ferric-core — molecular structure, basis sets, shells, elements, and the BSE-JSON / Gaussian-94 parsers. Also the home of shared infrastructure: configuration (ConfigVar), memory budgets (MemoryPlan), the BLAS-thread hazard model, and MPI context.

ferric-integrals — the libint2 FFI and its C++ shim. Coulomb, erf and erfc operators; 1-electron, 2-electron, 3-center and 2-center integrals; first derivatives. The shim wraps every libint2 call in try/catch and returns a sentinel, so a C++ exception never unwinds across the FFI boundary.

Method crates — ferric-scf, ferric-mp2, ferric-dft, ferric-rpa, ferric-gw, ferric-cc, ferric-ci, ferric-tddft, ferric-pcm, ferric-mm, ferric-xtb, ferric-d3 (Grimme D3(BJ) dispersion energy and gradient).

Support — ferric-tensors (the einsum! contraction macro), ferric-quadrature (Lebedev grids, minimax Laplace roots), ferric-export (cube files, NPZ, GTO evaluation).

Interfaces — ferric-cli (TOML-driven) and ferric-python (pyo3).

Cross-cutting conventions

Threading. rayon owns outer parallelism; BLAS is pinned to one thread inside any rayon worker, enforced at runtime rather than by convention. Throughput across many jobs comes from many single-threaded processes, not one wide job.

Determinism. Reductions fold in a fixed ascending order independent of thread count, so results are bit-identical across RAYON_NUM_THREADS. This is pinned by tests. A different-but-deterministic order — a tree-fold, say — would not be acceptable, because floating-point addition is not associative.

Memory. MemoryPlan expresses what a path will allocate and when. An allocation that exceeds the budget and has no spill-to-disk or recompute fallback is refused before allocating, with a breakdown naming the dominant term. The CLI installs one process-wide MemoryPool sized from the resolved budget, and the large, size-dependent allocations reserve their bytes from it, so two allocations alive at the same time cannot each claim the whole budget. Library and Python callers install no pool; each check then compares its own allocation with the whole budget. Allocations that have a fallback spill to disk or are recomputed instead of being refused, and basis-sized matrices, engines and scratch are not charged (see Sharp bits). Guards are tested in both directions: a starved budget must be refused (or must fall back), and an ample budget must still run.

Errors. Methods return Result; iterative solvers additionally carry a converged flag, since non-convergence is a result rather than an error. Callers are expected to check it.

Capabilities and validation

Implemented ≠ validated. Working code is not a checked number.

Read this page before trusting any result. It lists every CLI method.kind (which references it accepts, which tasks it supports, the matching Python function, an example input) and the Python entry points. Each graded CLI method.kind gets a grade (three kinds are ungraded, and Python entry points are not graded individually), and the page shows how far each grade's numbers have been checked against an independent reference and where they are known to fail.

How to read it: find the capability in the matrix or the Python entry points. A matrix Grade cell links to the Anchors table when a row there holds the evidence; the Known limits section records the measured negatives. Every matrix cell comes from the CLI dispatch in crates/ferric-cli/src/lib.rs or the bindings in crates/ferric-python/src/lib.rs. The page is maintained by hand, so check those two files if a cell looks wrong. For keyword details see the input reference; for the example files see the examples index.

Grades

The grades are the ones ferric's CLI uses. Every CLI method.kind except laplace-sos-mp2 has a grade. A kind's grade is for the EXACT method; the local approximation of rimp2, drpa and linlccd (the [local] section) is graded separately in that row's Caveat. The CLI prints the Smoke and Spike grades as a [warning] line on stderr at run time; Proven kinds and the ungraded kind print nothing. The warning text is the EPISTEMIC_WARNINGS table in crates/ferric-cli/src/lib.rs, and the "Caveat" column of the matrix condenses it.

Symbols: ✓ = supported. — = not supported (the CLI refuses it with an error). FD = finite differences of the analytic gradient (6N gradient evaluations). The analytic Hessian covers closed-shell RHF and UHF of any multiplicity with exact J/K, no ECP and a basis up to f functions; rhf and uhf frequencies use it there by default.

CLI method.kind matrix

"Reference" is the SCF reference the CLI actually builds for that kind. Open-shell support is listed only where the dispatch code handles it (see Open shells below).

task = "optimize" is accepted only for rhf, ksdft, uhf, rohf, pdep-rpa and rimp2. task = "frequencies" is accepted only for rhf, ksdft, uhf and rohf. On an open-shell molecule both tasks are accepted only for uhf, rohf and ksdft (UHF/UKS, ROHF/ROKS). Both tasks refuse [dft] grid_prune, and [scf] k_builder = "cosx" except for rhf, ksdft and uhf (closed-shell RHF/RKS and UHF; UKS with COSX is refused before the SCF); a COSX final pass is skipped for them. [dft] dispersion is refused for frequencies on an open-shell reference: the frequency driver threads the correction through the closed-shell SCF only, so it would report the Hessian of the uncorrected surface. The RKS, UKS and ROKS optimizers apply D3(BJ) and MBD@rsSCS (energy and exact gradient); an open-shell task = "energy" single point applies either. Any other combination exits with an error before the SCF runs.

Open shells in the CLI

• uhf and rohf read [molecule] multiplicity directly. With [dft] functional they run UKS and ROKS, for energies, optimizations and frequencies.
• ksdft with multiplicity > 1 runs UKS (the uhf route with the functional set). For ROKS use rohf with [dft] functional.
• rhf refuses multiplicity > 1 and points to uhf/rohf.
• rimp2 and oo-rimp2 run on the same plain UHF that kind = "uhf" runs, then take unrestricted RI-MP2 (UMP2, as PySCF mp.MP2(uhf)) and unrestricted OO-RI-MP2. [mp2] kappa is refused on an open shell. task = "energy" only: there is no unrestricted MP2 nuclear gradient.
• pdep-rpa, gw and mp2-v solve UHF with MOM after 5 iterations when multiplicity > 1, for task = "energy" only. For pdep-rpa and gw, setting [rpa] xc makes that reference UKS; gw with [gw] reference = "rohf" uses ROHF (ROKS with [rpa] xc) instead, semi-canonicalized per spin; U-G0W0 on a semi-canonicalized ROHF reference is in the Anchors, the ROKS reference is not compared. mp2-v does not read [rpa] xc and stays UHF.
• linlccd and wb97x-l-v refuse an open-shell molecule; their open-shell versions are library-only (ferric_cc::linlccd_u::u_linlccd, ferric_cc::double_hybrid::u_solve_wb97x_l_v).
• ccd, ccsd and ccsd(t) refuse an open-shell molecule, and no open-shell CCD, CCSD or CCSD(T) exists in the library either: every CC solver reads restricted orbitals.
• Every other kind refuses multiplicity > 1 with an error before any integral is computed.

Python entry points

These are Python functions without a method.kind of their own; some rows also have a CLI route through a TOML section or task, noted in the row. Scope is taken from the binding code and its docstrings. They are not graded individually.

Polarizable (Thole) embedding is Python and Rust only. IEF-PCM, UHF stability descent and terfc-attenuated MP2 are reachable from both: see [pcm], [scf] stability_descent and [mp2] att_operator in the Input reference.

Anchors

Where numbers are checked, they are checked against external references or exact limits, not against ferric's own earlier output. "Stated agreement" is the measured difference recorded in the repository; "test tolerance" is what the pinning test actually asserts, which is often looser. Proof links to the test file that asserts the row, and, where one exists, the script that generated its reference data; a row with nothing to link is not graded Proven. Proof files named validation_*.rs are #[ignore]d, so an ordinary cargo test skips them; they run in the weekly validation CI job and on demand (cargo nextest run --run-ignored only -E 'binary(/^validation_/)').

A "—" means the repository states no number for that cell. It is left empty on purpose rather than filled with an estimate.

Benchmark sweeps (GW100 and others) are kept in the project's working notes and are not reproduced here. Quoting a benchmark statistic from memory rather than from the record is the kind of unchecked claim this page exists to prevent.

Known limits and negatives

Reported rather than omitted:

• Analytic Hessians: closed-shell RHF and UHF only (exact J/K, no ECP, no embedding or solvent, basis up to f functions). ROHF has none (PySCF has no ROHF Hessian either). Every other method's frequencies take central finite differences of the analytic gradient.
• Local MP2: the integral-direct local MP2 (rimp2 with [local] integral_direct = true) has its correlation stage measured at about N^1.24 (erfc) and N^1.4 (Coulomb) on n-alkanes C20–C48 (6-31G / cc-pVDZ-RI, ε = 1e-4, frozen carbon cores, a calibrated pair gate; fitted to three points). That is one family of molecules in one basis: not shown to be linear, and not measured on 3-D or diffuse systems. Local MP2 without integral_direct still builds the global 3-index tensor and makes no scaling claim. See The MP2 family.
• Exact dRPA by the Riccati solve is a small-system path: its ring-product plan holds no³·nv² numbers, no times the amplitudes themselves (C12 thrashed, then was killed for memory). The CLI and run_drpa refuse a run that cannot fit before the SCF; full-rank pdep-rpa gives the same energy at far lower memory.
• Laplace SOS-MP2, AO-sparse variant: the domain truncation is accurate, and the radius it needs grows far more slowly than the molecule (chemical accuracy at 3 to 5 Bohr on n-alkanes C2 to C12, radius/diameter falling from 0.52 to 0.17; a 71-atom drug molecule is within 0.05% at 4 Bohr, about 13% of its diameter). The algebra is still dense, so no speedup is claimed. The regression test sos_ao_sparse_truncation_radius_is_transferable_across_sizes pins the STO-3G butane/octane comparison (12 Bohr exact on both; octane worse at 3 Bohr). The C2-C12 sweep and the drug-molecule figure are measurements, not regression tests.
• TDHF/RPAx C6: ~63% low regardless of gap. Do not use it for dispersion; its static α is not established either (see the matrix).
• TDDFT / TDA: closed-shell references only. Functionals without an f_xc kernel (meta-GGA, VV10, range-separated) are refused rather than run without it. B3LYP at aug-cc-pVDZ agrees with PySCF to 6.5e-4 eV, 40× worse than PBE in the same basis; the cause is not yet identified.
• COSX wins at high angular momentum, not at large system size. Measured on one thread at the flat (50,110) grid with the overlap fit:
  • Against LinK on n-alkanes at def2-SVP it is slower at every size measured: 1.59× (C20), 1.09× (C32) and 1.22× (C48), with no trend toward parity. At def2-TZVP on C20 it is faster (0.67×); that is the only triple-zeta size measured.
  • Against exact direct exchange on butane it is 3.7× slower at def2-TZVP (full SCF). At def2-QZVP its K build takes 90 s against about 400 s for the direct build's single J+K sweep, at a relative K error of 5.4e-5.
  • Its K build grows as N^1.29–N^1.32 between C20 and C48 at def2-SVP.
  • It applies to the Coulomb operator only: a range-separated functional takes its exchange from density-fitted short- and long-range fitters and ignores k_builder = "cosx", with a warning.
  • Its analytic gradient is exact for RHF, RKS and UHF with the overlap fit off, and for RHF and UHF with the default overlap fit, on flat and pruned (sgx) COSX grids. Fitted COSX with a KS functional, UKS and ROHF/ROKS are refused for gradient tasks, before the SCF runs.

Why the distinction is drawn so sharply

A quantum chemistry code can produce a plausible number in many ways that are wrong:

• a non-converged SCF returned as an ordinary result, because convergence is a flag rather than an error
• a method missing a physical term it does not mention (for example, TDDFT excitations computed without the f_xc kernel)
• a fallback model silently substituted for one atom in a molecule, changing a partitioning without changing the shape of the output
• a screening or truncation threshold that happens to be safe for the test system and not for yours

None of these look like failures, and every one of them has occurred in this codebase. The remedy is to grade each capability separately and say which ones are checked against ground truth.

Testing discipline

Some properties are pinned by tests rather than asserted in prose:

• Bit-identity across thread counts for reductions and permutations: results do not depend on RAYON_NUM_THREADS
• ERI 8-fold permutational symmetry: verified against the engine, not assumed, since the MP2 code relies on it to compute only ~1/8 of quartets
• Size-extensivity and rotational invariance of total energies
• Memory guards in both directions: a starved budget must be refused and an ample budget must still run, because an over-estimating guard is also a bug

New guards are mutation-tested: a deliberate defect is injected and the test confirmed to fail before the guard is trusted. This has caught guards that passed while proving nothing.

Related pages

• Input reference: every TOML key
• Examples: every shipped input file
• Python bindings

Input reference (TOML)

Every key the ferric CLI accepts, section by section.

Unknown keys are a hard error. The config structs in crates/ferric-cli/src/config.rs are #[serde(deny_unknown_fields)], so a misspelled key or section aborts the run before anything is computed. That includes [external_potential], its point charges, and each [[scf.ladder]] rung.

String values go through strict parsers, where an unknown value is an error at load time rather than a default. Most accept any capitalisation; the exceptions are noted per key.

This page is hand-maintained against crates/ferric-cli/src/config.rs at commit 4b64e6ce. Where a default is applied at the point of use rather than in config.rs, it was read from crates/ferric-cli/src/lib.rs at the same commit. If the code and this page disagree, the code wins. For which method.kind values exist and what each supports, see Capabilities and validation.

Units follow the code: [molecule] geometries are Å (XYZ), point charges and cutoffs named *_bohr are Bohr, keys named *_angstrom or documented as Å are Å, and energies are Hartree.

A minimal file:

[molecule]
xyz = "testdata/molecules/water.xyz"
[basis]
name = "cc-pvdz"
[method]
kind = "rimp2"
[mp2]
auxbasis = "cc-pvdz-ri"

[molecule] (required)

[basis] (required)

Auxiliary basis names used elsewhere (auxbasis, df_*_aux) come from the same table: cc-pvdz-ri, cc-pvtz-rifit, aug-cc-pv{d,t,q}z-rifit, def2-svp-rifit, def2-tzvp-rifit, def2-tzvpp-rifit, def2-qzvp-rifit, def2-qzvpp-rifit, def2-universal-jkfit, cc-pvdz-f12-optri. cc-pvdz-rifit is an alias of cc-pvdz-ri; both names load the same set.

[method] (required)

[scf]

Read by every kind, because every kind runs an SCF first.

[[scf.ladder]] rungs

Each rung overrides the flat [scf] settings. The rungs are walked in order, and the ladder stops at the first converged rung. Unknown keys are an error, as elsewhere.

[dft]

[mp2]

This section is shared by the whole MP2 family, ccsd, ccd, ccsd(t), linlccd, drpa, the double hybrids and tda/tddft, all of which read auxbasis and frozen_core from here. drpa, linlccd and a local rimp2 read nothing else from it (plus linlccd_variant for linlccd); any other [mp2] key on them is an error.

[local]

The local approximation of a correlated method. method.kind names the method (rimp2, drpa or linlccd); this section says whether and how its amplitudes are truncated. Without it (or with scheme = "none") the method is computed exactly. On any other kind the section is an error. The local runs are closed shell and task = "energy" only, and every printout and run-log result carries the model: "<method> (exact)" with "local": null, or the scheme, eps and kept fraction. See The MP2 family and Exact and local correlation.

[rpa]

Read by pdep-rpa, gw, bse-tda, tdhf-static-polarizability, and, for trunc_thresh only, rs-mp2-rpa.

[gw]

Read by gw. bse-tda and tdhf-static-polarizability read frozen_core, and tdhf-static-polarizability also reads scissor. The [rpa] section supplies the screened interaction.