# QTAIM & ELF/LOL

The [previous chapter](charges-bond-orders.md) 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](../../10-foundations/topological-analysis-theory.md).

### Usage

```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="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:

```python
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](charges-bond-orders.md), 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

```python
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](../visualization.md) draws the ELF isosurface.

## Worked example

```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="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}
:label: ex-qtaim

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](charges-bond-orders.md) 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} ex-qtaim
:class: dropdown

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](weak-interactions.md) turns to the
weak, non-covalent interactions *between* fragments — hydrogen bonds and van der Waals contacts — with NCI
and IGM.