Atomic charges & bond orders

Atomic charges & bond orders#

“How much negative charge is on the oxygen?” and “is this a single or double bond?” are two of the most common questions in chemistry — and two of the most subtle, because neither an atomic charge nor a bond order is a quantum-mechanical observable. They are partitioning schemes: recipes for dividing the molecule’s density among atoms and pairs. qc-rs implements many, and this chapter shows how to compute them and how to read the differences.

Theory: why there are so many charge schemes#

The electron density \(\rho(\mathbf r)\) is real and unambiguous, but “the charge on atom A” requires you to decide where atom A ends and atom B begins — and there is no unique answer. Different schemes make different choices:

Orbital-based (Mulliken, Löwdin): split the density by which atom’s basis functions carry it. Simple and fast, but basis-set dependent (Mulliken especially).
Density-based (Hirshfeld, MBIS, VDD): partition real space using atomic reference densities. More robust and basis-stable.
Electrostatic-potential fit (MK/ChelpG/RESP): choose charges that best reproduce the molecule’s ESP — the right choice for force-field parameterization.
Natural population (NPA): from the natural atomic orbitals; basis-stable and widely used.

No scheme is “correct”; each answers a slightly different question. The lesson is to pick one appropriate to your purpose and be consistent, and never over-interpret the last digit.

Atomic charges#

Every scheme is a leaf of qc.prop.chrg. The two orbital-based schemes have the simplest formulas — both reduce a per-AO gross population \(P_\mu\) to a per-atom charge \(q_A = Z_A - \sum_{\mu\in A} P_\mu\), and differ only in which population they use:

\[ \text{Mulliken: } P_\mu = (\mathbf{DS})_{\mu\mu} = \sum_\nu D_{\mu\nu}S_{\mu\nu}, \qquad \text{Löwdin: } P_\mu = \big(\mathbf S^{1/2}\mathbf D\,\mathbf S^{1/2}\big)_{\mu\mu}. \]

Mulliken’s \(\mathbf{DS}\) population is not symmetric in \(\mu,\nu\) and can even be negative for a diffuse, poorly-overlapping AO — the basis-set sensitivity the table below shows. Löwdin’s symmetric orthogonalization (\(\mathbf S^{1/2}\) on both sides) fixes that asymmetry, which is part of why it tends to be gentler than Mulliken. Density-based schemes (Hirshfeld, MBIS) instead partition the real-space density \(\rho(\mathbf r)\) directly using reference atomic densities, needing no AO population at all — a fundamentally different (and more basis-stable) recipe, detailed in Population analysis theory. Some leaves return a plain per-atom array, some a richer record with a "charges" field:

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

np.asarray(qc.prop.chrg.hirshfeld(m))          # array: [-0.3219, 0.161, 0.161]
qc.prop.chrg.mulliken(m)["charges"]            # record: [-0.3056, 0.1528, 0.1528]
qc.prop.chrg.npa(m)["charges"]                 # [-0.9136, 0.4568, 0.4568]

The oxygen of water is negative and the hydrogens positive in every scheme — the qualitative picture is robust — but the magnitude varies widely with the recipe:

scheme	`qc.prop.chrg.`	O charge	H charge	character
Löwdin	`lowdin`	−0.094	+0.047	orbital, symmetrized
Mulliken	`mulliken`	−0.306	+0.153	orbital, basis-sensitive
Hirshfeld	`hirshfeld`	−0.322	+0.161	density, “stockholder”
CM5	`cm5`	−0.658	+0.329	Hirshfeld + empirical correction
ADCH	`adch`	−0.731	+0.365	atomic-dipole-corrected Hirshfeld
MBIS	`mbis`	−0.862	+0.431	minimal-basis iterative stockholder
NPA	`npa`	−0.914	+0.457	natural population

The spread — from −0.09 to −0.91 on the same oxygen — is exactly why you must state which scheme you used. For comparing across a series of molecules, a density-based scheme (Hirshfeld/MBIS) or NPA is usually the safer choice than Mulliken.

