Conceptual DFT & reactivity

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):

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:

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?#

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 14

  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.

That completes the property tour’s reactivity leg. The final chapter, spectra & DOS, reads the orbital-energy spectrum — the density of states and the HOMO–LUMO gap.