# DFT exchange-correlation quadrature

Kohn-Sham DFT's exchange-correlation term $E_{\text{xc}}[\rho]$ has no closed form for any functional
beyond the trivial ones, so every practical KS-DFT code evaluates it — and its Fock contribution
$V_{\text{xc}}$ — by **numerical quadrature over real space**: a molecular grid of points and weights fine
enough that the integral is converged to microhartree-level accuracy, but coarse enough to be affordable
for a system with hundreds of AOs. This chapter derives that grid — the atom-centered radial/angular
product rule, the Becke space-partition that stitches per-atom grids into one molecular grid, the
NWChem-style angular pruning that keeps point counts down without losing accuracy, and how AO values,
densities, and the XC potential matrix are actually assembled block by block — grounded in
`.design/27qc.dftgrid.md` and qc-rs's current `qc-grid` implementation (per `AGENTS.md`'s qc-grid entry:
Treutler-Ahlrichs M4 radial in production, NWChem 5-region pruning, PySCF-Bragg Becke partitioning).

## Molecular quadrature as an atom-centered product rule

The exact XC energy is a 3-D integral over all space, $E_{\text{xc}}=\int\rho(\mathbf
r)\,\epsilon_{\text{xc}}(\mathbf r)\,d\mathbf r$. No single global quadrature handles this well because
the integrand is sharply peaked at every nucleus (from the AO cusp) and smooth everywhere else. The fix,
due to Becke (1988), is to attach one atom-centered spherical grid to every nucleus and glue them together
with a smooth space-partition:

$$
\int f(\mathbf r)\,d\mathbf r \;\approx\; \sum_A\sum_i\sum_k w_{Aik}\,f(\mathbf r_{Aik}), \qquad
\mathbf r_{Aik} = \mathbf R_A + r_{Ai}\,\mathbf n_k,
$$

where $A$ runs over atoms, $i$ over radial shells centered on $A$, and $k$ over angular directions on the
unit sphere at that shell. If the angular (Lebedev) weights are normalized to $\sum_k w_k^{\text{ang}}=1$,
the one-center spherical-shell weight is $\tilde w_{Aik}=4\pi r_{Ai}^2\,w_{Ai}^{\text{rad}}\,
w_k^{\text{ang}}$, and the final molecular weight folds in the atomic partition function $W_A$:

$$
w_{Aik} = \tilde w_{Aik}\,W_A(\mathbf r_{Aik}), \qquad \sum_A W_A(\mathbf r) = 1\ \text{everywhere}.
$$

The partition-of-unity condition is what makes gluing the atomic grids together *exact* in the limit of
converged 1-D radial/angular rules — no double-counting, no gaps — and is the first invariant qc-rs's test
suite checks numerically (target $10^{-12}$).

## The Becke space-partition

Becke's partition function is built from a **two-center hyperbolic coordinate**. For a point $\mathbf r$
and two nuclei $A,B$,

$$
\mu_{AB}(\mathbf r) = \frac{|\mathbf r-\mathbf R_A| - |\mathbf r-\mathbf R_B|}{|\mathbf R_A-\mathbf R_B|}
\in[-1,1],
$$

which is $-1$ at $A$, $+1$ at $B$, and $0$ on the perpendicular bisector plane. A raw switching function of
$\mu_{AB}$ alone would put the dividing surface exactly at the midpoint regardless of how different the
two atoms are (a heavy atom next to hydrogen would get an unreasonably small share of space near itself),
so Becke corrects it using the **Bragg-Slater atomic radii** $R_A^{\text{BS}}$:

$$
\chi_{AB} = \frac{R_A^{\text{BS}}}{R_B^{\text{BS}}}, \qquad
m_{AB} = \frac{\chi_{AB}-1}{\chi_{AB}+1}, \qquad
a_{AB} = \operatorname{clip}\!\left(\frac{m_{AB}}{m_{AB}^2-1},\,-\tfrac12,\,\tfrac12\right),
$$

$$
\nu_{AB} = \mu_{AB} + a_{AB}\left(1-\mu_{AB}^2\right).
$$

