# Conceptual DFT & reactivity

Where will a molecule react, and how readily? **Conceptual DFT** answers this by treating familiar chemical
ideas — electronegativity, hardness, electrophilicity — as **derivatives of the energy** with respect to
electron number and external potential. It turns "reactivity" from hand-waving into computable numbers, both
**global** (one per molecule) and **local** (one per atom). They live in `qc.prop.cdft`.

## Theory: reactivity as energy derivatives

Conceptual DFT expands the energy in the electron number $N$ and the external potential $v(\mathbf r)$. The
low-order derivatives *are* the classic reactivity descriptors:

- **Chemical potential** $\mu = (\partial E/\partial N)_v$ — the escaping tendency of the electrons.
  **Electronegativity** is $\chi = -\mu$.
- **Chemical hardness** $\eta = (\partial^2 E/\partial N^2)_v$ — resistance to changing the electron count; a
  hard molecule has a large HOMO–LUMO gap. **Softness** is $1/\eta$.
- **Electrophilicity index** $\omega = \mu^2/2\eta$ — how strongly the molecule attracts electrons.

Using a finite-difference (Koopmans) approximation with the ionization potential $I$ and electron affinity
$A$: $\mu \approx -(I+A)/2$ and $\eta \approx I - A$. The **local** counterpart is the **Fukui function**
$f(\mathbf r) = (\partial\rho(\mathbf r)/\partial N)$ — where the density changes most when you add or remove
an electron, i.e. **where the molecule reacts**.

## Global descriptors

`qc.prop.cdft.reactivity` returns the whole global set (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="cc-pvdz", unit="angstrom").scf(ref="r").run()

r = qc.prop.cdft.reactivity(m)
r["ip"], r["ea"]                 # 13.419, -5.049    ionization potential / electron affinity
r["chemical_potential"]          # -4.185            μ
r["electronegativity"]           # 4.185             χ = -μ
r["hardness"]                    # 18.468            η  (large -> hard, big gap)
r["electrophilicity"]            # 0.474             ω
```

Water comes out as a **hard** molecule (η ≈ 18.5 eV — a large gap, chemically inert) with a modest
electrophilicity — exactly the expected profile for a small, tightly-bound closed-shell molecule.

## Local reactivity: Fukui functions

`qc.prop.cdft.fukui` gives the **condensed** (per-atom) Fukui functions — the local sites of reactivity:

```python
f = qc.prop.cdft.fukui(m)
f["fukui_plus"]        # [0.261, 0.370, 0.370]   f⁺ : susceptibility to NUCLEOPHILIC attack (adding e⁻)
f["fukui_minus"]       # [0.664, 0.168, 0.168]   f⁻ : susceptibility to ELECTROPHILIC attack (removing e⁻)
f["fukui_zero"]        # [0.462, 0.269, 0.269]   f⁰ : radical attack
f["dual_descriptor"]   # [-0.403, 0.202, 0.202]  f⁺ − f⁻ : >0 electrophilic site, <0 nucleophilic site
```

Reading water: the **oxygen has the largest f⁻ (0.66)**, so it is the site most susceptible to
**electrophilic** attack — chemically correct, since the oxygen lone pairs are where an electrophile (or a
proton) attacks. The **dual descriptor** is negative on oxygen (−0.40, a nucleophilic/electron-rich site) and
positive on the hydrogens (+0.20, electron-poor), a single field that classifies each atom in one number.

:::{note} Fukui functions cost extra SCFs
$f(\mathbf r)$ is a finite difference in electron number, so a condensed-Fukui evaluation needs the $N$,
$N{-}1$, and $N{+}1$ electron densities — i.e. a couple of extra SCF-like calculations behind the one call.
`qc.prop.cdft` also has `dual`, `local_reactivity` (softness/philicity), and `superdelocalizability`.
:::

## Worked example: where does water react?

```python
import qc, numpy as np
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()

r = qc.prop.cdft.reactivity(m)
print(f"μ = {r['chemical_potential']:.2f} eV,  η = {r['hardness']:.2f} eV,  ω = {r['electrophilicity']:.3f}")
#  μ = -4.18 eV,  η = 18.47 eV,  ω = 0.474

f = qc.prop.cdft.fukui(m)
print("f⁻ (electrophilic-attack sites):", np.round(f["fukui_minus"], 3))   # O highest -> [0.664 0.168 0.168]
```

:::{exercise}
:label: ex-cdft

1. Molecule A has hardness η = 2 eV, molecule B has η = 10 eV. Which is more reactive toward a change in
   electron count, and which has the larger HOMO–LUMO gap?
2. On a carbonyl compound, which condensed Fukui function would you inspect to predict where a **nucleophile**
   (e.g. a hydride) attacks, and would you look for a high or low value?
3. The dual descriptor is +0.3 on one carbon and −0.4 on an adjacent oxygen. Interpret each site.
:::

:::{solution} ex-cdft
:class: dropdown

1. **Molecule A** (η = 2 eV) is **more reactive** — low hardness means it easily gains/loses electron density.
   **Molecule B** (η = 10 eV) is harder and has the **larger HOMO–LUMO gap** (hardness ≈ gap).
2. **f⁺** (`fukui_plus`) — it measures susceptibility to nucleophilic attack (the density response to *adding*
   an electron). The nucleophile attacks the atom with the **highest** f⁺ (typically the carbonyl carbon).
3. The carbon (dual +0.3) is an **electrophilic** site (electron-poor, attacked by nucleophiles); the oxygen
   (dual −0.4) is a **nucleophilic** site (electron-rich, attacked by electrophiles/protons).
:::

That completes the property tour's reactivity leg. The final chapter, [spectra & DOS](spectra-dos.md), reads
the orbital-energy spectrum — the density of states and the HOMO–LUMO gap.