Dispersion corrections (D3-BJ)

Standard density-functional approximations (LDA, PBE, even B3LYP) are local or semi-local: the exchange-correlation kernel at a point only sees the density around that point. That makes them blind to long-range dispersion (London / van-der-Waals) — the electron correlation that attracts non-overlapping fragments. For molecular crystals, layered materials, π-stacking, drug-receptor binding, and host-guest chemistry, the missing dispersion energy is often 1-10 kcal/mol per contact, which is everything.

Grimme’s D3(BJ) correction is a pairwise sum

\[ E_{\text{disp}} = -\sum_{A < B} \Bigl( \frac{s_6 \, C_6^{AB}}{R_{AB}^6 + f(R_0^{AB})^6} + \frac{s_8 \, C_8^{AB}}{R_{AB}^8 + f(R_0^{AB})^8} \Bigr) \]

parametrised against high-accuracy reference data, with functional-specific damping parameters (s6, s8, a1, a2). It’s additive to the SCF energy — no re-converging the density — and vibe-qc wires it through run_job with a single argument.

The one-liner

from vibeqc import Atom, Molecule, run_job

mol = Molecule([
    Atom(8, [ 0.0,  0.00,  0.00]),
    Atom(1, [ 0.0,  1.43, -0.98]),
    Atom(1, [ 0.0, -1.43, -0.98]),
])

run_job(
    mol,
    basis="6-31g*",
    method="rks",
    functional="PBE",
    dispersion="d3bj",         # <-- adds PBE-D3(BJ)
    output="water_pbe_d3",
)

The output file gets a new “Dispersion correction” block:

  Dispersion correction (D3-BJ)
  ----------------------------------------------------
          s6       1.000000
          s8       0.787500
          a1       0.428900
          a2       4.440700
      E_disp    -0.00070434 Ha  (-0.4420 kcal/mol)
       E_SCF   -76.31993075 Ha
     E_total   -76.32063509 Ha

The returned result object carries .energy (bare SCF), .e_dispersion (the correction), and .energy_total (sum):

r = run_job(mol, basis="6-31g*", method="rks", functional="PBE",
            dispersion="d3bj", output="water_pbe_d3",
            write_molden_file=False)
print(f"E_SCF    = {r.energy:.6f} Ha")
print(f"E_disp   = {r.e_dispersion:+.6e} Ha")
print(f"E_total  = {r.energy_total:.6f} Ha")

A worked example: the water dimer

Where dispersion matters is in relative energies — it barely changes the energy of a single water molecule, but it stabilises the dimer by a real amount. Compute both:

dimer = Molecule([
    Atom(8, [ 0.0,   0.00,   0.00]),
    Atom(1, [ 0.0,   0.00,   1.814]),    # 0.96 Å along z
    Atom(1, [ 0.0,   1.756, -0.454]),
    Atom(8, [ 0.0,   0.00,   5.575]),    # O-O 2.95 Å
    Atom(1, [ 0.0,   1.434,  6.686]),
    Atom(1, [ 0.0,  -1.434,  6.686]),
])
water = Molecule([
    Atom(8, [0.0, 0.0, 0.0]),
    Atom(1, [0.0, 0.0, 1.814]),
    Atom(1, [0.0, 1.756, -0.454]),
])

for mol, label in ((water, "monomer"), (dimer, "dimer")):
    r_bare = run_job(mol, basis="6-31g*", method="rks", functional="PBE",
                     output=f"{label}", write_molden_file=False)
    r_d3 = run_job(mol, basis="6-31g*", method="rks", functional="PBE",
                   dispersion="d3bj", output=f"{label}_d3",
                   write_molden_file=False)
    print(f"{label:10s}  PBE   {r_bare.energy:.6f}   "
          f"PBE-D3  {r_d3.energy_total:.6f}   "
          f"Δdisp = {r_d3.e_dispersion*627.509:+.2f} kcal/mol")

Results:

monomer     PBE   -76.319931   PBE-D3  -76.320636   Δdisp = -0.44 kcal/mol
dimer       PBE  -152.650077   PBE-D3 -152.652627   Δdisp = -1.60 kcal/mol