(qc-rs uses the **PySCF-Bragg** size-adjusted radii table, matching PySCF's default rather than the
original 1988 Becke radii — the derivation is identical, only the tabulated $R_A^{\text{BS}}$ values
differ.) The corrected coordinate $\nu_{AB}$ is then pushed through an iterated smoothing polynomial to
make the space-partition boundary have zero derivative to high order (avoiding a cusp in $W_A$ that would
hurt quadrature accuracy):

$$
p_0(x)=x, \qquad p_{n+1}(x) = \tfrac32p_n(x) - \tfrac12p_n(x)^3, \qquad
s_h(x) = \tfrac12\bigl(1-p_h(x)\bigr),
$$

with $h$ the **Becke hardness** (default 3 — iterating the smoothing polynomial 3 times). The unnormalized
cell function for atom $A$ multiplies the switching function over every other atom, and normalizing across
all atoms gives the partition function used above:

$$
P_A(\mathbf r) = \prod_{B\ne A} s_h\bigl(\nu_{AB}(\mathbf r)\bigr), \qquad
W_A(\mathbf r) = \frac{P_A(\mathbf r)}{\sum_C P_C(\mathbf r)}.
$$

By symmetry, at the midpoint between two identical atoms $W_A=W_B=\tfrac12$ exactly — one of the
regression tests qc-rs's grid module checks directly. Ghost atoms (basis functions with no nucleus) and
dummy atoms are never partition centers by default, since Bragg-Slater radii are a nuclear-center concept;
they still get their basis functions evaluated on the grid, just not their own dedicated radial/angular
shells.

## Radial quadrature: the Treutler-Ahlrichs M4 mapping

A radial grid discretizes $\int_0^\infty F(r)\,dr\approx\sum_iw_i^{\text{rad}}F(r_i)$. qc-rs's production
radial scheme is the **Treutler-Ahlrichs M4 mapping** (Treutler & Ahlrichs, *J. Chem. Phys.* **102**, 346
(1995)), a faithful port of PySCF's `dft.radi.treutler_ahlrichs` with $\xi=1$ (element-independent — the
atomic-size adjustment happens in the Becke partition via $R_A^{\text{BS}}$, not here). Starting from
Gauss-Chebyshev abscissae on $(-1,1)$, $x_i=\cos\theta_i$ with $\theta_i=\tfrac{i\pi}{n+1}$, and setting
$L=1/\ln2$, the radial node and weight are

$$
r_i = -L\,(1+x_i)^{0.6}\ln\!\left(\frac{1-x_i}{2}\right),
$$

$$
w_i = \frac{\pi}{n+1}\sin\theta_i \cdot L\,(1+x_i)^{0.6}\left[-\frac{0.6}{1+x_i}\ln\!\left(\frac{1-x_i}{2}\right)
+ \frac{1}{1-x_i}\right],
$$

the second term being the chain-rule Jacobian $dr/dx$ needed to convert the plain $x$-space Chebyshev rule
into an $r$-space one. This mapping is markedly more point-efficient at DFT accuracy than the older
Gauss-Chebyshev-with-Becke-rational-map scheme (`RadialGridKind::BeckeGaussChebyshev`, retained as a
legacy option) — fewer radial shells reach the same integrated-electron-count tolerance, which is why it
is qc-rs's production default. (The `qc-grid` crate verifies both schemes against the analytic Gaussian
radial moment $\int_0^\infty e^{-\alpha r^2}4\pi r^2\,dr=(\pi/\alpha)^{3/2}$.)

## Angular quadrature and NWChem pruning

Angular integration over the unit sphere at each radial shell uses **Lebedev-Laikov** grids, selected by
algebraic order (exact for spherical harmonics up to that order) rather than raw point count:

| Lebedev order | Points | Use |
|---:|---:|---|
| 9 | 38 | innermost pruned region |
| 11 | 50 | near-nucleus region |
| 15 | 86 | inner-shell region |
| 29 | 302 | `medium` preset peak |
| 35 | 434 | `fine` preset peak |
| 131 | 5810 | `ultrafine` (dense, unpruned) peak |

