# Dispersion correction theory: DFT-D3/D4

Neither Hartree-Fock nor most density functionals capture long-range van der Waals (London dispersion)
attraction — HF has no correlation at all, and semilocal/hybrid DFT functionals are built from the local
density and its gradient, which carry no information about the instantaneous-multipole correlation between
electrons on *distant* fragments. Grimme's **DFT-D3** and **DFT-D4** dispersion corrections add this back
as a geometry-dependent (D3) or geometry-and-charge-dependent (D4) empirical energy term, computed entirely
from the nuclear positions (and, for D4, self-consistently determined atomic charges) — no additional
integrals, no participation in the SCF loop at all. This chapter derives both models in full, grounded
directly in `crates/qc-disp/src/energy.rs`, `ncoord.rs`, `c6.rs`, `eeq.rs`, and `d4.rs` — a faithful,
independently-implemented Rust port of the Grimme group's `s-dftd3`/`dftd4`/`multicharge`/`mctc-lib`.

## Why geometry alone can approximate a correlation effect

London dispersion between two well-separated fragments is, physically, a correlated fluctuation effect —
an instantaneous dipole on one fragment induces a correlated instantaneous dipole on the other, and their
attraction falls off as $1/R^6$ at leading order (higher multipole orders contributing $1/R^8$,
$1/R^{10}$, …). The key empirical fact both D3 and D4 exploit is that the *leading coefficient* $C_6$ of
this attraction is not really a fixed per-element constant — it depends strongly on an atom's local
chemical environment (an sp³ carbon and an aromatic carbon have measurably different effective
polarizabilities, hence different $C_6$), and that environment is well captured by something purely
geometric: how many neighbors an atom has and how close they are. This is the entire justification for
using a **coordination number** as the interpolation variable, rather than needing to run any actual
electronic-structure calculation to estimate a dispersion coefficient.

## The D3 coordination number

For atom $A$, the D3 coordination number counts effective neighbors with a smooth, saturating function
rather than a hard cutoff:

$$
\text{CN}_A = \sum_{B\ne A}\frac{1}{1+\exp\bigl[-k_{\text{cn}}(r_{AB}^{\text{cov}}/r_{AB}-1)\bigr]}, \qquad
r_{AB}^{\text{cov}} = r_A^{\text{cov}}+r_B^{\text{cov}},\qquad k_{\text{cn}}=16,
$$

with $r_A^{\text{cov}}$ tabulated D3 covalent radii. This counting function is $\approx1$ for a neighbor
well inside the covalent-radius sum (a genuine bonded neighbor), $\approx0$ for one well outside it (too
far to matter), and exactly $\tfrac12$ right at the covalent-radius threshold — a smooth sigmoid, not a
hard step, which is what makes $\text{CN}_A$ (and hence the interpolated $C_6$ built from it) a
differentiable function of geometry, needed for analytic dispersion gradients.

## CN-interpolated $C_6$: Gaussian similarity to tabulated references

D3 does not compute $C_6$ from a physical polarizability formula at runtime — it interpolates between a
table of **reference** coordination-number/$C_6$ pairs computed once, offline, for small reference
molecules of each element. For element $Z$ with tabulated reference coordination numbers
$\{\text{cn}^{\text{ref}}_a\}$, the weight given to reference $a$ at the *actual* coordination number
$\text{CN}$ is a Gaussian similarity kernel, normalized across all references:

$$
w^a(\text{CN}) = \frac{\exp\bigl[-w_f(\text{CN}-\text{cn}^{\text{ref}}_a)^2\bigr]}
{\sum_k\exp\bigl[-w_f(\text{CN}-\text{cn}^{\text{ref}}_k)^2\bigr]},
$$

so a reference whose tabulated coordination number is close to the atom's actual coordination number gets
nearly all the weight, and references far away in coordination-number space contribute almost nothing —
exactly a smooth nearest-neighbor interpolation in one variable. (When every Gaussian underflows to zero —
an extreme, out-of-table coordination number — the implementation falls back to giving full weight to
whichever reference has the highest tabulated coordination number, rather than producing a `0/0`
indeterminate weight.) The atom-pair $C_6$ coefficient is then a weighted double sum over both atoms'
reference tables:

