DFT exchange-correlation quadrature

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:

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

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

Exercise 7

  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?

The block-local AO/density machinery here reappears verbatim in linear-response theory, where it evaluates transition densities instead of the ground-state one, and in analytic derivatives, 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.