Solvation & dispersion

Two physical effects are missing from a gas-phase HF or DFT calculation, and both have cheap, standard corrections in qc-rs. Solvation puts the molecule in a solvent instead of a vacuum. Dispersion adds the long-range van der Waals attraction that most functionals leave out. They are independent — use either, both, or neither — and both flow automatically into the gradient.

Implicit solvation (PCM)

Theory

Most chemistry happens in solution, not vacuum. Modelling every solvent molecule explicitly is expensive; an implicit (continuum) model replaces the solvent with a polarizable dielectric continuum characterized by one number, the dielectric constant \(\varepsilon\) (≈ 78.4 for water, ≈ 1 for vacuum). The solute sits in a molecule-shaped cavity carved out of that continuum; its charge distribution polarizes the dielectric, which in turn creates a reaction field that acts back on the solute. The polarizable continuum model (PCM) solves for that mutual polarization self-consistently. Because the reaction field depends on the solute density and vice versa, PCM folds directly into the SCF loop as an extra Fock term.

Concretely, the cavity surface is discretized (by GePol) into small area elements — tesserae — and the reaction field is captured by an apparent surface charge (ASC) \(\mathbf q\), one value per tessera, found by solving a boundary-integral equation against the solute’s electrostatic potential \(\mathbf v\) at each tessera. The two formulations qc-rs implements differ in how they enforce the dielectric boundary condition:

\[ \textbf{IEF-PCM (default):}\quad \mathbf T(\varepsilon)\,\mathbf q = -\mathbf R_\infty\,\mathbf v, \qquad \mathbf T(\varepsilon) = \big[2\pi f(\varepsilon)\,\mathbf I - \mathbf D\mathbf A\big]\mathbf S, \quad f(\varepsilon)=\frac{\varepsilon+1}{\varepsilon-1}, \quad \mathbf R_\infty = 2\pi\mathbf I - \mathbf D\mathbf A, \]

\[ \textbf{C-PCM (COSMO scaling):}\quad \mathbf S_{\text{scaled}}\,\mathbf q = -\mathbf v, \qquad \mathbf S_{\text{scaled}} = \mathbf S\Big/\frac{\varepsilon-1}{\varepsilon+k}, \]

where \(\mathbf S\) and \(\mathbf D\) are the single- and double-layer boundary-integral operators (built from the Green’s function of the dielectric — \(1/r\) for the default vacuum-outside case), and \(\mathbf A\) is the diagonal tessera-area matrix. Once \(\mathbf q\) is known, the polarization energy is simply \(E_{\text{pol}}=\tfrac12\,\mathbf v\cdot\mathbf q\). Because \(\mathbf v\) is built from the current SCF density, this linear system is solved every SCF cycle — which is exactly why PCM is implemented as an extra Fock term rather than a one-shot correction applied after convergence. The full cavity-construction and Green’s-function details are in Implicit solvation theory.

Usage

Pass pcm= to scf(...): a bare number is the dielectric constant, or a dict selects the model:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

gas = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
sol = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", pcm=78.39).run()

print(f"gas   E = {gas.scf.energy:.6f}")     # -76.026794
print(f"water E = {sol.scf.energy:.6f}")     # -76.036749
print(f"ΔG_solv = {(sol.scf.energy - gas.scf.energy)*627.509:.3f} kcal/mol")   # -6.247

The solvated energy is lower — the dielectric stabilizes the molecule by −6.25 kcal/mol here (the electrostatic solvation free energy). The dict form selects the formulation:

scf(..., pcm={"model": "iefpcm", "epsilon": 78.39})   # IEF-PCM (default)
scf(..., pcm={"model": "cpcm",   "epsilon": 78.39})   # C-PCM

• IEF-PCM (integral-equation-formalism, the default) and C-PCM (conductor-like) are the two standard formulations; they give very slightly different energies (here −76.036749 vs −76.036812).

• On a cuda build, add "device": "cuda" (or "auto") to run the AO↔cavity coupling on the GPU with an identical energy.

PCM composes with HF and KS-DFT, and its analytic gradient is included, so you can optimize a geometry in solvent.

Dispersion corrections (DFT-D3 / D4)

Theory

London dispersion — the weak, attractive van der Waals force between instantaneously-induced dipoles — is a pure long-range correlation effect. Hartree–Fock misses it entirely, and most semi-local DFT functionals capture little of it, so they underbind π-stacks, layered materials, host–guest complexes, and molecular crystals. The DFT-D corrections add it back as a cheap, geometry-only atom-pairwise sum.

The default, D3 with Becke–Johnson (rational) damping, is exactly

\[ E_{\text{disp}}^{\text{D3(BJ)}} = -\sum_{A<B} C_6^{AB} \left( \frac{s_6}{R_{AB}^6 + R_0^6} + \frac{s_8\,\big(C_8^{AB}/C_6^{AB}\big)}{R_{AB}^8 + R_0^8} \right), \qquad R_0 = a_1\sqrt{C_8^{AB}/C_6^{AB}} + a_2 . \]