$$
C_6^{AB}(\text{CN}_A,\text{CN}_B) = \sum_{a,b}w^a(\text{CN}_A)\,w^b(\text{CN}_B)\,C_{6,ab}^{\text{ref}}(Z_A,Z_B).
$$

## The Becke-Johnson damped two-body energy

Raw $1/R^6$ and $1/R^8$ terms diverge at short range, where the dispersion picture itself breaks down (the
fragments are no longer well-separated) and where the SCF's own exchange-repulsion should already be
handling the physics. **Becke-Johnson (BJ) rational damping** — qc-rs's default — smoothly attenuates both
terms using a *functional-independent* mathematical form with only two fitted parameters per functional
($a_1,a_2$; $s_6=1$ almost always, $s_8$ functional-specific):

$$
E_{\text{disp}}^{(2)} = -\sum_{A>B}C_6^{AB}\left[\frac{s_6}{R_{AB}^6+R_0^6} +
\frac{s_8\,\overline{R_{AB}^2}}{R_{AB}^8+R_0^8}\right], \qquad
R_0 = a_1\sqrt{3\,\langle R^4\rangle_A\langle R^4\rangle_B/\langle R^2\rangle_A\langle R^2\rangle_B} + a_2,
$$

where $\overline{R_{AB}^2}=3\langle R^4\rangle_A\langle R^4\rangle_B/(\langle R^2\rangle_A\langle
R^2\rangle_B)$ is the same tabulated per-element ratio ($\langle R^4\rangle/\langle R^2\rangle$, an
effective radius-squared measure) that also builds $R_0$, so both the damping radius and the $C_8/C_6$
ratio come from the same static, element-pair-specific quantity — no additional fitted parameter beyond
$a_1,a_2,s_8$ is needed to get the $R^{-8}$ term's relative weight right. **Rational damping's key
property**, unlike the older "zero damping" scheme, is that $E_{\text{disp}}^{(2)}$ approaches a finite,
nonzero constant (not zero) as $R_{AB}\to0$ — this reflects that at very short range the dispersion energy
should smoothly saturate rather than artificially vanish, since the underlying physical attraction does not
actually disappear, it is simply dominated by other (repulsive) terms in the total energy at that range.

## The three-body Axilrod-Teller-Muto term

Two-body dispersion is not always the whole story — for three atoms simultaneously close together, a
genuine three-body induced-dipole correlation effect (Axilrod-Teller-Muto, ATM) contributes an additional,
typically small but non-negligible correction:

$$
E_{\text{disp}}^{(3)} = \sum_{A>B>C}f_{\text{damp}}(\bar R)\;C_9^{ABC}\,\frac{3\cos\theta_A\cos\theta_B
\cos\theta_C+1}{(R_{AB}R_{BC}R_{CA})^3}, \qquad
C_9^{ABC} = -s_9\sqrt{|C_6^{AB}C_6^{AC}C_6^{BC}|},
$$

with $\theta_A,\theta_B,\theta_C$ the interior angles of the triangle formed by the three atoms (the
purely geometric angular factor that distinguishes a compact triangular arrangement from a nearly linear
one) and a zero-type damping function $f_{\text{damp}}(\bar R)=1/(1+6(R_0/\bar R)^{(\alpha+2)/3})$ built
from the same pairwise van der Waals radii $R_0$ the two-body term uses, with $\bar
R=\sqrt[3]{R_{AB}^2R_{BC}^2R_{CA}^2}$ and $s_9=1$ the standard three-body scaling (setting $s_9=0$ disables
ATM entirely). The $C_9$ coefficient's geometric-mean-of-products form — $-s_9\sqrt{|C_6^{AB}C_6^{AC}
C_6^{BC}|}$, not a separately-tabulated three-body reference — is itself an approximation reusing the
already-computed pairwise $C_6$ values rather than requiring a whole new set of three-body reference data.

Verified example — a water dimer, PBE/STO-3G, with and without D3(BJ):

