ferric

A Rust-native quantum chemistry engine, wrapping libint2 for electron integrals, with pyo3 Python bindings.

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

• Attenuated MP2 — MP2 builds dispersion from an uncoupled polarizability that over-polarizes, giving too-large \( C_6 \) and overestimated π-stacking. Attenuating the correlation operator tames that response error with a single tunable parameter.
• PDEP-RPA / GW — the dielectric matrix is the density–density response. PDEP keeps only its dominant low-rank eigenmodes, so RPA correlation and the GW screened interaction need no explicit sum over empty states.
• Constrained DFT — a constraint couples to the density and reads its response (\( \partial N / \partial \lambda \) is a susceptibility), building charge-localized diabatic states and their electron-transfer couplings.

The motivating claim is that response is local in real space and low-rank in its eigenspectrum, so organizing around it should make the computation cheaper.

Implemented ≠ validated

Working code is not a checked number. This documentation describes what exists; it is not a claim that every number is trustworthy.

For how strongly each capability's numbers are checked against ground truth — and where they are known to fail — see What is validated. Capability maturity varies a great deal between methods, and they are graded individually rather than presented as a flat list of equals.

Where to start

Source

github.com/mgoldey/ferric — dual-licensed MIT / Apache-2.0.

Electronic response

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

The object of interest appears under several names depending on which perturbation you apply:

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 motivating 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 matrix 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, then organizing the computation around response should make it cheaper: attenuate the operator, keep the dominant dielectric modes.

What is actually demonstrated

This is where honesty matters more than the pitch.

The low-rank half is demonstrated. PDEP's compression of the dielectric matrix works and is used in production paths — see RPA and GW. Keeping only the dominant eigenmodes removes the explicit sum over empty states that conventional RPA and GW require.

The real-space locality half remains a design premise, not a measured result. Several attempts to exploit it are implemented and measured negative:

• The AO-sparse Laplace SOS-MP2 variant's truncation radius tracks the molecular diameter instead of saturating — so it is not a reduced-scaling path.
• Local MP2 (amplitude-threshold) is implemented with localized virtuals and per-pair domain-local RI fits, but the J build is still dense-from-RI, so no scaling claim is made.
• RI-Laplace MP2 is dense; it serves as a correctness reference for the AO formulation, not as an O(N) path.

Those are reported as negative results rather than quietly omitted, because a locality claim that has not survived measurement is not a feature.

Why this framing is useful anyway

Even where the scaling payoff has not materialized, the response framing buys something concrete: it makes the error in one method diagnosable through another.

MP2's dispersion error is the clearest case. MP2 builds dispersion from an uncoupled polarizability, which over-polarizes — giving \( C_6 \) coefficients that are too large and overestimated π-stacking. That is not a mysterious failure of a wavefunction ansatz; it is a specific, identifiable defect in a response function, and it suggests a specific fix: attenuate the correlation operator so the over-polarizing long-range part is damped.

That is attenuated MP2, and it works for a reason the response picture predicts.

Where the methods come from

The three method families in ferric are not an arbitrary selection. Each one attacks the response function from a different direction.

Attenuated MP2 — fixing a response error

MP2 correlation is built from an uncoupled polarizability. Uncoupled means the density fluctuation does not feel the field it creates: there is no self-consistency in the response. The result over-polarizes, which shows up as:

• \( C_6 \) dispersion coefficients that are too large
• overestimated π-stacking energies
• basis-set superposition error that partly cancels the overestimate, disguising the problem in small basis sets

Attenuating the correlation operator — replacing \( 1/r \) with \( \mathrm{erfc}(\omega r)/r \) or a terfc form — damps the long-range part where the uncoupled approximation is worst, with a single tunable parameter.

This is Goldey & Head-Gordon (JPCL 2012); the dual-attenuated SCS variant is Goldey, Dutoi & Head-Gordon (PCCP 2013). See The MP2 family.

PDEP-RPA and GW — compressing the response

The dielectric matrix is the density–density response function. Conventional RPA and GW evaluate it through an explicit sum over empty orbital states, which is expensive and converges slowly with basis size.

