Density fitting & the resolution of identity

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

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), 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 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.

        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 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’s \(B^P_{ia}\) is not a separately-derived object — it is this section’s \(B\), transformed once more.

Exercise 9

  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?

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