NBO, IAO/IBO & orbital localization

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:

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.

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.

Algorithm 3 (The NBO Lewis-structure search)

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

Core. Take the core NAOs directly as core (CR) NBOs; deplete their occupancy from \(\mathbf P\).
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.
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.
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.

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#

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 17

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?
An IBO localization returns exactly 5 orbitals for water (10 electrons). Why 5, and what physical objects would you expect them to be?
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 to Exercise 17

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.
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.
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 turn from orbital-space analysis to a real-space electrostatic property you can map directly onto a molecular surface.