```python
import qc
water_dimer = """
O  -1.551007  -0.114520   0.000000
H  -1.934259   0.762503   0.000000
H  -0.599677   0.040712   0.000000
O   1.350625   0.111469   0.000000
H   1.680398  -0.373741  -0.758561
H   1.680398  -0.373741   0.758561
"""
m0 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe").run()
m0.scf.energy   # -150.46818414059214

m1 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe", dispersion="d3bj").run()
m1.scf.energy, m1.scf.dispersion   # (-150.46956323055147, -0.0013790899594410507)
```

The dispersion correction lowers (stabilizes) the total energy by $\sim1.4\ \text{m}E_h$ for this
hydrogen-bonded dimer — small in absolute terms at this crude basis, but the qualitatively correct sign and
order of magnitude for an intermolecular attraction PBE alone does not fully capture.

## D4: from a fixed coordination number to self-consistent atomic charges

D3's coordination number is a purely geometric proxy for chemical environment; it cannot, by construction,
distinguish two atoms that are geometrically identical but chemically different because of the *molecule's
overall charge state or electronegativity balance* — a genuine electronic effect. **D4** generalizes D3 by
first solving for **atomic partial charges** via **electronegativity equilibration (EEQ)**, then using
those charges (in addition to a coordination number) to interpolate $C_6$, folding a real (if simplified)
electronic-structure step into what is otherwise still a fast, geometry-driven correction.

The EEQ charge model itself solves a small, per-molecule linear system — not an SCF, no basis-set
integrals, just tabulated atomic parameters and the molecular geometry — for a set of atomic charges $q_A$
that minimize the electrostatic-equilibration energy subject to a fixed total charge $q_{\text{tot}}$
(usually the molecular charge). The linear system augments the natural $n\times n$ electronegativity
equations with one extra row/column enforcing the charge-conservation constraint:

$$
\begin{pmatrix}A & \mathbf 1\\ \mathbf 1^{\mathsf T} & 0\end{pmatrix}
\begin{pmatrix}\mathbf q\\ \lambda\end{pmatrix} = \begin{pmatrix}\mathbf x\\ q_{\text{tot}}\end{pmatrix},
\qquad
A_{ii} = \eta_i + \sqrt{2/\pi}/\rho_i, \qquad A_{ij} = \frac{\operatorname{erf}\bigl(r_{ij}\sqrt{\gamma_{ij}}
\bigr)}{r_{ij}}, \qquad \gamma_{ij} = \frac{1}{\rho_i^2+\rho_j^2},
$$

