QTAIM & ELF/LOL

The previous chapter partitioned the density with recipes. This one lets the density partition itself. Two topological analyses — QTAIM (the quantum theory of atoms in molecules) and ELF (the electron localization function) — find atoms, bonds, and lone pairs from the shape of a real-space field, with no arbitrary basis choice. They are the heart of “Multiwfn-class” analysis.

QTAIM: atoms from the density topology

Theory

Bader’s QTAIM analyzes the gradient field of the electron density \(\nabla\rho(\mathbf r)\). Its critical points are where \(\nabla\rho=0\), and each is classified by the signature of the local density Hessian — its number of negative eigenvalues, written \((3,\sigma)\) where \(\sigma\) is (positive count) minus (negative count):

Type	Signature	Hessian eigenvalues	Meaning
Nuclear (NCP)	\((3,-3)\)	all 3 negative	a density maximum — sits at a nucleus
Bond (BCP)	\((3,-1)\)	2 negative, 1 positive	a saddle along the bond, a maximum across it
Ring (RCP)	\((3,+1)\)	1 negative, 2 positive	a saddle inside a ring
Cage (CCP)	\((3,+3)\)	all 3 positive	a local minimum inside a cage

The bond path linking two nuclei through a BCP is QTAIM’s rigorous definition of a chemical bond, and a sanity check — the Poincaré–Hopf relation \(n_{\text{NCP}} - n_{\text{BCP}} + n_{\text{RCP}} - n_{\text{CCP}} = 1\) — must hold for a correct topology (it is a topological invariant of the density field, independent of the molecule). QTAIM also carves space into atomic basins — regions bounded by zero-flux surfaces \(\nabla\rho(\mathbf r)\cdot\mathbf n(\mathbf r)=0\), i.e. surfaces the density’s gradient never crosses; integrating \(\rho\) over a basin gives a basis-set-robust atomic charge (the “Bader charge”), robust precisely because the zero-flux boundary is a property of \(\rho\) itself, not of the AO basis. The full derivation, including how basins are found by steepest-ascent integration, is in Topological analysis theory.

Usage

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="cc-pvdz", unit="angstrom").scf(ref="r").run()

topo = qc.prop.qtaim.topology(m)
topo["counts"]         # {'nuclear': 3, 'bond': 2, 'ring': 0, 'cage': 0}
topo["poincare_hopf"]  # {'sum': 1, 'holds': True}
topo["critical_points"][0]   # {'type': 'nuclear', 'position': [...], 'rho': 297.19, 'atoms': [0]}

For water, QTAIM finds exactly 3 nuclear + 2 bond critical points (the two O–H bonds) and no rings or cages, and the Poincaré–Hopf sum is 1 — a topologically consistent result. Integrating the atomic basins gives the Bader populations:

bader = qc.prop.qtaim.basin_integrate(m)
bader["atoms"][0]      # {'atom': 0, 'label': 'O', 'integral': 9.244}  -> Bader charge 8 - 9.244 = -1.24
bader["total"]         # 10.023   (≈ the 10 electrons; small grid residual)

Oxygen’s basin holds ~9.24 electrons — a Bader charge of about −1.24, the most negative of all the schemes in the charges chapter, and the most physically grounded. qc.prop.qtaim also exposes bond-critical-point properties, basin multipoles, delocalization indices, and the source function.

Tip

QTAIM is grid-based — expect it to be slower Topology and basin integration walk a real-space grid, so they cost more than the algebraic charges of the last chapter. It is still a single call; just do not be surprised when it takes a few seconds on a bigger molecule.

ELF: seeing electron pairs

Theory

The electron localization function \(\text{ELF}(\mathbf r) \in [0,1]\) measures how strongly electrons are localized — how much the Pauli principle keeps a like-spin electron away from a reference point. Becke and Edgecombe’s construction compares the true kinetic energy density to what a uniform electron gas of the same local density would have:

\[ D(\mathbf r) = \tau(\mathbf r) - \frac{|\nabla\rho(\mathbf r)|^2}{8\rho(\mathbf r)}, \qquad D_h(\mathbf r) = \frac{3}{10}(3\pi^2)^{2/3}\rho(\mathbf r)^{5/3}, \qquad \text{ELF} = \frac{1}{1+\chi^2},\quad \chi = \frac{D}{D_h}. \]

\(D\) is the Pauli kinetic energy density — the excess kinetic energy beyond the lowest possible for that density, which vanishes only where a single orbital (or a same-spin electron pair) dominates. Where the true system localizes electrons more than the uniform gas would, \(D \ll D_h\), so \(\chi\to0\) and \(\text{ELF}\to1\); a region behaving exactly like a uniform gas gives \(\chi=1\), \(\text{ELF}=0.5\) (qc-rs’s implementation verifies both limits directly). This is why ELF reveals the electron pairs of a Lewis structure directly: cores, bonds, and lone pairs each show up as a separate region of high ELF (a basin), letting you see bonding and lone pairs without invoking orbitals at all.

Usage

basins = qc.prop.elf.basins(m)      # the ELF basins (water: 5 — an O core, two O–H bonds, two lone pairs)
domains = qc.prop.elf.domains(m)    # localization domains at an isovalue

Water’s ELF resolves into 5 basins — the oxygen core, the two O–H bonding basins, and the two oxygen lone pairs — recovering the familiar Lewis picture from the wavefunction alone. qc.prop.elf also provides LOL (an ELF-like localized-orbital locator) and core/valence bifurcation analysis. To see any of these, m.view3d("elf") from the visualization chapter draws the ELF isosurface.

Worked example

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="cc-pvdz", unit="angstrom").scf(ref="r").run()

topo = qc.prop.qtaim.topology(m)
print("critical points:", topo["counts"], "| Poincaré–Hopf holds:", topo["poincare_hopf"]["holds"])
print("O Bader population:", round(qc.prop.qtaim.basin_integrate(m)["atoms"][0]["integral"], 3))
print("ELF basins        :", len(qc.prop.elf.basins(m)))

Exercise 12

1. A QTAIM analysis of benzene should report how many ring critical points, and what must the Poincaré–Hopf sum equal for the topology to be consistent?

2. The charges chapter gave oxygen charges from −0.09 (Löwdin) to −0.91 (NPA). Where does the QTAIM (Bader) charge of ≈ −1.24 sit, and why is it considered especially well-defined?

3. You want to see the two lone pairs on water’s oxygen. Which ELF call, and which viewer call, do you use?

Solution to Exercise 12

1. One ring critical point (the center of the aromatic ring). Poincaré–Hopf must equal 1: \(n_{\text{NCP}} - n_{\text{BCP}} + n_{\text{RCP}} - n_{\text{CCP}} = 12 - 12 + 1 - 0 = 1\) for benzene.

2. It is the most negative — even beyond NPA. It is well-defined because the basin boundary is a zero-flux surface of the density (a property of \(\rho\) itself, not of the basis functions), so the Bader charge is essentially basis-set-independent.

3. Compute with qc.prop.elf.basins(m) / qc.prop.elf.domains(m), and visualize with m.view3d("elf") — the two lone-pair basins appear as distinct lobes above the oxygen.

QTAIM and ELF analyze a single molecule’s bonding. The next chapter turns to the weak, non-covalent interactions between fragments — hydrogen bonds and van der Waals contacts — with NCI and IGM.