# Conceptual DFT & reactivity theory

Every chapter so far in Part II has asked "what is the energy, and how is the density structured" for a
*fixed* molecule at a *fixed* electron count. This chapter asks a different kind of question: how would the
energy change if the electron count itself changed, and which atom in a molecule is most reactive toward a
nucleophile, an electrophile, or a radical? **Conceptual DFT** answers the first with a handful of global
scalar descriptors derived directly from density-functional perturbation theory in $N$; the **Fukui
function** answers the second by asking how the *density itself* responds locally to that same
perturbation. This chapter also derives the electron-delocalization family of **aromaticity indices** —
PDI, FLU, $I_{\text{ring}}$, MCI — a genuinely different, wavefunction-based complement to the purely
geometric HOMA index [already introduced in the guide chapter](../20-guide/properties/aromaticity.md).
Grounded directly in `crates/qc-prop/src/{cdft,fukui,aromaticity,homa,bird}.rs`.

## Conceptual DFT: reactivity from $\partial E/\partial N$

Kohn-Sham DFT's ground-state energy is, formally, a functional of the electron number $N$ as well as the
external potential — and its derivatives with respect to $N$ at fixed geometry are exactly the
chemically familiar concepts of electronegativity and hardness, made mathematically precise rather than
merely qualitative. The two energies that anchor everything here are the **vertical ionization potential**
$I$ and **vertical electron affinity** $A$ — energies of the cation and anion *at the neutral molecule's
own geometry*, not their own relaxed geometries — obtainable two ways: **Koopmans' theorem**
($I\approx-\varepsilon_{\text{HOMO}}$, $A\approx-\varepsilon_{\text{LUMO}}$, free from the neutral SCF's
own orbital energies, no extra calculation) or **$\Delta$SCF** ($I=E(N-1)-E(N)$, $A=E(N)-E(N+1)$, two
extra single-point energies, more accurate but requiring real electron-count perturbation).

From $(I,A)$ alone, the finite-difference approximation to $\partial E/\partial N$ and $\partial^2E/\partial
N^2$ gives the **chemical potential** and **hardness** directly:

$$
\mu = -\frac{I+A}{2}, \qquad \chi = -\mu = \frac{I+A}{2}\ \text{(Mulliken electronegativity)}, \qquad
\eta = I-A\ \Bigl(=\frac{\partial^2E}{\partial N^2}\Bigr), \qquad S=\frac1\eta.
$$

$\mu$ is literally the (finite-difference) slope of $E(N)$ at the neutral molecule — the energy cost per
electron added or removed, equalized across the whole molecule at equilibrium, which is exactly the
chemical meaning "electronegativity" is supposed to capture: a substance with a very negative $\mu$ (large
$\chi$) pulls electron density toward itself. $\eta$, the curvature, measures resistance to that
electron-transfer process — a large $\eta$ (a "hard" species) resists both gaining and losing electron
density, while a small $\eta$ (a "soft" species, large $S=1/\eta$) is polarizable and reactive.

Combining $\mu$ and $\eta$ gives several higher-order descriptors, each answering a more specific
reactivity question than electronegativity/hardness alone:

$$
\omega = \frac{\mu^2}{2\eta}\ \text{(electrophilicity index, Parr-Szentpaly-Liu)}, \qquad
\Delta N_{\max} = -\frac{\mu}{\eta}\ \text{(maximum electron flow to a perfect electron bath)},
$$

$$
\omega^- = \frac{(3I+A)^2}{16(I-A)}\ \text{(electrodonating power)}, \qquad
\omega^+ = \frac{(I+3A)^2}{16(I-A)}\ \text{(electroaccepting power, Gazquez-Cedillo-Vela)},
$$

$$
\Delta_{\text{nu}} = \frac{(\mu+\eta)^2}{2\eta}\ \text{(nucleofugality)}, \qquad
\Delta_{\text{el}} = \frac{(\mu-\eta)^2}{2\eta}\ \text{(electrofugality, Ayers-Anderson-Bartolotti)}, \qquad
\Delta\omega^\pm = \omega^++\omega^-.
$$