PDEP — projective dielectric eigenpotentials — builds a low-rank basis from the dominant eigenmodes of the dielectric matrix instead. Because the spectrum decays quickly, a modest number of modes captures the physics, and the sum over empty states disappears.

This is the part of the locality-and-low-rank claim that is actually demonstrated in this codebase. See RPA and GW.

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, and it yields charge-localized diabatic states whose electron-transfer couplings \( H_{ab} \) follow from non-orthogonal determinant overlaps.

See Constrained DFT.

What this buys

Three methods, 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 the measured negatives); that it pays off in diagnosis is already clear.

Installation

ferric links against libint2, a C++ integral library that must be built from source. That build takes roughly 30 minutes and is the main cost of getting started.

Prerequisites

• Rust 1.75+ — install via rustup
• libint2 2.7+ — from the mpqc4 tarball
• OpenBLAS and LAPACK
• Eigen3 headers
• Python 3.10+ and maturin — optional, for the Python bindings

What the mpqc4 export does and does not carry

These capabilities are fixed when the tarball is generated, so no cmake flag changes them. compiler.config inside the tarball records the exact settings.

Building against this tarball is correct for everything ferric currently validates. The G12-dependent tests detect its absence at run time and skip with an explicit message rather than failing.

Getting either missing capability requires re-generating libint2 from the upstream source repo with the corresponding --enable-* flags — a substantially longer build.

Building from source

# System dependencies (Ubuntu 22.04+)
sudo apt-get install -y build-essential cmake g++ gfortran wget \
    libeigen3-dev libopenblas-dev liblapack-dev pkg-config \
    python3-dev python3-pip python3-venv

# Build and install libint2 (~30 min)
wget https://github.com/evaleev/libint/releases/download/v2.7.2/libint-2.7.2-mpqc4.tgz
tar xzf libint-2.7.2-mpqc4.tgz
cd libint-2.7.2-mpqc4
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$HOME/.local -DCMAKE_POSITION_INDEPENDENT_CODE=ON
make -j$(nproc)
make install
cd ../..

# Build ferric
cargo build --release

# Run tests
OPENBLAS_NUM_THREADS=1 cargo test --workspace

Debug vs release

The difference is large enough to matter in practice:

• 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 permit.
• Release (cargo build --release) — slow to compile, fast to run. Use it for anything you will actually wait on: real molecules, benchmarks, and any RPA/GW/CC job.

Python bindings

# Set up the venv and install the extension in editable/develop mode
python3 -m venv .venv && source .venv/bin/activate
pip install maturin
maturin develop --release

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

See Python bindings.

Optional: MPI

