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. Some 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. Wiberg (computed in an orthogonalized/NAO basis) gives a slightly different 1.19. 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 10

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 10

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.