# Population analysis theory A converged wavefunction has no atoms in it — only a molecular density and a set of delocalized molecular orbitals. Every "atomic charge" a chemist reports is therefore the result of a **partitioning scheme**, a choice of how to divide an inherently molecular quantity back among nuclei, and different schemes make genuinely different choices, not just different numerical approximations to one "true" atomic charge (there is no such single truth once you leave a strictly atomic system). This chapter derives every population-analysis method qc-rs's `qc.prop.chrg`/`qc.prop.bond` namespace implements — Mulliken, Löwdin, the Hirshfeld stockholder family (plain, iterative, CM5, ADCH), MBIS, NPA/NAO, VDD, and the Mayer/Wiberg bond indices — grounded directly in the `qc-prop` crate's source (each a clean-room implementation from its cited paper, cross-validated against reference codes such as JANPA). ## AO-partition charges: Mulliken and Löwdin The simplest partitioning schemes work entirely in AO-index space, never touching real space at all. Both start from the same object, the AO-basis population matrix, but differ in *which* basis the density is expressed in before summing. **Mulliken** sums the density matrix contracted with the overlap, attributing each *off-diagonal* AO-pair population entirely to the atoms of the two AOs involved (splitting it evenly when they differ): $$ P_\mu = (DS)_{\mu\mu}, \qquad N_A = \sum_{\mu\in A}P_\mu, \qquad q_A = Z_A - N_A. $$ Mulliken's well-known weakness is basis-set sensitivity: adding more diffuse functions to *either* atom in a bonded pair shifts overlap population between them without changing the physical electron distribution at all, so $q_A$ can vary substantially, even in sign, between basis sets for the same molecule and geometry. **Löwdin** charges fix this by first symmetrically orthogonalizing the AO basis — transforming the density by $S^{1/2}$ on both sides before summing the diagonal, so there is no off-diagonal overlap contribution left to arbitrarily split: $$ P_\mu = \bigl(S^{1/2}DS^{1/2}\bigr)_{\mu\mu}, \qquad N_A = \sum_{\mu\in A}P_\mu, \qquad q_A = Z_A - N_A. $$ Löwdin charges are markedly more basis-set stable than Mulliken's for exactly this reason — there is simply no overlap-splitting decision left to make — but both schemes remain fundamentally AO-index partitions: they say nothing about *where in space* the electron density actually sits, only how to divide up a density matrix expressed in one basis or another. ## The Hirshfeld stockholder family **Hirshfeld** partitioning instead works in real space, using a **promolecule** — the superposition of spherically-averaged, isolated, neutral free-atom densities placed at the actual molecular geometry — as a reference for how "generously" each atom should be credited with the molecule's real density at every point: $$ w_A(\mathbf r) = \frac{\rho_A^0(\mathbf r)}{\sum_B\rho_B^0(\mathbf r)}, \qquad N_A = \int w_A(\mathbf r)\,\rho(\mathbf r)\,d\mathbf r, \qquad q_A = Z_A - N_A, $$ with $\rho_A^0$ the free neutral atom's density centered on $A$. Because $\sum_Aw_A(\mathbf r)=1$ everywhere by construction, this is a genuine partition of unity — every point in space contributes its full density weight to *some* combination of atoms — but the reference promolecule is always built from *neutral* free atoms regardless of the molecule's actual bonding, which is the well-documented reason plain Hirshfeld charges tend to be smaller in magnitude than chemical intuition or dipole-fitted charges would suggest: a genuinely charged atom in the real molecule is still compared against a neutral reference. **Hirshfeld-I** removes exactly this arbitrariness by making the promolecule **self-consistent**: instead of always comparing against the neutral atom, it interpolates the isolated-atom reference density between the two bracketing *integer* electron counts around the atom's own current fractional population $N_A$, $$ \rho_A^0[N_A](\mathbf r) = \bigl(\lceil N_A\rceil-N_A\bigr)\,\rho_A[\lfloor N_A\rfloor](\mathbf r) + \bigl(N_A-\lfloor N_A\rfloor\bigr)\,\rho_A[\lceil N_A\rceil](\mathbf r), $$ and iterates $N_A=\int w_A\rho\,d\mathbf r$ to self-consistency — each cycle's promolecule is built from the *previous* cycle's fractional populations, until they stop changing. This directly removes plain Hirshfeld's "always neutral reference" bias and gives systematically larger-magnitude, better-converged charges (Bultinck et al., 2007) at the cost of needing a small library of charged-atom reference densities (the bracketing integer electron counts) rather than only neutral-atom ones. **CM5** takes the *opposite* strategy from Hirshfeld-I: rather than fixing the underlying partitioning scheme's self-consistency, it applies a purely empirical, geometry-based **post-correction** directly to plain Hirshfeld charges, fitted to reproduce experimental/high-level-theory molecular dipole moments: $$ q_A^{\text{CM5}} = q_A^{\text{HPA}} + \sum_{B\ne A}T_{AB}\,B_{AB}(r_{AB}), \qquad B_{AB} = \exp\bigl[-\alpha(r_{AB}-R_A-R_B)\bigr], $$ with $\alpha=2.474\ \text{Å}^{-1}$, $R_A,R_B$ tabulated covalent radii, and $T_{AB}$ an atom-pair-specific (or, for most element pairs, a simple atomwise-difference) correction coefficient fitted against a large reference set of Hirshfeld-to-"true"-dipole mappings. CM5's correction decays exponentially with distance (so only *bonded or near-bonded* pairs correct each other significantly) and needs no self-consistent solve at all — it is a one-shot, closed-form correction applied after a single plain-Hirshfeld calculation. **ADCH** solves a related but structurally different problem: rather than correcting the *charges* empirically, it corrects for the fact that Hirshfeld charges alone cannot reproduce the molecular dipole moment, because Hirshfeld discards each atom's own **atomic dipole** — the first moment of how the real density deviates from spherical symmetry around that atom, $\boldsymbol\mu_A=-\int w_A(\mathbf r)(\mathbf r-\mathbf R_A)\,\rho(\mathbf r)\,d\mathbf r$. ADCH expands this discarded atomic dipole into a set of point-charge corrections placed on the atom's neighbors (a Thole-Duijnen-style redistribution), solving for correction charges $q_{AB}$ that satisfy two exact constraints simultaneously — charge neutrality of the redistribution and *exact* reproduction of the atomic dipole as a first moment: $$ \sum_Bq_{AB}=0, \qquad \sum_Bq_{AB}\,\mathbf r_B=\boldsymbol\mu_A, \qquad Q_A = q_A^{\text{HPA}}+\sum_Bq_{BA}. $$ Because these are exact linear constraints (solved via a local weighted covariance matrix $\boldsymbol\Lambda=\sum_Bv_{AB}(\mathbf r_B-\langle\mathbf r\rangle)(\mathbf r_B-\langle\mathbf r\rangle)^{\mathsf T}$, not a fit against external reference data), **the molecular dipole computed from ADCH's final point charges exactly reproduces the QM dipole moment** by construction — a hard mathematical guarantee, not a calibrated approximation, and the natural correctness test the qc-rs implementation uses. Unlike CM5, ADCH needs no element-specific fitted parameters at all. ## MBIS: a different route to self-consistency **MBIS** (Minimal Basis Iterative Stockholder) is a second self-consistent stockholder scheme, but reaches self-consistency by a different mechanism than Hirshfeld-I: instead of interpolating between tabulated isolated-atom reference densities, it represents each atom's pro-atom density directly as a **minimal expansion of atom-centered Slater (exponential) shells**, $$ \rho_A^0(\mathbf r) = \sum_iN_{Ai}\,f_{Ai}(\mathbf r), \qquad f_{Ai}(\mathbf r) = \frac{1}{8\pi\sigma_{Ai}^3}\exp\Bigl(-\frac{|\mathbf r-\mathbf R_A|}{\sigma_{Ai}}\Bigr), $$ with each shell individually normalized to one electron ($\int f_{Ai}=1$). The stockholder weight is built from these shells exactly as in Hirshfeld, $w_A=\rho_A^0/\sum_B\rho_B^0$, but now *both* the shell populations $N_{Ai}$ and widths $\sigma_{Ai}$ are refined to self-consistency directly against the *molecule's own* density — no external charged-atom reference library at all: $$ N_{Ai} = \int\rho(\mathbf r)\,\frac{\rho_{Ai}^0(\mathbf r)}{\rho_A^0(\mathbf r)}\,d\mathbf r, \qquad \sigma_{Ai} = \frac{1}{3N_{Ai}}\int\rho(\mathbf r)\,\frac{\rho_{Ai}^0(\mathbf r)}{\rho_A^0(\mathbf r)}\,|\mathbf r-\mathbf R_A|\,d\mathbf r, \qquad q_A = Z_A-\sum_iN_{Ai}. $$ This is the essential structural difference from Hirshfeld-I: Hirshfeld-I's self-consistency lives in *which pre-computed reference density* to compare against (an interpolation problem over a small, external table), while MBIS's self-consistency lives in *fitting the pro-atom's own shape* directly to the molecule at hand (an optimization problem with no external reference densities needed beyond the initial neutral-atom shell populations used to seed the iteration). ## NPA/NAO: basis-robust populations via natural orbitals Mulliken and Löwdin charges are AO-index partitions, and both remain sensitive — Löwdin less severely than Mulliken, but still measurably — to which particular AO basis was used, since neither construction distinguishes a chemically meaningful AO contribution from a numerically redundant, diffuse one that happens to overlap heavily with a neighboring atom's functions. **Natural Population Analysis (NPA)** addresses this at its root by first constructing an orthonormal, per-atom, symmetry-adapted basis — the **Natural Atomic Orbitals (NAOs)** — from the density itself, and only then reading off populations, rather than partitioning an arbitrary externally-chosen AO basis directly. The central object is the **occupation matrix** $\mathbf P=\mathbf S\mathbf D\mathbf S$ (not the bare density $\mathbf D$) — using the full $SDS$ sandwich rather than the density alone is precisely what makes the resulting populations basis-set robust, because it accounts for how the AO basis's own overlap structure inflates or deflates a naive density-matrix reading. Constructing the NAOs is a sequence of distinct algebraic stages (the JANPA re-derivation of Reed-Weinstock-Weinhold's original procedure), each solving a well-defined local problem: :::{prf:algorithm} Natural atomic orbital construction (sketch) :label: alg-nao **Input:** AO density $D$, AO overlap $S$, per-atom/angular-momentum AO labels. **Output:** the natural atomic orbitals and their occupancies. 1. **Intra-atomic naturalization**: within each $(atom,l)$ block, angular-average over $m$ and solve the generalized eigenproblem $\bar P\mathbf c=n\bar S\mathbf c$ — the eigenvectors are pre-NAOs, the eigenvalues $n$ their occupancy weights. 2. **Partition**: the highest-occupancy pre-NAOs per $(atom,l)$, up to the free-atom minimal-basis count, form the natural minimal basis (NMB); the remainder form the natural Rydberg basis (NRB). 3. **Occupancy-weighted symmetric orthogonalization** of the NMB — a weighted analogue of ordinary Löwdin orthogonalization that preserves the high-occupancy orbitals preferentially. 4. **Schmidt-orthogonalize** every NRB function against the (now-orthogonalized) NMB. 5. **Re-naturalize** the NRB set (intra-atomic naturalization again, now within the post-Schmidt NRB) to refine its weights. 6. **Occupancy-weighted orthogonalization** of the refined NRB set. 7. **Final intra-atomic naturalization** within the full NMB+NRB set, giving the final NAOs. ::: The NAO occupancies are the diagonal of $\mathbf C^{\mathsf T}(\mathbf S\mathbf D\mathbf S)\mathbf C$ in the final NAO basis $\mathbf C$; summing per atom gives $N_A$ and $q_A=Z_A-N_A$ exactly as in every other scheme here. What makes NPA basis-set robust is that steps 1–2 explicitly separate "genuinely occupied, chemically meaningful" orbitals (the NMB) from "formally present in the basis but essentially unoccupied, diffuse/polarization-function" orbitals (the NRB) *before* any orthogonalization mixes them — a bigger, more diffuse basis mostly adds NRB character with near-zero occupancy, rather than perturbing the NMB populations the way it would perturb a raw Mulliken/Löwdin partition. ## VDD: hard-partition real-space charges **Voronoi Deformation Density (VDD)** takes yet another real-space approach, avoiding both the AO-basis partition of Mulliken/Löwdin and the smooth stockholder weight of the Hirshfeld family. Instead, space is divided into hard **Voronoi cells** — every point belongs entirely to whichever nucleus is nearest, with no smooth weight function at all — and the charge is the integral of the **deformation density** (how much the real molecular density differs from the same neutral-atom promolecule Hirshfeld uses) over each cell: $$ q_A = -\int_{\text{Voronoi}(A)}\bigl[\rho(\mathbf r)-\rho^0(\mathbf r)\bigr]\,d\mathbf r, $$ with $\rho^0=\sum_B\rho_B^0$ the same promolecule superposition Hirshfeld builds. The physical interpretation is direct: VDD asks not "how much of the total density belongs to atom $A$" (Hirshfeld's question) but "how much *more or less* density is in atom $A$'s own region of space, relative to the neutral-promolecule baseline in that same hard-partitioned region" — a genuinely different question, using a genuinely different (hard rather than smooth) spatial partition, not merely a numerical variant of Hirshfeld. ## Mayer and Wiberg bond orders Charges answer "how much" density belongs to an atom; **bond orders** answer a related but distinct question — how much *shared, covalent character* exists between a pair of atoms. The **Mayer bond order** is built directly from the AO density and overlap, generalizing the textbook bond-order idea to a non-orthogonal AO basis: $$ B_{AB} = 2\sum_{\mu\in A,\nu\in B}\Bigl[(P^\alpha S)_{\mu\nu}(P^\alpha S)_{\nu\mu} + (P^\beta S)_{\mu\nu}(P^\beta S)_{\nu\mu}\Bigr], $$ reducing for a closed shell ($P^\alpha=P^\beta=D/2$) to $B_{AB}=\sum_{\mu\in A,\nu\in B}(DS)_{\mu\nu}(DS)_{\nu\mu}$. Each atom's **valence** $V_A=\sum_{B\ne A}B_{AB}$ sums its bond orders to every other atom (the diagonal $B_{AA}$, an intra-atomic term, is excluded from this sum, though it is retained in the full matrix for reference). The **Wiberg bond index** is exactly the Mayer construction evaluated in the Löwdin-symmetrically- orthogonalized basis — transform the density with $S^{1/2}$ on both sides first (so the overlap that would appear in the Mayer formula becomes the identity and drops out): $$ \tilde P^\sigma = S^{1/2}P^\sigma S^{1/2}, \qquad W_{AB} = 2\sum_{\mu\in A,\nu\in B} \Bigl[(\tilde P^\alpha_{\mu\nu})^2+(\tilde P^\beta_{\mu\nu})^2\Bigr]. $$ This is exactly the same relationship Löwdin charges bear to Mulliken charges — a basis change to the symmetrically-orthogonalized AO set before the same underlying formula is applied — and for the same reason: it removes the basis-set-arbitrary overlap-splitting decision that the plain (Mayer/Mulliken) non-orthogonal-basis formula must otherwise make. ## Verified example: every scheme, one molecule ```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() m.prop.chrg.mulliken()["charges"] # [-0.3471, 0.1735, 0.1735] m.prop.chrg.lowdin()["charges"] # [-0.1555, 0.0778, 0.0778] m.prop.chrg.hirshfeld() # [-0.3277, 0.1639, 0.1639] m.prop.chrg.hirshfeld_i() # [-0.9181, 0.4590, 0.4590] m.prop.chrg.cm5() # [-0.6642, 0.3321, 0.3321] m.prop.chrg.adch() # [-0.7576, 0.3788, 0.3788] m.prop.chrg.mbis() # [-0.8702, 0.4351, 0.4351] m.prop.chrg.npa()["charges"] # [-0.9115, 0.4557, 0.4557] m.prop.chrg.vdd() # [-0.3023, 0.1518, 0.1504] m.prop.bond.mayer()["valence"] # [1.9948, 1.0040, 1.0040] m.prop.bond.wiberg()["valence"] # [2.3054, 1.2081, 1.2081] ``` Every scheme agrees water's oxygen is negatively charged and both hydrogens equally (and oppositely) positive, as chemically expected — but the *magnitude* varies enormously, from Löwdin's modest $-0.16$ to Hirshfeld-I/NPA/MBIS's much larger $\sim-0.87$ to $-0.92$. This spread is not numerical disagreement about one true number; it is exactly the point made at the start of this chapter — these are genuinely different, self-consistent definitions of "how much charge belongs to an atom," and comparing charges *across* methods is far less meaningful than comparing them consistently *within* one method across a series of related molecules. The Mayer and Wiberg valences both correctly identify two O–H single bonds (valence $\approx1$ at each hydrogen, $\approx2$ at oxygen), with Wiberg's Löwdin-basis construction giving a somewhat larger absolute value than Mayer's raw non-orthogonal one — consistent with the same general basis-dependence pattern Löwdin/Mulliken charges show. :::{exercise} :label: ex-population-theory 1. Mulliken charges are known to be more basis-set-sensitive than Löwdin charges. Using the AO-partition derivation above, explain in one sentence *what specific step* Löwdin's construction removes that Mulliken's does not, and why that step is where basis-set sensitivity enters. 2. Hirshfeld-I and MBIS both reach self-consistency, but by different mechanisms. Explain in one sentence what varies from iteration to iteration in each scheme (what is the "unknown" being converged). 3. ADCH guarantees exact reproduction of the molecular dipole moment by construction, while CM5 only approximately reproduces it (via a fitted empirical correction). Both start from the same plain Hirshfeld charges. What is the structural difference in how each corrects those charges that explains this guarantee-vs-approximation distinction? ::: :::{solution} ex-population-theory :class: dropdown 1. Mulliken's formula keeps the off-diagonal overlap population $(DS)_{\mu\nu}$ ($\mu\ne\nu$, different atoms) and splits it evenly between the two atoms involved, with no principled way to decide that split other than "half each" — adding a diffuse function to either atom changes this off-diagonal overlap substantially without changing the physical density at all. Löwdin's construction symmetrically orthogonalizes the basis *before* summing ($S^{1/2}DS^{1/2}$), which makes the transformed basis's overlap matrix the identity — there is no off-diagonal overlap population left over to split arbitrarily, which is exactly the step (and exactly the reason) Mulliken's sensitivity disappears. 2. Hirshfeld-I converges the atom's *fractional population* $N_A$, which selects *which pre-tabulated reference density* (interpolated between the bracketing integer-electron-count isolated atoms) to compare the molecular density against — the unknown is a mixing weight between externally supplied reference densities. MBIS instead converges the pro-atom's own *shape parameters directly* — the shell populations $N_{Ai}$ and widths $\sigma_{Ai}$ of a Slater-shell expansion fitted to the molecule's own density — there are no external reference densities to interpolate between at all; the unknown is the pro-atom's own functional form. 3. ADCH solves an explicit system of *linear equality constraints* (zero net redistributed charge, and the redistributed charges' first moment exactly equal to the atomic dipole $\boldsymbol\mu_A$) — this is an exact algebraic condition that holds by construction for any molecule, with no fitting to external data at all. CM5 instead applies a *fixed, universal* set of empirical parameters ($\alpha$ and the pairwise/atomwise $T$ coefficients) fitted once, offline, across a large reference set of molecules to approximately reproduce dipole moments *on average* — for any single new molecule, the correction is whatever that pre-fitted formula gives, with no guarantee (only an empirical expectation, validated against the reference set) that it reproduces that particular molecule's own dipole exactly. ::: Charges and bond orders describe the density's distribution among atoms; [the topological-analysis chapter](topological-analysis-theory.md) instead studies the density's own critical-point structure directly (QTAIM, ELF/LOL) — a genuinely different, partition-free way of extracting chemical meaning from the same converged density.