# Density fitting & the resolution of identity

The two-electron integral tensor $(\mu\nu|\lambda\sigma)$ scales as $O(n_{\text{ao}}^4)$ — the single
biggest cost in most quantum-chemistry calculations. **Density fitting**, also called the **resolution of
the identity (RI)**, approximates every four-index integral by a contraction over a much smaller
three-index factor, turning an $O(n_{\text{ao}}^4)$ object into an $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$
one with controllable, small error. This chapter derives the approximation itself, the exact metric
whitening qc-rs uses to turn it into the "factor $B$" that every RI backend shares, and how that single
factor serves both mean-field J/K and post-HF correlation — grounded in `.design/54two-electron-integral-spine.md`
and `.design/68ri-jk-strategies.md`, and matching `AGENTS.md`'s unified `eri=` RI-strategy description.

## The resolution-of-identity approximation

The exact AO product density $\rho_{\mu\nu}(\mathbf r)=\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)$ cannot in
general be expanded exactly in any finite auxiliary basis, but it can be **approximated** by a linear
combination of auxiliary fitting functions $\{\varphi_P\}$,

$$
\rho_{\mu\nu}(\mathbf r) \approx \tilde\rho_{\mu\nu}(\mathbf r) = \sum_P c^{\mu\nu}_P\,\varphi_P(\mathbf r).
$$

