5.7.5 Many-Body Dispersion (MBD) Method

Unlike earlier DFT-D methods that were strictly (atomic) pairwise-additive, DFT-D3 includes three-body (triatomic) corrections. These terms are significant for non-covalent complexes assembled from large monomers,⁵³³ especially those that contain a large number of polarizable centers.²³⁵ The many-body dispersion (MBD) method of Tkatchenko et al.^{921, 34} represents a more general and less empirical approach that goes beyond the pairwise-additive treatment of dispersion. This is accomplished by including $n$-body contributions to the dispersion energy up to the number of atoms, and polarization screening contributions to infinite order. Even in small systems such as benzene dimer, the MBD approach consistently outperforms other popular vdW methods.⁹⁹

The essential idea behind MBD is to approximate the dynamic response of a system by that of dipole-coupled quantum harmonic oscillators (QHOs), each of which represents a fragment of the system of interest. The correlation energy of such a system can then be evaluated exactly by diagonalizing the corresponding Hamiltonian:³⁷⁶

$${\widehat{H}}_{\text{MBD}}=\frac{1}{2}\sum _A^{\text{atoms}}{\widehat{\nabla }}_{{\xi }_A}^2+\frac{1}{2}\sum _A^{\text{atoms}}{\omega }_A^2{\xi }_A^2+\frac{1}{2}\sum _{A,B}^{\text{atoms}}{\omega }_A{\omega }_B{({\alpha }_A^0{\alpha }_B^0)}^{1/2}{\xi }_A{T}_{AB}{\xi }_B.$$

(5.43)

Here, ${\xi }_A=m_A^{1/2}|{𝐫}_A-{𝐑}_A|$ is the mass-weight displacement of oscillator $A$ from its center ${𝐑}_A$, ${\omega }_A$ is the characteristic frequency, and ${\alpha }_A^0$ is the static polarizability. $T_{AB}$ is the dipole potential between the oscillators $A$ and $B$. The MBD Hamiltonian is obtained through coarse-graining of the long-range correlation (through the long-range dipole tensor ${𝐓}_{\mathrm{lr}}$) and approximating the short-range polarizability via the adiabatic connection fluctuation-dissipation formula:

$$E_{c,\text{lr}}^{\mathrm{MBD}}=-\sum _{n=2}^{\mathrm{\infty }}\frac{{(-1)}^n}{n}{\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}u}{2\pi }\sum _{AB}\text{tr}[⟨{({\alpha }_{\mathrm{eff}}{𝐓}_{\mathrm{lr}})}^n{⟩}_{AB}(iu)].$$

(5.44)

This approximation expresses the dynamic polarizability ${\alpha }_{\mathrm{eff}}$ (of a given fragment) in terms of the polarizability of the corresponding QHO,

$${\alpha }_A^{\mathrm{QHO}}(u)=\frac{q_A^2}{m_A({\omega }_A^2-u^2-i\delta u)}$$

(5.45)

in which $q_A$ is the charge, $m_A$ the mass, and ${\omega }_A$ the characteristic frequency of the oscillator. The integration in the frequency domain in Eq. (5.44) can be done analytically, leading to the so-called plasmon pole formula for the correlation energy,

$$E_c=\frac{1}{2}\sum _{p=1}^{3N}({\overline{\omega }}_p-{\omega }_p)$$

(5.46)

in which $N$ is the number of fragments, ${\overline{\omega }}_p$ are the frequencies of the interacting (dipole-coupled) system, and ${\omega }_p$ are the frequencies of the non-interacting system (i.e., the collection of independent QHOs). The sum runs over all $3N$ characteristic frequencies of the system.

A particular method within the MBD framework is defined by the models for the static polarizability (${\alpha }_{\text{eff},A}^0$), the non-interacting characteristic frequencies (${\omega }_A$), and the damping function [$f(R)$] used to define ${𝐓}_{\text{lr}}$. In Q-Chem, the MBD method is implemented following the “MBD@rsSCS” approach, where “rsSCS” stands for range-separated self-consistent screening.³⁴ In this approach, ${\alpha }_{\text{eff},A}^0$ is obtained in a two-step process:

1. 1.

   The free-volume scaling approach is applied to the free-atom polarizabilities, using the Hirshfeld-partitioned molecular electron density. This is the same procedure used in the TS-vdW method described in Section 5.7.4.

2. 2.

   The short-range atomic polarizabilities ${\alpha }_{\mathrm{sr},AB}(iu)$ are obtained by applying a Dyson-like screening on only the short range part of the polarizabilities. The same range-separation will later be used to define ${𝐓}_{\mathrm{lr}}$.

The short-range atomic polarizabilities are summed up along one fragment coordinate to obtain the local effective dynamic polarizability, i.e., ${\alpha }_{\mathrm{eff},A}^0={\sum }_B{\alpha }_{\mathrm{sr},AB}$, and are then spherically averaged. The range-separation (damping) function $f(R)$ used to construct ${\alpha }_{\mathrm{sr},AB}(iu)$ and ${𝐓}_{\mathrm{lr}}$ is the same as that in Eq. (5.38), except with $d=6$ instead of $d=20$, and again $s_R$ for a given functional obtained by fitting to interaction energies for non-bonded complexes. The MBD energy is than calculated by diagonalizing the Hamiltonian Eq. (5.43) and using the plasmon-pole formula, Eq. (5.46).

The MBD-vdW approach greatly improves the accuracy of the interaction energies for S66⁷⁸⁵ test set, even if a simple functional like PBE is used, with a mean absolute error of 0.3 kcal/mol and a maximum error of 1.3 kcal/mol, as compared to 2.3 kcal/mol (mean) and 7.2 kcal/mol (max) for plain PBE. In general, the MBD-vdW method is superior to pairwise a posteriori dispersion corrections.³⁷⁶

As mentioned above in the context of the TS-vdW method (Section 5.7.4), the FHI-aims or Quantum Espresso codes cannot perform exact unrestricted SCF calculations for the atoms and this leads to inconsistent free-atom volumes as compared to the spherical ones computed in Q-Chem, and thus inconsistent values for the vdW correction. Since the parameters of the TS-vdW and MBD-vdW models were fitted for use with FHI-aims and Quantum Espresso, results obtained using these codes are slightly closer to S66 benchmarks and thus scalar modifiers are available for the internally-computed Hirshfeld volume ratios in Q-Chem. For S66, the use of these modifiers leads to negligible differences between results obtained from all three codes.⁶⁰

The MBD-vdW correction is requested by setting MBDVDM = TRUE in the $rem section. Other job control variables, including the aforementioned modifiers for the free-atom volume ratios, are the same as those for the TS-vdW method and are described in Section 5.7.4.

$molecule
   0 1
   O
   H 1 0.95
   H 1 0.95 2 104.5
$end

$rem
   BASIS               6-31G*
   METHOD              PBE
   !vdw settings
   MBDVDW              TRUE
   !setting the SCF_CONVERGENCE for single
   !atom calculations to 6
   HIRSHFELD_CONV      6
   !Apply modifiers to the free-atom volumes
   !the elements H, C, N, and O
   HIRSHMOD            4

$end