5.7.2 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.1; 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).^{329, 324}

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:

(5.21)

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

$$f_{\text{damp}}^{\text{D2}}(R_{AB})={\left[1+e^{-d(R_{AB}/R_{0,AB}-1)}\right]}^{-1}$$

(5.22)

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.21) are atomic coefficients $C_{6,A}$, from which the pairwise parameters in Eq. (5.21) are obtained as geometric means, as is common in classical force fields:

$$C_{6,AB}={\left(C_{6,A}{C}_{6,B}\right)}^{1/2}$$

(5.23)

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.1 seems inconsistent and is not recommended. The global parameter $s_6$ in Eq. (5.21) 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.21) but with a slightly different damping function:

$f_{\text{damp}}^{\text{CHG}}(R_{AB})={\left[1+a{(R_{AB}/R_{0,AB})}^{-12}\right]}^{-1}$

(5.24)

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 more recent 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:

(5.25)

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

$$f_{\text{damp},n}^{\text{D3(0)}}(R_{AB})={\left[1+6{\left(\frac{R_{AB}}{s_{r,n}{R}_{0,AB}}\right)}^{-{\beta }_n}\right]}^{-1}$$

(5.26)

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

$$f_{\text{damp},n}^{\text{D3(BJ)}}(R_{AB})=\frac{R_{AB}^n}{R_{AB}^n+{\left({\alpha }_1{R}_{0,AB}+{\alpha }_2\right)}^n}$$

(5.27)

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.25) and uses a damping function with one parameter, ${\alpha }_1$:

$$f_{\text{damp},6}^{\text{D3(CSO)}}(R_{AB})=\frac{C_{AB}^6}{R_{AB}^6+{(2.5\text{Å})}^6}\left(s_6+\frac{{\alpha }_1}{1+\exp [R_{AB}-(2.5\text{Å})R_{0,AB}]}\right).$$

(5.28)

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

$$f_{\text{damp},n}^{\text{D3M(0)}}(R_{AB})={\left[1+6{\left(\frac{R_{AB}}{s_{r,n}{R}_{0,AB}}+{\alpha }_1{R}_{0,AB}\right)}^{-{\beta }_n}\right]}^{-1}.$$

(5.29)

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

$$f_{\text{damp},n}^{\text{D3(op)}}(R_{AB})=\frac{R_{AB}^{{\beta }_n}}{R_{AB}^{{\beta }_n}+{({\alpha }_1{R}_{0,AB}+{\alpha }_2)}^{{\beta }_n}}.$$

(5.30)

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.26).

• •

  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.27).

• •

  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.28).

• •

  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.29).

• •

  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.27).

• •

  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.27).

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.

$comment
   Geometry optimization, followed by single-point calculations using a larger
   basis set.
$end

$molecule
   0 1
   C       0.000000    -0.000323     1.755803
   H      -0.887097     0.510784     1.390695
   H       0.887097     0.510784     1.390695
   H       0.000000    -1.024959     1.393014
   H       0.000000     0.001084     2.842908
   C       0.000000     0.000323    -1.755803
   H       0.000000    -0.001084    -2.842908
   H      -0.887097    -0.510784    -1.390695
   H       0.887097    -0.510784    -1.390695
   H       0.000000     1.024959    -1.393014
$end

$rem
   JOBTYPE         opt
   EXCHANGE        B3LYP
   BASIS           6-31G*
   DFT_D           D3_ZERO
   DFT_D3_S6       100000
   DFT_D3_RS6      126100
   DFT_D3_S8       170300
   DFT_D3_3BODY    FALSE
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE         sp
   EXCHANGE        B3LYP
   BASIS           6-311++G**
   DFT_D           D3_ZERO
   DFT_D3_S6       100000
   DFT_D3_RS6      126100
   DFT_D3_S8       170300
   DFT_D3_3BODY    FALSE
$end