Bond orders#

A bond order quantifies how many electron pairs are shared between two atoms. The two workhorses are in qc.prop.bond, each returning a full atom–atom "matrix" plus a per-atom "valence":

mayer = qc.prop.bond.mayer(m)
mayer["matrix"][0, 1]     # 1.0207   the O–H Mayer bond order
mayer["valence"][0]       # 2.041    oxygen's total valence (≈ 2 bonds)

qc.prop.bond.wiberg(m)["matrix"][0, 1]   # 1.1886   the O–H Wiberg bond order

The O–H Mayer bond order is ~1.02 — a single bond, as chemistry expects, and oxygen’s valence ~2.04 correctly reflects its two O–H bonds. Both bond orders come from squaring the same density-times-overlap matrix product that drove Mulliken charges above, just summed differently:

\[ B_{AB}^{\text{Mayer}} = 2\!\!\sum_{\mu\in A,\,\nu\in B}\!\!\Big[(\mathbf P^\alpha\mathbf S)_{\mu\nu}(\mathbf P^\alpha\mathbf S)_{\nu\mu} + (\mathbf P^\beta\mathbf S)_{\mu\nu}(\mathbf P^\beta\mathbf S)_{\nu\mu}\Big], \qquad B_{AB}^{\text{Wiberg}} = 2\!\!\sum_{\mu\in A,\,\nu\in B}\!\!\Big[(\tilde{\mathbf P}^\alpha)_{\mu\nu}^2 + (\tilde{\mathbf P}^\beta)_{\mu\nu}^2\Big], \]

with \(\mathbf P^\alpha,\mathbf P^\beta\) the spin density matrices and \(\tilde{\mathbf P}=\mathbf S^{1/2}\mathbf P\,\mathbf S^{1/2}\) the Löwdin-orthogonalized density. Wiberg is exactly Mayer evaluated in the orthogonalized basis — where the overlap is the identity, so the formula simplifies to a sum of squared density-matrix elements — which is why it gives a slightly different number (1.19 here vs. Mayer’s 1.02): the two are the same idea (how much density is shared between \(A\) and \(B\)) measured in two different bases. Each atom’s total valence, \(\sum_{B\ne A}B_{AB}\), is the row sum excluding the diagonal. qc.prop.bond also holds fuzzy-atom, multicenter (for delocalized/aromatic bonding), IBSI, and delocalization-index bond measures.

Worked example: charges and bonds together#

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

print("Hirshfeld charges:", np.round(np.asarray(m.prop.chrg.hirshfeld()), 3))   # [-0.322 0.161 0.161]
print("NPA charges      :", np.round(np.asarray(m.prop.chrg.npa()["charges"]), 3))
print("O–H Mayer BO     :", round(m.prop.bond.mayer()["matrix"][0, 1], 3))       # 1.021

Exercise 11

Two papers report the oxygen charge in water as −0.31 and −0.91. Can both be right? What single piece of information reconciles them?
You are fitting a classical force field and need point charges. Which family of charge scheme should you use, and why not Mulliken?
Oxygen’s Mayer valence in water comes out ≈ 2.04. What does that number tell you, and why is it not exactly 2?

Solution to Exercise 11

Yes — atomic charge is not an observable. −0.31 is Mulliken-like, −0.91 is NPA-like; stating the partitioning scheme reconciles them. Comparisons are only meaningful within one scheme.
An electrostatic-potential-fit scheme (MK / ChelpG / RESP, in qc.prop.chrg) — the charges are chosen to reproduce the ESP the force field must model. Mulliken charges are not fit to the ESP and are basis-sensitive, so they transfer poorly.
It says oxygen participates in ~2 bonds’ worth of shared electron pairs (its two O–H bonds) — consistent with its divalency. It is not exactly 2 because real bonds are slightly polarized/ionic, which Mayer’s covalent bond order does not count as full sharing.

Charges and bond orders partition the density numerically. The next chapter partitions it topologically — QTAIM and ELF find the atoms, bonds, and lone pairs from the shape of the density itself.