Solvation & dispersion

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.

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,

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

with tabulated dispersion coefficients \(C_n^{AB}\) and a damping function that switches the correction off at short range (where the functional already describes the interaction). The parameters depend on the functional.

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 7

  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.

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