Intra-molecular dispersion contributes ~0.44 kcal/mol to each isolated water; the dimer picks up 1.60 kcal/mol total. The extra 1.60 - 2·0.44 = 0.72 kcal/mol is the genuine inter-molecular dispersion — the part that shows up in binding energies:

Quantity

Value

Binding energy (PBE)

−6.41 kcal/mol

Binding energy (PBE-D3)

−7.13 kcal/mol

D3 contribution

−0.72 kcal/mol

Both values over-bind the experimental (CCSD(T)/CBS) reference of ~5.0 kcal/mol — that’s a known PBE issue unrelated to dispersion. D3 is additive, not corrective.

Picking damping parameters

dispersion="d3bj" auto-looks up the damping parameters for the functional= you passed. For a handful of common functionals vibe-qc ships a builtin parameter table (the “D1a starter set”); for the full ~80-functional Grimme set, install the [dispersion] extra:

pip install 'vibe-qc[dispersion]'

which pulls in the reference dftd3 package — dispersion="d3bj" then transparently resolves any functional it knows.

You can also pass an explicit functional-name string to pick a different parametrisation than the one the SCF ran on (useful for HF + D3 or for cross-comparison):

from vibeqc import D3BJParams, d3bj_params_for

params = d3bj_params_for("b3lyp")          # look up B3LYP's damping
print(params)                              # D3BJParams(s6=..., s8=..., a1=..., a2=...)

# Pass an explicit D3BJParams instance instead of a string:
custom = D3BJParams(s6=1.0, s8=1.5, a1=0.4, a2=4.5)
run_job(mol, basis="6-31g*", method="rks", functional="pbe",
        dispersion=custom, output="custom_d3")

D3 with Hartree-Fock

HF has no correlation at all, so dispersion is a bigger miss there than for DFT. HF + D3(BJ) uses Grimme’s HF-specific damping set, which lives in the full dftd3 reference parameter file rather than vibe-qc’s builtin starter set. Install the extra first:

pip install 'vibe-qc[dispersion]'

then pass the functional lookup name Grimme uses for HF:

run_job(mol, basis="6-31g*", method="rhf",
        dispersion="hf",           # damping params for HF, via dftd3
        output="hf_d3")

Any basis-set-superposition-error concerns that apply to plain HF still apply — counterpoise correction is a separate topic, not baked into vibe-qc yet.

D3 in geometry optimisation

D3 has analytic gradients, so it rides through optimize=True without extra configuration:

run_job(mol, basis="6-31g*", method="rks", functional="pbe",
        dispersion="d3bj",
        output="dimer_d3_opt",
        optimize=True, fmax=0.2)

Each BFGS step adds the dispersion gradient to the SCF gradient; the final .out file reports E_SCF, E_disp, and E_total at the optimised geometry.

Caveats

  • D3 is pairwise by design. It doesn’t see the electron density at all — just atomic positions and element types. It’s fast (adds microseconds to a DFT calculation) but it can’t fix an already-wrong density (e.g. PBE’s tendency to over-polarise H-bonds).

  • No replacement for a correlated wavefunction. If you need chemical accuracy for weak-interaction energies, use CCSD(T)-F12 or RPA, not DFT-D3. D3 gets you most of the way with small-basis DFT for a lot of work.

  • Parameters are functional-specific. Using PBE’s damping with B3LYP double-counts short-range correlation. Always match functional= and the dispersion= lookup (the default behaviour does this automatically).

  • BSSE is a separate topic. Pair binding energies still benefit from counterpoise correction when the basis is moderate. Run the monomer in the dimer basis for a CP-corrected number — vibe-qc doesn’t automate this yet.

Next

  • Bigger systems where dispersion dominates: π-stacked aromatics (benzene dimer, DNA bases), layered materials, rare-gas crystals. Standalone examples will land as the relevant geometries make it into examples/.

  • D3 in periodic calculations needs a cutoff treatment for the pairwise sum — not yet wired into run_rhf_periodic / run_rks_periodic. Track that in the roadmap.