5.7 DFT Methods for van der Waals Interactions

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,533 especially those that contain a large number of polarizable centers.235 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.99

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:376

H^MBD=12Aatoms^ξA2+12AatomsωA2ξA2+12A,BatomsωAωB(αA0αB0)1/2ξATABξB. (5.43)

Here, ξA=mA1/2|𝐫A-𝐑A| is the mass-weight displacement of oscillator A from its center 𝐑A, ωA is the characteristic frequency, and αA0 is the static polarizability. TAB 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 𝐓lr) and approximating the short-range polarizability via the adiabatic connection fluctuation-dissipation formula:

Ec,lrMBD=-n=2(-1)nn0du2πABtr[(αeff𝐓lr)nAB(iu)]. (5.44)

This approximation expresses the dynamic polarizability αeff (of a given fragment) in terms of the polarizability of the corresponding QHO,

αAQHO(u)=qA2mA(ωA2-u2-iδu) (5.45)

in which qA is the charge, mA the mass, and ω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,

Ec=12p=13N(ω¯p-ωp) (5.46)

in which N is the number of fragments, ω¯p are the frequencies of the interacting (dipole-coupled) system, and ω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 (αeff,A0), the non-interacting characteristic frequencies (ωA), and the damping function [f(R)] used to define 𝐓lr. In Q-Chem, the MBD method is implemented following the “MBD@rsSCS” approach, where “rsSCS” stands for range-separated self-consistent screening.34 In this approach, αeff,A0 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 α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 𝐓lr.

The short-range atomic polarizabilities are summed up along one fragment coordinate to obtain the local effective dynamic polarizability, i.e., αeff,A0=Bαsr,AB, and are then spherically averaged. The range-separation (damping) function f(R) used to construct αsr,AB(iu) and 𝐓lr is the same as that in Eq. (5.38), except with d=6 instead of d=20, and again sR 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 S66785 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.376

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.60

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.

MBDVDW
       Flag to switch on the MBD-vdW method
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 Do not calculate MBD. 1 Calculate the MBD-vdW contribution to the energy. 2 Calculate the MBD-vdW contribution to the energy and the gradient.
RECOMMENDATION:
       NONE

Example 5.12  Sample input illustrating a calculation of a water molecule, including the MBD-vdW energy.

$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