Understanding DFT Calculations of Weak Interactions: Density-Corrected Density Functional Theory

Abstract

In this work, we discuss where the failure of Kohn-Sham Density Functional Theory (DFT) occurs in weak interactions. We have adopted density-corrected density functional calculations and dispersion correction separately to find out whether the failure is due to density-driven error or functional error. The results of Benzene Ar complex, one of the most common examples of van der Waals interactions, show that DFT calculations of van der Waals interaction suffer from functional error, rather than density-driven error. In addition, errors in DFT calculations of the S22 dataset, which contains small to relatively large (30 atoms) complexes with non-covalent interactions, are governed by functional errors.

INTRODUCTION

Kohn Sham Density Functional Theory (KS-DFT), widely used in quantum chemistry, is known to produce relatively accurate results using local density functionals, but to make large errors with semilocal density functionals.¹⁻⁵ Accurate calculations are important for semilocal DFT because weak intermolecular forces, such as London dispersion force and hydrogen bond, play an important role in determining the structure and energy of chemical species.⁶⁻⁸ Because of its importance, several attempts have been made to correct the error in the semilocal density functionals: (i) Density-Corrected Density Functional Theory (DC-DFT),⁹⁻¹¹ (ii) Density Functional Theory Dispersion Correction (DFTD),¹² (iii) Coulomb-Attenuate method (CAM).¹³ In this study, we use the following methods to determine how to correct DFT calculations of weak intermolecular interactions. At first, DC-DFT is used to see if improved electron density can produce better results. In addition, we also investigate whether DFT-D, developed to ensure the accuracy of KS-DFT, can also work for DC-DFT without further modification. Finally, we employed CAM to examine the effect of longrange interactions, since it is characterized in considering range-separated exchange-correlation potential.

THEORY AND COMPUTATIONAL METHOD

KS-DFT fails to calculate accurate dissociation curves of hetero-diatomic species such as NaCl and CH− : the potential energy curve at the dissociation limit does not converge to zero.¹⁴ This problem is due to the delocalization error of KS-DFT.¹⁵ When Mülliken charge for atoms in a stretched molecule is monitored at the dissociation limit, a fractional charge on the atoms is observed rather than an integer charge. The artificial charge transfer attributes to the delocalization error in KS-DFT and causes errors in the selfconsistent electron density. In recent works, we have shown that the electron density from Hartree-Fock (HF) calculation can correct the delocalization error, i.e., self-interaction induced density error.^16,17 And that the use of HF density instead of DFT approximate density, so-called HFDFT, can well represent the convergence of the binding energy curve at the dissociation limit. Therefore, HF-DFT is shown to agree excellently with CCSD(T) calculations, which is generally used as a reference in quantum calculations.¹⁴

There is another problem that KS-DFT with semilocal density functionals cannot predict London dispersion force accurately. There have been various attempts to solve this problem of dispersion interaction, which is one of the representative functional errors. In this work, a semiclassical correction method is adopted, namely DFT-D.12 DFT-D effectively corrects for the dispersion interaction by adding silico-empirical dispersion energy to the total energy from KS-DFT. In particular, the recent version, DFT-D3, has been amended from its predecessors. In this work, therefore, we employed the D3 for all dispersion corrected calculations.¹⁸ DFT-D3 has the following advantages over its predecessors:^19,20

(1) It is less empirical than previous works, the most important parameters are computed from first principles by standard KS-DFT.

(2) The approach is asymptotically correct with all density functionals for finite systems (molecules) or nonmetallic infinite systems. It gives nearly accurate dispersion energy for a gas of weakly interacting neutral atoms and smoothly interpolates to molecular (bulk) regions.

(3) It provides a consistent description of all chemically relevant elements of the periodic system (nuclear charge Z=1-94).

(4) Atom pair-specific dispersion coefficients and cutoff radii are explicitly computed.

(5) Coordination number (geometry) dependent dispersion coefficients are used that do not rely on atom connectivity information (differentiable energy expression).

