5.7.3 Empirical Dispersion Corrections: DFT-D

A major development in DFT during the mid-2000s was the recognition that, first of all, semi-local density functionals do not properly capture dispersion (van der Waals) interactions, a problem that has been addressed only much more recently by the non-local correlation functionals discussed in Section 5.7.2; and second, that a cheap and simple solution to this problem is to incorporate empirical potentials of the form $-C_6/R^6$, where the $C_6$ coefficients are pairwise atomic parameters. This approach, which is an alternative to the use of a non-local correlation functional, is known as dispersion-corrected DFT (DFT-D).^{400, 395}

There are currently three unique DFT-D methods in Q-Chem. These are requested via the $rem variable DFT_D and are discussed below.

The oldest of these approaches is DFT-D2,³⁹⁹ in which the empirical dispersion potential has the aforementioned form, namely, pairwise atomic $-C/R^6$ terms:

This function is damped at short range, where $R_{AB}^{-6}$ diverges, via

which also helps to avoid double-counting of electron correlation effects, since short- to medium-range correlation is included via the density functional. (The quantity $R_{0,AB}$ is the sum of the van der Waals radii for atoms $A$ and $B$, and $d$ is an additional parameter.) The primary parameters in Eq. (5.23) are atomic coefficients $C_{6,A}$, from which the pairwise parameters in Eq. (5.23) are obtained as geometric means, as is common in classical force fields:

The total energy in DFT-D2 is of course $E_{\text{DFT-D2}}=E_{\text{KS-DFT}}+E_{\mathrm{disp}}^{\text{D2}}$.

DFT-D2 is available in Q-Chem including analytic gradients and frequencies, thanks to the efforts of David Sherrill’s group. The D2 correction can be used with any density functional that is available in Q-Chem, although its use with the non-local correlation functionals discussed in Section 5.7.2 seems inconsistent and is not recommended. The global parameter $s_6$ in Eq. (5.23) was optimized by Grimme for four different functionals,³⁹⁹ and Q-Chem uses these as the default values: $s_6=0.75$ for PBE, $s_6=1.2$ for BLYP, $s_6=1.05$ for BP86, and $s_6=1.05$ for B3LYP. For all other functionals, $s_6=1$ by default. The D2 parameters, including the $C_{6,A}$ coefficients and the atomic van der Waals radii, can be modified using a $empirical_dispersion input section. For example:

$empirical_dispersion
   S6 1.1
   D 10.0
   C6 Ar 4.60 Ne 0.60
   VDW_RADII Ar 1.60 Ne 1.20
$end

Values not specified explicitly default to the values optimized by Grimme.

Note:  1. DFT-D2 is only defined for elements up to Xe. 2. B97-D is an exchange-correlation functional that automatically employs the DFT-D2 dispersion correction when used via METHOD = B97-D.

An alternative to Grimme’s DFT-D2 is the empirical dispersion correction of Chai and Head-Gordon,¹⁸³ which uses the same form as Eq. (5.23) but with a slightly different damping function:

This version is activated by setting DFT_D = EMPIRICAL_CHG, and the damping parameter $a$ is controlled by the keyword DFT_D_A.

Note:  1. DFT-CHG is only defined for elements up to Xe. 2. The $\omega$B97X-D and $\omega$M05-D functionals automatically employ the DFT-CHG dispersion correction when used via METHOD = wB97X-D or wM05-D.

Grimme’s DFT-D3 method³⁹² constitutes an improvement on his D2 approach, and is also available along with analytic first and second derivatives, for any density functional that is available in Q-Chem. The D3 correction includes a potential akin to that in D2 but including atomic $C_8$ terms as well:

The total D3 dispersion correction consists of this plus a three-body term of the Axilrod-Teller-Muto (ATM) triple-dipole variety, so that the total D3 energy is $E_{\text{DFT-D3}}=E_{\text{KS-DFT}}+E_{\text{D3,2-body}}+E_{\text{ATM,3-body}}$

Several versions of DFT-D3 are available as of Q-Chem 5.0, which differ in the choice of the two damping functions. Grimme’s formulation,³⁹² which is now known as the “zero-damping” version [DFT-D3(0)], uses damping functions of the form

for $n=6$ or 8, ${\beta }_6=12$, and ${\beta }_8=14$. The parameters $R_{0,AB}$ come from atomic van der Waals radii, $s_{r,6}$ is a functional-dependent parameter, and $s_{r,8}=1$. Typically $s_6$ is set to unity and $s_8$ is optimized for the functional in question.

The more recent Becke–Johnson-damping version of DFT-D3,³⁹⁴ DFT-D3(BJ), is designed to be finite (but non-zero) as $R_{AB}\to 0$. The damping functions used in DFT-D3(BJ) are

where ${\alpha }_1$ and ${\alpha }_2$ are adjustable parameters fit for each density functional. As in DFT-D3(0), $s_6$ is generally fixed to unity and $s_8$ is optimized for each functional. DFT-D3(BJ) generally outperforms the original DFT-D3(0) version.³⁹⁴