with $x_i=-\chi_i+k_{\text{cn}\chi,i}\sqrt{\text{CN}_i}$ (a per-element electronegativity $\chi_i$, coupled
to a *different* coordination number — the D4 model's own error-function CN, distinct in its exact
functional form from D3's sigmoid CN above, though the same underlying physical idea), $\eta_i$ a chemical
hardness parameter, and $\rho_i$ a charge-density radius that sets how "smeared out" each atom's charge
density is taken to be (the $\operatorname{erf}(\cdot)/r$ off-diagonal form is exactly the Coulomb
interaction between two such smeared Gaussian charge distributions, not a bare point-charge $1/r$ — this is
what keeps the equations well-conditioned even at short interatomic distance). The Lagrange multiplier
$\lambda$ enforces $\sum_iq_i=q_{\text{tot}}$ exactly. Once solved, D4 uses these charges — alongside its
own coordination number — as a *second* interpolation axis for the reference $C_6$ table, giving each
atom's effective dispersion coefficient sensitivity to both its local geometric environment and its
electronic charge state.

Verified example — the same water dimer with D4 instead of D3(BJ):

```python
m2 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe", dispersion="d4").run()
m2.scf.energy, m2.scf.dispersion   # (-150.46918114898963, -0.0009970083974916537)
```

D4's correction here ($\sim-1.0\ \text{m}E_h$) is smaller in magnitude than D3(BJ)'s ($\sim-1.4\
\text{m}E_h$) for this system — the two models are calibrated independently against different (though
overlapping) reference data sets and are not expected to agree exactly, only to each capture the same
qualitative physics: dispersion stabilizes a hydrogen-bonded dimer relative to the same geometry with no
dispersion correction at all.

## Where dispersion sits relative to the SCF

Both D3 and D4 are evaluated **entirely outside the SCF self-consistency loop** — they depend only on
nuclear positions (D3) or nuclear positions plus a one-time EEQ solve (D4), never on the converged electron
density itself, and contribute no term to the Fock matrix. This is the structural opposite of [the PCM
solvation term](solvation-theory.md), which genuinely is density-dependent and must participate in every
SCF cycle and (per [the linear-response chapter](linear-response-cphf.md)) in the response kernel as well.
Dispersion's energy is simply added once, after (or independent of) SCF convergence — which is also why an
[analytic gradient](analytic-derivatives.md) of the dispersion term needs no orbital response of its own:
it is a function of nuclear coordinates alone, differentiated directly, with no Hellmann-Feynman subtlety
and no coupled-perturbed solve involved at all.

:::{exercise}
:label: ex-dispersion-theory

1. D3's coordination number uses a smooth sigmoid counting function rather than a hard "bonded/not bonded"
   cutoff. Explain in one sentence why a hard cutoff would make the resulting dispersion energy
   discontinuous as a function of geometry, and why that would be a serious problem specifically for
   geometry optimization and analytic gradients.
2. D4 uses EEQ charges as a *second* interpolation axis for $C_6$, in addition to (not instead of) a
   coordination number. Give one concrete scenario — two atoms with identical local geometric environments
   but different molecular contexts — where D3's coordination-number-only interpolation would be unable to
   distinguish them, but D4's added charge axis could.
3. Why does dispersion's analytic gradient need no CPHF/coupled-perturbed solve, unlike the analytic
   Hessian ([derived earlier](analytic-hessian-thermo.md)) or the PCM reaction-field response?
:::

:::{solution} ex-dispersion-theory
:class: dropdown

1. A hard cutoff would make the coordination number (and hence the interpolated $C_6$, and hence the
   dispersion energy) jump discontinuously the instant an atom pair's distance crosses the cutoff radius,
   even for an infinitesimally small change in geometry — the energy surface would have a genuine
   discontinuity (or at least a kink) at that distance. This is fatal for geometry optimization (which
   relies on the energy and its gradient varying smoothly to take sensible steps) and for analytic
   gradients specifically (which require the energy to be differentiable at all — a discontinuous function
   has no well-defined derivative exactly at the jump).
2. Two structurally identical fragments (say, two chemically equivalent aromatic carbons in geometrically
   similar local environments) that end up with different formal or partial charges because they sit in
   different overall molecular charge states, or different regions of a strongly polarized molecule, would
   receive *identical* D3 coordination numbers and therefore identical interpolated $C_6$ values — D3 has
   no variable through which to distinguish them. D4's EEQ-derived charges, computed from the whole
   molecule's electronegativity-equilibration balance rather than purely local geometry, would in general
   differ between the two, giving D4 a genuine second axis along which to assign them different, more
   chemically realistic dispersion coefficients.
3. Both D3 and D4 are functions purely of nuclear geometry (D4 additionally solving a one-time,
   density-independent EEQ linear system) — neither term ever depends on the SCF electron density or
   contributes to the Fock matrix. Differentiating a function of nuclear coordinates alone with respect to
   those same coordinates never needs an orbital-response term at all, because there are no orbitals to
   respond in the first place; CPHF/coupled-perturbed machinery exists specifically to handle *orbital*
   response to a perturbation, which is simply not a question dispersion's purely geometric energy
   expression ever raises.
:::

This chapter completes qc-rs's three Part II "Hamiltonian modification" chapters — [ECP](ecp-theory.md),
[PCM solvation](solvation-theory.md), and dispersion — each independently toggleable and composable with
any SCF reference and exchange-correlation functional. The next chapters turn from energetics to
**analysis**: [population/bond-order theory](population-analysis-theory.md) derives how to read chemical
meaning back out of a converged density and MOs.