The mpi feature additionally needs an MPI implementation (OpenMPI or MPICH) and libclang (libclang-dev, for mpi-sys's bindgen step).

Threading

Set OPENBLAS_NUM_THREADS=1. ferric uses rayon for outer parallelism and pins BLAS to a single thread inside rayon workers; letting OpenBLAS thread on top of that oversubscribes the machine and produces unstable timings.

Quick start

Both interfaces cover most methods: a TOML-driven CLI and Python bindings.

Build first — see Installation. ferric needs libint2 built and on the linker path.

CLI

Calculations are described by a TOML file. The examples/ directory has one per method.

# RHF on water with STO-3G
cargo run --release -- examples/water-rhf.toml

# RI-MP2 on water with cc-pVDZ / cc-pVDZ-RI
cargo run --release -- examples/water-rimp2.toml

# Attenuated RI-MP2 (short-range correlation only, r0 = 1.05 Å)
cargo run --release -- examples/water-attmp2.toml

# SCS-MP2 (Grimme spin-component scaling)
cargo run --release -- examples/water-scs-mp2.toml

# SCS-MP2(2terfc) (dual-attenuated, Goldey/Head-Gordon 2013)
cargo run --release -- examples/water-scs-mp2-2terfc.toml

# CCSD (H2/STO-3G)
cargo run --release -- examples/water-ccsd.toml

# LinLCCD(hh) — linearized hole-hole ladder CCD (closed-shell only)
cargo run --release -- examples/water-linlccd.toml

# wB97X-L-V — a double hybrid built on LinLCCD(hh) instead of MP2
cargo run --release -- examples/water-wb97xlv.toml

CLI coverage is not complete. Only method.kind = "ccsd" is wired for coupled cluster; CCD and CCSD(T) are library/Python-only. Use ferric.run_ccd / ferric.run_ccsd_t from Python until a CLI arm is added.

Python

import ferric

mol = ferric.Molecule.from_xyz("testdata/molecules/water.xyz")
bs  = ferric.BasisSet.bundled("cc-pvdz")
aux = ferric.BasisSet.bundled("cc-pvdz-ri")

# Standard RI-MP2
mp2 = ferric.run_rimp2(mol, bs, aux)
print(f"RI-MP2 total: {mp2.total_energy:.10f} Ha")

# Attenuated RI-MP2 (omega in Å⁻¹)
att = ferric.run_attenuated_rimp2(mol, bs, aux, omega=0.420)
print(f"Att-MP2 total: {att.total_energy:.10f} Ha "
      f"(E_OS={att.e_os:.6f}, E_SS={att.e_ss:.6f})")

# SCS-MP2 (Grimme defaults)
scs = ferric.run_scs_mp2(mol, bs, aux)

# SCS-MP2(2terfc) — thesis defaults r0_1=0.75Å, r0_2=1.05Å, c_OS=1.27, c_SS=4.05
terfc = ferric.run_scs_mp2_2terfc(mol, bs, aux)

# Coupled cluster — RI-CCSD(T)
cc = ferric.run_ccsd_t(mol, bs, aux)
print(f"CCSD(T) total: {cc.correlation_energy + cc.t_correction:.10f} Ha")

See Python bindings for the full surface and threading notes.

Threading

Set OPENBLAS_NUM_THREADS=1 when running tests or benchmarks. ferric uses rayon for outer parallelism and pins BLAS to one thread inside rayon workers; letting OpenBLAS thread on top of that oversubscribes the box and can produce unstable timings.

OPENBLAS_NUM_THREADS=1 cargo test --workspace

For throughput across many independent jobs, prefer many single-threaded processes over one multi-threaded job.

Python bindings

The pyo3 bindings expose most of the library. Build them with maturin develop --release — see Installation.

Basics

import ferric

mol = ferric.Molecule.from_xyz("testdata/molecules/water.xyz")
bs  = ferric.BasisSet.bundled("cc-pvdz")
aux = ferric.BasisSet.bundled("cc-pvdz-ri")

Bundled orbital bases include STO-3G, 6-31G, cc-pVDZ and def2-SVP; bundled auxiliary bases include cc-pVDZ-RI, the def2-*-RIFIT family, and def2-universal-jkfit. Both BSE-JSON and Gaussian-94 basis files can be parsed from disk.

Ground state

rhf = ferric.run_rhf(mol, bs)
print(f"RHF: {rhf.energy:.10f} Ha, converged={rhf.converged}")

uhf  = ferric.run_uhf(mol, bs, multiplicity=3)
rohf = ferric.run_rohf(mol, bs, multiplicity=3)
dft  = ferric.run_dft(mol, bs, functional="b3lyp")

Always check .converged. These functions return a result whether or not the SCF converged — a non-converged result is not an error, it is a result with converged = False. Treating it as success is a common way to get a plausible, wrong number.

Correlation

mp2   = ferric.run_rimp2(mol, bs, aux)
att   = ferric.run_attenuated_rimp2(mol, bs, aux, omega=0.420)   # Å⁻¹
scs   = ferric.run_scs_mp2(mol, bs, aux)
terfc = ferric.run_scs_mp2_2terfc(mol, bs, aux)

ccd    = ferric.run_ccd(mol, bs, aux)
ccsd   = ferric.run_ccsd(mol, bs, aux)
ccsd_t = ferric.run_ccsd_t(mol, bs, aux)

run_ccd and run_ccsd_t are Python-only — they are not yet CLI-wired.

Response and excited states

gw    = ferric.run_gw(mol, bs, aux)      # G0W0 quasiparticle energies
tddft = ferric.run_tddft(mol, bs, aux)   # TDA or Casida

run_tddft warns on stderr when the reference is not pure Hartree–Fock: the (ia|f_xc|jb) XC-kernel response is unimplemented, so with a DFT reference the excitation energies omit a physical term and are approximate. Only c_hf = 1.0 (CIS/TDHF) is exact.

Memory budgets

Most drivers accept memory_budget_gb. The budget is enforced, not advisory: a job whose predicted peak exceeds it is refused with a breakdown naming the dominant term, rather than being OOM-killed partway through.

mp2 = ferric.run_rimp2(mol, bs, aux, memory_budget_gb=8.0)

Threading and concurrency

The compute drivers release the GIL, so independent jobs submitted from a ThreadPoolExecutor genuinely run in parallel rather than serializing at the FFI boundary.

Set OPENBLAS_NUM_THREADS=1. ferric uses rayon for outer parallelism and pins BLAS to one thread inside rayon workers; for throughput across many jobs, prefer many single-threaded processes over one wide job.

Property export

ESP at nuclei, electric fields, static and atom-partitioned polarizabilities, Hirshfeld and Löwdin charges, and density matrices are available, with NPZ export of ML-ready features (MO coefficients, orbital energies, PDEP eigenvectors, ESP, polarizability tensors, charges) for downstream model conditioning.

For the full signature list, see the API documentation.

Methods overview

What exists, grouped by family. Maturity varies a great deal between these — see What is validated before trusting any particular number.

Infrastructure

Shared machinery underneath all of the above:

• QQR screening (Maurer, Lambrecht & Ochsenfeld 2012) and LinK exchange (Ochsenfeld, White & Head-Gordon 1998) for the Fock build
• Spherical and Cartesian basis support, with BSE-JSON and Gaussian-94 parsers
• einsum! — a tensor-contraction macro routing contractions through BLAS3 GEMMs, used throughout the CC and MP3 code
• Memory budgets — enforced allocation ceilings that refuse an oversized job with a named breakdown rather than letting it be OOM-killed
• Python bindings (pyo3) and a TOML-driven CLI

Properties

ESP at nuclei, electric field, static and atom-partitioned polarizabilities, Hirshfeld and Löwdin charges, density matrices, and NPZ export of ML-ready features for downstream generative-model conditioning.

A note on negative results

Several reduced-scaling approaches in this codebase are implemented and measured negative — they are documented as such rather than omitted:

• AO-sparse Laplace SOS-MP2 — truncation radius tracks the molecular diameter instead of saturating
• Local MP2 — the J build is still dense-from-RI, so no scaling claim is made
• TDHF/RPAx \( C_6 \) — stays ~60% low regardless of gap

Knowing which ideas did not work is part of the documentation.

SCF and DFT

Ground-state self-consistent field methods, plus analytical nuclear gradients.

Hartree–Fock

• RHF (closed-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
• UHF / ROHF (open-shell) with per-spin DIIS, virtual-space level shifting, augmented-Hessian Newton, and Maximum-Overlap-Method (MOM) orbital tracking for near-degenerate cases

The convergence machinery is not decoration. Heavy atoms and near-degenerate frontier orbitals genuinely break plain DIIS; the virtual-block level shift and MOM exist because specific systems failed without them.

Kohn–Sham DFT

Closed- and open-shell (RKS / UKS / ROKS) via libxc:

• LDA, GGA, hybrid, and range-separated-hybrid functionals — LDA, PBE, B3LYP, ωB97X-V
• Becke–Lebedev grids, with pruning
• VV10 nonlocal correlation

Gradients

Analytical nuclear gradients for RHF, UHF, ROHF and KS-DFT — including grid response — validated against finite differences.

Hessians are not implemented. The mpqc4 libint2 export does not carry second-derivative integrals, so the CPKS machinery is stubbed with explicit TODOs rather than silently absent. See Installation.

Convergence

A few things worth knowing before debugging a stubborn SCF:

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

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

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

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

Screening

• Schwarz bounds on every 4-centre path, built so they can never underestimate (a zero-valued table entry once cost 1.5e-4 Ha; it is now floored)
• LinK (Ochsenfeld, White & Head-Gordon 1998) — exchange via significant-pair and density-pair lists; the lists were corrected in #50 and 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 before the #50 list fix)
• COSX shell-pair screening uses a primitive-level Hölder bound that provably never underestimates; an earlier overlap-based bound did, and silently corrupted K

Choosing how exchange is built

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

* full SCF at TZ was 358 s for COSX against 98 s direct — COSX is slower at TZ.

Start with density fitting. RI-JK is two to three orders of magnitude faster than anything else here whenever its three-index tensor fits in memory (n_aux × n_bf² × 8 bytes — 1 GB for butane/QZVP, 7 GB for octane/QZVP), and it spills to disk when it does not. Its error with a JK-fitting auxiliary basis is a few µHa. If your system fits, this is the answer and the rest of this section is about when it does not.

Need exact exchange? Direct or LinK; both are exact to the screening threshold as K builders. LinK's pair lists were fixed in #50 (three pair-list defects; butane/def2-SVP link == direct to 9e-12 Ha). Every LinK timing taken before that fix was against a kernel that skipped quartets, so none is repeated here; its cost against the corrected kernel is being re-measured. k_builder is honoured by UHF and ROHF as well as RHF (it was silently ignored for open-shell runs before 2026-09-08): 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. It is skipped with a warning, never silently, whenever density-fitted J/K is active, the functional uses no exact exchange, or the functional is range-separated (exchange then comes from the SR/LR fitters).

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 reaches analytic exchange only at quadruple-zeta: on butane/def2-QZVP a COSX K build is 137 s against 400 s for the default direct J+K build (parity; J and K share that sweep), while at TZ the full COSX SCF is 3.7× slower than direct (358 s vs 98 s). Below QZ it is the wrong tool. Ratios against LinK are withheld until LinK is re-measured with the #50 lists.

Its integral work is sub-quadratic in system size — a density-driven pair screen (on the product of the integral bound and the local half-transformed density) gives an A-build tail exponent of N^1.5 on C12–C20 alkanes at def2-SVP, with a K error below 2e-6 Ha at the default threshold. The half-transforms D·X and X·Gᵀ are still dense GEMMs, which grow faster and are a third of the build by C20; until they are made sparse (the standard next step), expect the full build to scale roughly N^2 past a dozen heavy atoms even though the integrals do not.

COSX's error is a grid error, and it is not µHa-small: 5e-6 Ha on water/cc-pVDZ and 1.2e-4 Ha on butane/def2-TZVP at the default grid. Reaction energies cancel most of it (0.02 kcal/mol on an isodesmic alkane reaction at the same grid); absolute energies do not. Three knobs, all optional:

• cosx_grid = { radial = 50, angular = 110 } is the default and the coarsest grid that meets a 0.1 kcal/mol reaction-energy bar. Coarser grids fail it. Finer grids reduce the error roughly tenfold per step and cost proportionally.
• cosx_overlap_fit = true (default) applies the Izsák–Neese overlap correction. At the default grid it helps; on coarser grids it makes things worse, and its benefit is strongly molecule-dependent — large on water, nil to negative on ethane — so do not expect the factor quoted in the literature.
• cosx_backend = "md3c1e" (default) is the batched McMurchie–Davidson integral kernel. "cosx-a" is the per-point libint2 path, about three times slower and kept only as the cross-check the kernel is anchored against.
• cosx_screen_thresh = 1e-7 (default) is the density-driven pair-screening threshold. 0.0 disables screening bit-identically; 1e-6 already fails a 1e-6 Ha K-error bar on butane. Not available with the "cosx-a" backend.

Setting any cosx_* key without k_builder = "cosx", or k_builder together with df_k_aux, is refused or warned about rather than silently ignored.

Determinism

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

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

The MP2 family

The largest method family here, and the one most directly tied to the response framing.

Why attenuation

MP2 builds dispersion from an uncoupled polarizability — the density fluctuation does not feel the field it creates. That over-polarizes, producing:

• \( C_6 \) coefficients that are too large
• overestimated π-stacking
• an error partly masked by BSSE in small basis sets, which is why the problem is easy to miss

Attenuating the correlation operator damps the long-range part where the uncoupled approximation is worst. One parameter, and the error it targets is identifiable rather than empirical.

Variants

RI-MP2 — density-fitted via 3-center/2-center integrals. Canonical MP2 is also implemented, for cross-validation rather than production.

Attenuated RI-MP2 — \( \mathrm{erfc}(\omega r)/r \) and terfc operators (Goldey & Head-Gordon, JPCL 2012).

SCS-MP2 — Grimme spin-component scaling (JCP 2003), and SCS-MP2(2terfc), dual-attenuated (Goldey, Dutoi & Head-Gordon, PCCP 2013).

OO-RI-MP2 — orbital-optimized, with level-shifted Newton, orbital DIIS, Cayley rotations and backtracking.

MP3 — spin-orbital third-order Møller–Plesset via the einsum! framework.

MP2-V — attenuated MP2 combined with VV10 nonlocal correlation.

Laplace formulations

RI-Laplace MP2 — AO-Laplace via pseudo-density matrices. The implementation is dense; it is the correctness reference for that formulation, not a reduced-scaling path. No O(N) has been measured.

Laplace SOS-MP2 — opposite-spin-only with a Laplace-factorized denominator (c_os scaling, minimax quadrature). The MO and AO formulations agree to machine precision.

The AO-sparse variant is measured negative: its truncation radius tracks the molecular diameter instead of saturating, so it does not deliver reduced scaling. That is reported rather than omitted.

Local MP2

Amplitude-threshold LMP2 — WSHG23 single-threshold, with localized virtuals and per-pair domain-local RI fits.

Counters only — no scaling claim is made. The J build is still dense-from-RI, so while the amplitude machinery is in place and anchored, the end-to-end cost is not reduced. The assembly step, not the solve, is the measured wall.

Size-extensivity

RI-MP2's total energy is size-extensive to 2e-12 Ha for a well-separated dimer versus twice the monomer — pinned by a test, not asserted. That is five orders inside the test's tolerance, so it passes on physics rather than on a loose bound.

Robust fitting

When the RI metric differs from the physical kernel — as it does for attenuated operators — robust (Dunlap) density fitting is required, not optional. Domain-local \( V^{-1} \) plus robust fitting gives size-extensive µHa-level MP2 error; the non-robust form collapses in a way driven by the metric, not by the domain size.

Coupled cluster

RI-based coupled cluster, validated against exact-integral and PySCF references.

What exists

• RI-CCD — doubles only
• RI-CCSD — singles and doubles
• (T) — the perturbative triples correction

H2O/cc-pVDZ CCSD(T) matches PySCF to roughly 1e-6 Ha.

CLI coverage is partial

Only method.kind = "ccsd" is wired into the CLI — see examples/water-ccsd.toml. CCD and CCSD(T) are library/Python-only:

ccd    = ferric.run_ccd(mol, bs, aux)
ccsd_t = ferric.run_ccsd_t(mol, bs, aux)

LinLCCD

Linearized hole-hole ladder CCD, closed-shell only. Its main use here is as the correlation component of wB97X-L-V — a double hybrid that converges its own KS reference and then adds a short-range LinLCCD(hh) correction, rather than the MP2 correction a conventional double hybrid uses.

The [dft] lambda and omega keys override the published 0.6 / 0.1 Bohr⁻¹ values; omitting them gives the published parameters.

Implementation

All contractions route through einsum!, a macro that maps tensor contractions onto BLAS3 GEMMs. That matters for performance: the alternative — many small explicit loops — leaves most of the machine's throughput unused.

The permutation copies that feed those GEMMs are parallelized. This is less trivial than it sounds: for a strided permutation the copy can dominate the contraction it feeds — measured at 47% at nv=40 and 70% at nv=80, since it is memory-bandwidth-bound and gets relatively worse with size.

Those copies are bit-identical regardless of thread count. A permutation is pure data movement — every output element written exactly once — so unlike a reduction there is no summation order to perturb. That property is pinned by a test that has been verified to fail when deliberately broken.

Memory

The amplitude tensors dominate, and they 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 rather than being OOM-killed midway.

The DLPNO family additionally reads its own budget, which it previously ignored.

RPA and GW

The part of the response claim that is actually demonstrated.

PDEP

The dielectric matrix is the density–density response function. Conventional RPA and GW evaluate it through an explicit sum over empty orbital states — expensive, and slow to converge with basis size.

PDEP — projective dielectric eigenpotentials — instead builds a low-rank basis from the dominant eigenmodes of the dielectric matrix. The spectrum decays quickly, so a modest number of modes captures the physics and the empty-state sum disappears.

This low-rank compression works and is used in production paths. It is the demonstrated half of the design premise; the real-space locality half is not (see Electronic response).

RPA

• PDEP-RPA — RPA correlation via a low-rank W basis in Gaussians
• U-PDEP-RPA — open-shell, over a spin-summed dielectric
• Attenuated RPA — short-range correlation via erfc

The solver defaults to Lanczos; a dense path is used for small problems. Note that the eigensolve is serial by design — that is a deliberate choice, not an oversight, and RPA here is already faster than the PySCF reference.

GW

• G0W0, COHSEX, evGW0, evGW
• U-GW — unrestricted

G0W0@HF matches MOLGW to roughly 5 meV.

The quasiparticle solve runs a Newton root-find on the self-energy, which is fragile near \( \Sigma_c \) poles. That fragility has a consequence worth recording: the frequency-quadrature loop inside the self-energy is a sequential floating-point accumulation, so it cannot be parallelized without changing summation order — and reordering would perturb quasiparticle energies in a thread-count-dependent way. The loop is deliberately left serial.

TDDFT

Linear response in both the Tamm–Dancoff approximation (TDA/CIS) and the full Casida equations, closed-shell references.

Important limitation. The \( (ia|f_{xc}|jb) \) XC-kernel response is not implemented. With a pure Hartree–Fock reference (\( c_{HF} = 1 \)) that term is identically zero and the result is exactly CIS/TDHF. With any DFT reference it is not zero, and the excitation energies omit it — they are approximate.

The code warns on stderr when \( c_{HF} \neq 1 \) rather than returning silently incomplete numbers.

Double hybrids

B2PLYP and DSD-PBEP86, plus wB97X-L-V — see Coupled cluster.

Dispersion and polarizability

Static and atom-partitioned polarizabilities, Casimir–Polder \( C_6 \) coefficients, and many-body dispersion.

TDHF/RPAx is a measured negative for dispersion: the static α it produces is reasonable, but the \( C_6 \) stays roughly 60% low regardless of gap. It is a polarizability tool, not a dispersion one.

Dynamic dRPA@PBE α, by contrast, gives \( C_6 \) roughly 3× better than the static Tkatchenko–Scheffler (TS) model.

Constrained DFT

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

The response connection

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

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

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

Implementation

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

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

Electron-transfer coupling

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

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

A caveat

The cdft_lambda_tol convergence tolerance interacts with the coupling calculation in a way worth checking: a loosely converged \( \lambda \) produces diabatic states that are not quite the ones you asked for, and \( H_{ab} \) inherits that error. Tighten it before trusting a coupling.

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

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

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

Cross-cutting conventions

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

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

Memory. MemoryPlan expresses what a path will allocate and when, so an oversized job is refused before allocating, with a breakdown naming the dominant term. Guards are tested in both directions: a starved budget must be refused, 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.

What is validated

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

This is the most important page in this documentation, and the one to read before trusting any result.

The authority is the wiki

The project wiki's VALIDATION.md is the authority on what each capability's numbers are checked against, and where they are known to fail. It grades each capability — proven / smoke / stub — rather than presenting them as a flat list of equals.

This documentation describes what exists. It is not a claim that everything here is equally trustworthy, and the feature list on these pages should not be read as a validation claim.

Why the distinction is drawn so sharply

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

• a non-converged SCF returned as an ordinary result, because convergence is a flag rather than an error
• a method missing a physical term it does not mention — TDDFT's XC-kernel response is exactly this case, which is why it now warns
• a fallback model silently substituted for one atom in a molecule, changing a partitioning without changing the shape of the output
• a screening or truncation threshold that happens to be safe for the test system and not for yours

None of these look like failures. All of them have occurred and been fixed in this codebase. The remedy is not optimism — it is grading each capability separately and saying which ones are checked against ground truth.

Some anchors

Where numbers are checked, they are checked against external references rather than self-consistency:

These are the figures stated in the repository itself. The wiki's VALIDATION.md carries the full per-capability grading, including benchmark sweeps (GW100 and others) whose numbers are not reproduced here — quoting a benchmark MAE from memory rather than from the record is exactly the kind of unchecked claim this page exists to discourage.

Known negatives

Reported rather than omitted:

• AO-sparse Laplace SOS-MP2 — truncation radius tracks molecular diameter instead of saturating; not a reduced-scaling path
• Local MP2 — J build still dense-from-RI; no scaling claim made
• TDHF/RPAx \( C_6 \) — ~60% low regardless of gap
• TDDFT with a DFT reference — omits the \( (ia|f_{xc}|jb) \) kernel response; warns at run time
• Hessians / frequencies — not implemented; the mpqc4 libint2 export lacks second-derivative integrals

Testing discipline

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

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

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

API documentation

The crate-level API docs are generated by rustdoc and published alongside this book.

Generating locally

cargo doc --workspace --no-deps --open

Drop --no-deps to include dependency documentation as well (much slower, and much larger).

If that fails with "Only one may be documented at once since they output to the same path", add --exclude ferric-python. The pyo3 crate's lib is deliberately named ferric — that is what makes Python's import ferric work — which collides with the ferric facade crate. Excluding it costs nothing: it is a cdylib, and its surface is documented in Python bindings.

Entry points

The most useful starting points, by crate:

On reading the docs

The doc comments in this codebase carry more than signatures. Where a design decision was hard-won — a convergence hazard, a memory-accounting subtlety, a parallelization that would be unsafe — the reasoning is recorded at the point of use, including cases where an apparently obvious optimization was rejected and why.

Those notes are often the most useful part of the documentation for anyone modifying the code, and they are deliberately kept next to the code rather than in prose docs that drift.

References

Dependencies

• libint2 — Obara–Saika integral engine
• pyo3 — Rust/Python interop
• ndarray — N-dimensional arrays for Rust
• ndarray-linalg — LAPACK bindings for ndarray
• libxc — exchange–correlation functionals

Methods

General

• Szabo & Ostlund, Modern Quantum Chemistry (1996)
• Pulay, Chem. Phys. Lett. 73, 393 (1980) — DIIS convergence acceleration

MP2 family

• Weigend, Phys. Chem. Chem. Phys. 4, 4285 (2002) — RI-MP2 auxiliary basis sets
• Grimme, J. Chem. Phys. 118, 9095 (2003) — SCS-MP2
• Bozkaya & Sherrill, J. Chem. Phys. 135, 104103 (2011) — orbital-optimized MP2
• Goldey & Head-Gordon, J. Phys. Chem. Lett. 3, 3592 (2012) — attenuated MP2
• Goldey, Dutoi & Head-Gordon, Phys. Chem. Chem. Phys. 15, 15869 (2013) — SCS-MP2(2terfc)

Coupled cluster

• Scuseria, Janssen & Schaefer, J. Chem. Phys. 89, 7382 (1988) — CCSD
• Raghavachari et al., Chem. Phys. Lett. 157, 479 (1989) — CCSD(T) triples
• Bartlett & Musiał, Rev. Mod. Phys. 79, 291 (2007) — coupled-cluster theory

Screening and scaling

• Ochsenfeld, White & Head-Gordon, J. Chem. Phys. 109, 1663 (1998) — LinK exchange
• Maurer, Lambrecht & Ochsenfeld, J. Chem. Phys. 136, 144107 (2012) — QQR screening

Constrained DFT and nonlocal correlation

• Wu & Van Voorhis, J. Chem. Phys. 125, 164105 (2006) — cDFT electron-transfer coupling \( H_{ab} \)
• Vydrov & Van Voorhis, J. Chem. Phys. 133, 244103 (2010) — VV10 nonlocal correlation

License

Dual-licensed under either

• Apache License, Version 2.0 (LICENSE-APACHE)
• MIT License (LICENSE-MIT)

at your option.