# Topological analysis theory: QTAIM & ELF [The population-analysis chapter](population-analysis-theory.md) derived a family of ways to divide the density among atoms — every one of them a *partitioning choice*. This chapter instead studies the electron density (and closely related scalar fields) directly, as a genuine mathematical object with its own critical points, basins, and localization structure — extracting chemical meaning (bonds, lone pairs, rings, cages) from the topology of $\rho(\mathbf r)$ itself, with no partitioning decision to make at all. This is Bader's **Quantum Theory of Atoms in Molecules (QTAIM)** and Becke-Edgecombe's **Electron Localization Function (ELF)**, grounded directly in `crates/qc-prop/src/{qtaim,basin,realspace}.rs`. ## QTAIM: the density as a topological object The electron density is a smooth, everywhere-defined scalar field with a single, chemically decisive structural fact: it has local maxima almost exactly at every nucleus (the nuclear cusp), and *saddle points* connecting pairs of those maxima along what turns out to be exactly the chemical bond paths a structural chemist would draw by hand. QTAIM makes this rigorous by studying $\nabla\rho$'s **critical points** — points where $\nabla\rho=\mathbf 0$ — and classifying each by the eigenvalue signature of the density Hessian there. ### Classifying critical points by curvature signature At a critical point, the (real, symmetric) Hessian of $\rho$ has three real eigenvalues; every chemically meaningful critical point is *nondegenerate* (all three eigenvalues nonzero — a genuine local extremum or saddle, not a flat/degenerate direction), and its type is decided purely by how many of the three eigenvalues are negative: | Signature $(3,\sigma)$ | Negative eigenvalues | Type | Chemical meaning | |---|---:|---|---| | $(3,-3)$ | 3 | **Nuclear (NCP)** | a density maximum — sits on (almost exactly) a nucleus | | $(3,-1)$ | 2 | **Bond (BCP)** | a saddle between two bonded atoms — density decreasing in the two directions perpendicular to the bond, increasing along it | | $(3,+1)$ | 1 | **Ring (RCP)** | the density minimum *within* a ring of bonds | | $(3,+3)$ | 0 | **Cage (CCP)** | a genuine local density minimum, enclosed by rings | $\sigma$ (the trace's sign count: $+1$ per positive eigenvalue, $-1$ per negative) together with the rank (always 3 for these nondegenerate types) is the traditional $(3,\sigma)$ notation; a Bond CP's defining property — two negative curvatures perpendicular to the bond path and one positive curvature along it — is precisely the mathematical statement "this is a mountain pass between two peaks," which is exactly the topological signature of a genuine chemical bond in this framework, not merely a suggestive analogy. ### The Poincaré-Hopf invariant: a free, whole-molecule consistency check A remarkable topological fact ties every critical point count together for a *complete*, isolated molecular system, independent of which molecule it is: $$ n_{\text{NCP}} - n_{\text{BCP}} + n_{\text{RCP}} - n_{\text{CCP}} = 1. $$ This is the Poincaré-Hopf index theorem applied to the density's gradient vector field, and it holds for *any* correctly and completely characterized molecular topology — water's 3 nuclear + 2 bond + 0 ring + 0 cage critical points give $3-2+0-0=1$; benzene's 12 nuclear + 12 bond + 1 ring + 0 cage give $12-12+1-0=1$; a tetrahedral $P_4$ cage's 4 nuclear + 6 bond + 4 ring + 1 cage give $4-6+4-1=1$. Any other value means a critical point was missed (an expected bond or ring point the search failed to locate) or a spurious one was found — this identity is therefore a genuinely free correctness check on a numerical critical-point search, needing no reference calculation to compare against, only internal consistency. ## Bond-critical-point energetics: the local virial theorem A bond critical point's position and density value alone already carry qualitative information, but Bader's **local virial theorem** extracts genuine *energy densities* at that single point, turning a purely topological classification into an energetic one. The theorem relates the density's Laplacian to the local kinetic and potential energy densities: $$ \tfrac14\nabla^2\rho = 2G+V, \qquad G(\mathbf r) = \tfrac12\sum_{\mu\nu}D_{\mu\nu}\nabla\varphi_\mu\cdot \nabla\varphi_\nu\ (\ge0), $$ with $G$ the (always positive) Lagrangian kinetic energy density. This one relation, together with the already-known $\rho$, $\nabla^2\rho$, and $G$ at the critical point, is enough to solve directly for the potential energy density and the total local energy density: $$ V = \tfrac14\nabla^2\rho - 2G\ (\le0\text{ in bonding regions}), \qquad H = G+V = \tfrac14\nabla^2\rho - G. $$ Two dimensionless ratios built from these classify the bond's character directly from local energetics rather than from the Hessian signature alone: the **bond-character ratio** $|V|/G$ (empirically, $>2$ signals a covalent/shared-shell interaction, $<1$ a closed-shell/ionic-or-van-der-Waals one, with a partially-covalent regime in between) and the **bond degree** $H/\rho$ (negative for covalent bonds — the total local energy density is stabilizing — positive for closed-shell interactions). Verified example — water's full QTAIM topology, computed and classified automatically: ```python import qc water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469" m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r").run() topo = m.prop.qtaim.topology() topo["counts"] # {'nuclear': 3, 'bond': 2, 'ring': 0, 'cage': 0} topo["poincare_hopf"] # {'sum': 1, 'holds': True} bcp = topo["critical_points"][3] # one O-H bond critical point bcp["rho"], bcp["laplacian"] # (0.3659441002679062, -2.5397077300947384) bcp["v_over_g"], bcp["h_over_rho"] # (11.152930010504011, -1.9245989766373959) ``` Both O–H bond CPs (which are, by the molecule's $C_{2v}$ symmetry, equivalent) give $|V|/G\approx11.2\gg2$ and $H/\rho\approx-1.92<0$ — an unambiguous covalent classification, exactly as chemical intuition demands for an O–H single bond, but now derived from a rigorous local-energy-density calculation at a single, mathematically well-defined point rather than from any partitioning scheme. ### Bond ellipticity The two *negative* eigenvalues at a bond CP ($\lambda_1,\lambda_2$, with $|\lambda_1|\ge|\lambda_2|$ by convention) need not be equal in magnitude — real bonds are not perfectly cylindrically symmetric — and their ratio defines the **bond ellipticity**, $$ \varepsilon = \frac{\lambda_1}{\lambda_2} - 1, $$ which measures how much more the density is depleted in one perpendicular direction than the other. A perfectly cylindrical bond (a genuine single bond, or a triple bond with full axial symmetry) has $\varepsilon\approx0$; a double bond with a preferential $\pi$-accumulation direction shows a measurably larger ellipticity. The water O–H bond CP above has $\varepsilon\approx0.024$ — small, consistent with an ordinary, close-to-cylindrical single bond with no special $\pi$ character. ## Basin integration: Bader charges from the density's own watershed The critical points alone classify bonds; **basin integration** goes further, using the topology to *partition all of space* without any external reference density (unlike every scheme in [the population-analysis chapter](population-analysis-theory.md)) — each atomic nucleus is an attractor of $\nabla\rho$'s gradient flow, and its **basin** $\Omega_A$ is the set of every point whose steepest-ascent path in $\rho$ terminates there. The basins are bounded by **zero-flux surfaces** ($\nabla\rho(\mathbf r)\cdot\hat{\mathbf n}(\mathbf r)=0$ for every point $\mathbf r$ on the surface) — surfaces the gradient field never crosses, which is exactly what makes the resulting atomic populations $N_A=\int_{\Omega_A}\rho\,d\mathbf r$ a genuine, unambiguous *partition* of the total electron count, with no smooth weight function and no arbitrary reference density anywhere. Numerically, on a cubic real-space grid, qc-rs assigns every grid point to a basin by **on-grid steepest ascent** (a Henkelman-style algorithm): walk from each point toward its distance-weighted, density-increasing neighbor until a local maximum (a nucleus) is reached, with path compression so repeated walks from nearby points are cheap. Verified example — water's Bader (basin-integrated) atomic populations and charges: ```python b = m.prop.qtaim.basin_integrate() b["atoms"] # O: 9.2826916251811 electrons -> q = 8 - 9.283 = -1.283 # H: 0.37835478988956356 -> q = 1 - 0.378 = +0.622 # H: 0.3626332357145994 -> q = 1 - 0.363 = +0.637 ``` Bader charges here ($q_O\approx-1.28$) are noticeably *larger in magnitude* than every scheme derived in [the population-analysis chapter](population-analysis-theory.md) — even NPA/MBIS/Hirshfeld-I's $\sim-0.87$ to $-0.92$. This is not a bug or an outlier: QTAIM's zero-flux partitioning genuinely divides space by a different, topology-driven criterion than any stockholder-weight or AO-index scheme, and (as that chapter's own closing point emphasized) there is no single "true" atomic charge to converge toward — QTAIM is simply one more, independently well-motivated definition, not a more-correct version of the others. ## ELF and LOL: localization without any critical-point search A structurally different real-space tool asks not "where are the topological features of $\rho$ itself" but "how spatially localized is the electron pair density at each point" — the **Electron Localization Function (ELF)**, due to Becke and Edgecombe. Its construction starts from the **Pauli kinetic energy density**, the excess of the true (positive-definite, Kohn-Sham-orbital) kinetic energy density $\tau$ over the von Weizsäcker kinetic energy $\tau_W$ a single-orbital (bosonic, fully delocalized) description of the same density would have: $$ D(\mathbf r) = \tau(\mathbf r) - \tau_W(\mathbf r), \qquad \tau_W = \frac{|\nabla\rho|^2}{8\rho}. $$ $D$ vanishes exactly wherever the local density is genuinely describable by a single orbital (no Pauli repulsion penalty at all — the theoretical floor), which is why regions of high localization (covalent bonds, lone pairs, core shells) show small $D$. To turn this into a bounded, chemically interpretable number, $D$ is compared against the same quantity for a **uniform electron gas** at the local density (the Thomas-Fermi kinetic energy density $D_h=\tfrac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}$, a purely local function of $\rho$ alone): $$ \chi(\mathbf r) = \frac{D(\mathbf r)}{D_h(\mathbf r)}, \qquad \text{ELF}(\mathbf r) = \frac{1}{1+\chi^2} \in(0,1]. $$ The normalization is deliberately chosen so that a region genuinely described by one orbital gives $D=0\Rightarrow\chi=0\Rightarrow\text{ELF}=1$ exactly (an exact identity qc-rs's test suite checks directly, not merely a limiting approximation), and a region with exactly the uniform-electron-gas kinetic energy density gives $D=D_h\Rightarrow\chi=1\Rightarrow\text{ELF}=\tfrac12$ exactly — the canonical reference point that gives the function its interpretive scale (ELF above $\tfrac12$ means "more localized than a uniform electron gas at this density," below $\tfrac12$ means less). In practice, beyond a molecule's outer boundary $D$ can decay *faster* than $D_h$, spuriously driving ELF toward 1 in the vacuum tail (a "completely localized vacuum," clearly unphysical); qc-rs follows Multiwfn's convention of adding a small floor to $D$ specifically to tame this tail without perturbing the chemically interesting bonding/lone-pair regions. **LOL** (the Localized Orbital Locator, Schmider-Becke) is built from the *same* two ingredients, $\tau$ and $D_h$, but combined differently — a ratio rather than a normalized-square form: $$ t(\mathbf r) = \frac{D_h(\mathbf r)}{\tau(\mathbf r)}, \qquad \text{LOL}(\mathbf r) = \frac{t}{1+t} \in[0,1), $$ giving a bounded function with the same qualitative localization-mapping purpose as ELF but a different numerical scale and (in practice) somewhat different visual emphasis on core-shell structure — the two are complementary presentations of closely related underlying physics, not competing theories. A related quantity, the **Electron Localizability Indicator (ELI-D)**, $\Upsilon_D=\rho_\sigma[12/g^\sigma]^{3/8}$ with $g^\sigma=\rho_\sigma D_\sigma$ (the same Pauli kinetic energy density numerator ELF uses), shares ELF's exact topology for a single-determinant wavefunction — it marks the same shells, bonds, and lone pairs at the same locations — but needs **no uniform-electron-gas reference at all**, and its value has an absolute physical meaning (electrons per fixed fraction of a same-spin electron pair, with the uniform gas giving $\approx1.10$ as a universal constant rather than an arbitrary normalization point) rather than ELF's relative, gas-normalized scale. :::{exercise} :label: ex-topology-theory 1. The Poincaré-Hopf sum for a complete molecular topology must equal exactly $+1$. If a numerical critical-point search for a moderately large, ring-containing molecule reports a sum of $0$, what does that tell you has almost certainly gone wrong, without needing any reference calculation to compare against? 2. Water's Bader charge on oxygen ($\approx-1.28$) is noticeably larger in magnitude than its NPA charge ($\approx-0.91$, [from the previous chapter](population-analysis-theory.md)). Explain in one sentence why this difference does not mean one of the two methods is simply wrong. 3. ELF and ELI-D share the exact same topology (same critical points, same qualitative bonding/lone-pair picture) for a single-determinant wavefunction, yet ELI-D needs no uniform-electron-gas reference while ELF's entire construction is built around comparing against one. What does this tell you about whether the uniform-electron-gas comparison is *essential physics* or merely a *convenient normalization choice* for capturing the same underlying localization information? ::: :::{solution} ex-topology-theory :class: dropdown 1. A sum of $0$ instead of $1$ means the search is missing exactly one critical point relative to a complete, self-consistent topology (or, less commonly, found one spurious extra point) — for a ring-containing molecule specifically, the most likely culprit is a missed **ring critical point** (RCPs are notoriously easy for a naive numerical search to skip, since they sit in the geometric middle of a ring, away from any bond path, at a comparatively flat and low-density region of the density surface). This diagnosis requires no comparison to any other calculation — it follows purely from the topological identity itself always needing to equal 1 for a complete, correctly characterized topology. 2. QTAIM's zero-flux basin partition and NPA's natural-atomic-orbital partition are answering the *same qualitative question* (how much electron density/population belongs to this atom) using *genuinely different mathematical partitioning criteria* — one a real-space topological watershed of $\nabla\rho$, the other an orthonormal natural-orbital occupation analysis in AO-index space. As emphasized throughout [the population-analysis chapter](population-analysis-theory.md), there is no single "true" atomic charge that both methods are independently trying (and one failing) to approximate; a difference in magnitude between two self-consistent, well-defined partitioning schemes is expected, not evidence that either is in error. 3. It tells you the uniform-electron-gas comparison is a **convenient normalization choice**, not essential physics — the underlying localization information (captured by the Pauli kinetic energy density $D$, or equivalently the same-spin Fermi-hole curvature $g^\sigma$) is exactly the same input both functions build from, and the topology (where $D$ or $g^\sigma$ is small vs. large, i.e. where electrons are localized vs. delocalized) is entirely determined by that shared input alone. Comparing $D$ against a uniform-gas reference (ELF's $\chi=D/D_h$) versus reporting a value with its own absolute physical scale (ELI-D's $\Upsilon_D$) are two different *presentations* of the identical underlying physical content — neither presentation choice changes where the localized/delocalized regions actually are. ::: Topology and localization complete this manual's real-space toolkit for reading chemical structure directly out of a converged density. [The reactivity-theory chapter](reactivity-theory.md) turns from *structure* to *reactivity* — conceptual DFT descriptors (electronegativity, hardness, Fukui functions) and aromaticity indices that predict how a molecule *will* behave, not just how its density is currently arranged.