The electrophilicity index $\omega$ is the most widely used of these — it measures the energy a species
*stabilizes by* upon accepting the maximum charge flow a perfect electron reservoir would deliver, so a
large $\omega$ identifies a good electrophile independent of whether it is $I$ or $A$ individually driving
that behavior. Nucleofugality and electrofugality quantify a related but distinct question — not "how
reactive is this whole molecule" but "how readily would this species depart as an anion (nucleofugality)
or a cation (electrofugality)" once bonded elsewhere, the leaving-group question rather than the
attacking-nucleophile/electrophile question.

Verified example — water's global reactivity descriptors (from $\Delta$SCF-quality $I,A$, in eV):

```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()
r = m.prop.cdft.reactivity()
r["ip"], r["ea"]                    # (13.554, -4.795)  eV
r["electronegativity"], r["hardness"]  # (4.379, 18.349)  eV
r["electrophilicity"]               # 0.523  eV -- small: water is a poor electrophile
```

Water's large hardness ($\eta\approx18.3\ \text{eV}$) and small electrophilicity ($\omega\approx0.52\
\text{eV}$) are exactly what chemical intuition predicts for a small, closed-shell, kinetically inert
molecule with a large HOMO-LUMO gap — it strongly resists both gaining and losing electron density, and
correspondingly has little tendency to stabilize by accepting charge from a reservoir.

## Fukui functions: where global reactivity lands locally

Global descriptors answer "how reactive is this molecule overall," but not "which atom reacts first." The
**Fukui function** $f(\mathbf r)=\partial\rho(\mathbf r)/\partial N$ answers exactly that — how the density
at a specific point responds to adding or removing an electron — and its condensed-to-atoms form gives a
per-atom reactivity map directly comparable across a molecule's different sites:

$$
f_A^+ = q_A(N) - q_A(N+1)\ \text{(nucleophilic attack site — favors gaining electron density)}, \qquad
f_A^- = q_A(N-1) - q_A(N)\ \text{(electrophilic attack site)},
$$

$$
f_A^0 = \tfrac12(f_A^++f_A^-)\ \text{(radical attack site)}, \qquad
\Delta f_A = f_A^+-f_A^-\ \text{(dual descriptor — sign identifies the dominant character)}.
$$

qc-rs implements the **frozen-molecular-orbital (FMO)** approximation to these — a parameter-free,
essentially free-to-compute shortcut that needs only the *neutral* molecule's converged wavefunction, no
extra $N\pm1$ single-point calculations at all. The approximation replaces the true density-difference
Fukui functions with the frontier-orbital densities directly, $f^+(\mathbf
r)\approx|\varphi_{\text{LUMO}}(\mathbf r)|^2$ and $f^-(\mathbf r)\approx|\varphi_{\text{HOMO}}(\mathbf
r)|^2$ — physically, "where would the next electron go" is approximated by "where the LUMO already has
amplitude," and symmetrically for removal and the HOMO. Condensed to atoms via a Mulliken-style orbital
composition, this becomes exactly the frontier rows of that composition matrix:

$$
f_A^+ = \text{comp}_{\text{LUMO},A}, \qquad f_A^- = \text{comp}_{\text{HOMO},A},
$$

and because every row of a Mulliken orbital-composition matrix sums to exactly 1 by construction, the sum
rules $\sum_Af_A^+=\sum_Af_A^-=1$ and $\sum_A\Delta f_A=0$ hold *exactly*, not approximately — a genuinely
free correctness check on the condensed values with no reference calculation needed, in the same spirit as
[QTAIM's Poincaré-Hopf identity](topological-analysis-theory.md). (A more accurate, but more expensive,
finite-difference Fukui variant — differencing real $N,N\pm1$ SCF atomic charges rather than approximating
with frozen frontier orbitals — is a natural future addition; the FMO form is what qc-rs implements today.)

Verified example — water's condensed Fukui functions:

```python
f = m.prop.cdft.fukui()
f["fukui_minus"]      # [0.665, 0.167, 0.167]  -- HOMO composition: O dominates (electrophilic attack site)
f["fukui_plus"]       # [0.237, 0.381, 0.381]  -- LUMO composition: H dominates (nucleophilic attack site)
f["dual_descriptor"]  # [-0.428, 0.214, 0.214] -- negative at O (electrophilic), positive at H (nucleophilic)
```

