# 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](../../10-foundations/population-analysis-theory.md). Some leaves return a plain
per-atom array, some a richer record with a `"charges"` field:

```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="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"`:

```python
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

```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="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}
:label: ex-charges

1. 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?
2. You are fitting a classical force field and need point charges. Which *family* of charge scheme should you
   use, and why not Mulliken?
3. Oxygen's Mayer valence in water comes out ≈ 2.04. What does that number tell you, and why is it not exactly
   2?
:::

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

1. 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.
2. 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.
3. 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](qtaim-elf.md) partitions it
*topologically* — QTAIM and ELF find the atoms, bonds, and lone pairs from the shape of the density itself.