Three physically distinct ingredients feed this:

• \(C_6^{AB}\) is not a fixed atomic constant — it is interpolated from each atom’s coordination number (a continuous count of its neighbors, from a smooth distance-based switching function), so a carbon in a \(\text{sp}^3\) environment gets a different \(C_6\) than an aromatic carbon. This local-environment sensitivity is what the “3” in D3 refers to.

• \(C_8^{AB}/C_6^{AB} = 3\sqrt{\langle r^4\rangle_A/\langle r^2\rangle_A}\cdot\sqrt{\langle r^4\rangle_B/\langle r^2\rangle_B}\) from tabulated per-element radii.

• \(s_6,s_8,a_1,a_2\) are the four Becke–Johnson parameters, refit for each DFT functional (looked up by the functional name you pass to dispersion(...)).

Rational (BJ) damping’s whole virtue is in that \(+R_0^6\) inside the denominator: unlike an exponential switching function, it never diverges as \(R_{AB}\to0\), so the correction stays finite (and small, since \(C_6/R_0^6\) is bounded) even for a short, strongly-bonded contact — no separate short-range cutoff logic is needed. d3zero uses an older, exponential switching form instead (hence its near-zero contribution for the compact water molecule in the table below); d4 replaces the coordination-number-only \(C_6\) with geometry-dependent EEQ partial charges, refining the local-environment sensitivity further.

Usage

mychk.dispersion(functional, method=...) returns the dispersion energy to add to the SCF energy (it is geometry-only, so it needs no SCF first):

mol = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom")

mol.dispersion("pbe")                      # -0.00035947   D3(BJ), the default
mol.dispersion("pbe", method="d3zero")     # -0.00000363   D3 zero-damping
mol.dispersion("pbe", method="d3bjatm")    # -0.00035947   D3(BJ) + ATM 3-body term
mol.dispersion("pbe", method="d4")         # -0.00019523   D4 (EEQ charges + ATM)

total = mol.scf(ref="r", xc="pbe").run().scf.energy + mol.dispersion("pbe")   # dispersion-corrected total

method=	correction
d3bj (default)	D3 with Becke–Johnson damping — the standard choice
d3zero	D3 with zero damping (an older damping form)
d3bjatm	D3(BJ) plus the Axilrod–Teller–Muto 3-body term
d4	D4 — geometry-dependent charges (EEQ) + ATM

Only real atoms contribute — ghost and dummy atoms are excluded. The dispersion energy for water is tiny (it has no dispersion-bound contacts); the correction matters for larger, weakly-bound assemblies, where it can be several kcal/mol. The d3bjatm here equals d3bj to the printed digits because water’s 3-body (ATM) term is negligible — it grows in dense, polarizable systems.

Tip

Dispersion inside the SCF workflow Passing dispersion=... to scf(...) (e.g. scf(ref="r", xc="pbe", dispersion="d3bj")) folds the correction into the reported scf.energy and the gradient automatically, which is what you want for a dispersion-corrected geometry optimization. The standalone mychk.dispersion(...) above is the geometry-only accessor for the energy term itself.

Worked example: both corrections

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

gas   = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe").run().scf.energy
solv  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe", pcm=78.39).run().scf.energy
disp  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").dispersion("pbe")

print(f"gas-phase PBE      : {gas:.6f}")
print(f"+ solvation (water): {solv:.6f}   (ΔG_solv = {(solv-gas)*627.509:.2f} kcal/mol)")
print(f"+ dispersion (D3BJ): {gas + disp:.6f}")

Exercise 8

1. You compute a gas-phase reaction energy and it disagrees with an aqueous experiment. Which correction in this chapter is the natural first thing to add, and how do you request it?

2. Why does mychk.dispersion(...) need no .scf().run() first, whereas pcm= must go inside scf(...)?

3. For water the dispersion correction is ~0.0004 Ha (~0.2 kcal/mol). Name a system where you would expect it to be far larger, and why.

Solution to Exercise 8

1. Implicit solvation — the experiment is in water, your calculation is in vacuum. Add pcm=78.39 (or pcm={"model":"iefpcm","epsilon":78.39}) to the scf(...) call.

2. Dispersion (D3/D4) is a geometry-only function of the atomic positions and the functional — it does not depend on the electron density — so it needs no wavefunction. PCM’s reaction field does depend on the density and is solved self-consistently with the SCF, so it must be part of scf(...).

3. Anything bound by van der Waals forces: a π-stacked dimer (e.g. two benzenes), a layered material (graphite), a host–guest complex, or a molecular crystal. There the dispersion attraction is the bonding, so the correction can reach several kcal/mol and change whether the system is bound at all.

Next, visualization turns these densities and orbitals into pictures you can inspect.