Using the peak angular order at *every* radial shell would waste enormous numbers of points near the
nucleus, where the density is nearly spherically symmetric and a low-order rule already integrates it
exactly, and in the far tail, where the density is exponentially small. qc-rs's production pruning is the
**NWChem 5-region scheme** (PySCF/gpu4pyscf-compatible; the earlier continuous ORZ-style ramp is
deprecated legacy, retained only for API compatibility). For an atom of Bragg radius $R_A^{\text{BS}}$, the
radial coordinate is expressed as a fraction $r/R_A^{\text{BS}}$ and split into five regions by four
element-band-dependent cut points $\alpha_1<\alpha_2<\alpha_3<\alpha_4$ (three bands: H/He, Li–Ne, else —
`.design/27qc.dftgrid.md`'s `NWCHEM_ALPHAS` table), and the five regions take angular orders

$$
\bigl[\,11,\ 15,\ \text{ladder}[j-1],\ n_{\max},\ \text{ladder}[j-1]\,\bigr],
$$

where $j$ indexes how far $n_{\max}$ sits along the fixed order ladder — dense only in the physically
important mid-radius bonding region, and light near the nucleus and in the tail. This is a discrete,
element-band/region table rather than a continuous formula, in contrast to the older deprecated ramp
$n_{\text{ang}}(r_i)=\min\{n\in\mathcal L\mid n\ge\min(\lfloor k_sr_in_{\max}/R_A^{\text{BS}}\rfloor,\,
n_{\max})\}$ — both express the same idea (dense angular sampling only where it is needed), but the
NWChem table is calibrated to match PySCF/ORCA reference grids exactly rather than following a smooth
ramp.

### The preset ladder: coarse → ultrafine

qc-rs exposes four named presets, an ORCA-`DefGrid`-like ladder trading cost for accuracy:

| Preset | Radial shells | Peak angular order (points) | Pruned? | Worst-case error vs `ultrafine` |
|---|---:|---:|:---:|---:|
| `coarse` | 50 | 23 (194) | yes | $\lesssim 8\ \mu E_h$/atom |
| `medium` (default) | 60 | 29 (302) | yes | $\lesssim 0.4\ \mu E_h$/atom |
| `fine` | 75 | 35 (434) | yes | $\lesssim 0.2\ \mu E_h$/atom |
| `ultrafine` | 200 | 131 (5810) | **no** (dense reference) | — |

Only `ultrafine` is unpruned — it is the dense reference tier the other three are calibrated against, not
a tier meant for routine production work. Additionally, 3d transition metals (Sc–Zn) get denser
radial shells and a higher angular peak on the pruned tiers (`qc-grid::tuning`, per-element grid tuning,
#390) — the presets above are the *general* recipe; the tuning table refines it for elements whose compact
core density is harder to integrate accurately with the generic point counts.

Verified example — water/def2-SVP, B3LYP, across the ladder:

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
for g in ("coarse", "medium", "fine"):
    m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r", xc="b3lyp", grid=g).run()
    print(g, m.scf.energy, m.scf.converged)
# coarse -76.35812209974867 True
# medium -76.35812511536525 True
# fine   -76.35812496750043 True
```

The `coarse`→`medium` energy shift here is $\approx3\times10^{-6}\ E_h$ across 3 atoms — consistent with
the calibrated $\lesssim8\to\lesssim0.4\ \mu E_h$/atom tolerances above — and `medium`→`fine` changes the
9th digit, well past chemical accuracy; `medium` is qc-rs's default for exactly this reason.

## AO values and density on the grid

Each AO must be evaluated (and, for GGA/meta-GGA, its gradient) at every retained grid point. A contracted
Cartesian Gaussian centered on atom $A$ is

$$
\chi_{\kappa,lmn}(\mathbf r) = \sum_p c_{\kappa p}\, x_A^l y_A^m z_A^n\, e^{-\alpha_p r_A^2}, \qquad
\mathbf r_A = \mathbf r - \mathbf R_A,
$$

with the gradient obtained by the ordinary product rule, e.g.
$\partial_x(x^ly^mz^ne^{-\alpha r^2}) = (lx^{l-1}y^mz^n - 2\alpha x^{l+1}y^mz^n)e^{-\alpha r^2}$; spherical
AOs apply the same solid-harmonic transform $\chi^{\text{sph}}_s=\sum_cC^{(l)}_{sc}\chi^{\text{cart}}_c$
used everywhere else in qc-rs (both to the value and to each gradient component). The restricted and
unrestricted densities and their gradients follow directly from the AO values and the density matrix:

$$
\rho(\mathbf r) = \sum_{\mu\nu}D_{\mu\nu}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r), \qquad
\nabla\rho^\sigma(\mathbf r) = \sum_{\mu\nu}D^\sigma_{\mu\nu}\bigl[(\nabla\chi_\mu)\chi_\nu +
\chi_\mu(\nabla\chi_\nu)\bigr],
$$

and libxc's GGA contracted-gradient invariants are $\sigma_{\alpha\alpha}=\nabla\rho^\alpha\cdot
\nabla\rho^\alpha$, $\sigma_{\alpha\beta}=\nabla\rho^\alpha\cdot\nabla\rho^\beta$,
$\sigma_{\beta\beta}=\nabla\rho^\beta\cdot\nabla\rho^\beta$ (restricted GGA needs only
$\sigma=\nabla\rho\cdot\nabla\rho$). A crucial regression test integrates the density back against the
overlap-consistent electron count, $\sum_gw_g\rho_g\approx\operatorname{Tr}[DS]$ — a discrepancy here means
either the grid or the AO evaluator has a bug, independent of any XC functional at all.

AO evaluation is never done densely over the whole molecule at once: each grid **block** screens for
*active* AOs — those whose magnitude exceeds a cutoff $\tau_\phi$ somewhere in the block,
$\max_g|\chi_\mu(\mathbf r_g)|>\tau_\phi$ — and only the active subset is evaluated, contracted into the
density, and later contracted back into $V_{\text{xc}}$. For a spatially extended molecule this keeps the
per-block AO buffer small regardless of total system size, the same screening philosophy as the Schwarz
shell-quartet screening used for the two-electron integrals.

## Assembling $E_{\text{xc}}$ and $V_{\text{xc}}$

Given libxc's per-particle energy density $\epsilon_{\text{xc}}$ and its functional derivatives (evaluated
pointwise, [the SCF chapter](../20-guide/scf.md) covers the wrapper), the XC energy is a straightforward
weighted sum,

$$
E_{\text{xc}} = \int\rho(\mathbf r)\,\epsilon_{\text{xc}}(\mathbf r)\,d\mathbf r \approx
\sum_g w_g\,\rho_g\,\epsilon_{\text{xc},g}.
$$

The Fock contribution follows from differentiating the discretized energy with respect to the density
matrix. For restricted LDA,

$$
V^{\text{xc}}_{\mu\nu} = \sum_g w_g\,v_{\rho,g}\,\chi_{\mu g}\chi_{\nu g},
$$

and restricted GGA adds the gradient term, using the product-gradient identity
$\nabla(\chi_\mu\chi_\nu)=(\nabla\chi_\mu)\chi_\nu+\chi_\mu(\nabla\chi_\nu)$:

$$
V^{\text{xc}}_{\mu\nu} = \sum_g w_g\Bigl[v_{\rho,g}\chi_{\mu g}\chi_{\nu g} +
2v_{\sigma,g}\,\nabla\rho_g\cdot\nabla(\chi_{\mu g}\chi_{\nu g})\Bigr].
$$

The unrestricted GGA form couples both spin channels through the cross gradient invariant
$\sigma_{\alpha\beta}$:

$$
V^{\alpha}_{\mu\nu} = \sum_g w_g\Bigl[v_{\rho^\alpha,g}\chi_{\mu g}\chi_{\nu g} +
\bigl(2v_{\sigma_{\alpha\alpha},g}\nabla\rho^\alpha_g + v_{\sigma_{\alpha\beta},g}\nabla\rho^\beta_g\bigr)
\cdot\nabla(\chi_{\mu g}\chi_{\nu g})\Bigr],
$$

with $\beta$ obtained by exchanging $\alpha\leftrightarrow\beta$ throughout. This slots into the composable
Fock builder alongside the Coulomb/exchange terms exactly as [the SCF chapter](../20-guide/scf.md)
describes: pure RKS is $F=H+J[D]+V_{\text{xc}}[D]$, global hybrids add $-a_xK[D^\sigma]$, and
range-separated hybrids further split $K$ into short/long-range pieces — the grid and libxc machinery in
this chapter is exactly what supplies $V_{\text{xc}}$ and $E_{\text{xc}}$ to that Fock expression,
regardless of which reference/hybrid combination is in play.

```{mermaid}
flowchart TD
    SPEC["grid= preset<br/>(coarse/medium/fine/ultrafine)"] --> RAD["Treutler-Ahlrichs<br/>radial shells per atom"]
    RAD --> PRUNE["NWChem 5-region<br/>angular pruning"]
    PRUNE --> PART["Becke space-partition<br/>W_A(r), sum_A W_A = 1"]
    PART --> BLOCK["grid blocks<br/>(workspace-lane sized)"]
    BLOCK --> AOSCREEN["AO screening<br/>+ AO/grad evaluation"]
    AOSCREEN --> DENS["rho, grad-rho<br/>(+ tau for meta-GGA)"]
    DENS --> LIBXC["libxc: eps_xc, v_rho, v_sigma, ..."]
    LIBXC --> ACC["accumulate E_xc<br/>+ V_xc into FockAccumulator"]
```

## Performance model: Workspace, screening, and threading

Every quantity in this pipeline scales with the number of grid points, which in turn scales linearly with
the number of atoms — so, unlike the two-electron integral tensor, the grid workload never explodes
combinatorially, but it still dominates wall-clock time for mid-size DFT jobs if handled naively. qc-rs's
performance discipline mirrors the [`Workspace` rule](../30-hpc/threads-and-blas.md) used everywhere else
in the codebase:

- **No hot-loop allocation.** The grid engine never owns its points/weights outside of debug/export modes
  (`store_points=True`); production code fills a run-scoped `Workspace` once and reborrows
  `MolecularGridView`/`GridBlock` slices, so iterating blocks allocates nothing.
- **Block-local AO screening** keeps the per-block active-AO buffer proportional to the *local* AO density
  near that block, not the whole-molecule AO count — the same locality argument that makes Schwarz
  screening effective for ERIs.
  With block size $B$ and $N_{\text{AO,active}}$ active AOs, an LDA worker lane needs roughly
  $M_{\text{lane}}^{\text{LDA}}\sim8B(N_{\text{AO,active}}+c_{\text{LDA}})$ bytes, and a GGA lane —
  needing AO gradients too — roughly $M_{\text{lane}}^{\text{GGA}}\sim8B(4N_{\text{AO,active}}+
  c_{\text{GGA}})$; the workspace planner picks the largest block size $B$ that keeps every worker lane and
  the reduction buffer inside the process memory budget.
- **Threading**: the grid-block loop is the one part of the SCF cycle that is genuinely thread-parallel
  today (LDA, GGA, and meta-GGA including $\tau$ and the Laplacian, both RKS and UKS) — each worker thread
  owns disjoint blocks and a private/tiled partial $V_{\text{xc}}$, reduced in a fixed thread-id order at
  the end for bit-reproducibility. Per-block contraction itself is **single-threaded** GEMM
  (`qc_array::blas::gemm_st`), so there is no nested oversubscription to coordinate between the outer
  thread-parallel grid loop and an inner multi-threaded BLAS call — one thread budget, no
  `BlasThreadGuard` needed (unlike, hypothetically, a design that ran multi-threaded BLAS *inside* each
  worker).

## The two energy formulas: RKS vs UKS

Once $E_{\text{xc}}$ and $V_{\text{xc}}^\sigma$ are known, the electronic energy for pure (non-hybrid) KS-DFT
follows the same trace structure as HF, replacing exact exchange with the XC functional:

$$
\text{RKS:}\quad E_{\text{elec}} = \operatorname{Tr}[DH] + \tfrac12\operatorname{Tr}[DJ] + E_{\text{xc}},
\qquad F = H + J[D] + V_{\text{xc}}[D],
$$

$$
\text{UKS:}\quad E_{\text{elec}} = \operatorname{Tr}\bigl[(D^\alpha{+}D^\beta)H\bigr] +
\tfrac12\operatorname{Tr}\bigl[(D^\alpha{+}D^\beta)J\bigr] + E_{\text{xc}}, \qquad
F^\sigma = H + J[D^\alpha{+}D^\beta] + V_{\text{xc}}^\sigma.
$$

Neither formula has an exchange term at all for a pure functional — this is the fundamental structural
difference from Hartree-Fock, not merely a different numeric coefficient, and it is why a pure GGA/LDA run
costs no two-electron exchange contraction while a hybrid (`xc="b3lyp"`, `xc="pbe0"`) needs both this grid
machinery *and* the exact-exchange $K$ build from [the SCF chapter](../20-guide/scf.md).

:::{exercise}
:label: ex-dft-grid-theory

1. Water/def2-SVP/B3LYP gives energies $-76.358122100$, $-76.358125115$, and $-76.358124968\ E_h$ on
   `coarse`, `medium`, `fine` respectively. Using the calibrated worst-case tolerances
   ($\lesssim8$, $\lesssim0.4$, $\lesssim0.2\ \mu E_h$/atom, 3 atoms), is the observed
   `coarse`$\to$`medium` shift ($\approx3\times10^{-6}\ E_h$) consistent with those tolerances, and what
   does the much smaller `medium`$\to$`fine` shift tell you about where `medium` sits on the accuracy
   ladder?
2. The Becke partition uses Bragg-Slater radii to correct $\mu_{AB}$ into $\nu_{AB}$. Explain in one
   sentence why using the *raw* $\mu_{AB}$ (the naive perpendicular-bisector boundary) would be a poor
   choice for a heavy atom bonded to hydrogen.
3. Ghost atoms carry basis functions but are never partition centers. What would go wrong with the Becke
   partition's normalization $\sum_AW_A=1$ if a ghost atom *were* given a Bragg-Slater radius and treated as
   a center?
:::

:::{solution} ex-dft-grid-theory
:class: dropdown

1. Yes: the observed $\approx3\times10^{-6}\ E_h$ (3 $\mu E_h$) `coarse`$\to$`medium` shift is comfortably
   inside the $3\times8=24\ \mu E_h$ worst-case ceiling for `coarse`, and well inside the corresponding
   `medium` ceiling ($3\times0.4=1.2\ \mu E_h$) once `medium` has itself converged — consistent with
   `medium` already being close to its own accuracy floor for this small, well-behaved molecule. The tiny
   `medium`$\to$`fine` shift (sub-$\mu E_h$, changing only the 9th digit) confirms `medium` is already deep
   in the asymptotic convergence region for this system — exactly why it is the sensible default rather
   than `fine`.
2. A heavy atom and a hydrogen have very different Bragg-Slater radii, so the physically reasonable
   dividing surface between "belongs to the heavy atom" and "belongs to hydrogen" is much closer to the
   hydrogen than the geometric midpoint — using the raw, radius-blind $\mu_{AB}$ would assign hydrogen an
   unreasonably large share of space (and therefore radial/angular points) around a bond it barely spans,
   wasting quadrature effort where the density is dominated by the heavy atom's core.
3. $\sum_AW_A(\mathbf r)=1$ only holds because *every* atom that could compete for a point's partition
   weight is included in the same normalizing sum $\sum_CP_C(\mathbf r)$ used to build every $W_A$. If a
   ghost atom were silently given a competing cell function $P_{\text{ghost}}$ without also being included
   consistently everywhere the partition is built and used, some fraction of real-atom weight would be
   siphoned into a center with no nucleus and no density contribution of its own there, breaking the
   completeness identity and systematically undercounting $\int\rho\,d\mathbf r$ near the ghost center —
   which is exactly why qc-rs's default policy excludes ghost/dummy atoms from the partition centers.
:::

The block-local AO/density machinery here reappears verbatim in [linear-response
theory](linear-response-cphf.md), where it evaluates *transition* densities instead of the ground-state
one, and in [analytic derivatives](analytic-derivatives.md), where the Becke weight derivative
$\partial W_A/\partial R_C$ derived from the same $s_h(\nu_{AB})$ switching function becomes part of the
DFT nuclear gradient.