# NBO, IAO/IBO & orbital localization

The canonical molecular orbitals an SCF returns are delocalized — the HOMO of water is not "the O–H bond," it
is some spread-out combination of AOs across the whole molecule. **Localization** methods find an equivalent,
*chemically interpretable* set of orbitals — Lewis-structure bonds, lone pairs, core orbitals — that
represent the *same* density and energy. This chapter covers qc-rs's two localization families: **NBO**
(natural bond orbitals, the classic Weinhold recipe) and **IAO/IBO** (intrinsic atomic/bond orbitals, a
newer, simpler construction), plus general-purpose orbital localization.

## Theory: minimal atom-like orbitals — IAO

The starting point for both IAO and IBO is building a small, **orthonormal** basis that spans exactly the
occupied space while looking as much as possible like **free-atom** orbitals — the **intrinsic atomic
orbitals (IAOs)**. The construction (Knizia, 2013) projects the occupied MOs $\mathbf C$ onto a **minimal**
reference basis (`minao`) sitting on the same atoms, depolarizes that projection by Löwdin-orthogonalizing
it, and combines the result with $\mathbf C$ itself so the final IAOs are exactly as many as the minimal
basis has functions, yet span the occupied space exactly. Practically: a merged overlap between the target
and minimal bases is built once, the depolarized projector is formed, and the whole construction is driven
entirely by the **density** (via its natural orbitals) rather than needing molecular-orbital bookkeeping.

**IAO charges** — a Mulliken-style population taken in this minimal, atom-like basis — give **NPA-quality**
atomic charges from a single linear-algebra recipe, without NPA's own iterative machinery:

```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()

qc.prop.nbo.iao(m)["charges"]   # [-0.696, 0.348, 0.348]  -- matches pyscf.lo.iao to ~1e-6
```

## Theory: localizing into bonds — IBO

**Intrinsic bond orbitals (IBOs)** take the occupied canonical MOs and **rotate** them, within the occupied
space, to maximize how sharply each localized orbital sits on as few atoms as possible — a **Pipek–Mezey**
localization, but carried out in the IAO basis (so "how much charge sits on atom $A$" is well-defined even
though IAOs, unlike full AOs, span only the occupied space):

$$
L = \sum_i \sum_A \big(Q^i_A\big)^p ,
$$

maximized over $2\times2$ orbital rotations (a Jacobi sweep), where $Q^i_A$ is orbital $i$'s IAO-basis
Mulliken population on atom $A$ and $p$ (default 4) controls how aggressively the functional rewards
concentration on one or two atoms. The result is a set of orbitals you can *look at* and recognize: a core
orbital sitting on one atom, a bonding orbital split between two, a lone pair almost entirely on one.

```python
ibo = qc.prop.nbo.ibo(m)
len(ibo["orbitals"])                 # 5 -- an O core + 2 lone pairs + 2 O-H bonds, for water
ibo["orbitals"][0]["centers"]        # which atoms this orbital sits on (>10% weight)
```

`qc.prop.orb.localize` provides the same Pipek–Mezey machinery as a general-purpose tool (not tied to the IAO
basis), useful whenever you just want *a* set of localized orbitals rather than the specific IAO/IBO
construction.

## Theory: NBO — a Lewis structure from the density

**Natural bond orbitals (NBO)** take a different, older (Weinhold) route: rather than a smooth optimization,
they **search** for the best Lewis structure directly from the natural-orbital occupancy matrix
$\mathbf P = \mathbf C^\top(\mathbf{SDS})\mathbf C$ in the natural-atomic-orbital (NAO) basis.

:::{prf:algorithm} The NBO Lewis-structure search
:label: alg-nbo

**Input:** the NAO-basis occupancy matrix $\mathbf P$; the target pair count $n_{\text{pairs}}=\text{round}(N/2)$.
**Output:** a depleted, orthonormal set of Lewis-structure orbitals (core, lone pair, bond).

