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.
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\) 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:
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):
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:
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:
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. (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:
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 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) and the atomic overlap matrices \(S(A)\) in the occupied-MO basis:
\(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:
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 18
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.
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?
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 to Exercise 18
\(\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\).
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).
“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), forces and vibrations (analytic derivatives, the Hessian), environment (ECP, solvation, dispersion), and now charge distribution, topology, and reactivity. The multireference outlook closes Part II by looking at what qc-rs does not yet implement — CASSCF-based correlation methods — and why.