(6) It offers similar or better accuracy for “light” molecules and has greatly improved the description of metallic and “heavier” systems.

The energy of DFT-D3 is

\(E_{DFT-D3} = E_{ KS-DFT} + E_{disp}\)       (1)

where EKS-DFT is the energy calculated from KS-DFT. Edisp is originated from London dispersion forces and decomposed into two- (E⁽²⁾) and three-body (E⁽³⁾) terms as

\(E_{disp} = E^{(2)} + E^{(3)}\)       (2)

The two-body energy is

\(E^{(2)} = \sum_{AB_n = 6, 8, 10, \cdots } \sum - s_n \frac{ C^{AB}_n}{r^n_{AB}} f_{d,n} (r_{AB})\)       (3)

where C_n^AB is the n^th-order dispersion coefficient for consisting atoms. The global scaling factor, sn is parameterized for the datasets and is functional dependent. The damping function, f_d,n is given as Eq. (4).

\(f_{d,n} (r_{AB} ) = \frac{ 1 }{ 1+6( {r_{AB}\over{ (s_{r,n} R^{AB}_0)}} )^{- \alpha _0} }\)       (4)

where the cutoff radii, R0 AB reduces the cost of numerous calculations.¹⁴ αn represents the degree of decaying of London dispersion forces and is defined by the following recurrence formula:

\(\alpha_{n+2} = \alpha_n + 2 \), \(\alpha_6 = 14\)       (5)

Another way to consider weak interactions is to apply a long-range correction approximation. CAM-B3LYP functional can correct what B3LYP poorly predicts; the polarizability of the long chain, excitations using time-dependent DFT, and charge transfer excitations.²¹⁻²⁵ In general, the reason why standard DFT cannot predict the long-range interaction is that, at the long-distance limit where the exact exchange potential is given by −0.2r⁻¹ , the approximate exchange potential of DFT is expressed as −r⁻¹.¹³ To solve this problem a range-separated approximation divides the exchange potential into two parts: short-range or long-range.²⁶ In such range-separated exchange functionals, the exchange potential at short range uses DFT exchange, while the HF exchange is used at long-distance.²⁷⁻³⁰

We calculated Benzene·Ar complex (Scheme 1), one of the simplest van der Waals (or London dispersion) systems, and the S22 dataset,³¹ which includes hydrogen bonds and London dispersion interactions. All calculations are performed self-consistently (indicated by SC-DFT) and DCDFT as well as employing dispersion correction, respectively, i.e., DFT-D3 and DC-DFT-D3. Density functionals, used in this research, are Perdew-Burke-Ernzerhof (PBE),³² Becke Three parameter and Lee-Yang-Parr (B3LYP),³³ and TaoPerdew-Staroverov-Scuseria (TPSS),³⁴ which are widely adapted functionals from generalized gradient approximation (GGA), hybrid and meta-GGA, respectively, in addition to Coulomb-Attenuating Method-Becke Three parameter and Lee-Yang-Parr (CAM-B3LYP).¹³

RESULTS AND DISCUSSION

The results of Benzene·Ar complex are shown in Fig. 1. According to Fig. 1(a), there is no potential well in KSDFT curve, indicating that KS-DFT cannot estimate the binding state of the complex. Moreover, the result with B3LYP functional is even repulsive in the binding region. This trend may be due to the exponential decay of the exchange-correlation potential energy of DFT, because London dispersion forces depend on R^-6 terms in real systems. As with other density functionals, CAM-B3LYP, which has rangeseparated exchange potential, does not predict accurate weak interactions either. Therefore, for this system, the rangeseparation of exchange-correlation potential is not effective to reduce the dispersion functional error. The results of Fig. 2(a), which were performed by the density-correction calculation, show the same tendency as in Fig. 1(a). Thus, the failure of DFT in vdW interactions is due to reasons other than inaccurate density or long-range exchange-correlation potential.

Figure 1. Dissociation curve of Benzene·Ar complex using (a) KS-DFT and (b) DFT-D3.

