Dispersion correction theory: DFT-D3/D4

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

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{split} \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}, \end{split}\]

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

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, which genuinely is density-dependent and must participate in every SCF cycle and (per the linear-response chapter) 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 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 15

  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) or the PCM reaction-field response?

This chapter completes qc-rs’s three Part II “Hamiltonian modification” chapters — ECP, PCM solvation, 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 derives how to read chemical meaning back out of a converged density and MOs.