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):
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:
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:
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\),
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:
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:
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,
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:
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:
Algorithm 7 (Natural atomic orbital construction (sketch))
Input: AO density \(D\), AO overlap \(S\), per-atom/angular-momentum AO labels. Output: the natural atomic orbitals and their occupancies.
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.
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).
Occupancy-weighted symmetric orthogonalization of the NMB — a weighted analogue of ordinary Löwdin orthogonalization that preserves the high-occupancy orbitals preferentially.
Schmidt-orthogonalize every NRB function against the (now-orthogonalized) NMB.
Re-naturalize the NRB set (intra-atomic naturalization again, now within the post-Schmidt NRB) to refine its weights.
Occupancy-weighted orthogonalization of the refined NRB set.
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:
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:
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):
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#
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 16
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.
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).
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 to Exercise 16
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.
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.
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 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.