# 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](gradients-geomopt.md).

## 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](../10-foundations/solvation-theory.md).

### Usage

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

```python
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:

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

```python
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

```python
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}
:label: ex-solv-disp

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} ex-solv-disp
:class: dropdown

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](visualization.md) turns these densities and orbitals into pictures you can inspect.