Water's oxygen dominates the HOMO (its lone pairs are the highest-energy occupied density, exactly where
an electrophile would attack) while the hydrogens dominate the LUMO (the $\sigma^*$ character a
nucleophile's electron pair would flow into) — the dual descriptor's sign flip between O (negative,
electrophilic) and H (positive, nucleophilic) makes this qualitative picture a single, unambiguous number
per atom.

## Aromaticity beyond geometry: electron-delocalization indices

[The aromaticity chapter in the guide](../20-guide/properties/aromaticity.md) already introduced **HOMA**
(Harmonic Oscillator Model of Aromaticity), a purely geometric index built from how close each ring bond's
length sits to a tabulated "ideal aromatic" reference length. HOMA's structural limitation is exactly that
it *is* purely geometric — it says nothing about the actual electronic delocalization the geometry is
supposed to be a proxy for. The **electron-delocalization** aromaticity family instead measures
delocalization directly from the wavefunction, built on the same **delocalization index**
$\delta(A,B)$ (a measure of shared electron pairs between atoms $A,B$, closely related in spirit to [the
Mayer bond order](population-analysis-theory.md)) and the **atomic overlap matrices** $S(A)$ in the
occupied-MO basis:

$$
\text{PDI} = \text{mean of the three }\textit{para}\text{ }\delta\text{'s in a 6-ring}\ \text{(higher =
more aromatic)}, \qquad
\text{FLU} = \Bigl\langle\Bigl(\frac{\delta_i-\delta_i^{\text{ref}}}{\delta_i^{\text{ref}}}\Bigr)^2
\Bigr\rangle\ \text{(0 = aromatic)},
$$

$$
I_{\text{ring}} = 2^n\operatorname{Tr}\bigl[S(A_1)S(A_2)\cdots S(A_n)\bigr]\ \text{(the genuine
}n\text{-center delocalization around one ring path)}, \qquad
\text{MCI} = \text{mean of }I_{\text{ring}}\text{ over all ring permutations}.
$$

$I_{\text{ring}}$ is the conceptually central object here — a genuine $n$-center generalization of the
2-center delocalization index (the 2-center case of $I_{\text{ring}}$ reduces exactly to $\delta(A,B)$
itself), measuring how much electron density is delocalized *around the entire ring path at once*, not
just pairwise between adjacent atoms. Because $I_{\text{ring}}$'s value depends on which direction you
traverse the ring and where you start, **MCI** (multicenter index) averages it over every distinct
permutation of the ring atoms to give a single, permutation-invariant number — the price of that invariance
being a combinatorially larger number of ring-path evaluations for bigger rings. PDI and FLU are cheaper,
more targeted quantities built directly from pairwise $\delta$ values rather than the full multicenter
trace: PDI specifically exploits the *para* relationship unique to 6-membered rings (opposite atoms across
the ring), while FLU measures how uniformly delocalization is spread around the ring relative to a
reference aromatic system at the same level of theory (FLU needs that external reference $\delta^{\text{ref}}$
supplied by the caller, since "how much delocalization is normal" is inherently a comparison, not an
absolute scale).

Verified example — benzene, HOMA/Bird (geometric) alongside the electron-delocalization indices:

```python
benzene = """
C  0.0000  1.3970  0.0000
C  1.2098  0.6985  0.0000
C  1.2098 -0.6985  0.0000
C  0.0000 -1.3970  0.0000
C -1.2098 -0.6985  0.0000
C -1.2098  0.6985  0.0000
H  0.0000  2.4810  0.0000
H  2.1486  1.2405  0.0000
H  2.1486 -1.2405  0.0000
H  0.0000 -2.4810  0.0000
H -2.1486 -1.2405  0.0000
H -2.1486  1.2405  0.0000
"""
mb = qc.chk.new(atom=benzene, ao="sto-3g", unit="angstrom").scf(ref="r").run()
ring = [0, 1, 2, 3, 4, 5]
mb.prop.arom.homa(ring=ring)      # 0.9792  -- near 1: bonds close to the ideal aromatic length
mb.prop.arom.bird(ring=ring)      # 99.987  -- near 100: bond orders essentially uniform
mb.prop.arom.indices(ring=ring)
# {'pdi': 0.1029, 'flu': 0.000259, 'iring': 0.0434, 'mci': 0.0662, ...}
```

