Dispersion interactions

Long-range interaction energy¶

The interaction energy between two molecular systems $A$ and $B$ is given by

For charge neutral, nonpolar, systems separated by a sufficiently large distance $R$ this interaction energy is known as dispersion (or van der Waals) energy.

Dispersion and correlation¶

The prototypical system exemplifying dispersion interactions is the He dimer for which the interaction energy equals

We will determine this interaction energy with respect to the interatomic separation distance at the Hartree–Fock and MP2 levels of theory.

Determine the atomic energies the HF and MP2 levels of theory.

Determine the dimer energies over a range of interatomic separation distance.

Plot the interaction energies in units of micro-Hartree.

The failure of the HF method to describe the minimum of the potential energy curve for the helium dimer is due to the lack of correlation. More specifically, let us consider a system setup where the $z$-axis is chosen as the internuclear axis and the helium atoms are placed at $z = 0$ and $z=R_0$, respectively (see the inset in the figure below). We will then be concerned with the two-particle density, $n(\mathbf{r}_1, \mathbf{r}_2)$, at coordinates

At the HF level of theory, the two-particle density will be independent of the angle $\theta$ whereas this is not the case when electron correlation is accounted for. The MP2 density shows a maximum when the two electrons are positioned as far apart as possible, or in other words as $\theta = 180^\circ$.

This asymmetry in the two-particle density is referred to as fluctuating induced dipoles, or instantaneous dipole–dipole interactions, and it is the source origin of the weakly attractive van der Waals energies.

Dispersion by perturbation theory¶

The long-range interaction energy between two randomly oriented molecules $A$ and $B$ is given by the Casimir–Polder potential

where $R$ is the intermolecular separation, $c$ is the speed of light, and $\overline{\alpha}_A(i\omega^I)$ is the isotropic average of the electric dipole polarizability tensor of molecule $A$ evaluated at a purely imaginary frequency. The molecules are here considered to be polarizable entities. Additional contributions to the energy will arise from higher-order electric multipole interactions as well as magnetic interactions.

The interaction energy is often expressed in terms of dispersion coefficients $C_m$, where $m = 6, 8, 10, \ldots$, and long-range coefficients $K_n$, where $n = 7, 9, 11, \ldots$, so that, in the van der Waals region, we have

and, at very large intermolecular separation, we have

The reason for the deviation of the interaction energy at large $R$ from a $1/R^6$ dependence, as given by London–van der Waals dispersion theory, is the finite speed of the photons mediating the electromagnetic interaction, or, equivalently, retardation.

The leading dispersion coefficient is given by

Based on Hartree–Fock data and with access to the $C_6$ dispersion coefficient, a much more realistic potential energy curve can be constructed as follows

where the first and second term contribute the repulsive and attractive parts of the potential, respectively.

A FCI reference value for the binding energy, 29 $\mu E_h$, and equilibrium distance, 2.96 Å Vos et al. (1990), is depicted in the figure as the crossing point of the vertical and horizontal lines. Both the MP2 and the dispersion-corrected HF minima are quite close to this FCI reference point.

It has become a standard procedure to adopt a similar approach to add dispersion corrections to DFT potentials as standard functionals do not offer an inclusion of dispersion interactions. This can be very important when performing molecular structure optimizations of systems with nonbonded interactions.

Dispersion by symmetry-adapted perturbation theory (SAPT)¶

In the supermolecular approach above, the weak interaction energy is obtained as a small difference between large total energies---and nothing in the calculation tells us which physical effects make up the result. Symmetry-adapted perturbation theory (SAPT) Jeziorski et al. (1994) instead treats the interaction between the monomers as a perturbation---the zeroth-order problem being the two isolated monomers---and the interaction energy is obtained directly, order by order, as a sum of physically distinct contributions. At the lowest level, known as SAPT0, the monomers are described at the HF level and the interaction operator is treated to second order

where the first superscript gives the order in the intermolecular interaction and the second the order in the intramonomer correlation, which we keep at zerohere. The terms are

• Electrostatics, the classical Coulomb interaction between the unperturbed HF charge distributions of the monomers, nuclei included.

• Exchange, the Pauli repulsion that arises when the dimer wavefunction is antisymmetrized with respect to exchange of electrons between the monomers. Note that it is here evaluated in the single-exchange, or $S^2$, approximation.

• Induction, the energy lowering due to the polarization of each monomer in the static field of the other, obtained from the CPHF equations.

• Exchange--induction, the exchange counterpart of induction, which quenches it at short range.

• Dispersion, the attraction due to correlated fluctuations of the two charge distributions. This is exactly the effect that is absent at the HF level and that we discussed above in terms of the two-particle density.