The $C_6$-only (CSO) approach of Schröder et al.¹⁰⁰⁶ discards the $C_8$ term in Eq. (5.27) and uses a damping function with one parameter, ${\alpha }_1$:

The DFT-D3(BJ) approach was re-parameterized by Smith et al.¹⁰⁵¹ to yield the “modified” DFT-D3(BJ) approach, DFT-D3M(BJ), whose parameterization relied heavily on non-equilibrium geometries. The same authors also introduces a modification DFT-D3M(0) of the original zero-damping correction, which introduces one additional parameter (${\alpha }_1$) as compared to DFT-D3(0):

Finally, optimized power approach of Witte et al.¹²²⁶ treats the exponent, ${\beta }_6$, as an optimizable parameter, given by

Note that ${\beta }_8={\beta }_6+2$.

To summarize this bewildering array of D3 damping functions:

• •

  DFT-D3(0) is requested by setting DFT_D = D3_ZERO. The model depends on four scaling parameters ($s_6$, $s_{r,6}$, $s_8$, and $s_{r,8}$), as defined in Eq. (5.28).

• •

  DFT-D3(BJ) is requested by setting DFT_D = D3_BJ. The model depends on four scaling parameters ($s_6$, $s_8$, ${\alpha }_1$, and ${\alpha }_2$), as defined in Eq. (5.29).

• •

  DFT-D3(CSO) is requested by setting DFT_D = D3_CSO. The model depends on two scaling parameters ($s_6$ and ${\alpha }_1$), as defined in Eq. (5.30).

• •

  DFT-D3M(0) is requested by setting DFT_D = D3_ZEROM. The model depends on five scaling parameters ($s_6$, $s_8$, $s_{r,6}$, $s_{r,8}$, and ${\alpha }_1$), as defined in Eq. (5.31).

• •

  DFT-D3M(BJ) is requested by setting DFT_D = D3_BJM. The model depends on four scaling parameters ($s_6$, $s_8$, ${\alpha }_1$, and ${\alpha }_2$), as defined in Eq. (5.29).

• •

  DFT-D3(op) is requested by setting DFT_D = D3_OP. The model depends on four scaling parameters ($s_6$, $s_8$, ${\alpha }_1$, ${\alpha }_2$, and ${\beta }_6$), as defined in Eq. (5.29).

The scaling parameters in these damping functions can be modified using the $rem variables described below. Alternatively, one may simply set DFT_D = D3, and a D3 dispersion correction will be selected automatically, if one is available for the selected functional.

Note:  1. DFT-D3(0) is defined for elements up to Pu ($Z\mathrm{=}\mathrm{94}$). 2. The B97-D3(0), $\omega$B97X-D3, $\omega$M06-D3 functionals automatically employ the DFT-D3(0) dispersion correction when invoked by setting METHOD equal to B97-D3, wB97X-D3, or wM06-D3. 3. The old way of invoking DFT-D3, namely through the use of EMPIRICAL_GRIMME3, is still supported, though its use is discouraged since D3_ZERO accomplishes the same thing but with additional precision for the relevant parameters. 4. When DFT_D = D3, parameters may not be overwritten, with the exception of DFT_D3_3BODY; this is intended as a user-friendly option. This is also the case when EMPIRICAL_GRIMME3 is employed for a functional parameterized in Q-Chem. When any of D3_ZERO, D3_BJ, etc. are chosen, Q-Chem will automatically populate the parameters with their default values, if they available for the desired functional, but these defaults can still be overwritten by the user.

The three-body interaction term, $E^{(3)}$,³⁹² must be explicitly turned on, if desired.

More recently, Grimme published an extended D3 model, D4.^{146, 147, 148} The main feature of D4 is that the coefficients are generated through Casimir-Polder integration of the dynamic atomic polarizabilities $\alpha (i\omega )$ where electronic density information is employed via atomic partial charges. Benchmark results show that the proposed D4 model yields significantly lower error bars. The DFT-D4 dispersion energy similar to D3 model is given by

The Becke-Johnson damping is utilized as default. The coordination number dependent $C_6^{\text{AB}}$ coefficients are obtained on-the-fly via Casimir-Polder integration

where

and

${\alpha }^{A_m{H}_n}$ denotes the reference polarizailities which represents the molecular polarizability of symmetric hydride systems $A_m{H}_n$. $W_{A/B}^{A,\mathrm{ref}/B,\mathrm{ref}}$ are weighting factors determining the contributions of all element specific reference systems $N^{A,\mathrm{ref}/B,\mathrm{ref}}$. $z^{H_A}$ describes the effective charge of hydrogen connectd to atom A in the reference system $A_m{H}_n$. The effective charge $z^A$ is computed self-consistently via Mulliken charge $q^A$,

The coefficients $a$ and $b$ are parametrized to match cationic static polarzizabilies and TD-DFT derived molecular dispersion coefficients, respectively.