Dunlap's variational RI chooses the fitting coefficients $c^{\mu\nu}_P$ to minimize the **Coulomb
self-energy of the residual** $\rho_{\mu\nu}-\tilde\rho_{\mu\nu}$ — not the residual's real-space norm —
because it is the Coulomb energy, not the pointwise density error, that determines the accuracy of $J$ and
$K$. Minimizing $\iint\frac{(\rho-\tilde\rho)(\mathbf r)(\rho-\tilde\rho)(\mathbf r')}{|\mathbf r-\mathbf
r'|}\,d\mathbf r\,d\mathbf r'$ over $c^{\mu\nu}_P$ gives the stationary condition
$\sum_Q(P|Q)c^{\mu\nu}_Q=(P|\mu\nu)$, i.e.

$$
c^{\mu\nu}_P = \sum_Q (V^{-1})_{PQ}\,(Q|\mu\nu), \qquad V_{PQ} = (P|Q) = \iint
\frac{\varphi_P(\mathbf r)\varphi_Q(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf r\,d\mathbf r',
$$

where $(Q|\mu\nu)$ is the 3-center Coulomb integral between one auxiliary function and one AO pair, and
$V$ is the 2-center Coulomb metric over the auxiliary basis alone. Substituting the fit into any
4-index integral gives the approximation every RI backend actually uses:

$$
(\mu\nu|\lambda\sigma) \approx \sum_{PQ}(\mu\nu|P)\,(V^{-1})_{PQ}\,(Q|\lambda\sigma).
$$

Because this comes from a genuine variational minimization (not an ad hoc truncation), the RI error is
smooth, systematically improvable by enlarging the auxiliary basis, and empirically tiny for a
properly matched fitting basis (`cc-pvdz-jkfit`, `def2-svp-jkfit`, …, tabulated per orbital basis) — this
is why RI is qc-rs's *recommended default* mean-field path, not merely a cheaper approximation reached for
only under memory pressure.

## Whitening: turning the metric into a single 3-index factor

Working directly with $(V^{-1})_{PQ}$ inside every J/K contraction would mean re-applying an
$n_{\text{fit}}\times n_{\text{fit}}$ operator on every Fock build. qc-rs instead **whitens** the 3-center
integral once, folding the entire metric into a single 3-index factor:

$$
B^P_{\mu\nu} = \sum_Q\bigl(V^{-1/2}\bigr)_{PQ}\,(Q|\mu\nu), \qquad
(\mu\nu|\lambda\sigma) \approx \sum_P B^P_{\mu\nu}\,B^P_{\lambda\sigma},
$$

which is exactly $\sum_{PQ}(\mu\nu|P)(V^{-1})_{PQ}(Q|\lambda\sigma)$ rewritten using
$V^{-1}=V^{-1/2}V^{-1/2}$ — a symmetric square-root factorization rather than an asymmetric one, chosen
so that $B$ carries the whole metric and every downstream contraction becomes a *plain* inner product over
$P$ with no metric left in the loop. qc-rs computes $V^{-1/2}$ two equivalent ways: an eigendecomposition
of $V$ (dropping near-zero eigenvalues below the `ri.metric_cutoff` IOP — the metric can be mildly
rank-deficient for large/diffuse auxiliary bases) or a Cholesky solve $V=LL^{\mathsf T}$, $B=L^{-1}\cdot(\cdot)$
(no eigendecomposition, numerically more stable, done in place) — both give the *same* contracted
$(\mu\nu|\lambda\sigma)$ for the same operator, so which one is used is a numerical-stability choice, not
a physics choice.

Verified example — the RI approximation error on water, and the exact agreement between the two whitening
implementations (`ri-ram` = store-$B$, `ri-recomp` = recompute-$B$ per cycle, both built from the same
math):

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

m4c = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="4c-incore").scf(ref="r").run()
m4c.scf.energy   # -76.02679364497408 -- exact 4-center reference

mri = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-ram").scf(ref="r").run()
mri.scf.energy   # -76.02677274798874
mri.scf.energy - m4c.scf.energy   # 2.09e-05 Ha -- the RI fitting error, with the auto-derived jkfit aux

mri2 = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-recomp").scf(ref="r").run()
mri2.scf.energy - mri.scf.energy   # -1.4e-14 -- ri-ram and ri-recomp agree to the FP-reduction floor
```

A $2\times10^{-5}\ E_h$ fitting error against the exact 4-center reference, for a default aux basis with no
manual tuning, is typical of RI-JK — well below the accuracy of the underlying orbital basis itself. When
no explicit `rijk=` auxiliary basis is given, qc-rs auto-derives the default JK-fitting basis from the AO
basis (`default_jk_aux`) rather than leaving it unset — there is no "auto-RI" that silently *promotes* a
requested exact method to RI, but there is always a sensible aux-basis default once RI is explicitly
requested.

## $J$ and $K$ from the whitened factor

Given $B$, both Coulomb and exchange become simple contractions with the density matrix, entirely local
once $B$ is resident (in RAM, or streamed from disk/RMA — the storage backend is a policy, not a change to
this math):

$$
J_{\mu\nu} = \sum_P B^P_{\mu\nu}\,q_P, \qquad q_P = \sum_{\lambda\sigma}B^P_{\lambda\sigma}D_{\lambda\sigma},
$$

$$
K_{\mu\nu} = \sum_P\sum_i B^P_{\mu i}\,B^P_{\nu i}, \qquad B^P_{\mu i}=\sum_\nu B^P_{\mu\nu}C_{\nu i}
\ \text{(half-transformed to occupied MOs)}.
$$

$J$ needs only the small aux-space vector $q_P$ (one contraction over the density, one over $B$ again);
$K$ needs a half-transform of $B$ to the occupied-MO index first (exactly analogous to the AO→MO
half-transform correlation methods use, [next chapter](post-hf-correlation.md)), then an outer-product
sum. Neither ever touches the full 4-index tensor. This is also why RI is markedly cheaper for $K$ than
conventional 4-center exchange at scale — the half-transform shrinks one AO index to $n_{\text{occ}}$
before the outer product, rather than contracting $n_{\text{ao}}^2$ against $n_{\text{ao}}^2$ directly.

## Two ways to get $B$: recompute vs. store

qc-rs's production RI machinery is **one $\mu$-distributed mechanism with two $L$-source modes**
(`MuRiJk`/`MuRiJkDistributed`), selected by the `eri=` value:

- **`ri-recomp`** stores nothing: it recomputes the 3-center $L=(Q|\mu\nu)$ from libcint *every SCF
  cycle* and contracts with the resident whitening $W=V^{-1/2}$ on the fly — it never materializes the
  full $B$ tensor at all, trading CPU (a re-evaluation each cycle) for memory (only the small metric factor
  is resident). `ri-recomp-disk` is its out-of-core sibling for even tighter memory budgets, spooling the
  per-cycle intermediate to the checkpoint's `scratch` directory.
- **`ri-ram`** builds the whitened factor $\tilde B=W\cdot L$ **once**, in place, and keeps it resident for
  every subsequent SCF cycle — the per-cycle J and K both read the *same* resident $\tilde B$, with no
  per-cycle 3-center recompute at all. This is faster per cycle at the cost of holding the full factor in
  memory.

Both are mathematically identical — bit-identical to each other up to floating-point reduction order, as
the verified example above shows directly (`ri-ram` and `ri-recomp` differ by $\sim10^{-14}$, the ordinary
floating-point noise floor, not a physics difference) — the choice between them is a pure memory/CPU
trade-off, not an accuracy trade-off. `ri-ram` is the incore-speed default at scale; `ri-recomp` is the
memory-frugal choice for the largest bases, and is also what `MoTransformSource` reuses directly for
correlation, since correlation methods need the same 3-center building block regardless of whether SCF
itself stored or recomputed it.

## Distributing $B$ without a distributed metric

The one thing that makes distributed-memory RI hard in general is that a naive scheme needs the metric
$V^{-1/2}$ applied *across* whatever distribution scheme holds the 3-center tensor — an all-to-all or a
remote-memory matrix-vector product on every SCF cycle. qc-rs's design avoids this entirely by exploiting
that **$B$ is SCF-invariant** (built once from integrals only, reused unchanged across every SCF
iteration), so any redistribution cost is paid exactly once, at construction:

1. **The metric is a duplicated array, not a distributed one.** $V$ and its factor
   ($V^{-1}$/$V^{-1/2}$/the Cholesky solve operator) are $O(n_{\text{fit}}^2)$ — tiny next to the
   $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$ 3-center tensor — so every rank simply holds its own full copy
   (the [`Workspace` RMA pool's](../30-hpc/threads-and-blas.md) `dup` array kind: `get`/`put` are local
   memcpy, zero MPI). Duplicating something this small is free; it is what makes every subsequent step
   local.
2. **Whitening is applied once, locally, during a one-time transpose.** Each rank's slab of the raw
   3-center tensor is whitened by a local GEMM against the duplicated metric factor — no distributed
   metric matvec anywhere in the *iterative* SCF loop, only (at most) during this one-time build pass.
3. **The whitened factor is distributed by the auxiliary index $P$, full AO pair per node.** Every rank
   owns a contiguous row-block of $P$ and holds the *complete* $\mu\nu$ range for its own rows. This is
   exactly what makes both $J$ (owned-row local GEMM against $D$, then a tiny $O(n_{\text{ao}}^2)$
   `Allreduce`) and $K$ (owned-row half-transform, then a local outer product) fully local per rank —
   the only inter-rank communication per SCF cycle is one small `Allreduce` of the Fock matrix, never a
   distributed tensor contraction.

Memory aggregates linearly: each rank holds only its own $P$-row share of $B$, so a molecule whose $B$
tensor is too large for one node's RAM can fit comfortably across a handful of nodes — the ≤4-node,
memory-aggregation regime `AGENTS.md`'s qc-hpc/RMA rules describe as the primary distributed-memory target
(not massive HPC scale-out). The multi-rank $K$ build itself has a further choice, the `ri.ram.kdist` IOP
key: `"p-transpose"` (the default — a one-time `MPI_Alltoallv` transposes $\tilde B$ to a
fit-index-distributed layout so every SCF cycle's $J$/$K$ is rank-local plus one small
$\operatorname{Allreduce}(n_{\text{ao}}^2)$) or `"mu-stream"` ($\mu$-distributed, Gram-streaming the
half-transform each cycle instead) — a per-cycle communication/one-time-transpose trade-off, not an
accuracy one.

```{mermaid}
flowchart TD
    L["raw 3-center L = (Q|mu nu)<br/>(libcint, per aux shell)"] --> WHITEN["whiten: B = V^-1/2 . L<br/>(duplicated metric, local GEMM)"]
    WHITEN --> DIST["distribute B by aux index P<br/>(full mu-nu per rank)"]
    DIST --> J["J: q_P = sum B^P.D (owned rows)<br/>J = sum_P B^P q_P<br/>+ Allreduce(nao^2)"]
    DIST --> K["K: half-transform B to occ MOs<br/>K = sum_P sum_i B^P_mu-i B^P_nu-i<br/>+ Allreduce(nao^2)"]
    DIST --> MO["correlation: AO->MO half-transform<br/>B^P_pq (arbitrary MO space pair)"]
```

## The same factor serves correlation

The AO→MO half-transform that [RI-MP2](post-hf-correlation.md) needs is built from the *same* whitened
factor, just contracted against MO coefficients instead of the density matrix:

$$
B^P_{pq} = \sum_{\mu\nu}C^L_{\mu p}\,B^P_{\mu\nu}\,C^R_{\nu q}
$$

for an arbitrary left/right MO-space pair (occupied×virtual for RI-MP2, but the same machinery serves
occupied×occupied or virtual×virtual for any future correlation method). One subtlety qc-rs's design
exploits deliberately: the *order* of whitening and transforming can be swapped — "whiten-then-transform"
(the J/K form above) and "transform-then-whiten" give the *exact same* final $B^P_{pq}$, because whitening
is a linear operation on the auxiliary index alone and commutes with the AO→MO contraction on the orbital
indices. Correlation uses **transform-then-whiten**: shrink $n_{\text{ao}}^2\to n_{\text{left}}\cdot
n_{\text{right}}$ first (a much smaller object once occupied/virtual spaces are far from $n_{\text{ao}}$),
*then* apply the metric — $5$–$20\times$ fewer whitening FLOPs than whitening the full AO-pair tensor
first, and a correspondingly smaller distributed reduction. This is why [the RI-MP2
chapter](post-hf-correlation.md)'s $B^P_{ia}$ is not a separately-derived object — it is this section's $B$,
transformed once more.

:::{exercise}
:label: ex-ri-theory

1. The verified example shows `ri-ram` and `ri-recomp` agreeing to $\sim10^{-14}$ but both differing from
   exact 4-center by $\sim2\times10^{-5}$. Explain in one sentence why these are two completely different
   kinds of "error" — one a genuine physics approximation, the other floating-point noise — and why only
   one of them would shrink with a larger auxiliary basis.
2. Why can the metric $V$ always be treated as a `dup` (fully duplicated, zero-MPI) array regardless of
   how large the molecule or AO basis gets, while the 3-center tensor $B$ cannot?
3. Correlation methods apply the metric *after* the AO→MO transform ("transform-then-whiten"), while
   mean-field J/K applies it *before* ("whiten-then-transform"). Both are mathematically exact
   reorderings of the same operation. What is the practical (not mathematical) reason correlation prefers
   the opposite order from mean-field?
:::

:::{solution} ex-ri-theory
:class: dropdown

1. The `ri-ram` vs `ri-recomp` difference ($\sim10^{-14}$) is pure floating-point rounding from
   accumulating the same exact sum in a different order (store-once-and-reuse vs recompute-every-cycle) —
   it is bounded by machine epsilon and would not shrink with *any* basis change, because both methods are
   computing the identical mathematical quantity. The `ri-ram`/`ri-recomp` vs 4-center difference
   ($\sim2\times10^{-5}$) is the genuine RI approximation error from truncating the resolution of the
   identity to a finite auxiliary basis — it is a real physical/mathematical approximation, and it
   *would* shrink (systematically, and eventually to numerical noise) if the auxiliary basis were enlarged
   toward completeness.
2. The metric $V=(P|Q)$ scales as $O(n_{\text{fit}}^2)$ — quadratic only in the *auxiliary* basis size —
   while the 3-center tensor $B$ scales as $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$, quadratic in the
   (typically much larger and faster-growing) AO basis size as well. For any realistically sized molecule
   the metric stays small enough to duplicate cheaply on every rank, while $B$ is exactly the object whose
   size forces memory-aggregation across nodes in the first place — duplicating it instead of
   distributing it would defeat the entire point of going distributed.
3. Because the AO→MO transform itself is far more expensive to redo than a cheap metric application, and
   it *shrinks* the tensor from $n_{\text{ao}}^2$ down to $n_{\text{left}}\cdot n_{\text{right}}$ (occupied
   × virtual is far smaller than AO × AO once you're past a small basis). Applying the metric on the
   smaller, already-transformed object is a strictly cheaper GEMM than applying it on the larger AO-pair
   object, and it also means any distributed reduction during whitening operates on the smaller tensor.
   Mean-field J/K, by contrast, needs the *AO*-basis result every SCF cycle (Fock matrices are AO
   objects), so there is no shrinking transform to defer the metric past — whitening once at the AO level
   up front is already the cheapest option there.
:::

The next two chapters build directly on this factor: [analytic derivatives](analytic-derivatives.md) and
[linear-response theory](linear-response-cphf.md) both need RI-consistent derivative integrals of exactly
this $B$, and [RI-MP2](post-hf-correlation.md) is, as shown above, nothing more than one further
contraction of it.