1. **Core.** Take the core NAOs directly as core (CR) NBOs; deplete their occupancy from $\mathbf P$.
2. **Lewis search.** Walk a descending occupancy threshold from ${\approx}2.0$ down to a floor. At each
   level, accept the highest-occupancy candidate eigenvector of $\mathbf P$ — but *prefer* a one-center lone
   pair over a two-center bond of similar occupancy (this keeps symmetric lone-pair combinations from being
   misread as bonds, and lets a strongly delocalized lone pair, e.g. an amide nitrogen's, still be found
   rather than absorbed into a spurious bond). A candidate counts as a genuine two-center **bond** only if it
   has significant weight on *both* centers; otherwise it is a lone-pair tail. Deplete the accepted orbital,
   $\mathbf P \leftarrow \mathbf P - n\,|v\rangle\langle v|$, and stop at $n_{\text{pairs}}$ orbitals.
3. **Hybrids.** Each NBO's amplitude on a centre, resolved by angular momentum, gives its natural hybrid
   composition ($\%s/\%p/\%d$) and — for a bond — its polarization between the two atoms.
4. **Completion.** Every bond gets its orthogonal complement, the antibond (BD*); the remaining directions
   complete an orthonormal transform from the NAO basis to the full NBO basis (bonds, lone pairs, and the
   non-Lewis Rydberg orbitals).
:::

The payoff of building a full Lewis structure this way is the **second-order donor–acceptor analysis**:
transform the Fock matrix into the NBO basis, $F = T^\top(\mathbf C^\top F_{\text{AO}}\mathbf C)T$, and every
filled orbital $i$ interacting with an empty one $j$ (typically a lone pair donating into an antibond)
stabilizes the energy by

$$
E^{(2)}_{i\to j} = -2\,\frac{F_{ij}^2}{\varepsilon_j-\varepsilon_i} .
$$

This is exactly the textbook second-order perturbation-theory stabilization energy, evaluated in the basis
where "donor" and "acceptor" orbitals are already chemically labelled — which is why NBO's $E^{(2)}$ is the
standard way to quantify hyperconjugation, anomeric effects, and other donor–acceptor interactions.

```python
nbo = qc.prop.nbo.nbo(m)
nbo["e2"][0]   # {'donor': ..., 'acceptor': ..., 'donor_desc': 'LP O1', 'acceptor_desc': 'RY H2', 'energy_kcal': 2.07}
```

## Worked example

```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="def2-svp", unit="angstrom").scf(ref="r").run()

print("IAO charges :", np.round(qc.prop.nbo.iao(m)["charges"], 3))     # [-0.696  0.348  0.348]

ibo = qc.prop.nbo.ibo(m)
print("IBO orbital count:", len(ibo["orbitals"]))                      # 5

nbo = qc.prop.nbo.nbo(m)
top = nbo["e2"][0]
print(f"largest E(2): {top['donor_desc']} -> {top['acceptor_desc']}, {top['energy_kcal']:.2f} kcal/mol")
```

:::{exercise}
:label: ex-nbo

1. IAO charges and NPA charges are both meant to be "NPA-quality," yet they are built completely
   differently. What is the one thing both methods have in common that Mulliken lacks?
2. An IBO localization returns exactly 5 orbitals for water (10 electrons). Why 5, and what physical
   objects would you expect them to be?
3. You find a large NBO $E^{(2)}$ from a lone pair into a neighbouring antibond. Name one real chemical
   phenomenon this quantifies, and explain using the $E^{(2)}$ formula why a *small* orbital-energy gap
   $\varepsilon_j-\varepsilon_i$ makes the effect stronger.
:::

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

1. Both are evaluated in a **minimal, atom-centred reference basis** (IAO's depolarized projection onto
   `minao`; NPA's natural atomic orbitals) rather than the raw calculation basis — which is exactly why both
   are far less basis-set-sensitive than plain Mulliken.
2. **5** because water has 5 occupied orbitals (10 electrons ÷ 2 per orbital): the O 1s core, two O lone
   pairs, and the two O–H bonds — precisely the textbook Lewis structure of water.
3. This is the signature of **hyperconjugation** (or, for a heteroatom lone pair into an antibonding orbital
   across a ring, the **anomeric effect**). From $E^{(2)}=-2F_{ij}^2/(\varepsilon_j-\varepsilon_i)$, a
   smaller energy denominator directly *inflates* $E^{(2)}$ for the same coupling $F_{ij}$ — donor and
   acceptor orbitals closer in energy interact more strongly, exactly as ordinary perturbation theory
   predicts.
:::

Next, [ESP surfaces](esp-surfaces.md) turn from orbital-space analysis to a real-space electrostatic
property you can map directly onto a molecular surface.