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, and for no other platform:

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. Which build am I running shows how to check.

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 and the testdata/ molecules 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.

The large-libint wheel (GitHub releases)

Starting with 0.1.0rc8, each GitHub release also carries a second wheel for each CPython. It has the same version and commit as the PyPI wheel, but it is built against conda-forge's libint2 2.13.1 (the build scripts/install-libint.sh installs), which it bundles as a shared library. Its file name carries the wheel build tag 1largeam:

ferric-0.1.0rc8-1largeam-cp311-cp311-manylinux_2_28_x86_64.whl

Compared with the PyPI wheel it adds l = 7 shells, analytic 3- and 2-centre second derivatives, multipoles above order 2 and G12 (F12) integrals; see What the libint2 build carries. It does not extend analytic Hessians: 4-centre second derivatives stop at f functions in both wheels. RI-J/RI-K Hessians already work with the PyPI wheel, so they are no reason to switch. Nor is it a strict superset: its one-electron integrals reach l = 7 / 6 / 3 (energy / 1st / 2nd derivative), while the PyPI wheel's reach 6 / 6 / 6, so one-electron second derivatives of g, h or i functions are available only in the PyPI wheel.

Install the PyPI wheel first, so its dependencies are present, then swap in the file for your Python from the release (check it against the release's SHA256SUMS):

pip install ferric==0.1.0rc8
pip install --force-reinstall --no-deps \
  https://github.com/mgoldey/ferric/releases/download/v0.1.0rc8/ferric-0.1.0rc8-1largeam-cp311-cp311-manylinux_2_28_x86_64.whl

Both wheels are the distribution ferric at the same version, so pip list shows no difference and only one can be installed at a time:

• pip install ferric==0.1.0rc8 then reports "already satisfied" and keeps the large wheel.
• pip install -U ferric at the next release replaces it with that release's PyPI wheel; swap in the large one again from the new release.
• With both files on one --find-links index, pip picks the build-tagged (large) wheel.

To see which one is loaded:

python -c "import ferric; i = ferric.build_info(); print(i['libint_version'], i['libint_linkage'], i['libint_max_am'])"

The large wheel prints 2.13.1 shared and 'eri3': [7, 7, 4]; the PyPI wheel prints 2.7.2 static and 'eri3': [6, 6]. The libint: line of ferric --version reports the same limits.

The large wheel contains a libint2 binary (LGPL-3.0-only) built by conda-forge. Its source is the loriab/libint v0.1 generated tarball, built by the conda-forge libint feedstock; each release's notes link both and give the tarball's sha256. The licence texts ship with the release as libint2-LICENSE.txt and libint2-COPYING.LESSER.txt, covered by its SHA256SUMS.

Building from source

Build from source if you are changing ferric, need the MPI build, or are on a platform without a wheel (only Linux x86_64 is tested; see the table above). 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 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.

Building from the sdist

Each release publishes an sdist next to the wheels, so pip download ferric==<version> --no-binary :all: succeeds. The sdist holds the Rust workspace crates the extension depends on, Cargo.lock and the vendored libecpint. It does not hold libint2, so it is not a self-contained build: pip install from it needs everything in Prerequisites already present, with LIBINT2_PREFIX set to a libint2 install. The libint2 source is generated and too large to bundle (the export below is 183 MiB against PyPI's 100 MiB file limit).

Two ways to get a libint2 for an sdist build:

• Linux x86-64: scripts/install-libint.sh (conda-forge 2.13.1). The script is not in the sdist; take it from the repository.
• Any platform with CMake and Eigen3: build the export the wheels use, libint-2.7.2-ferric-small-1.tgz from the libint-2.7.2-ferric-small-1 release (verify its sha256 against LIBINT_SHA256 in .github/workflows/wheels.yml). Configure it with cmake .. -DCMAKE_INSTALL_PREFIX=$PREFIX -DCMAKE_PREFIX_PATH=$PREFIX -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DCMAKE_BUILD_TYPE=Release, then make and make install, and set LIBINT2_PREFIX=$PREFIX. The build takes about 30 minutes. It has the wheel's integral classes (see the table below), so F12 methods still need the 2.13.1 build. The recipe is tested only as part of the Linux x86-64 wheel build.

On a platform with no wheel, plain pip install ferric falls back to this sdist and fails at the libint2 step unless those prerequisites are present. An sdist build carries no git metadata, so the build identity reports commit: "unknown".

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):

The 3-centre limit is for the auxiliary shell; libint caps the two orbital shells of a 3-centre integral at its default limit for that derivative order (the wheel: l ≤ 6 for energies, l ≤ 4 for first derivatives; 2.13.1: 6 / 6 / 3).

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.

Which build am I running

python -c "import ferric; print(ferric.__version__); print(ferric.build_info()); print(ferric.__build__)"
ferric --version

The wheel and a conda-forge 2.13.1 build share a version and commit; they differ in the libint2 limits in the table above, which ferric --version (the libint: line) and ferric.build_info()["libint_max_am"] report.

ferric.build_info() and ferric --version report the version, the full git commit, whether tracked files differed from that commit (dirty), the cargo profile and the libint2 version. ferric.__build__ is the compact stamp {"git_sha": ..., "dirty": ...}. All of it is fixed when the extension is compiled, so it describes the loaded .so even if that is a symlink into another checkout.

• A released wheel carries the plain version (0.1.0), the tagged commit and dirty: False. To release, push a vX.Y.Z tag; there is nothing to commit, and the build fails rather than ship from a dirty or unidentified tree.
• A build from a checkout carries a .devN version (0.1.0.dev0) and the commit it was built from, with dirty: True if tracked files had uncommitted changes. Untracked files do not count.
• A build with no git metadata (an sdist or source tarball) reports commit: "unknown", and ferric.__build__ is None.

MPI

MPI is a source build only. The wheels workflow publishes no MPI wheel to PyPI.

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.

GPU (CUDA)

The CUDA backend is a source build only and is off by default even when built in.

cargo build --release -p ferric-cli --features ferric-cli/gpu
FERRIC_GPU=auto ./target/release/ferric examples/water-ccsd.toml

Building needs no CUDA toolkit: the driver and cuBLAS libraries are loaded at run time. Running needs an NVIDIA driver and the CUDA 12 runtime libraries (libcuda.so.1, libcublas.so.12) on the loader path; without them mode = "auto" prints a notice and runs on the CPU, and mode = "on" stops the ferric binary with the reason. Library and Python callers that set FERRIC_GPU=on with no usable device get the same notice and a CPU run, not an error; check ferric.gpu_status(). The CUDA 12 libraries support compute capability 5.0 and above; CUDA 13 does not support Pascal (compute capability 6.x), so a Pascal card needs the CUDA 12 libraries.

What runs on the device, under mode = "auto" or "on":

• Dense f64 contractions issued through einsum! (the coupled-cluster and MP2-family drivers) above a size threshold.
• The closed-shell RI-MP2 energy, with B_ov resident on the card: f64 by default ([precision] mode = "f64"), or with [precision] mode = "mixed" B_ov stored as f32 and accumulated in f64 (see the validation page for its measured error).
• The unrestricted RI-MP2 energy (open-shell reference, amplitude-free path), with one resident B_ov per spin: f64 by default, or with precision = "mixed" each spin's B_ov stored as f32 and the αα, ββ and αβ blocks accumulated in f64 by the same kernel.
• Mixed applies to the RI-MP2 energy under the Coulomb, erfc and terfc operators, closed-shell or unrestricted (plain RI-MP2, attenuated RI-MP2, SCS-MP2 with two terfc operators, the MP2 half of att-MP2+VV10); the erf operator, composite fitted operators, SR-MP2 in RS-MP2+RPA, OO-MP2 and kappa-regularised runs use the f64 device kernel.
• The density-fitted Coulomb matrix (RI-J): the raw three-index tensor is stored on the card once per geometry (packed lower triangle, f64) and each J build moves only the density weights and the result.
• The density-fitted exchange matrix (DF-K, occupied-orbital path): the dressed three-index tensor is stored on the card once per geometry and each K build moves only the occupied coefficients and the result, in f64.
• The Kohn-Sham exchange-correlation grid (closed- and open-shell, LDA, GGA and meta-GGA), only under mode = "on" (or a preset that sets it): auto keeps the XC grid on the CPU, because the device path was measured slower than the CPU. Measured on a GTX 1080 (Pascal, sm_61), PBE/cc-pVDZ SCF wall time: benzene 2.0 s on the CPU, 3.6 s with the device XC grid in f64 and 5.5-6.8 s in mixed; n-octane 8.5-8.9 s, 18-23 s and 17 s; each device XC build took 2.2-3x the host time (6-core host, RAYON_NUM_THREADS=3). With mode = "on" the screened per-batch basis-function values are stored on the card once per grid, and each Fock build runs the two per-batch matrix products (density Φ = D·χ and potential V = A·χᵀ) on the device in f64; naming xc-vxc in mixed_kernels (with precision = "mixed") forms the potential product from f32 operands with f64 accumulation while the SCF density is still changing (opt-in; it is not in the default mixed set). Grid construction, screening, the libxc functional evaluation, the density reductions and the fixed-order V accumulation stay on the CPU, so the result is reproducible and independent of the thread count. The VV10 nonlocal grid and XC gradients run on the CPU.

RI-J and DF-K run on the CPU when their tensor does not fit the device pool (or, for DF-K, is spilled or recomputed on the host) and under multi-rank MPI; the notice names the reason. When one SCF has both and the pool cannot hold both, DF-K gets the device. The XC grid runs on the CPU when its basis-function blocks do not fit what is left of the device pool, or after any device error (the notice names the reason; the failed pass is recomputed on the CPU). Integrals, SCF diagonalisation, DIIS and grid construction run on the CPU. Device results are reproducible run to run on one device and agree with the CPU to the accuracy of a different summation order, not bit for bit.

The shortest way to use the device is one line in the input, gpu = "on" (f64) or gpu = "mixed", or ferric --gpu mixed input.toml on the command line (the flag overrides the input file and FERRIC_GPU_PRESET). The [gpu] table in the input reference lists the keys and what each preset turns on.

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

Behaviour that differs from what you might expect. 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 an SCF run tens of iterations past a settled density, or never converge.

Why: convergence is driven by the density criterion (density_conv). The energy change between iterations has a noise floor, and a tight energy bound can sit below it, so it is met only by chance. For a closed-shell run with exact J and K the floor comes from the incremental Fock update and its integral screening (integral_thresh). ferric estimates it as 0.0125 * integral_thresh * nbf^2.25 + 0.5 * eps * |E| * sqrt(nbf) Ha. The estimate against the measured RMS of the energy change once the density has settled:

The estimate is good to a factor of a few and errs high for large systems; it was fitted on systems up to about 340 basis functions, and beyond that it is an extrapolation. With density-fitted J (and in the open-shell solvers, which rebuild J and K every iteration) only rounding remains, about 0.5 * eps * |E| * sqrt(nbf): measured 3e-13 for benzene/cc-pVDZ with RI-JK and 1.3e-12 for alkane C20/6-31G with RI-JK. The CLI default is loose (1e-3) by design.

What ferric tells you: when you set energy_conv below 0.3 times the estimated floor, the CLI prints warning: [scf] energy_conv ... is below the estimated noise floor of dE ... to stderr and the Python run_rhf, run_uhf and run_rohf emit a RuntimeWarning, once, before the SCF starts. A default energy_conv is never warned about. Nothing else changes: the threshold is not clamped and the SCF stops exactly where it otherwise would. A run that ends without converging on the dE criterion alone repeats the statement in RhfResult.convergence_summary. No estimate is made for functionals (grid noise), the LinK and COSX exchange builders, PCM / COSMO, polarizable embedding or non-default screening: those runs get no warning, which is not a statement that the tolerance is reachable. The |E| in the rounding term is estimated from the nuclear charges (ECP cores removed, about 5 % for all-electron molecules, 3.7x high for an iodine ECP case); the point-charge or field energy of an embedding is not included.

A symptom of the same floor: with a large point-charge field (point_charges= on run_rhf's closed-shell direct-J/K path) the number of SCF iterations can differ between two orderings of the same charges, because a change of a few ulps in the one-electron potential decides which iteration the dE test happens to pass. On rc7 a 6458-charge field took 64 iterations nearest-first and 111 in file order at energy_conv=1e-10. The potential itself changes only by roundoff (3e-15 to 1e-14 relative) and the converged energy by about 1e-12 Ha. ferric turns the incremental update off when energy_conv < 1e4 * integral_thresh or density_conv < 1e3 * integral_thresh (1e-8 and 1e-9 at the default integral_thresh), where all measured orders took 20 iterations; at the default energy_conv and density_conv every measured order took the same number (15 iterations on the field and orders tried; the tests do not assert iteration counts, so this is a measurement, not a guarantee for other fields). run_uhf and run_rohf use a separate open-shell incremental path and run_dft's default RI-J does not use the direct incremental path: for those the order changes the potential by roundoff too, but iteration counts were not measured. ferric does not reorder the charges.

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. 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/. 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 spilled RI-JK path sizes its dressing and DF-K read-back widths to the budget. Below a floor of about 200 aux rows for the dressing (raw stream, one unpacked block and one GEMM result) and two DF-K chunks for the read-back, a budget overshoots, with a warning naming the phase. The read-back is sized for 8 workers; a run with more can overshoot by the extra per-worker scratch, and says so.

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

• The SCF's three-index tensor is stored whole when n_aux × n_bf² × 8 bytes fit. When only the symmetric half (n_aux × n_bf(n_bf+1)/2 × 8 bytes, plus 64 unpacked aux rows of scratch) fits, it is kept in memory as that packed triangle and the RI-J build reads the packed rows directly. When even that does not fit, it is either spilled to a file in $TMPDIR (/tmp by default) and reread on every iteration, or recomputed on every iteration without writing anything. auto chooses between those two from a rate measured on the machine; [scf] jk_storage forces memory, disk or direct. Holding the K dressing to a small budget takes several passes over the raw tensor, and a run making many says so when a pass costs IO or integrals. The dressing and K are designed to be bitwise the same spilled or in core, and tests confirm that at OPENBLAS_NUM_THREADS=1 on one OpenBLAS kernel (Haswell). Threaded BLAS moves bits, and other kernels are untested. J is a separate matter: which storage raw gets can follow the shared ledger, so J can differ by a few ulp between runs that charge the ledger differently.
• A dressed tensor that fits the budget but not beside the resident raw tensor spills instead of stopping the job, unless jk_storage = "memory".
• A DF-guess pre-stage that runs out of memory is abandoned with a warning, and the SCF starts from its usual guess. That changes iterations and can select a different SCF solution, a larger effect than any bit-level one. An explicit jk_storage = "memory" that does not fit, and any memory failure under MPI, stop the job instead. A spill directory that is full also stops it, and df_increments has no fallback.
• The DFT grid AO cache is recomputed at every Fock build instead of stored. The energy is bit-identical. Under an explicit budget, XC takes its memory after the RI-J/K tensors, so it never pushes them out. An auto-detected budget first gives up 0.8 of the AO blocks XC expects to keep, so the raw tensor can drop to a packed or spilled tier while XC stays resident. If what is left cannot hold its per-build working set either, that set runs uncharged, a warning gives its size, and the run log records it. The job can exceed its budget by more than that figure, because other scratch is not charged either: danuglipron/def2-SVP warned of 0.46 GB and peaked 0.92 GB over a 15.3 GB budget. The working set scales with the thread count. The VV10 NLC cache has no fallback and is charged before RI-J/K, so a budget that cannot hold it alone stops the job before any integrals are built. It also leaves the MINAO guess's free-atom solves less room. If one runs out, that element's block of the guess is left zero, so the SCF starts from a partial MINAO guess rather than from hcore. If that was every element in the molecule, the density is zero and the first Fock matrix is hcore. A warning names the element and says which, the run log records it, and the SCF may take more iterations. Very near that limit, repeat runs of one input have differed in whether this happens. jk_storage = "memory" with a VV10 functional is held to the budget less the cache, so it refuses budgets that cannot hold the cache and the tensors together.
• 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.

GPU offload is opt-in and may not help

The CUDA backend is not in the PyPI wheel; it needs a source build with the gpu feature (see GPU (CUDA)), and even then it is off by default (FERRIC_GPU=auto or [gpu] mode turns it on). Consumer GeForce cards run FP64 at 1/32 of their FP32 rate or less, so offloading f64 contractions may be slower than the CPU once transfer time is counted; whether it helps depends on the card and the problem, and has to be measured on your own hardware.

The Kohn-Sham XC grid runs on the device only under mode = "on" (or gpu = "on" / "mixed"), never under auto: on a GTX 1080 every device XC build took 2.2-3x the CPU time and whole PBE SCFs ran about 2x slower, so auto keeps it on the CPU rather than slowing every DFT run on a GPU box.

Device memory is a separate pool ([gpu] memory_gb); the host [memory] budget_gb is not reduced by it. A contraction that does not fit the device pool runs on the CPU and is counted, not refused. FERRIC_GPU_TRACE=1 prints each fallback.

[precision] mode = "mixed" (or its deprecated alias [gpu] precision = "mixed") changes the numbers as well as the speed: the audit line, ferric.gpu_status()["precision"] and each result's precision provenance record (run_end.extra.precision in the run log, result.precision in Python) say which precision a run used, and the per-kernel error map is on the validation page. The one shipped mixed kernel is rimp2-energy (f32-resident B_ov; closed-shell and unrestricted RI-MP2 under the Coulomb, erfc and terfc operators); its measured energy errors are 3e-11 to 7e-9 Eh closed-shell Coulomb (total, with the OS and SS parts up to 5e-9 Eh each), 5e-10 to 1.4e-8 Eh unrestricted Coulomb (total, with the αα, ββ and αβ blocks up to 1.1e-8 Eh each) and 7e-11 to 6e-9 Eh for erfc and terfc on the systems listed there, graded as a measured error map, not Proven; a mixed run prints precision mixed in its [ferric] gpu: line, and its energies are not interchangeable with f64 results in a convergence or reference comparison. A tiny problem still pays the device upload when the mode is on, whatever the precision. gpu = "mixed" (or [gpu] preset = "mixed") is mode = "on" plus precision = "mixed", so it carries the same change in the numbers; a mode or precision setting that contradicts the preset is an error, not an override. A mixed_kernels list given next to a preset is kept as given instead of the preset's set (it may narrow it); a name this build does not ship is refused, and the opt-in xc-vxc and dfk-occ are admitted only when the list names them and mode resolves to "on" (so gpu = "mixed" with the list works, auto-mixed with it is an error).

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".

Outputs for recipes 1 and 3 come from running them as written (2026-09-23); recipe 2's output is from an earlier run. MEASURED marks numbers that came out of the program.

Prerequisite: a working ferric. The fast path is the prebuilt wheel (about a minute). A source build takes about 10 minutes, mostly compiling ferric:

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/ and scripts/ 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, so a job's resident memory can exceed it. Only a cgroup puts a hard ceiling on the process; without one, an overshoot can trigger the system-wide OOM killer, which may kill unrelated processes. 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

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 running 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:

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" (the analytic Hessian where it applies, otherwise 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.
• 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. 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 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:

Each RIFIT set in the table also loads under ORCA's <orbital>/c name (for example def2-tzvp/c), each JKFIT set under <orbital>/jk (for example aug-cc-pvtz/jk), and def2-universal-jkfit as def2/jk. STO-3G, the Pople sets and the pcseg family have no matching RI set; pair them with a def2 or cc RIFIT set of similar quality. Any other orbital 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, UHF and ROHF (exact J/K, or RI-J / RI-K with no external potential; an ECP is fine; basis up to f functions; RHF and UHF with exact J/K also allow an external potential, ROHF allows none on any route), and closed-shell KS LDA/GGA (pure or global hybrid, moving-grid Hessian). Other KS-DFT (UKS, ROKS, meta-GGA, range-separated hybrids) 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

Integrals and building blocks

With ints = ferric.AoIntegrals(mol, bs) (see Building blocks). AO matrices follow libint2's basis-function conventions, not PySCF's, so raw AO matrices differ between the two codes; energies and other invariant quantities agree.

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.

Only some drivers raise on an unconverged SCF. 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 64 sets plus aliases (see Bundled basis sets). BasisSet.from_bse_json(path) loads a Basis Set Exchange JSON file; there is no per-element basis dictionary.

The Pople sets use Gaussian's d convention, not PySCF's default. PySCF builds every molecule spherical (5D) unless you pass cart=True, and the flag covers the whole molecule. ferric types each shell:

The 6D/5D difference is not small: water/6-31G* differs by about 1.4 mHa between the two. For the 6-311G family ferric departs from the Basis Set Exchange, which marks the Na-Ar d shells Cartesian; ferric makes them spherical like every other element, so a molecule such as CH3Cl is 5D throughout and matches PySCF's default. The f shells the 6-31G* family carries on Sc-Zn are spherical (7F), as in Gaussian; cart=True would make them Cartesian too, so PySCF's single flag cannot reproduce 6-31G* on Sc-Zn.

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.

What is missing: no AMBER prmtop reader (go through OpenMM), and no periodic boundary conditions.

Before you start: relax the geometry with xtb

DFT on an unrelaxed structure wastes hours and often fails.

Measured on 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

An SCF has no reason to converge on a structure with a 0.9 Å C–H bond. xtb fixes it in seconds; DFT spends hours failing on it.

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.

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:

# 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: 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.

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=: without a force field, nothing holds the moving MM atoms' geometry together.

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: the 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 QM–MM Lennard-Jones term across a covalent cut is checked against OpenMM's own exclusion and 1-4 handling (ethanol cut at C–C, realistic LJ): energy and gradient agree to ~1e-18 Ha. With no bond path between the regions it is bit-identical to the plain every-pair sum.

Not validated: periodic boundary conditions (absent).

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 contributes nothing. [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 count a contribution twice.

Scope of the CLI section: it reads only a PQR, because it needs charges and geometry together, and does electrostatic embedding only. Polarizable sites, smeared charges, the MM force field and MM relaxation are Python-only.

Known limits

• No AMBER prmtop reader. Go through OpenMM, or build the arrays.
• No periodic boundary conditions.
• 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.
• QM/MM dispersion is the force field's Lennard-Jones only. D3/D4/XDM/VV10 are QM-atom-pairwise and never see the MM atoms. With mm_topology= set, the MM force field supplies Lennard-Jones between every QM real atom and every MM atom (link atoms and boundary charges carry none). Without it there is no QM–MM van der Waals term at all.
• QM–MM Lennard-Jones follows the force field's exclusion rules across the cut. Pairs 1-2 or 1-3 apart through the topology's bond list (the cut bond included) get no LJ, and 1-4 pairs are scaled by the topology's scale_lj_14 (AMBER default 0.5). The MM–MM terms use the same bond graph, so the cut bond has to be in the topology's bonds. QM–MM electrostatics do not use these exclusions or scale_coul_14. The embedding puts every MM charge in the QM Hamiltonian, and the boundary scheme (keep/delete-host/rc/rcd) handles the charges next to the cut.
• The MM crate is AMBER-form (harmonic bonds and angles, periodic torsions, Lennard-Jones, Coulomb), validated against OpenMM.

Export: cube files and NPZ bundles

ferric writes volumetric data as Gaussian cube files and per-molecule property bundles as compressed NumPy .npz archives. The Python functions evaluate basis functions, the density and orbitals at arbitrary points; the CLI [rpa] section writes the NPZ bundle and PDEP eigenpotential cubes. The array schema is in crates/ferric-export/theory.md.

Units

These differ from Molecule.from_xyz_string, which reads Å, and from the periodic and QM/MM Python entry points, which take Å. The functions below take Bohr throughout.

Cube files from Python

write_cube writes the total density or one restricted MO on an axis-aligned grid. Without origin and shape the grid is the atoms' bounding box plus padding. It returns the (origin, shape) used. what="mo" needs mo_index (0-based, energy order) and a restricted RHF or RKS result. Orbital signs are the SCF's arbitrary phase.

import ferric

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
""")
bs = ferric.BasisSet.bundled("sto-3g")
rhf = ferric.run_rhf(water, bs)

origin, shape = ferric.write_cube("density.cube", water, bs, rhf, spacing=0.3)
print(shape)
origin, shape = ferric.write_cube("homo.cube", water, bs, rhf,
                                  what="mo", mo_index=4, spacing=0.3)
print(shape)

[28, 38, 32]
[28, 38, 32]

Values at points

eval_ao_on_points, eval_density_on_points and eval_mo_on_points take an (N, 3) array of points in Bohr.

result is an RhfResult or DftResult. eval_mo_on_points refuses open-shell results. All three check the allocation against the memory budget before evaluating.

import numpy as np

pts = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])    # Bohr
print(np.round(ferric.eval_density_on_points(water, bs, rhf, pts), 4))
print(ferric.eval_ao_on_points(water, bs, pts).shape)
print(ferric.eval_mo_on_points(water, bs, rhf, [4], pts).shape)

[10.2392  0.8103]
(2, 7)
(2, 1)

NPZ bundles from the CLI

A pdep-rpa run with [rpa] export_npz = "<path>" writes the bundle. The compute_* keys then default to true; the full key list and defaults are in Input file (TOML).

[molecule]
xyz = "testdata/molecules/water.xyz"

[basis]
name = "cc-pvdz"

[method]
kind = "pdep-rpa"
task = "energy"

[rpa]
auxbasis = "cc-pvdz-ri"
n_quad = 20
trunc_thresh = 1e-4
export_npz = "h2o_pdep_rpa.npz"
compute_esp = true
compute_polarizability = true

A requested property that cannot be computed fails the run unless allow_partial_npz = true. Read the archive with numpy.load:

import numpy as np
d = np.load("h2o_pdep_rpa.npz")
print(d.files)

Keys written for water/cc-pVDZ with the input above (N atoms, nao = 24):

mo_coeffs, orbital_energies, orbital_centers and orbital_spreads are written for restricted references only. c6_iso and c6_aniso are per-atom pair tensors under a partition convention. They do not sum to the molecular C6; read c6_molecular_iso. alpha_ct + alpha_atomic.sum(0) == alpha_tensor. The definitions of alpha_atomic and alpha_ct are in the compute_alpha_atomic row of the input reference.

PDEP eigenpotential cubes

Setting export_eigpot_prefix in [rpa] writes <prefix>_eigpot_NNN.cube for the leading export_eigpot_count eigenpotentials (default 10; the count is capped at the number retained). cube_spacing (default 0.2 Bohr) and cube_margin (default 4.0 Bohr) set the grid. examples/benzene-pdep-rpa-export.toml is a complete input.

Golden paths: ferric in applications

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

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

Golden path A — Reaction energetics of a cation cascade

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

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

Step 0 — gate before you spend hours

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

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

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

Step 1 — build the species, and check parity

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

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

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

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

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

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

Step 1b: relax with xtb before any DFT

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

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

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

Measured on the same molecule, method and basis:

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

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

Step 2 — optimize each intermediate

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

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

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

[dft]
functional = "PBE"

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

[optimize]
max_steps = 60

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Step 3 — confirm the ordering survives the functional

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

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

Step 4 — properties at the key intermediate

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

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

Step 5 (optional) — barriers

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

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

What this cannot tell you

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

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

Reading the output

Before publishing any number

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

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 about 10 minutes, mostly compiling ferric.

pip install ferric
git clone https://github.com/mgoldey/ferric && cd ferric   # for examples/, testdata/, 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/ and scripts/ 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, instead of building from source.

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 §4 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, so the setting affects correctness as well as speed.
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 slow each other by more than half. 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.

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 run. 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 the same as "the job produced nothing". Test the filter against known-good output first.

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 QM is the product.

Its analytic Hessian covers RHF, UHF, ROHF and closed-shell KS LDA/GGA (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 matter when picking 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; a value below the estimated noise floor of the energy change draws a warning (see sharp bits).

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 / TPSS 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, up to 75 × 110 points per atom, pruned by default (grid_prune = "nwchem-radial": NWChem angular regions, and 50 radial shells for H and He, 66 for Li–Ne, 75 for heavier atoms) with a 1e-14 weight cutoff (grid_weight_cut, ORCA's WeightCut). That is 42-47% fewer points than the flat 75 × 110 grid on organic molecules, with energies within 2e-7 to 4e-6 Ha of it on PBE/def2-SVP (the flat grid's own error against a 99 × 590 grid is 4.5e-4 to 3.3e-3 Ha on C10–C32 alkanes).
• The flat 75 × 110 grid is grid_prune = "none" with grid_weight_cut = 0. "nwchem" prunes the angular grid only (about 23% fewer points, examples/water-pbe-pruned-grid.toml); "defgrid1"–"defgrid3" are ORCA 6.1's DEFGRIDs, with their own sizes: fewer radial shells than the default (H 30/35/44, Li–Ne 38/44/53, Na–Ar 46/54/63 against 50/66/75) but more angular points (peak orders 194/302/590 for heavier atoms and 110/194/434 for H and He, against 110). defgrid2, ORCA's default, moved PBE/def2-SVP energies by 1e-3 to 4e-3 Ha from the flat grid. Setting only grid_radial / grid_angular keeps the default scheme at that size. Optimization, frequency (including the analytic moving-grid Hessian) and TDDFT runs use the grid the SCF ran on.
• SCAN is the exception to "pruning costs nothing". On water/6-31G its gradient error against a 99 × 590 grid is 2.1e-4 Ha/Bohr on the flat grid and 5.6e-4 on the default, because the hydrogens lose radial shells. r2SCAN (3.7e-4 flat, 3.6e-4 default) and PBE (unchanged) are not affected. For SCAN, use grid_prune = "nwchem" or a larger 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/TPSS energies are pinned against PySCF (ferric-scf/tests/dft_scan.rs); the SCAN/r2SCAN closed-shell gradients against PySCF with grid response (dft_gradient_mgga.rs). TPSS gradients have no PySCF pin. Other meta-GGA variants (deorbitalized, +VV10, hybrid meta-GGA) are refused, not approximated. The \( V_{xc} \), \( f_{xc} \), VV10, Becke-weight and gradient formulas are derived in DFT derivations.

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, plus the Axilrod–Teller–Muto three-body term when its scale s9 is nonzero. s9 defaults to 0 (off), the reference default for a named functional (simple-dftd3 without --atm).

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, s9=0.0) for the correction alone; s9=1.0 adds the three-body term. The CLI and run_dft(dispersion=...) use s9 = 0.

Scope. The damping parameters are fitted per functional, so the key is refused on any method other than ksdft. Tabulated functionals, each equal to simple-dftd3 1.6.0's rational-damping entry: BLYP, BP86, B97-D, PBE, PBE0, revPBE, B3LYP, TPSS, TPSS0, B2PLYP, BHLYP, SCAN, r2SCAN, ωB97X. M06-2X has no D3(BJ) fit (only zero-damping D3(0)) and is an error. task = "optimize" works (the D3 gradient is implemented; ferric-d3/tests/gradient_vs_fd.rs), and so does task = "frequencies" on RKS, UKS and ROKS: see Frequencies with dispersion.

Accuracy. Energy and analytic gradient agree with simple-dftd3 1.6.0 with the three-body term off and on, on nine molecules from water to benzene and CHBr3: 1.2e-13 Ha and 5.5e-13 Ha/Bohr on the totals, 2.6e-15 Ha on the three-body part (ferric-d3/tests/vs_reference_dftd3_atm.rs). The three-body term is 0.1% of the two-body energy at benzene and grows with size. Unlike simple-dftd3, no real-space cutoff is applied (it drops triples beyond 40 Bohr), and the triple sum costs O(N³).

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. For a linear radical whose state breaks the cylindrical symmetry (²Π: OH, NO, CH), a rotation of the electrons about the molecular axis is an orbital rotation that leaves the energy unchanged, so the orbital Hessian has a null mode; only the anisotropy of the XC grid gives it a small eigenvalue of either sign. The UKS and ROKS Z-vectors detect it on linear molecules and solve on its complement, which is exact for a property that is symmetric about the axis. On OH, NO and CH (STO-3G and 6-31G, PBE and PBE0) the gradient agrees with finite differences of the full pipeline to 7e-12–5e-11 Hartree/Bohr. A state with any other negative or near-null orbital-Hessian direction is still refused (the UKS-PBE CH state is a saddle point and is refused). 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). With COSX exchange the overlap fit must be off (cosx_overlap_fit = false); the UKS-PBE0 + COSX gradient on NH2/STO-3G agrees with finite differences of the full pipeline to 2.3e-11. Other setups are refused rather than given an approximate gradient. Python's run_dft is closed-shell only; open-shell frequencies with MBD@rsSCS run through the CLI or run_frequencies.

Frequencies with dispersion

What it is. Harmonic frequencies on the dispersion-corrected surface E(KS) + E(disp), on an RKS, UKS or ROKS reference. The Hessian is the central difference of the corrected analytic gradient: at each of the 6N displaced geometries the SCF of that reference 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, with the Z-vector of the reference solved (RKS, UKS or ROKS), 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 KS SCF: ksdft, or rhf/uhf/rohf with [dft] functional. 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; reference="uhf" or "rohf" selects UKS or ROKS.

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.

On open shells (PBE/STO-3G, bent NH2 doublet on UKS, HCO doublet on ROKS), the D3(BJ) part agrees with 4-point second differences of the D3(BJ) energy to 3.7e-9 and 6.5e-9 Hartree/Bohr² and the frequencies to 9.4e-6 and 3.0e-5 cm⁻¹; the D3 Hessian obeys the translational sum rule to 1e-11 and the rotational one (with the D3 gradient) to 3e-9. The MBD@rsSCS part agrees along three fixed directions with 5-point second differences of the full SCF + MBD energy to ≤ 4.8e-8 Hartree/Bohr² for UKS NH2 and CH2 (triplet) and ROKS HCO and CH2, against curvatures of 6e-5–2e-4; without the orbital-relaxation term the misses are 8.3e-7–1.6e-5. A correction that does nothing reproduces the plain UKS or ROKS finite-difference Hessian bit for bit.

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. IEF-PCM has an analytic nuclear gradient (ferric_pcm::pcm_gradient: the tessera-charge potential term plus the S/D/area response term, switching-function derivatives included; equations in crates/ferric-pcm/theory.md §8), added to the RHF, UHF, ROHF and KS gradients whenever RhfConfig.pcm is set. It matches PySCF IEF-PCM nuc_grad_method() to 7e-11 Ha/Bohr on PySCF's 302-point cavity and finite differences of the solvated SCF energy to 3e-7 Ha/Bohr. It needs the cavity built from the molecule (an injected PcmConfig::cavity is refused). There is no PCM Hessian: frequencies with [pcm] are refused. COSMO has no gradient term and no Hessian: every analytic-gradient entry (the SCF gradient dispatchers, optimize_geometry*, optimize_qmmm, finite-difference frequencies, and CLI task = "optimize" or "frequencies" with [cosmo]) refuses a COSMO SCF with "COSMO has no analytic gradient; use [pcm] ... or finite differences". COSMO is energy-only; use [pcm] for geometry work. 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, examples/nh2-rohf-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 or with density-fitted J and/or K (df_j_aux / df_k_aux, Coulomb operator, no external potential on the fitted routes, a well-conditioned fitting metric), no solvent or polarizable embedding, and a basis up to f functions (RHF and UHF also with an ECP, on the fitted routes too, and with exact J/K also with fixed external point or smeared charges, or a uniform field), and for ROHF with exact four-centre J/K or RI-J / RI-K, an ECP allowed, no external potential 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; for ROHF, one solve over the closed-virtual, open-virtual and closed-open rotations). This needs a libint2 with second derivatives (the build scripts/install-libint.sh installs).
• Analytic on the moving grid for closed-shell Kohn-Sham (ksdft, or rhf with a functional): LDA and GGA functionals, pure or global hybrid, exact or density-fitted J and K, Becke weight response included, so it is the Hessian of the surface the KS gradient differentiates. One SCF plus a coupled-perturbed KS solve.
• Central finite differences of the analytic gradient everywhere else (UKS, ROKS, meta-GGA, VV10, range-separated hybrids, KS with an ECP or an external potential, pruned grids, [dft] dispersion, RI J/K with an external potential or an ill-conditioned fitting metric, solvent, polarizable embedding, g functions, ROHF with an external potential): 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, the moving-grid Hessian for closed-shell KS), "analytic" (the Hartree-Fock analytic Hessian; an error where it does not apply, including every KS reference), "analytic-moving-grid" (the closed-shell KS one; an error where it does not apply) or "fd". The output names the one that ran (Hessian = analytic, analytic-moving-grid or finite-difference; Python .hessian_source).

The moving-grid KS Hessian agrees with "fd" up to the finite-difference error: ≤ 5e-7 cm⁻¹ against Richardson-extrapolated FD at a tight SCF (density_conv 1e-11) on water / STO-3G, for one SCF plus a coupled-perturbed KS solve instead of 6N gradients. Against the default "fd" (step 5e-3 Bohr, default SCF), which "auto" ran for KS before, the wavenumbers move by ≤ 0.071 cm⁻¹ (water / STO-3G, LDA, PBE, B3LYP). Under "auto", open-shell KS, meta-GGA, VV10, range-separated hybrids, ECPs, external potentials, pruned grids, [dft] dispersion, COSX, implicit solvent, fractional occupations, MOM, and systems whose response AO table (basis functions × grid points × 4 for a GGA, 8 bytes each) exceeds 2 GiB (on the default 75×110 grid, a GGA above about 8,100 basis functions × atoms, e.g. about 30 atoms at def2-SVP) run finite differences; "analytic-moving-grid" refuses them before the SCF.

The choice is made before the SCF. An error the analytic path raises after the SCF (for ROHF: the stored spin Fock matrices differing from those rebuilt from the final densities by more than 1e-4 Ha, which a loose density_conv can cause; a response solve that does not converge; non-positive orbital curvature; a first-order rotation above 10 rad/Bohr) is returned as an error under "auto" too, never retried with finite differences; the same holds for the moving-grid KS Hessian (a non-converged SCF included).

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. hessian="analytic" (or "auto", which falls back to finite differences where the analytic Hessian does not apply) swaps in the analytic RHF Hessian (with xc, "auto" takes the moving-grid KS Hessian), which needs no displaced gradients; "analytic" raises for KS xc, point charges / an external field, ECPs or a basis above the libint second-derivative limit, delta= is only valid with "fd", and result.hessian_source says which ran. The default is "fd". 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; half that when stored as the packed symmetric triangle). When it does not fit, J is built by Split-RI-J, which stores no tensor and on danuglipron PBE/def2-SVP took 81 s for the whole SCF against 46 s with the tensor in memory and about 590 s spilling or recomputing it. A tensor K still needs is spilled to disk or recomputed each iteration, whichever a measurement on the machine says is cheaper ([scf] jk_storage). 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. The 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_grid_schedule = false (default). true runs the early SCF iterations on a coarse pruned { radial = 25, angular = 110, prune = "sgx" } grid and switches to cosx_grid once the largest density-matrix change between iterations drops below 1e-3, restarting DIIS at the switch. Convergence is only accepted after a full iteration on cosx_grid, so the converged energy is the cosx_grid energy; the final pass is unchanged. On water/cc-pVDZ the coarse grid has 4062 points against 10866, and the scheduled SCF reaches the unscheduled energy to 1e-13 Ha in 15 iterations (8 coarse, 7 on cosx_grid) against 12. Its wall time has not been measured. RHF/UHF and their Kohn–Sham variants; ROHF refuses 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 depends strongly on the molecule (large on water, nil to negative on ethane), so the factor quoted in the literature does not carry over.
• 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). Geometry tasks (optimize, frequencies) use this grid, 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, UHF, UKS, ROHF and ROKS 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 open-shell exchange term is differentiated per spin; an exact-exchange gradient paired with the same COSX energy misses those finite differences by 4e-5–2e-4 Ha/Bohr (HO2 and NH2, B3LYP and ROHF, STO-3G). The MBD@rsSCS gradient's orbital-response term on a UKS or ROKS reference uses the same COSX derivative. 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). Refused for gradient tasks before the SCF runs: fitted COSX with a Kohn–Sham functional (RKS, UKS, ROKS; the fitted energy's response needs the nuclear derivative of the XC Fock matrix) and fitted COSX for ROHF (no restricted-open-shell response).

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. For systems like alkane chains, different initial guesses converge to different SCF solutions, each in its own 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.

Secondary-basis guess (Q-Chem BASIS2)

[scf] guess_basis = "6-31g" (Python guess_basis="6-31g" on run_rhf, run_uhf, run_rohf, run_dft, run_ksdft) converges the same SCF (method, functional, charge, multiplicity and settings) in the named basis first, projects the result into the target basis, and starts the target SCF from it. guess_projection picks how the projection is done:

• "fock" (default; Q-Chem FOPPROJECTION): the small-basis spin densities are projected, D = P D_small Pᵀ with P = S_bb⁻¹ S_bs, and the target Fock built at that density is diagonalised. Closed shell, that Fock build is the target SCF's first iteration. Open shell, the guess Fock is h + J − K_σ with no XC term, which is the form ferric's MINAO open-shell guess also uses.
• "overlap" (Q-Chem OVPROJECTION): the occupied MOs are mapped, C = S_bb⁻¹ S_bs C_small, and re-orthonormalised in the target metric. For ROHF, the doubly occupied block is orthonormalised first and the singly occupied block is made orthogonal to it.

S_bb⁻¹ is the lindep-filtered pseudo-inverse from the canonical orthogonalizer. A secondary basis missing an element, a mismatched ECP core count, or a projected occupied space that is near-singular in the target basis are errors. A secondary basis larger than the target is allowed, with a warning. The guess is refused with symmetry = true and on periodic systems. An explicit starting density or MOs take precedence, as they do over guess = "minao". With guess_basis set, guess = "minao" | "sad" | "hcore" selects the starting point of the secondary SCF only; the target starts from the projection. The dual-basis energy correction (Q-Chem DUAL_BASIS_ENERGY) is not implemented.

A secondary SCF that does not converge never stops the run. The target SCF then starts from its default guess, a warning is printed, and the report says converged = false, why, and used_as_guess = false. Projecting the unconverged last iterate was tried and rejected: on CN UHF/aug-cc-pVDZ with an STO-3G secondary (which oscillates for 200 iterations), the target started from it failed to converge as well, for both projections, while the default guess converges in 22 iterations. The result reports the secondary SCF's iterations, energy and convergence (Python: result.guess_basis). The report is kept when a later stage supplies the final result: a stability descent, a later ladder rung, or the exact stage after a DF pre-stage.

Cost on re-solves. Every solve without an explicit starting point runs the secondary SCF again: each geometry of an optimization and each displaced geometry of a finite-difference Hessian. A ladder's later rungs, the stability descent's re-solves and the cDFT λ loop start from a carried density or MOs, so they do not repeat it (cDFT builds the guess once, before the λ loop).

The guess can select the state. The converged energy is not guaranteed to equal the default-guess energy. Like any guess, it decides which stationary state the SCF falls into, in either direction:

• NO UKS-PBE/aug-cc-pVDZ: MINAO −129.7833808133, BASIS2 −129.7833798480 Ha. These are two states 0.97 µHa apart, and BASIS2 lands on the higher.
• He₂ at 0.3 Å RHF/aug-cc-pVTZ: BASIS2 and the core guess give −3.2458864766 Ha, while MINAO gives −3.2388590827 Ha. Here BASIS2 lands 7 mHa lower.

Use check_stability (and stability_descent) when the state matters.

What it buys, measured in ferric (build 352eabab9, density_conv 1e-8; target-SCF iterations, with the secondary SCF's own iterations listed separately; converged energies agree to ≤ 3e-11 Ha in every row; wall times were taken on a loaded machine and are not quoted):

On closed shells the saving over MINAO is 0 to 1 target iterations, so there is no wall-time win. On open-shell radicals it saves up to 7 target iterations (CN), but the secondary SCF costs about as many iterations as it saves, or more; they are cheap iterations in a small basis. The value is a better starting point on hard open-shell cases, not a faster closed-shell SCF. Fock projection beats overlap projection by 0 to 1 iterations.

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.

A tree-fold reduction would be deterministic but would produce different bits, since floating-point addition is not associative, so the sum must run in ascending order.

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.

DFT derivations

The exchange-correlation pieces of the Kohn–Sham solver, each stated in the code's conventions and proved by a sympy script that exits nonzero on a mismatch. The full equations, function locations and test cross-references are in crates/ferric-dft/theory.md. Each script also checks mutated formulas (a dropped factor of 2, a wrong libxc packing slot) and requires them to fail.

Conventions

Closed shell, total density matrix \( D \), real AOs \( \chi_\mu \):

\[ \rho = \sum_{\mu\nu} D_{\mu\nu}\chi_\mu\chi_\nu,\qquad \sigma = |\nabla\rho|^2,\qquad \tau = \tfrac12 \sum_b \sum_{\mu\nu} D_{\mu\nu},\partial_b\chi_\mu,\partial_b\chi_\nu . \]

Spin-polarized variables use libxc's packing: \( (\rho_\alpha,\rho_\beta) \), \( (\sigma_{\alpha\alpha},\sigma_{\alpha\beta},\sigma_{\beta\beta}) \) with \( \sigma_{\alpha\beta} = \nabla\rho_\alpha\cdot\nabla\rho_\beta \), and \( (\tau_\alpha,\tau_\beta) \). \( v_\sigma = \partial(\rho\epsilon)/\partial\sigma \), not twice it. Meta-GGAs are \( \tau \)-only (no density Laplacian).

What is proved

The numerical checks of the same identities (finite differences of the energy, PySCF references) are listed per section in theory.md.

Known open issue

The VV10 AO matrix \( V_{nl} \) is not exactly \( \partial E_{nl}/\partial D \): on NH₂/6-31G the directional derivative misses by \( 3.6\times10^{-6} \) relative (\( -1.31\times10^{-7} \) absolute), independent of the step size, while ωB97X without VV10 agrees to \( 8.6\times10^{-14} \). The per-point partials are proved exact, so the residual lies elsewhere on the path from \( D \) to \( V_{nl} \); the cause is not identified. The test is a strict xfail (crates/ferric-python/tests/test_xc_blocks.py).

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. The working equations of every variant on this page (spin factors, frozen core, operator and unit conventions), with the code path behind each, are in crates/ferric-mp2/theory.md; the identities it relies on are checked by scripts/mp2_variant_identities.py.

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.

The reference SCF: exact, RI-JK or RIJCOSX. The auxbasis above fits the correlation step only. The reference SCF is a separate choice, scf_jk on every MP2-family Python driver ([scf] jk in the CLI, task energy): "exact" (four-centre J and K, the default), "rijk" (RI-J + RI-K) or "rijcosx" (RI-J + COSX K). The fitted modes use jk_auxbasis ([scf] jk_aux), a JK-fit basis with default def2-universal-jkfit. It is never the MP2 auxbasis, because an MP2-fit set used for K costs mHa. rs-mp2-rpa is the exception: its default is "rijk". Every result echoes scf_jk and jk_auxbasis. An unknown mode, jk_auxbasis with "exact", or a k_builder the mode would override is an error (so is any k_builder on rs-mp2-rpa's RI-JK default, in the CLI as in Python). A jk_auxbasis whose name lacks "jkfit" draws a warning. oo-rimp2 accepts only "exact": orbital optimisation rebuilds exact four-centre J/K at every orbital step, so a fitted reference would change only the starting orbitals; threading RI-JK through the OO Fock builds is not implemented. In the CLI a fitted reference cannot be combined with [scf] df_increments (an exact-J/K scheme) or an explicit df_guess = true, and it skips the default DF-guess pre-stage. At aug-cc-pVTZ the four-centre SCF is most of an MP2 single point, so a fitted reference is the usual choice for interaction-energy benchmarks. The fitting error is small and largely cancels in a counterpoise-corrected interaction energy. Measured with PySCF DF-MP2 (frozen core, aug-cc-pvXz-ri for MP2) on the S66 water dimer:

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, so attenuated MP2 has no long-range dispersion: 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, 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's implementation is Proven (narrow): the damped VV10 half, the terfc-attenuated MP2 half and the erfc control each match an independent PySCF/numpy reference. Its parameterization is unbenchmarked: there is no published MP2-V total energy to compare with, only an aug-cc-pVTZ S66 subset run (n=5, H-bonded/mixed, S66x8 geometry, no CP-anchored reference), 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.

Fitting metric. Each attenuated energy fits its integrals in the metric of the same operator, \( (P|\hat w|Q) \), so the fit is variational for that operator. Robust (Dunlap) fitting is required only when the metric operator differs from the integral operator.

The range split is not additive. MP2 is quadratic in the integrals, so the erfc and erf energies at the same ω do not sum to the Coulomb MP2 energy; the difference is the short-range × long-range cross term. RS-MP2 + LR-RPA (below) is built on it.

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; 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 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 amplitude equations are solved by preconditioned conjugate gradients to a relative residual of 1e-6 by default. The energy is the stationary Hylleraas form, so its solver error is second order in that residual. Using that form with a loose stop comes from the method's authors: Wang et al. (JCTC 19, 7577, 2023; Eq. 10 and Sec. 3.1) stop at an absolute ‖r‖ < 2.5·√eps, and the 2025 follow-up (ChemRxiv preprint of JCTC 21, 10910) at 5·√eps. ferric uses a fixed relative tolerance instead and gates it on a per-run bound, 3‖PJ‖²·relres²/λ_min with λ_min ≥ 2(ε_LUMO − ε_HOMO). The bound holds for the true residual, up to rounding, so it is evaluated on ‖J + At‖ from the energy's own matvec, not on the solver's recursive residual. If the bound at the solver's stop exceeds 1e-10 Ha, the solve continues from where it stopped to a tighter tolerance, never past 1e-11; with no positive orbital gap the bound is undefined and the solve runs at 1e-11 from the start. Running out of iterations before the tolerance being sought, first or tighter, stops the run with an error; cg_max_iter and an explicit cg_rtol are fields of the Rust AmplitudeLmp2Config, not CLI or Python options. The certificate (result.local in Python, "local" in the CLI JSON) records cg_rtol, cg_rtol_source, cg_tightened, cg_relres, cg_true_relres, cg_energy_error_bound, cg_energy_error_target and cg_energy_target_met, and solver_warnings says if a converged solve did not meet the target, because the target needs more than the 1e-11 floor or the true residual drifted from the recursive one. The bound limits solver error only, not the truncation error of eps.

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.

Schwarz triple cut (schwarz_skip, default 1e-5) applies to the Coulomb, erfc and terfc kernels alike. A shell triple is skipped when sqrt(max (P|P)) * sqrt(max (ab|ab)) < schwarz_skip, both factors from the run's own operator; schwarz_skip = 0 evaluates every triple. For terfc the factors come from the table engine's own diagonals and carry a 1e-5 relative inflation for table error, and the cut is only available inside the checked box: orbital shells up to g, r0 between 0.3 and 2.5 Å and r0·ω ≤ 8.9 (outside it the run errors and names schwarz_skip = 0.0). Measured for terfc (n-alkanes, 6-31G with cc-pVDZ-RI, eps = 5e-5, default maps, plain Schwarz, 2026-10-11, a loaded machine so ratios only, order alternated, two repeats):

At 1e-5 the cut skips 55% of the shell triples at C12 and 67% at C16 (51% and 63% at 1e-6). The pair gate (gate_cal = 0.7) leaves the C12 result unchanged. At 1e-6 the energy changes by 1e-10 (C12 linked) to 1.8e-9 (C16 NCP).

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. Choose the aux basis for the accuracy you need.

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.

Particle-particle ladder: dense or RI

The particle–particle ladder Σ_cd W_abcd τ_ij^cd of closed-shell CCSD can be evaluated two ways, chosen by [cc] ladder (Python run_ccsd(ladder=) / run_ccsd_t(ladder=)):

• dense holds the chemist (ab|cd) block in core and builds W_abcd every iteration: n_v⁴ memory, a few times over.
• ri holds neither the (ab|cd) block nor W_abcd. It rebuilds the (ac|bd) entries from the RI 3-index tensors every iteration in a_blk·n_v³ panels: at most one n_v⁴-sized panel at full width (dense holds four n_v⁴ buffers at its peak), and narrower panels when the budget requires. It is the same algebra over the same tensors, so it is exact: on water/cc-pVDZ |ΔE| and max|Δt| are ≲ 3e-17 (round-off), and the alkane energies below agree to all 10 printed digits.
• auto (the default) runs dense whenever its memory plan fits the budget and ri otherwise, so it is never refused where ri would run. When it falls back to ri it says so on stderr. The route taken is printed (ladder = dense (auto)), logged, and returned (CcResult::ladder_route; Python CcResult.ladder_route).

Reproducibility. Under auto, a closed-shell CCSD job whose dense plan exceeds the memory budget runs on the RI ladder instead of being refused. Because auto decides against the budget, a run without an explicit [memory] budget_gb (auto-detected as 0.8 × the RAM available at start-up) can take different routes on the same input from one run to the next. The energies agree to round-off either way; the time does not. Set ladder explicitly, or fix the budget, for reproducible timings.

Measured 2026-10-26, alkane/cc-pVDZ (cc-pvdz-ri), OPENBLAS_NUM_THREADS=1, RAYON_NUM_THREADS=4, routes interleaved one per process, load average 3–6 from other jobs (benchmarks/harness/examples/ri_ccsd_memory.rs). RI ran at full block width (a_blk = n_v, the default-budget choice). Per-iteration time is (time for k iterations − setup time) / k, repeated 4× per route for water, 3× for C3H8 and 2× for C4H10/C5H12; ratios are of the medians. Peak RSS is one converged run per route.

RI is 1.3–1.5× slower per iteration from n_v ≈ 70 up and uses 2.1–2.3× less memory, which is why auto takes it only when dense does not fit. The RI route still holds the n_o·n_v³ (ia|bc) block in core, so it lowers the n_v⁴ ceiling but does not remove every large tensor.

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

The equations each solver implements, with file and line, are in crates/ferric-cc/theory.md.

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 once; a test pins this.

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.

CAS-CI

ferric runs complete-active-space configuration interaction (CAS-CI) from a closed-shell RHF reference. It is a Rust library routine only: the CLI and the Python bindings do not expose it.

The Hamiltonian, the determinant phase convention and the Slater-Condon rules are derived in crates/ferric-ci/theory.md. The phase conventions are checked by scripts/casci_slater_condon_phases_sympy.py.

Scope

With no inactive orbitals the method is FCI in the basis. A call that fails the bookkeeping check, 2 * active_start + n_alpha + n_beta == nelec, returns an error. So does a non-restricted reference, n_active == 0, or n_active > 63 (the u64 string limit).

Energy

$$ E_\text{CASCI} = E_\text{core} + \lambda_\text{min}(H_\text{act}), $$

where $E_\text{core}$ is the nuclear repulsion plus the energy of the active_start doubly occupied inactive orbitals, and $H_\text{act}$ carries the inactive-core Fock potential in its one-electron part. e_active is e_total - e_core, the quantity PySCF reports as CASCI.e_cas.

Library use

use ferric_ci::{run_cas_ci, CasCiConfig};

// `mol`, `prep` and `rhf` come from a converged restricted solve_rhf
// with exact J/K (no density fitting), as in the validation tests.
let cfg = CasCiConfig {
    n_active: 6,
    n_elec_active: (3, 3),
    active_start: 4,          // MOs 0..4 are doubly occupied and inactive
    ..Default::default()
};
let res = run_cas_ci(&mol, &prep, &rhf, cfg)?;
println!("{res}");            // energy, determinant count, iterations

CasCiResult holds e_total, e_active, e_core, n_determinants, ci_vector (index ia * n_beta + ib), converged and iterations. Defaults: conv_thresh = 1e-9, max_subspace = 50, max_iter = 200.

Limits

• The active-space integrals come from a dense AO ERI tensor, O(nbasis⁴) memory. Use it on minimal and double-zeta basis sets.
• The sigma build is an O(N_det²) loop over determinant pairs. It does not use the string-based Knowles-Handy/Olsen algorithm.
• The Davidson subspace needs 2 * max_subspace * N_det * 8 bytes. memory_budget_bytes bounds the integral transform and the subspace; None resolves through FERRIC_MEM_BUDGET_GB, then 0.8 of available RAM, then a 2 GiB floor.

Validation

Compared with PySCF mcscf.CASCI on an exact-integral RHF with the same orbital window, and with a dense diagonalization of the full active-space Hamiltonian:

Each system is also run with the window moved down one orbital. The test tolerance is 3e-9 Ha. The stretched N2 reference is the symmetric RHF, which is internally unstable toward a symmetry-broken RHF 0.19 Ha lower; both codes use the symmetric state.

The sign convention is checked three ways, because energy comparisons cannot detect a diagonal $\pm1$ similarity of $H$:

1. tests/slater_condon_brute_force.rs compares general_element, sigma and the diagonal with an independent bit-loop second-quantized $H$ on random asymmetric integrals, CAS(4,4) to CAS(6,6), odd and even $N_\alpha$. The largest difference is 8.9e-15.
2. scripts/casci_slater_condon_phases_sympy.py builds the Jordan-Wigner operators symbolically, verifies the anticommutation relations, and shows that a line-for-line port of general_element equals the explicit $H$ as an exact linear form in $h_{pq}$ and $(pq|rs)$ for four $(n, N_\alpha, N_\beta)$ sectors. Four sign mutations and the interleaved orbital order each fail.
3. The PySCF energies above fix the spectrum.

The grade is on Capabilities and validation. Sources: crates/ferric-ci/tests/validation_casci.rs and scripts/validation/gen_casci.py.

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). It runs the Lanczos (default) and Davidson eigensolvers, PDEP truncation (trunc_thresh), the Dense and Laplace χ₀ backends and run_diagnostics (full-basis U-RI-dRPA check). It refuses chi0_sparsity = "boys", and "auto" when that resolves to Boys for the molecule: the Boys-screened χ₀ is closed shell only.
• 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

The equations and conventions of this section, the BSE-TDA and the RPAx polarizability, with the code paths and tests behind each, are in crates/ferric-gw/theory.md; the identities it relies on are checked by scripts/gw_self_energy_sympy.py.

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. τ-only meta-GGAs (SCAN, r2SCAN, TPSS) take the τ kernel; 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 8.1e-4 eV on the default pruned grid (2.8e-6 eV for B3LYP/aug-cc-pVDZ on the reference's flat grid), test bar 1e-3 eV, and water/6-31G r2SCAN matches on the flat 75 × 110 grid (ferric-tddft/tests/validation_tddft.rs).

Meta-GGA response needs a finer grid. A TDA/TDDFT run on a meta-GGA (SCAN, r2SCAN, TPSS) with no [dft] grid key (Python: no grid_* kwarg) runs its reference SCF and its f_xc kernel on 150 × 434 with the default nwchem-radial pruning and 1e-14 weight cut. Every other functional keeps the production 75 × 110 grid. Any grid key replaces the meta-GGA default. The run prints the grid it used (CLI: the TDDFT grid line; Python: TddftResult.grid).

On the same grid, ferric and PySCF agree: fed ferric's own coordinates and weights, PySCF reproduces ferric's r2SCAN excitation energies to 6e-8 eV and its SCF energy to 1e-12 Ha. The grid itself is the error. Largest deviation of the five lowest TDA roots from a flat 150 × 590 grid (Casida roots follow the same pattern):

All values in eV. 150 × 434 is the smallest grid tested with every r2SCAN case at or below 1e-4 eV (H2CO: 9.96e-5 TDA, 1.02e-4 Casida); the margin is comparable to the reference's own distance from flat 200 × 590 (2-4e-5 eV).

PBE on the production grid is 4.4e-4 eV from flat 150 × 590, and radial pruning moves it by 1.6e-7 eV. For r2SCAN the radial shells the production grid removes from H (50) and O (66) shift the roots; the radial count matters more than the angular order (99 × 434 is worse than 99 × 302). SCAN is not grid-converged even at flat 150 × 590 (1.2e-3 eV from 200 × 590), so no grid here gets its roots to 1e-4 eV; the 150 × 434 default is 3.6e-3 eV from flat 150 × 590. Compare codes only on the same grid.

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,
    scf_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.
• A constraint is applied at least once before it counts as met. The loop starts at λ = 0, where the state is the unconstrained one. If its population is within lambda_tol of the target but more than 1e-3 electrons away, the loop still takes a Newton step rather than returning the unconstrained state as if it were constrained. Within 1e-3 the unconstrained state is returned with λ = 0: that is the measured spread of one state's population across grids, level shifts and guesses (up to 8.9e-5), with a margin. A loose tolerance can therefore return a state that is only partly constrained, so check r.max_constraint_error (or the per-constraint r.residuals) to see how far from the target it is. The stability-descent restarts follow the same rule.
• 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 no comparison for energies along a target scan or for couplings: 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. Constrained UHF state selection on HeNe⁺ and LiH⁺ is compared through that DFT module, with Hartree–Fock exchange and a vanishing Slater term extrapolated to zero (ferric-scf/tests/validation_cdft_state.rs).
• 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 apart from the HeNe⁺ and LiH⁺ states in validation_cdft_state.rs. On HeNe⁺/def2-SVP at R = 2.0 Å the hole-on-Ne state (N_He = 2) is NWChem's lower state once NWChem's Becke partition is given ferric's He/Ne radius ratio (0.667; NWChem's own is 0.700). Against NWChem's native partition it sits 0.563 eV higher. The offset comes from the partition: with NWChem's radii patched into ferric's radius table (a one-off measurement, not a test), ferric reaches −130.447404 Ha against NWChem's native −130.447401 Ha. 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). On HeNe⁺, an asymmetric pair, the same file checks that both diabats carry the requested charge, that the hole-on-He state is reached through either fragment (3.9e-6 Ha apart), that the hole prefers Ne, the lower-IP atom (gap 1.7646 eV with the default stability descent; 1.0975 eV is a saddle point), and that H_ab is invariant under W → W − cS. These run at R = 2.0 Å only: from the hcore start, at 2.5, 3.5 and 4.0 Å the outer loop does not converge or returns the unconstrained state, and from MINAO at 2.5 Å the two routes to the hole-on-Ne state disagree by 6e-4 Ha. |H_ab| is not asserted on HeNe⁺; cdft_coupling_hene.rs probes it over a distance series. The tests also 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. Apart from the HeNe⁺ and LiH⁺ states above, 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} \]

is the charge moved per unit constraint potential, a susceptibility. cDFT therefore 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

Each outer step is a full converged inner solve, so cDFT costs several plain SCFs.

Electron-transfer coupling

The coupling \( H_{ab} \) between two charge-localized diabatic states follows from a non-orthogonal determinant overlap, computed via Löwdin biorthogonalization. It is the matrix element governing electron-transfer rates in Marcus theory.

Tolerance and coupling error

The λ convergence tolerance (lambda_tol in Python, cdft_lambda_tol in Rust, default 1e-5 electrons) affects the coupling: a loosely converged \( \lambda \) produces diabatic states that miss their targets, 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 with a UHF reference: on He₂⁺/def2-SVP UHF at 2.5–3.5 Å the outer loop does not converge in 40 iterations at 1e-5, while at 1e-2 (what the tests use) it stops after 3 with λ ≈ 0.91 and residuals of 1.6e-3 to 7.8e-3 electrons, so check r.max_constraint_error before trusting the state. With a UKS/PBE reference the He₂⁺ diabats do converge from λ = 0 at 1e-10 (validation_cdft_et.rs).

Cite

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

Periodic systems

ferric runs 3-D periodic Gaussian-basis calculations through the ferric-pbc crate. Conventions follow PySCF pyscf.pbc: Ewald nuclear repulsion with a neutralising background, the dropped-G = 0 Coulomb convention, and the Madelung exchange correction (exxdiv = "ewald"). Energies are Hartree per cell. The equations, with code and test cross-references, are in crates/ferric-pbc/theory.md; the identities are checked by scripts/pbc_ewald_pairft_stress_sympy.py.

What runs

Meshes are Gamma-centred (centring = "gamma") or shifted Monkhorst-Pack ("mp"). The exchange divergence uses the Gamma Madelung constant, or for a mesh that of the diag(N) supercell.

Two-electron builders

• jk = "dense" (default): a dense pure-AFT nao⁴ ERI tensor. It is a correctness oracle for toy cells, hard-capped by max_eri_gb (default 0.5 GiB; range-separated hybrids need 40 nao⁴ bytes).
• jk = "rsgdf": range-separated Gaussian density fitting (PySCF RS-GDF conventions). Requires an auxiliary basis (auxbasis, e.g. "cc-pvdz-ri") and is bounded by the unified memory budget ([memory] budget_gb / memory_budget_gb). Options: gdf_omega (the fit's Ewald split), range_split (energy only), sr_column_rotation.

The nuclear attraction is split at its own omega (default min(2.5 √π/Ω^{1/3}, ω_cap)); any value gives the same energy up to truncation.

Entry points

CLI. Add a [cell] section (lattice, unit, optional kmesh, exxdiv, jk, …) to an input file; method.kind is rhf, uhf, rohf, ksdft, rimp2 (periodic MP2) or pdep-rpa (periodic dRPA). Full key list: Input file, [cell]. Examples: examples/h2-cell-rhf.toml, examples/h2-cell-kpts-rhf.toml, examples/h2-cell-kpts-rhf-stress.toml, examples/h2-cell-rks-opt.toml.

Python. Coordinates and lattice rows in Å, omega/gdf_omega in Å⁻¹:

• Gamma point: run_rhf_gamma, run_uhf_gamma, run_rohf_gamma, run_rks_gamma, run_uks_gamma, run_roks_gamma, run_mp2_gamma, run_drpa_gamma.
• k-point mesh: run_rhf_kpts, run_uhf_kpts, run_rohf_kpts, run_rks_kpts, run_uks_kpts, run_roks_kpts, run_mp2_kpts, run_drpa_kpts.
• auto_gdf_omega chooses an RS-GDF split.

SCF drivers take with_gradient / with_stress; results expose energy, e_nuc (Ewald), madelung, gradient() (Hartree/Bohr per cell) and stress() (σ = (1/Ω) dE/dε, Hartree/Bohr³, not symmetrised; pressure −tr σ/3). Signatures and refusals are in ferric.pyi.

import ferric
mol = ferric.Molecule.from_xyz("testdata/molecules/h2.xyz")
a = 4.0 * 0.529177210903  # 4 Bohr cube, in Angstrom
res = ferric.run_rhf_gamma(
    mol, [[a, 0, 0], [0, a, 0], [0, 0, a]],
    ferric.BasisSet.bundled("sto-3g"), exxdiv="ewald", with_stress=True,
)
print(res.energy, res.stress())

Limits

• 3-D periodicity only; no slabs, wires or 2-D Coulomb truncation.
• Charged cells are an error.
• No lattice-vector relaxation: task = "optimize" moves atoms at a fixed lattice. task = "frequencies" is an error.
• Gamma-point stress is Python-only; the CLI stress = true requires kmesh.
• The dense AFT builder is toy-scale by construction (memory 8 nao⁴ bytes, G count growing as Ω p_max^{3/2}); production-size cells need jk = "rsgdf".
• VV10 is refused on every periodic derivative path; per-functional and per-route support for forces and stress is listed in the input reference and Capabilities and validation.
• Open-shell forces and stress refuse a non-stationary SCF state rather than return a number: ROHF/ROKS orbital gradient above 1e-5 Ha, UKS spin commutator above 1e-5 Ha (k-point), a stored Fock that differs from the rebuilt one by more than 1e-8 Ha, or SCF k-points that are not the mesh's. Converge the SCF tighter (density_conv 1e-10 is FD-grade). Derivation and the full list of tolerances: crates/ferric-pbc/theory.md section 6.
• Gamma-point local MP2 is library-only and its default pair gate keeps every pair below a volume onset.

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)). ferric uses response as the organizing principle instead.

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 demonstrated

Low rank: demonstrated, with a narrower scope. 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 approaches 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 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

Each method family in ferric 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, optional ATM three-body term).

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. Another deterministic order, such as a tree-fold, would change the bits, 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. The grades themselves are METHOD_KIND_GRADES in crates/ferric-core/src/certificate.rs, a copy of the Grade column below that a test keeps equal to it.

rhf, uhf, rohf and ksdft energy runs print a certificate block after the energy and write it to the JSON run log (run_end.extra.certificate; a certificate record for the open-shell runs); Python RhfResult, UhfResult and DftResult carry it as .certificate. It states the kind's grade, the build, the approximations in effect (RI-J / RI-K, LinK / COSX, the XC grid, a level shift, ECPs, embedding, solvation) and warnings: a Plateau exit, a non-Stable stability verdict, and an unchecked UHF/UKS stability. An unconverged SCF is capped to not graded.

Post-SCF kinds print a certificate block after their result: the MP2 family (rimp2 exact, local and open-shell, mp3, oo-rimp2, att-rimp2, mp2-v, scs-mp2, scs-mp2-2terfc, laplace-mp2, laplace-sos-mp2), rs-mp2-rpa, drpa, pdep-rpa, ccsd, ccd, ccsd(t), linlccd, wb97x-l-v, b2plyp, dsd-pbep86, tda, tddft, gw, bse-tda and tdhf-static-polarizability. Kinds that write a result record carry it as result.components.certificate; pdep-rpa, wb97x-l-v, b2plyp, dsd-pbep86, tda, tddft, gw, bse-tda and tdhf-static-polarizability write a certificate record. Python RiMp2Result, Mp3Result, OoRiMp2Result, AttenuatedMp2Result, ScsMp2Result, Mp2VResult, LaplaceMp2Result, SosMp2Result, RsMp2RpaResult, CcResult, PdepRpaResult, GwResult, UGwResult, BseResult, TdhfStaticPolarizabilityResult, DoubleHybridResult and TddftResult carry it as .certificate; the run_drpa and run_linlccd dicts carry a "certificate" key. The grade is the kind's row, capped at the grade of the SCF reference's certificate, whose approximations and warnings it repeats as scf <name> and scf reference: .... It records the auxiliary basis, frozen core, a non-Coulomb correlation operator and each method's settings (SCS scales, Laplace points, frequency quadrature, PDEP truncation, GW variant, double-hybrid scaling or λ and ω, the TDDFT response grid). Caps to not graded: an [local] amplitude threshold at finite eps, an unconverged orbital optimization (oo-rimp2), drCCD fixed point (drpa) or evGW eigenvalue loop. Warnings: GW and BSE-TDA runs at settings other than n_quad = 100, trunc_thresh = 0, unconverged quasiparticle roots, suspect quasiparticle energies, an unconverged PDEP eigensolver. Periodic runs (the CLI [cell] section and the Python Gamma-point and k-point results) carry a certificate with the grade not graded.

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 (both also with RI-J and/or RI-K, without an external potential; all with ECPs), and ROHF with exact J/K or RI and no external potential, with a basis up to f functions; rhf, uhf and rohf frequencies use it there by default.

CLI method.kind matrix

"Reference" is the SCF reference the CLI 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, uhf and rohf (the overlap fit with a functional, and for rohf, is refused before the SCF); a COSX final pass is skipped for them. The RKS, UKS and ROKS optimizers and frequency drivers 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 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 / validation-* CI jobs and on demand (cargo nextest run --run-ignored only -E 'binary(/^validation_/)').