The exchange--dispersion term and the so-called $\delta$HF correction, which collects higher-order induction effects from a supermolecular HF calculation, are left out in the implementation below.

There is no SAPT driver in VeloxChem, but all terms can be assembled from the integrals supplies together with the monomer HF solutions. Each monomer is computed in the full dimer basis set by placing ghost basis functions (Bq) on the partner atom---the standard practice in SAPT. The equations are collected below, and the implementation is found in the code cell that follows, and are based on Xie et al. (2022).

With ${\mathbf{C}^X_\mathrm{occ}}$ and ${\mathbf{C}^X_\mathrm{vir}}$ denoting the occupied and virtual MO coefficients of monomer $X$, we define the AO density matrices ${\mathbf{P}^X = \mathbf{C}^X_\mathrm{occ} \mathbf{C}^{X\dagger}_\mathrm{occ}}$ and ${\mathbf{P}^X_\mathrm{vir} = \mathbf{C}^X_\mathrm{vir} \mathbf{C}^{X\dagger}_\mathrm{vir}}$, the Coulomb and exchange matrices ${\mathbf{J}[\mathbf{D}]_{\mu\nu} = \sum_{\lambda\sigma} D_{\lambda\sigma} (\mu\nu|\lambda\sigma)}$ and ${\mathbf{K}[\mathbf{D}]_{\mu\nu} = \sum_{\lambda\sigma} D_{\lambda\sigma} (\mu\lambda|\nu\sigma)}$, the electrostatic potential of monomer $X$, ${\boldsymbol{\omega}^X = \mathbf{V}^X + 2\mathbf{J}[\mathbf{P}^X]}$, with $\mathbf{V}^X$ the nuclear attraction matrix of $X$, as well as ${\mathbf{h}^X = \boldsymbol{\omega}^X - \mathbf{K}[\mathbf{P}^X]}$ and ${\mathbf{O} = \mathbf{P}^A \mathbf{S} \mathbf{P}^B}$, where $\mathbf{S}$ is the AO overlap matrix. With ${\mathbf{A} \cdot \mathbf{B} = \sum_{\mu\nu} A_{\mu\nu} B_{\mu\nu} }$ and $V_\mathrm{nuc}$ the nuclear repulsion between the monomers, the first-order terms read

The induction energy requires the response of monomer $A$ to the field of $B$, with $\tilde{\boldsymbol{\omega}}^B = \mathbf{C}^{A\dagger}_\mathrm{occ} \boldsymbol{\omega}^B \mathbf{C}^A_\mathrm{vir}$, the CPHF amplitudes $\mathbf{x}^A$ are obtained from

and correspondingly for $\mathbf{x}^B$. The induction and exchange--induction energies are then given by

with the exchange--induction potential

and $\mathbf{W}^{B \leftarrow A}$ obtained by swappin the monomer labels, and replacing $\mathbf{O}$ by $\mathbf{O}^T$. Finally, the uncoupled dispersion energy is a sum over the occupied ($i$, $j$) and virtual ($a$, $b$) orbitals of $A$ and $B$, respectively,

We can now compute the SAPT0 terms for the same set of separations as above.

Let us add the SAPT0 interaction energies to the previous figure:

The minimum of the SAPT0 potential is close to the MP2 and the dispersion-corrected HF results. The similarity to MP2 is not a coincidence, as the uncoupled dispersion term is essentially the intermolecular part of the second-order correlation energy, and SAPT0 and MP2 interaction energies are therefore typically of similar quality for dispersion-bound systems. What is missing compared to FCI is mainly the intramonomer correlation, i.e., the zero in the second superscript of the SAPT0 terms, which affects the monomer charge distributions, and so all contributions.

The real strength of SAPT, however, lies in the decomposition of the interaction energy, providing insight into how each term contributes. The decomposition is plotted below, together with the supermolecular results from above.

The decomposition makes clear what we inferred above from the comparison of the HF and MP2 curves. The electrostatic, exchange, induction, and exchange--induction terms add up to a curve that is nearly indistinguishable from the supermolecular HF interaction energy, the difference being the $\delta$HF correction and the $S^2$ approximation. In other wrods, the purely repulsive HF potential results from the exchange term dominating over the (weak) electrostatic attraction. The dispersion term, in turn, closely follows the correlation contribution to the MP2 interaction energy, $\Delta E^\mathrm{MP2} - \Delta E^\mathrm{HF}$. As such, SAPT not only provides interaction energies without the subtraction of large total energies, but also an interpretation of the result in terms of well-defined physical contributions. The dispersion corrections to DFT mentioned above can in this light be seen as adding the one term in this decomposition that approximate exchange--correlation functionals lack.