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(no·nv³)-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.

QuantitySystem / basisReferenceAgreementPinned by
CCSD correlation energyH2 / STO-3Gexact-integral numpy−0.02052453 Haferric-cc/src/ccsd.rs::test_ccsd_h2_sto3g
Closed-shell (T)H2O / cc-pVDZPySCF ccsd_t()~1e-6 Ha (test asserts 1e-4)ferric-cc/src/ccsd_t_closed_shell.rs::closed_shell_t_h2o_ccpvdz_matches_pyscf
Streaming vs dense spin-orbital (T)H2O / cc-pVDZferric's former dense path5e-16 Haferric-cc/src/ccsd_t.rs::streaming_matches_dense_h2o_ccpvdz

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 Wabcd τijcd 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 Wabcd every iteration: nv⁴ memory, a few times over.
  • ri holds neither the (ab|cd) block nor Wabcd. It rebuilds the (ac|bd) entries from the RI 3-index tensors every iteration in ablk·nv³ panels: at most one nv⁴-sized panel at full width (dense holds four nv⁴ 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 (ablk = nv, 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.

SystemnvPer iteration, densePer iteration, RIRI / densePeak RSS, densePeak RSS, RI
water197.8 ms8.3 ms1.06——
C3H8691.27 s1.64 s1.301.50 GB0.70 GB
C4H10893.2–3.7 s5.2–5.3 s1.503.79 GB1.67 GB
C5H121098.8–10.8 s14.5–14.6 s1.487.90 GB3.40 GB

RI is 1.3–1.5× slower per iteration from nv ≈ 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 no·nv³ (ia|bc) block in core, so it lowers the nv⁴ 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 nv = 40 and 70% at nv = 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.