Every index independently confirms benzene's textbook aromaticity — HOMA and Bird near their respective
ideal values (1 and 100), PDI near the reference "aromatic" value of $\approx0.1$ the delocalization-index
literature reports for benzene specifically, and FLU essentially zero (uniform delocalization, no
significant bond-alternation-driven deviation). The geometric (HOMA/Bird) and electron-delocalization
(PDI/FLU/$I_{\text{ring}}$/MCI) families are answering related but genuinely distinct questions — one about
bond-length/bond-order uniformity, the other about the wavefunction's own multicenter electron sharing —
and their simultaneous agreement here is a strong, multi-method confirmation, not a redundant repetition of
the same calculation.

:::{exercise}
:label: ex-reactivity-theory

1. Water's electrophilicity index $\omega\approx0.52\ \text{eV}$ is small, and its hardness
   $\eta\approx18.3\ \text{eV}$ is large. Explain in one sentence why a *large* $\eta$ in the denominator
   of $\omega=\mu^2/(2\eta)$ is exactly what you would expect to produce a *small* electrophilicity for a
   hard, unreactive closed-shell molecule, independent of what $\mu$ happens to be.
2. The FMO approximation to the Fukui function needs no extra $N\pm1$ calculations at all, unlike the more
   accurate finite-difference variant. What is the one specific physical assumption FMO makes that a true
   $N\pm1$ calculation would not need to make, and under what circumstance (think: near-degenerate frontier
   orbitals, or a system with strong orbital relaxation upon ionization) would you expect FMO to become
   noticeably less reliable?
3. PDI is defined only for 6-membered rings (via the *para* relationship), while $I_{\text{ring}}$/MCI have
   no such restriction. Why does the *para* concept specifically require a 6-membered ring, and what does
   this tell you about why $I_{\text{ring}}$/MCI are the more generally applicable members of this index
   family?
:::

:::{solution} ex-reactivity-theory
:class: dropdown

1. $\omega=\mu^2/(2\eta)$ has $\eta$ in the denominator, so for any fixed numerator $\mu^2$, a larger
   $\eta$ mechanically produces a smaller $\omega$ — and $\eta=I-A$ being large means the HOMO-LUMO-like
   gap driving both ionization and electron attachment is wide, i.e. exactly the textbook definition of a
   "hard," resistant-to-electron-transfer species. A hard molecule that strongly resists both gaining and
   losing electron density should, by the electrophilicity index's own physical motivation (how much
   energy is gained by accepting the maximum charge flow from a reservoir), stabilize very little upon
   that hypothetical charge transfer — which is exactly what a large denominator enforces regardless of
   the specific value of $\mu$.
2. FMO assumes the frontier orbitals themselves do not change shape upon adding or removing an electron —
   it uses the *neutral* molecule's HOMO/LUMO densities as stand-ins for where the actual $N\pm1$ density
   difference would appear, rather than allowing the orbitals to relax in response to the changed electron
   count. This assumption becomes noticeably worse whenever ionization/attachment causes substantial
   **orbital relaxation** (a system where the cation's or anion's true orbitals differ significantly in
   shape from the neutral's frozen ones) or whenever the **frontier orbitals are near-degenerate** (so
   which orbital is "the" HOMO or LUMO, and hence which single orbital FMO picks to approximate the Fukui
   function, becomes an arbitrary or unstable choice rather than a clearly dominant single orbital).
3. "Para" specifically means "directly across the ring" — a well-defined, unique relationship that only
   exists for an *even*-membered ring where every atom has exactly one atom diametrically opposite it; a
   6-membered ring is the smallest (and chemically most common) ring where this geometric relationship is
   unambiguous. $I_{\text{ring}}$ and MCI, by contrast, are built from the full $n$-center trace around the
   *entire* ring path in sequence — a construction that needs no special "opposite atom" relationship at
   all and generalizes immediately to rings of any size (5-membered, 7-membered, or larger), which is
   exactly why they are the more general-purpose members of this aromaticity-index family, at the cost of
   the extra permutation-averaging MCI needs for permutation invariance.
:::

This chapter completes Part II's tour of what a converged density and wavefunction can tell you: energies
and structure ([Parts I-II's foundational chapters](hartree-fock.md)), forces and vibrations ([analytic
derivatives](analytic-derivatives.md), [the Hessian](analytic-hessian-thermo.md)), environment
([ECP](ecp-theory.md), [solvation](solvation-theory.md), [dispersion](dispersion-theory.md)), and now
charge distribution, topology, and reactivity. [The multireference outlook](multireference-outlook.md)
closes Part II by looking at what qc-rs does *not yet* implement — CASSCF-based correlation methods — and
why.