Topological analysis theory: QTAIM & ELF

The population-analysis chapter 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:

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

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 — 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 17

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). 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 to Exercise 17

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