Figure 2. Dissociation curve of Benzene·Ar complex using (a) HF-DFT and (b) HF-DFT-D3.​​​​​​​

Interestingly, the calculation results using DFT-D3 shows the binding state, as shown in Fig. 1(b). The experimentally known distance between the center of benzene and argon atom is 3.592Å.^35,36 In Table 1, DFT-D3 has an equilibrium distance close to the experimental value, regardless of the approximate functional. There is no noticeable difference between Figs. 1(b) and 2(b), showing that DFTD3 methods which are parametrized for self-consistent DFT calculations, is applicable to HF-DFT. This result indicates that dispersion correction is required for vdW interactions whether we use DFT or DC-DFT. That is, standard DFT approximations have large functional errors due to missing dispersion interactions and are not density sensitive.

Table 1. Equilibrium distance and binding energy using DFT-D3​​​​​​​

DFT-D3 is parameterized for KS-DFT, but the correction effect on HF-DFT is remarkable for Benzene·Ar complex. We extend our analysis to the S22 dataset consisting of various long-range interactions such as (1) Hydrogen-bonded complexes, (2) complexes with predominant dispersion contribution, and (3) mixed complexes. Among them, the second category is the most corrected case with DFT-D3. For instance, the error of adenine-thymine stacked (Table 2) system from the calculation without using DFT-D3 is 10~13 kcal/mol, but the correction by DFT-D3 reduces the error to 1 kcal/mol or less, which is accepted as chemical accuracy. On the other hand, in the first category, CAM-B3LYP results are within 1 kcal/mol for several systems.

Table 2. Errors of energies for model complexes (S22 dataset³¹) (A: B3LYP, B: CAM-B3LYP, ND3 means result is without DFT-D3.) Reference Energies are from CCSD(T)/Complete Basis Set(CBS)³⁷

Throughout the entire dataset, CAM-B3LYP results show that the error is reduced when compared with B3LYP (in the ΔE_ND3^SCF-A and ΔE_ND3^SCF-B column in Table 2). However, it seems that the correction is not enough. Fortunately, by adding a D3 correction, mean absolute error (MAE) of B3LYP and CAM-B3LYP reduces to within chemical accuracy. This is due to the fact that the S22 dataset has been included when the silico-empirical parameters of D3 are optimized. Density-corrected DFT energies were not considered in the parameter optimization process, nevertheless, D3 works on HF-DFT and performs even better on HF-CAM-B3LYP.

These results of Benzene·Ar complex and the S22 dataset demonstrate that D3 correction is applicable to HF-DFT without additional parameter optimization for these vdW and hydrogen bond systems.

CONCLUSION

Although various calculation methods have been used, benzene and argon atomic bonds cannot be predicted when using conventional DFT and DC-DFT. Nevertheless, D3 corrected DFT or HF-DFT produces a relatively accurate binding energy compared to the reference CCSD (T) value. In addition, for the S22 dataset including hydrogen bonding and London dispersion, DFT-D3 or HF-DFT-D3 reduces errors to chemical accuracy. Therefore, for weak dispersion such as vdW interaction, optimized parameters of DFTD3 for KS-DFT can be used for HF-DFT without additional optimization. Furthermore, based on these results, it can be seen that the self-consistent DFT density-error of the S22 dataset is relatively smaller than the functional error. Future studies should investigate molecular systems that include stronger dispersive forces or dipole attraction that are not included in the S22 dataset to better understand weak interactions. In particular, extensive studies using DC-DFT are needed to investigate molecular systems with high density-sensitivity.

COMPUTATION DETAILS

We used Ahlrichs’s def2-QZVP basis set for all calculations because it provides moderate computational cost and accuracy and DFT-D3 was parameterized with def2- QZVP basis set. TURBOMOLE 7.0.2^38,39 was used for PBE, B3LYP, and TPSS but Gaussian 16⁴⁰ was used for CAMB3LYP functional. The dissociation curve of Benzene·Ar complex was obtained by single point calculation at the interval of 0.1Å from 3 to 5Å.