Article pubs.acs.org/JPCA
Dispersion Correction Derived from First Principles for Density Functional Theory and Hartree−Fock Theory Emilie B. Guidez and Mark S. Gordon* Department of Chemistry, Iowa State University, Ames, Iowa 50011, United States ABSTRACT: The modeling of dispersion interactions in density functional theory (DFT) is commonly performed using an energy correction that involves empirically fitted parameters for all atom pairs of the system investigated. In this study, the first-principles-derived dispersion energy from the effective fragment potential (EFP) method is implemented for the density functional theory (DFT-D(EFP)) and Hartree−Fock (HF-D(EFP)) energies. Overall, DFT-D(EFP) performs similarly to the semiempirical DFT-D corrections for the test cases investigated in this work. HF-D(EFP) tends to underestimate binding energies and overestimate intermolecular equilibrium distances, relative to coupled cluster theory, most likely due to incomplete accounting for electron correlation. Overall, this firstprinciples dispersion correction yields results that are in good agreement with coupledcluster calculations at a low computational cost.
■
developed by Grimme et al. The expression for the R−6 term in DFT-Dn is given by9
INTRODUCTION Intermolecular dispersion forces arise from the interaction between induced multipoles.1 These forces are at the origin of many chemical and biological processes such as protein folding,2−5 molecular recognition,6 and DNA base pair stacking.7 The modeling of dispersion interactions has been the subject of many investigations.8−11 Density functional theory (DFT) is frequently used to model molecular systems that contain on the order of hundreds of atoms.12,13 However, most commonly used density functionals (as well as Hartree− Fock (HF) theory) cannot account for dispersion interactions. Certain density functionals such as MPWB1K,14 M06-2X,15 M08-HX,16 M08-SO,16 and M1117 developed by Truhlar and co-workers correctly capture attractive noncovalent interactions where intermolecular overlap is nonnegligible (at the van der Waals minima). In addition, other functionals such as the van der Waals nonlocal correlation functionals (vdW-DF)18−22 include dispersion. In order to enable popular density functionals (GGAs, hybrids, ...) to account for dispersive interactions, Grimme and co-workers introduced a series of empirical corrections,23−25 collectively referred to here as “-D”. In DFT-D, the -D dispersion interaction energy correction is added to the Kohn−Sham energy. In general, the dispersion interaction between two molecules A and B can be expressed as26 C C C Edisp = 66 + 77 + 88 + ... (1) R R R R is the intermolecular distance and Cn are coefficients. The terms with an odd power of R are orientation-dependent,27,28 and average to zero in freely rotating systems. The dispersion energy expression is often truncated at the R−6 term, which corresponds to the interaction between induced dipoles of two species A and B. Three empirically parametrized dispersion corrections, called DFT-Dn (n = 1, 2, 3),23−25 have been © 2015 American Chemical Society
Nat − 1
Edisp = −S6
Nat
∑ ∑ I=1 J=I+1
C6IJ RIJ 6
fdmp (RIJ )
(2)
S6 is a scaling parameter that depends on the functional. CIJ6 , the dispersion coefficient for the atom pair IJ, is a fitted parameter, Nat is the number of atoms, and RIJ is the distance between atoms I and J. The repulsive interaction between the nuclei at very small RIJ is taken into account by using the damping function fdmp. In the DFT-Dn implementations by Grimme, the damping function also involves a fitted parameter. While DFT-Dn (n = 1, 2) only include the R−6 term,23,24 DFT-D3 also includes the R−8 term in a recursive manner.25 In DFT-D1 and DFT-D2, dispersion coefficients are predetermined using atomic ionization potentials and static polarizabilitites. There is no dependence on the atomic environment, which can lead to substantial errors. On the other hand, DFT-D3 considers the effective volume of the atoms by taking into account the number of neighbors in the environment to calculate the dispersion coefficients. Other methods such as the Tkatchenko−Scheffler29−31 and the Becke−Johnson32−34 models also consider the chemical environment of the atoms to calculate dispersion coefficients. The former method relies on the Hirshfeld partitioning of the total electron density between the atoms (which is obtained from electronic structure calculations) to compute the dispersion coefficients.29 The latter computes the dipole moment generated by the exchangecorrelation hole.34 Both methods are minimally empirical. The dDsC method, based on the generalized gradient approximation of the Becke−Johnson model, was developed later on.35,36 Received: January 13, 2015 Published: February 4, 2015 2161
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
α̅ represents one-third of the trace of the dynamic polarizability tensor. The use of the isotropic approximation is justified since anisotropic effects are minimal in this distributed approach.41 In addition, it permits the direct comparison of the C6 coefficients with experimental and other theoretical data and reduces computational cost.41 The integral in eq 7 is evaluated using a 12-point Gauss−Legendre quadrature formula.47 With a change of variable
One drawback of DFT-D can be the double counting of some of the electron correlation that is already included in the correlation functional. The parametrization attempts to minimize the effect of double counting. On the other hand, the Hartree−Fock method does not include any electron correlation except for the Fermi hole.37 Therefore, there is no possibility of double counting when one adds the -D correction to HF theory. However, other correlation effects are not included in HF-D. HF-D3 calculations using the Grimme corrections have been performed on several systems.38 Both DFT and HF have similar computational costs, which is one reason for the popularity of DFT. The effective fragment potential (EFP) method has been developed to treat intermolecular interactions.39,40 In the EFP method, all interaction energy terms, including the dispersion energy, are derived from first principles.41 In the EFP method, the interaction energy between the molecules or fragments is divided into five contributions: the Coulomb, polarization, exchange−repulsion, charge transfer, and dispersion interaction energies:40,42,43 E = Ecoul + Epol + Eex−rep + Ect + Edisp
ω = ω0
1+t , 1−t
dω =
−2ω0 dt (1 + t )2
(8)
Equation 7 can be rewritten as AB Edisp =−
3ℏ π
12
∑ ∑ ∑ Wk i∈A j∈B k=1
α̅ i(iω) α̅ j(iω) (1 − tk) R ij 6 2ω0
2
(9)
Wk and tk are the Gauss−Legendre weighting factor and abscissa.41 ω0 has an optimal value of 0.3.48 In order to avoid a singularity near R = 0, each term in eq 9 is multiplied by a damping function.49 An overlap-based damping function is used here: 6
(3)
fdmp (i , j) = 1 − Sij
2
∑ n=0
The objective of this work is to add the dispersion energy derived from first principles in the EFP method to the DFT Kohn−Sham energy (DFT-D(EFP)) and Hartree−Fock energy (HF-D(EFP)).
( −2 ln|Sij|)n /2 n!
(10)
Sij in eq 10 is the overlap between the localized molecular orbitals i and j. This damping function differs slightly from the one used originally in EFP49 in that it is now equally applicable to even and odd powers of n. An overlap-based damping function is a logical choice since the intermolecular overlap integral depends on the intermolecular distance, so the dispersion interactions are diminished as the overlap increases. In addition, unlike other damping functions, the overlap damping function contains no empirically fitted parameters.49 The present approach is derived from first principles and does not involve any empirically fitted parameters. The Tkatchenko− Scheffler method determines atomic polarizabilities based on the partitioning of the electron density as well as the effective atomic volumes. In contrast, the method described here solves the time-dependent Hartree−Fock equations using localized
■
METHOD In the DFT-D(EFP) and HF-D(EFP) implementations, each molecule represents a fragment within the EFP scheme. The dispersion energy between the molecules (fragments) is calculated as it would be in an EFP calculation. The dispersion interaction energy calculated in this manner is then added to the quantum mechanical (QM) energy of the system (Kohn−Sham or Hartree−Fock). The total energy of the system is given by E = Edisp + EQM (4) The expression for the dispersion energy derived from Rayleigh−Schrödinger perturbation theory is1,41,44 AB =− Edisp
ℏ 2π
x ,y,z
∑ ∑ ∑ Tαβij Tγδij ∫
0
i ∈ A j ∈ B αβγδ
∞
i ααγ (iω) αβδj (iω) dω
(5)
In eq 5, i and j label the localized molecular orbitals (LMOs) of molecules A and B, respectively.45 α(iω) represents the dynamic dipole polarizability tensor. The dynamic polarizability tensor is calculated by solving the time-dependent Hartree− Fock equations.41,46 T is a second order electrostatic tensor given by1,41,44 ij Tαβ
3R ijαR ijβ − R ij 2δαβ 1 1 = ∇α ∇β = R ij 4πε0R ij 4πε0R ij 5
(6)
Rij = Ri − Rj, where Ri and Rj are the coordinates of the LMO centroids i and j of fragments A and B, respectively. One can 1 ij ij 6 substitute ∑x,y,z αβγσ TαβTγσ = 6/Rij into eq 5. Within the isotropic 41 approximation, eq 5 further reduces to
Figure 1. Systems with a dominant dispersion interaction (first row): (A) benzene (sandwich); (B) methane dimer; (C) hydrogen dimer; (D) π-stacked adenine−thymine. Hydrogen-bonded systems (second row): (E) water dimer; (F) ammonia dimer; (G) methanol dimer; (H) adenine−thymine base pair (WC). Mixed systems (third row): (I) ethene−ethyne; (J) T-shape benzene dimer; (K) benzene−water dimer; (L) benzene−ammonia dimer. Color coding: black = carbon; white = hydrogen; red = oxygen; blue = nitrogen.
∞
AB Edisp
3ℏ =− π
∑∑ i∈A j∈B
∫0 α̅ i(iω) α̅ j(iω) dω R ij 6
(7) 2162
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 2. Potential energy curves of a set of test dimer cases with dominant dispersion interactions. (A) Benzene sandwich; (B) methane dimer; (C) hydrogen dimer with a 6-311++G(d,p) basis set; (D) hydrogen dimer with a 6-311++G(3df,3p) basis set; (E) adenine−thymine DNA stack. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set except for (D), where a 6-311++G(3df,3p) basis set is used. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations except for (D) where a 6-311++G(3df,3p) basis set is used.
molecular orbitals to calculate molecular polarizabilities. The DFT-D(EFP) binding energy between two molecules is calculated using the formula DFT‐D(EFP) ΔE = Edimer −
DFT ∑ Emonomer
was used with the 6-311++G(d,p) basis set. The HF-D(EFP) calculations were performed with the same basis set. For the hydrogen dimer, an additional set of calculations was done with the 6-311++G(3df,3p) basis set. The dimerization energies (eqs 11 and 12) were calculated at various intermolecular distances to generate a potential energy surface for each dimer. The internal monomer geometries were held fixed in the dimer calculations. All potential energy surfaces were generated by moving the molecules along the axis that connects their centers of mass at the equilibrium dimer orientation. The intermolecular equilibrium distance R0 and the equilibrium binding energy ΔE0, corresponding to the minimum on the potential energy surface, are reported for all methods. The equilibrium geometries of the benzene dimers (sandwich and T-shape) were optimized at the CCSD(T)/aug-cc-pVQZ* level of theory with frozen monomers.54 Equilibrium geometries of the methane, hydrogen, water, ammonia, and methanol dimers were taken from ref 49. The equilibrium geometries for the Watson−Crick DNA pair adenine−thymine (AD-WC) and the adenine−thymine stack
(11)
An analogous formula is used for HF-D(EFP) calculations: HF‐D(EFP) ΔE = Edimer −
HF ∑ Emonomer
(12)
In this work, the potential energy surfaces (PESs) of several test dimers are generated and compared with the PESs obtained with the Grimme DFT-Dn methods and with CCSD(T) and MP2 calculations. Note that, for DFT-Dn calculations, intramolecular dispersion forces are included in the energy of the monomer but they are not included in this new DFT-D(EFP) method.
■
COMPUTATIONAL DETAILS Calculations reported in this work were performed using the GAMESS software package.50,51 The DFT functional B3LYP52,53 2163
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
2164
a The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized MP2/cc-pVTZ CP geometry, except for the methane dimer which was optimized at the CCSD(T)/cc-pVTZ level of theory.56 cThese values are obtained at the CCSD(T)/aug-cc-pVQZ* level of theory from ref 54. dThese values are from the CCSD(T)/CBS potential energy surface of the T-shape hydrogen dimer by Johnson and Diep.57
(−3.08) (−0.38) (−0.052) (−0.087) 3.88 4.08 3.68 3.44 − (−4.55) (−0.34) (−0.042) (−0.073) (−18.21) benzene (sandwich) (CH4)2 (H2)2 (6-311++G(d,p)) (H2)2 (6-311++G(3df,3p)) A−T (stack)
4.02 3.90 3.16 3.16 3.28
(−0.55) (−0.21) (−0.049) (−0.055) (−10.45)
3.82 3.76 3.22 3.20 3.10
(−1.21) (−0.38) (−0.035) (−0.040) (−13.97)
3.86 3.96 3.30 3.30 3.22
(−1.66) (−0.38) (−0.068) (−0.074) (−12.51)
3.82 (−1.11) 3.88 (−0.24) N/A 3.30 (−0.027) 3.08 (−12.27)
4.24 4.10 3.92 3.60 3.26
(−0.60) (−0.28) (−0.022) (−0.060) (−10.98)
3.68 4.12 3.72 3.50 3.10
MP2 HF-D(EFP) DFT-D(EFP) DFT-D3 DFT-D2 DFT-D1
RESULTS AND DISCUSSION The three classes of systems that were investigated are shown in Figure 1. The first class corresponds to systems with dominant dispersion interactions, which include the sandwich configuration of the benzene dimer, the methane dimer, the hydrogen dimer, and the π-stacked adenine−thymine (Figure 1A−D). The second class corresponds to hydrogen-bonded systems, which include the water dimer, the ammonia dimer, the methanol dimer, and the adenine−thymine Watson−Crick DNA base pair (Figure 1E−H). The last class of complexes investigated is the mixed systems: the ethene−ethyne complex, the T-shape benzene dimer, and the benzene−water and benzene−ammonia complexes (Figure 1I−L). This classification was taken from ref 9. The potential energy surfaces of the first class of systems are displayed in Figure 2. The intermolecular equilibrium distances R0 and equilibrium binding energies ΔE0 of all of these systems are shown in Table 1. The benzene dimer is first investigated in a sandwich configuration. Like the semiempirical DFT-Dn methods (n = 1, 2, 3), DFT-D(EFP) and HF-D(EFP) both tend to underestimate the binding energy of the benzene sandwich in comparison to CCSD(T). MP2 overestimates the binding energy for this system. DFT-D(EFP) gives distances and binding energies that are similar to those of both DFT-D2 and DFT-D3, with an intermolecular equilibrium distance R0 of 3.82 Å and a binding energy ΔE0 of −1.11 kcal/mol. HF-D(EFP) overestimates the intermolecular equilibrium distance and underestimates the binding energy compared to CCSD(T), DFT-D2, DFT-D3, and DFT-D(EFP). In fact, the HF-D(EFP) results are very similar to those predicted by DFT-D1. For the methane dimer, all of the DFT-Dn methods underestimate the intermolecular equilibrium distance, while the HF-D(EFP) intermolecular distance is in excellent agreement with the CCSD(T) value. The binding energy at R0 for all of the -D methods is within 0.2 kcal/mol of the CCSD(T) value. The CCSD(T) interaction between two hydrogen molecules is very weak (smaller than 0.1 kcal/mol); it is therefore a challenge to model the PES accurately. As may be seen in Table 1, two basis sets were used for the H2 dimer: 6-311++G(d,p) and 6-311++G(3df,3p). The larger basis set increases the MP2 and CCSD(T) binding energies by ∼0.03 kcal/mol. The effect on DFT-Dn is small because the dependence on the basis set of the empirically determined -D methods is minimal. This is not the case for the -D(EFP) method, since it arises from first principles and therefore has an explicit basis set dependence. For the smaller basis set, DFT-D(EFP) does not have a minimum on the H2 dimer potential energy surface, while the HF-D(EFP) binding energy is smaller than that predicted by CCSD(T) by ∼0.03 kcal/mol. The binding energies predicted by both -D(EFP) methods are similar to those predicted by MP2 and CCSD(T). The π-stacking interaction between the DNA bases adenine and thymine (A−T) for the methods considered here are
dimer
Table 1. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Systems with Dominant Dispersion Contributionsa
■
CCSD(T)
CCSD(T)/CBSb
optimized at the MP2/cc-pVTZ level of theory with counterpoise correction (CP) were taken from ref 55. The equilibrium geometries of the benzene−water and benzene−ammonia complexes optimized at the MP2/cc-pVTZ level of theory with counterpoise correction were also taken from ref 55. The ethene−ethyne complex optimized at the CCSD(T)/cc-pVQZ level of theory was adapted from ref 55. The CCSD(T)/CBS energies were taken from refs 56, 54, and 57.
3.9c (−1.70) 3.71 (−0.53) 3.36d (−0.11) 3.36d (−0.11) 3.17 (−11.66)
The Journal of Physical Chemistry A
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 3. Potential energy curves of a set of hydrogen-bonded complexes. (A) Water dimer; (B) ammonia dimer; (C) methanol dimer; (D) adenine−thymine Watson−Crick pair. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations.
Table 2. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Hydrogen-Bonded Systemsa dimer (H2O)2 (NH3)2 (CH3OH)2 A···T (WC)
DFT-D1 2.90 3.28 3.58 5.96
(−6.25) (−3.89) (−6.44) (−16.55)
DFT-D2 2.88 3.28 3.54 5.94
(−6.48) (−4.19) (−6.88) (−18.17)
DFT-D3 2.90 3.28 3.56 5.96
(−6.41) (−4.06) (−6.82) (−17.90)
DFT-D(EFP) 2.84 3.20 3.50 5.84
(−6.52) (−4.24) (−6.94) (−20.11)
HF-D(EFP) 2.94 3.36 3.60 5.92
(−5.50) (−3.13) (−5.90) (−17.07)
MP2 2.92 3.32 3.58 5.98
(−5.93) (−3.75) (−6.47) (−16.03)
CCSD(T)
CCSD(T)/CBSb
2.92 (−5.85) 3.34 (−3.67) 3.58 (−6.46) −
2.91 (−5.02) 3.21 (−3.17) − 5.97 (−16.74)
a The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized CCSD(T)/cc-pVQZ geometry, except for the adenine−thymine DNA pair which was optimized with MP2/cc-pVTZ CP.56
only slightly smaller than the CCSD(T) value of 2.92 Å. The intermolecular equilibrium distance of 2.94 Å obtained with HF-D(EFP) is very close to the CCSD(T) and MP2 values. Similar trends are observed for the ammonia and methanol dimers, for which the DFT-D(EFP) method underestimates the intermolecular equilibrium distance and overestimates the binding energy. HF-D(EFP) overestimates R0 by only 0.02 Å for these two systems in comparison to CCSD(T) but underestimates the binding energy by about 0.55 kcal/mol. The adenine−thymine Watson−Crick pair contains two hydrogen bonds. The PESs obtained with all DFT-D methods for this complex are similar to the ones obtained with MP2. The values of R0 obtained with DFT-Dn (n = 1, 2, 3), and with DFT-D(EFP) and HF-D(EFP) are in good agreement with the MP2/cc-pVTZ CP optimized value. The DFT-Dn and HF-D(EFP) binding energies are all within 2 kcal/mol of the CCSD(T)/CBS value, while the DFT-D(EFP) method overestimates the binding energy by nearly 3.4 kcal/mol. The last class of systems investigated is the mixed systems. The PESs are shown in Figure 4. The values of R0 and ΔE0 are displayed in Table 3.
summarized in Table 1. CCSD(T) calculations were not performed for this system, but the value of R0 optimized at the MP2/cc-pVTZ CP level of theory55 and the value of ΔE0 obtained at the CCSD(T)/CBS level56 are reported in Table 1. MP2 greatly overestimates the binding energy of A−T: 18.2 kcal/mol vs 11.66 kcal/mol for CCSD(T)/CBS. In contrast, all of the DFT-D methods predict binding energies that are within 2 kcal/mol of the CCSD(T)/CBS value. Specifically, DFT-D(EFP) predicts a binding energy that is 0.61 kcal/mol higher than the CCSD(T)/CBS value. The HF-D(EFP) energy is just 0.68 kcal/mol smaller than the CCSD(T)/CBS value. All of the -D methods predict intermolecular distances that are in good agreement with the MP2/cc-pVTZ CP optimized value. The second class of systems is the hydrogen-bonded complexes. The PESs of these systems are shown in Figure 3, and the values of R0 and ΔE0 are reported in Table 2. For the water dimer, all DFT-D methods overestimate the binding energy in comparison to CCSD(T) by 0.5−0.7 kcal/mol, while the HF-D(EFP) binding energy is only 0.35 kcal/mol too low. The intermolecular equilibrium distances obtained with the DFT-D methods are between 2.84 and 2.90 Å, which is 2165
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 4. Potential energy curves of a set of mixed complexes. (A) Ethene−ethyne complex; (B) T-shape benzene dimer; (C) water−benzene complex; (D) ammonia−benzene complex. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations.
Table 3. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Mixed Systemsa dimer ethene−ethyne benzene (T-shape) benzene−H2O benzene−NH3
DFT-D1 4.44 5.06 3.42 3.62
(−1.45) (−2.10) (−3.31) (−1.98)
DFT-D2 4.30 4.86 3.28 3.42
(−1.85) (−3.25) (−4.43) (−2.96)
DFT-D3 4.40 4.98 3.36 3.54
(−1.81) (−3.14) (−4.21) (−2.81)
DFT-D(EFP)
HF-D(EFP)
4.36 4.88 3.32 3.50
4.60 5.10 3.52 3.74
(−1.45) (−2.51) (−3.26) (−2.04)
(−1.23) (−2.16) (−2.93) (−1.80)
MP2 4.44 4.80 3.32 3.48
(−1.68) (−5.35) (−4.34) (−3.42)
CCSD(T) 4.48 4.88 3.36 3.52
(−1.53) (−4.38) (−4.04) (−3.09)
CCSD(T)/CBSb 4.42 (−1.51) 5.0c (−2.61) 3.38 (−3.29) 3.56 (−2.32)
a
The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized MP2/cc-pVTZ geometry, except for the ethene−ethyne complex which was optimized at the CCSD(T)/cc-pVQZ level of theory.56 c These values are obtained at the CCSD(T)/aug-cc-pVQZ* level of theory from ref 54.
by ∼2 kcal/mol. The DFT-D2 and DFT-D3 methods underestimate the binding energy by ∼1 kcal/mol, while the DFT-D1 method is very similar to HF-D(EFP). For the water−benzene and ammonia−benzene complexes, DFT-D(EFP) slightly underestimates the intermolecular equilibrium distances and underestimates the equilibrium binding energies. On the other hand, HF-D(EFP) overestimates R0 for both the benzene−ammonia complex and the water− benzene complex. Overall, DFT-D(EFP) and HF-D(EFP) yield PESs that are similar to those of DFT-D1, while the DFT-D3 method is in better agreement with CCSD(T). In general, for the mixed species, the -D(EFP) method tends to underestimate binding energies slightly more than the DFT-D3 method, but the differences are typically 1 kcal/mol or less.
For the ethene−ethyne complex, the R0 values predicted by all of the DFT-D methods are smaller than those predicted by CCSD(T). DFT-D(EFP) and HF-D(EFP) slightly underestimate the binding energy (by 0.08 and 0.3 kcal/mol respectively), while the DFT-D2 and DFT-D3 methods overestimate the binding energy by 0.3 kcal/mol. For this system, the DFT-D1 method is in very good agreement with CCSD(T). For the T-shape benzene dimer, all DFT-Dn methods as well as HF-D(EFP) and DFT-D(EFP) underestimate the binding energy in comparison to CCSD(T). MP2 on the other hand overestimates the interaction between the monomers. DFT-D(EFP) predicts an equilibrium intermolecular distance that is identical to that of CCSD(T) and a binding energy that is nearly 2 kcal/mol too small. HF-D(EFP) overestimates R0 and also underestimates the equilibrium binding energy 2166
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
■
Models: The Benzene Dimer Case. J. Phys. Chem. A 2014, 118, 9561− 9567. (9) Antony, J.; Grimme, S. Density Functional Theory Including Dispersion Corrections for Intermolecular Interactions in a Large Benchmark Set of Biologically Relevant Molecules. Phys. Chem. Chem. Phys. 2006, 8, 5287−5293. (10) Grimme, S. Density Functional Theory with London Dispersion Corrections. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 211−228. (11) Klimeš, J.; Michaelides, A. Perspective: Advances and Challenges in Treating Van Der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (12) Levine, I. N. Quantum Chemistry; Pearson Prentice Hall: Upper Saddle River, NJ, 2009; Vol. 6. (13) Dreizler, R. M.; Engel, E. Density Functional Theory; Springer: Berlin, 2011. (14) Zhao, Y.; Truhlar, D. G. Hybrid Meta Density Functional Theory Methods for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions: The MPW1B95 and MPWB1K Models and Comparative Assessments for Hydrogen Bonding and Van Der Waals Interactions. J. Phys. Chem. A 2004, 108, 6908−6918. (15) Zhao, Y.; Truhlar, D. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (16) Zhao, Y.; Truhlar, D. G. Exploring the Limit of Accuracy of the Global Hybrid Meta Density Functional for Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2008, 4, 1849−1868. (17) Peverati, R.; Truhlar, D. G. Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810−2817. (18) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van Der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (19) Vydrov, O. A.; Van Voorhis, T. Nonlocal Van Der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103, 063004. (20) Vydrov, O. A.; Van Voorhis, T. Improving the Accuracy of the Nonlocal Van Der Waals Density Functional with Minimal Empiricism. J. Chem. Phys. 2009, 130, 104105. (21) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. Van Der Waals Density Functional: Self-Consistent Potential and the Nature of the Van Der Waals Bond. Phys. Rev. B 2007, 76, 125112. (22) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy Van Der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (23) Grimme, S. Accurate Description of Van Der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (24) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (25) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (26) London, F. The General Theory of Molecular Forces. Trans. Faraday Soc. 1937, 33, 8b−26. (27) Xu, P.; Zahariev, F.; Gordon, M. S. The R−7 Dispersion Interaction in the General Effective Fragment Potential Method. J. Chem. Theory Comput. 2014, 10, 1576−1587. (28) Buckingham, A. D. Theory of Long-Range Dispersion Forces. Discuss. Faraday Soc. 1965, 40, 232−238. (29) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and FreeAtom Reference Data. Phys. Rev. Lett. 2009, 102, 073005.
CONCLUSIONS The main conclusion of the present work is that the DFT-D(EFP) and HF-D(EFP) methods predict potential energy surfaces for a variety of types of intermolecular complexes with an accuracy that is comparable to that of the -D methods of Grimme and co-workers. The cost of computing the dispersion correction with this new method is slightly higher than the empirical -D methods but still small in comparison to the DFT part of the calculation. The importance of the -D(EFP) approach is that the dispersion correction is derived from first principles, without the need for empirically fitted parameters. In principle, this means that the EFP-based -D method is more broadly applicable. In the -D(EFP) methods, the intermolecular dispersion energy is calculated using the approach of the effective fragment potential method. The DFT-D(EFP) and HF-D(EFP) methods were applied to multiple test sets: complexes with dominant dispersion interactions, hydrogen-bonded complexes, and mixed complexes. For all of these systems, DFT-D(EFP) performs similarly to the semiempirical DFT-Dn (n = 1, 2, 3) corrections by Grimme. The HF-D(EFP) method also performs surprisingly well. Unlike the DFT methods, HF contains no other source of electron correlation (no correlation functional). Therefore, one cannot expect the HF-D method to perform consistently as well as DFT-D. Nonetheless, the HF-D method is appealing as a low-cost approximate alternative to MP2. The dispersion correction in DFT-D(EFP) and HF-D(EFP) only contains the R−6 term. Results could potentially be improved by including higher order terms.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Tel.: 515-294-0452. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by a U.S. National Science Foundation Software Infrastructure (SI2) grant, ACI-1047772.
■
REFERENCES
(1) Stone, A. The Theory of Intermolecular Forces, 2nd ed.; Oxford University Press: Oxford, 2013. (2) Vondrasek, J.; Kubar, T.; Jenney, F. E.; Adams, M. W. W.; Kozisek, M.; Cerny, J.; Sklenar, V.; Hobza, P. Dispersion Interactions Govern the Strong Thermal Stability of a Protein. Chem.Eur. J. 2007, 13, 9022−9027. (3) Serrano, A. L.; Waegele, M. M.; Gai, F. Spectroscopic Studies of Protein Folding: Linear and Nonlinear Methods. Protein Sci. 2012, 21, 157−170. (4) Ruiz-Blanco, Y. B.; Marrero-Ponce, Y.; Paz, W.; García, Y.; Salgado, J. Global Stability of Protein Folding from an Empirical Free Energy Function. J. Theor. Biol. 2013, 321, 44−53. (5) Pace, C. N.; Scholtz, J. M.; Grimsley, G. R. Forces Stabilizing Proteins. FEBS Lett. 2014, 588, 2177−2184. (6) Tewari, A. K.; Srivastava, P.; Singh, V. P.; Singh, P.; Khanna, R. S. Molecular Recognition Phenomenon in Aromatic Compounds. Res. Chem. Intermed. 2013, 39, 2925−2944. (7) Riley, K. E.; Hobza, P. On the Importance and Origin of Aromatic Interactions in Chemistry and Biodisciplines. Acc. Chem. Res. 2013, 46, 927−936. (8) Strutyński, K.; Gomes, J. A. N. F.; Melle-Franco, M. Accuracy of Dispersion Interactions in Semiempirical and Molecular Mechanics 2167
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
(53) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (54) Sinnokrot, M. O.; Sherrill, C. D. Highly Accurate Coupled Cluster Potential Energy Curves for the Benzene Dimer: Sandwich, TShaped, and Parallel-Displaced Configurations. J. Phys. Chem. A 2004, 108, 10200−10207. (55) Jurečka, P.; Šponer, J.; Č erny, J.; Hobza, P. Benchmark Database of Accurate (MP2 and CCSD(T) Complete Basis Set Limit) Interaction Energies of Small Model Complexes, DNA Base Pairs, and Amino Acid Pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985−1993. (56) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. Basis Set Consistent Revision of the S22 Test Set of Noncovalent Interaction Energies. J. Chem. Phys. 2010, 132, 144104. (57) Diep, P.; Johnson, J. K. An Accurate H2-H2 Interaction Potential from First Principles. J. Chem. Phys. 2000, 112, 4465−4473.
(30) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body Van Der Waals Interactions. Phys. Rev. Lett. 2012, 108, 236402. (31) Kronik, L.; Tkatchenko, A. Understanding Molecular Crystals with Dispersion-Inclusive Density Functional Theory: Pairwise Corrections and Beyond. Acc. Chem. Res. 2014, 47, 3208−3216. (32) Becke, A. D.; Johnson, E. R. Exchange-Hole Dipole Moment and the Dispersion Interaction. J. Chem. Phys. 2005, 122, 154104. (33) Becke, A. D.; Johnson, E. R. A Density-Functional Model of the Dispersion Interaction. J. Chem. Phys. 2005, 123, 154101. (34) Johnson, E. R.; Becke, A. D. A Post-Hartree−Fock Model of Intermolecular Interactions. J. Chem. Phys. 2005, 123, 024101. (35) Steinmann, S. N.; Corminboeuf, C. A Generalized-Gradient Approximation Exchange Hole Model for Dispersion Coefficients. J. Chem. Phys. 2011, 134, 044117. (36) Steinmann, S. N.; Corminboeuf, C. Comprehensive Benchmarking of a Density-Dependent Dispersion Correction. J. Chem. Theory Comput. 2011, 7, 3567−3577. (37) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Courier Dover Publications: Mineola, NY, 2012. (38) Conrad, J. A.; Gordon, M. S. Modeling Systems with π-π Interactions Using the Hartree−Fock Method with an Empirical Dispersion Correction. J. Phys. Chem. A 2014, DOI: 10.1021/ jp510288k. (39) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968−1986. (40) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646−9663. (41) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103, 379−387. (42) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632−672. (43) Pruitt, S. R.; Bertoni, C.; Brorsen, K. R.; Gordon, M. S. Efficient and Accurate Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2786−2794. (44) Amos, R. D.; Handy, N. C.; Knowles, P. J.; Rice, J. E.; Stone, A. J. Ab-Initio Prediction of Properties of Carbon Dioxide, Ammonia, and Carbon Dioxide···Ammonia. J. Phys. Chem. 1985, 89, 2186−2192. (45) Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. Rev. Mod. Phys. 1963, 35, 457−464. (46) Yamaguchi, Y.; Osamura, Y.; Goddard, J.; Schaefer, H., III. A New Dimension to Quantum Chemistry; Oxford University Press: New York, 1994. (47) Press, W.; Flannery, B.; Teukolsky, S.; Vetterling, W. Numerical Recipes: The Art of Scientific Computing (Fortran Version); Cambridge University Press: Cambridge, 1989. (48) Gross, E. K. U.; Ullrich, C. A.; Grossman, U. J. Density Functional Theory of Time-Dependent Systems. NATO ASI Ser., Ser. B: Phys. 1995, 337, 149−171. (49) Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107, 999− 1016. (50) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (51) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: GAMESS a Decade Later; Elsevier: Amsterdam, 2005. (52) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. 2168
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Dispersion Correction Derived from First Principles for Density Functional Theory and Hartree−Fock Theory Emilie B. Guidez and Mark S. Gordon* Department of Chemistry, Iowa State University, Ames, Iowa 50011, United States ABSTRACT: The modeling of dispersion interactions in density functional theory (DFT) is commonly performed using an energy correction that involves empirically fitted parameters for all atom pairs of the system investigated. In this study, the first-principles-derived dispersion energy from the effective fragment potential (EFP) method is implemented for the density functional theory (DFT-D(EFP)) and Hartree−Fock (HF-D(EFP)) energies. Overall, DFT-D(EFP) performs similarly to the semiempirical DFT-D corrections for the test cases investigated in this work. HF-D(EFP) tends to underestimate binding energies and overestimate intermolecular equilibrium distances, relative to coupled cluster theory, most likely due to incomplete accounting for electron correlation. Overall, this firstprinciples dispersion correction yields results that are in good agreement with coupledcluster calculations at a low computational cost.
■
developed by Grimme et al. The expression for the R−6 term in DFT-Dn is given by9
INTRODUCTION Intermolecular dispersion forces arise from the interaction between induced multipoles.1 These forces are at the origin of many chemical and biological processes such as protein folding,2−5 molecular recognition,6 and DNA base pair stacking.7 The modeling of dispersion interactions has been the subject of many investigations.8−11 Density functional theory (DFT) is frequently used to model molecular systems that contain on the order of hundreds of atoms.12,13 However, most commonly used density functionals (as well as Hartree− Fock (HF) theory) cannot account for dispersion interactions. Certain density functionals such as MPWB1K,14 M06-2X,15 M08-HX,16 M08-SO,16 and M1117 developed by Truhlar and co-workers correctly capture attractive noncovalent interactions where intermolecular overlap is nonnegligible (at the van der Waals minima). In addition, other functionals such as the van der Waals nonlocal correlation functionals (vdW-DF)18−22 include dispersion. In order to enable popular density functionals (GGAs, hybrids, ...) to account for dispersive interactions, Grimme and co-workers introduced a series of empirical corrections,23−25 collectively referred to here as “-D”. In DFT-D, the -D dispersion interaction energy correction is added to the Kohn−Sham energy. In general, the dispersion interaction between two molecules A and B can be expressed as26 C C C Edisp = 66 + 77 + 88 + ... (1) R R R R is the intermolecular distance and Cn are coefficients. The terms with an odd power of R are orientation-dependent,27,28 and average to zero in freely rotating systems. The dispersion energy expression is often truncated at the R−6 term, which corresponds to the interaction between induced dipoles of two species A and B. Three empirically parametrized dispersion corrections, called DFT-Dn (n = 1, 2, 3),23−25 have been © 2015 American Chemical Society
Nat − 1
Edisp = −S6
Nat
∑ ∑ I=1 J=I+1
C6IJ RIJ 6
fdmp (RIJ )
(2)
S6 is a scaling parameter that depends on the functional. CIJ6 , the dispersion coefficient for the atom pair IJ, is a fitted parameter, Nat is the number of atoms, and RIJ is the distance between atoms I and J. The repulsive interaction between the nuclei at very small RIJ is taken into account by using the damping function fdmp. In the DFT-Dn implementations by Grimme, the damping function also involves a fitted parameter. While DFT-Dn (n = 1, 2) only include the R−6 term,23,24 DFT-D3 also includes the R−8 term in a recursive manner.25 In DFT-D1 and DFT-D2, dispersion coefficients are predetermined using atomic ionization potentials and static polarizabilitites. There is no dependence on the atomic environment, which can lead to substantial errors. On the other hand, DFT-D3 considers the effective volume of the atoms by taking into account the number of neighbors in the environment to calculate the dispersion coefficients. Other methods such as the Tkatchenko−Scheffler29−31 and the Becke−Johnson32−34 models also consider the chemical environment of the atoms to calculate dispersion coefficients. The former method relies on the Hirshfeld partitioning of the total electron density between the atoms (which is obtained from electronic structure calculations) to compute the dispersion coefficients.29 The latter computes the dipole moment generated by the exchangecorrelation hole.34 Both methods are minimally empirical. The dDsC method, based on the generalized gradient approximation of the Becke−Johnson model, was developed later on.35,36 Received: January 13, 2015 Published: February 4, 2015 2161
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
α̅ represents one-third of the trace of the dynamic polarizability tensor. The use of the isotropic approximation is justified since anisotropic effects are minimal in this distributed approach.41 In addition, it permits the direct comparison of the C6 coefficients with experimental and other theoretical data and reduces computational cost.41 The integral in eq 7 is evaluated using a 12-point Gauss−Legendre quadrature formula.47 With a change of variable
One drawback of DFT-D can be the double counting of some of the electron correlation that is already included in the correlation functional. The parametrization attempts to minimize the effect of double counting. On the other hand, the Hartree−Fock method does not include any electron correlation except for the Fermi hole.37 Therefore, there is no possibility of double counting when one adds the -D correction to HF theory. However, other correlation effects are not included in HF-D. HF-D3 calculations using the Grimme corrections have been performed on several systems.38 Both DFT and HF have similar computational costs, which is one reason for the popularity of DFT. The effective fragment potential (EFP) method has been developed to treat intermolecular interactions.39,40 In the EFP method, all interaction energy terms, including the dispersion energy, are derived from first principles.41 In the EFP method, the interaction energy between the molecules or fragments is divided into five contributions: the Coulomb, polarization, exchange−repulsion, charge transfer, and dispersion interaction energies:40,42,43 E = Ecoul + Epol + Eex−rep + Ect + Edisp
ω = ω0
1+t , 1−t
dω =
−2ω0 dt (1 + t )2
(8)
Equation 7 can be rewritten as AB Edisp =−
3ℏ π
12
∑ ∑ ∑ Wk i∈A j∈B k=1
α̅ i(iω) α̅ j(iω) (1 − tk) R ij 6 2ω0
2
(9)
Wk and tk are the Gauss−Legendre weighting factor and abscissa.41 ω0 has an optimal value of 0.3.48 In order to avoid a singularity near R = 0, each term in eq 9 is multiplied by a damping function.49 An overlap-based damping function is used here: 6
(3)
fdmp (i , j) = 1 − Sij
2
∑ n=0
The objective of this work is to add the dispersion energy derived from first principles in the EFP method to the DFT Kohn−Sham energy (DFT-D(EFP)) and Hartree−Fock energy (HF-D(EFP)).
( −2 ln|Sij|)n /2 n!
(10)
Sij in eq 10 is the overlap between the localized molecular orbitals i and j. This damping function differs slightly from the one used originally in EFP49 in that it is now equally applicable to even and odd powers of n. An overlap-based damping function is a logical choice since the intermolecular overlap integral depends on the intermolecular distance, so the dispersion interactions are diminished as the overlap increases. In addition, unlike other damping functions, the overlap damping function contains no empirically fitted parameters.49 The present approach is derived from first principles and does not involve any empirically fitted parameters. The Tkatchenko− Scheffler method determines atomic polarizabilities based on the partitioning of the electron density as well as the effective atomic volumes. In contrast, the method described here solves the time-dependent Hartree−Fock equations using localized
■
METHOD In the DFT-D(EFP) and HF-D(EFP) implementations, each molecule represents a fragment within the EFP scheme. The dispersion energy between the molecules (fragments) is calculated as it would be in an EFP calculation. The dispersion interaction energy calculated in this manner is then added to the quantum mechanical (QM) energy of the system (Kohn−Sham or Hartree−Fock). The total energy of the system is given by E = Edisp + EQM (4) The expression for the dispersion energy derived from Rayleigh−Schrödinger perturbation theory is1,41,44 AB =− Edisp
ℏ 2π
x ,y,z
∑ ∑ ∑ Tαβij Tγδij ∫
0
i ∈ A j ∈ B αβγδ
∞
i ααγ (iω) αβδj (iω) dω
(5)
In eq 5, i and j label the localized molecular orbitals (LMOs) of molecules A and B, respectively.45 α(iω) represents the dynamic dipole polarizability tensor. The dynamic polarizability tensor is calculated by solving the time-dependent Hartree− Fock equations.41,46 T is a second order electrostatic tensor given by1,41,44 ij Tαβ
3R ijαR ijβ − R ij 2δαβ 1 1 = ∇α ∇β = R ij 4πε0R ij 4πε0R ij 5
(6)
Rij = Ri − Rj, where Ri and Rj are the coordinates of the LMO centroids i and j of fragments A and B, respectively. One can 1 ij ij 6 substitute ∑x,y,z αβγσ TαβTγσ = 6/Rij into eq 5. Within the isotropic 41 approximation, eq 5 further reduces to
Figure 1. Systems with a dominant dispersion interaction (first row): (A) benzene (sandwich); (B) methane dimer; (C) hydrogen dimer; (D) π-stacked adenine−thymine. Hydrogen-bonded systems (second row): (E) water dimer; (F) ammonia dimer; (G) methanol dimer; (H) adenine−thymine base pair (WC). Mixed systems (third row): (I) ethene−ethyne; (J) T-shape benzene dimer; (K) benzene−water dimer; (L) benzene−ammonia dimer. Color coding: black = carbon; white = hydrogen; red = oxygen; blue = nitrogen.
∞
AB Edisp
3ℏ =− π
∑∑ i∈A j∈B
∫0 α̅ i(iω) α̅ j(iω) dω R ij 6
(7) 2162
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 2. Potential energy curves of a set of test dimer cases with dominant dispersion interactions. (A) Benzene sandwich; (B) methane dimer; (C) hydrogen dimer with a 6-311++G(d,p) basis set; (D) hydrogen dimer with a 6-311++G(3df,3p) basis set; (E) adenine−thymine DNA stack. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set except for (D), where a 6-311++G(3df,3p) basis set is used. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations except for (D) where a 6-311++G(3df,3p) basis set is used.
molecular orbitals to calculate molecular polarizabilities. The DFT-D(EFP) binding energy between two molecules is calculated using the formula DFT‐D(EFP) ΔE = Edimer −
DFT ∑ Emonomer
was used with the 6-311++G(d,p) basis set. The HF-D(EFP) calculations were performed with the same basis set. For the hydrogen dimer, an additional set of calculations was done with the 6-311++G(3df,3p) basis set. The dimerization energies (eqs 11 and 12) were calculated at various intermolecular distances to generate a potential energy surface for each dimer. The internal monomer geometries were held fixed in the dimer calculations. All potential energy surfaces were generated by moving the molecules along the axis that connects their centers of mass at the equilibrium dimer orientation. The intermolecular equilibrium distance R0 and the equilibrium binding energy ΔE0, corresponding to the minimum on the potential energy surface, are reported for all methods. The equilibrium geometries of the benzene dimers (sandwich and T-shape) were optimized at the CCSD(T)/aug-cc-pVQZ* level of theory with frozen monomers.54 Equilibrium geometries of the methane, hydrogen, water, ammonia, and methanol dimers were taken from ref 49. The equilibrium geometries for the Watson−Crick DNA pair adenine−thymine (AD-WC) and the adenine−thymine stack
(11)
An analogous formula is used for HF-D(EFP) calculations: HF‐D(EFP) ΔE = Edimer −
HF ∑ Emonomer
(12)
In this work, the potential energy surfaces (PESs) of several test dimers are generated and compared with the PESs obtained with the Grimme DFT-Dn methods and with CCSD(T) and MP2 calculations. Note that, for DFT-Dn calculations, intramolecular dispersion forces are included in the energy of the monomer but they are not included in this new DFT-D(EFP) method.
■
COMPUTATIONAL DETAILS Calculations reported in this work were performed using the GAMESS software package.50,51 The DFT functional B3LYP52,53 2163
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
2164
a The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized MP2/cc-pVTZ CP geometry, except for the methane dimer which was optimized at the CCSD(T)/cc-pVTZ level of theory.56 cThese values are obtained at the CCSD(T)/aug-cc-pVQZ* level of theory from ref 54. dThese values are from the CCSD(T)/CBS potential energy surface of the T-shape hydrogen dimer by Johnson and Diep.57
(−3.08) (−0.38) (−0.052) (−0.087) 3.88 4.08 3.68 3.44 − (−4.55) (−0.34) (−0.042) (−0.073) (−18.21) benzene (sandwich) (CH4)2 (H2)2 (6-311++G(d,p)) (H2)2 (6-311++G(3df,3p)) A−T (stack)
4.02 3.90 3.16 3.16 3.28
(−0.55) (−0.21) (−0.049) (−0.055) (−10.45)
3.82 3.76 3.22 3.20 3.10
(−1.21) (−0.38) (−0.035) (−0.040) (−13.97)
3.86 3.96 3.30 3.30 3.22
(−1.66) (−0.38) (−0.068) (−0.074) (−12.51)
3.82 (−1.11) 3.88 (−0.24) N/A 3.30 (−0.027) 3.08 (−12.27)
4.24 4.10 3.92 3.60 3.26
(−0.60) (−0.28) (−0.022) (−0.060) (−10.98)
3.68 4.12 3.72 3.50 3.10
MP2 HF-D(EFP) DFT-D(EFP) DFT-D3 DFT-D2 DFT-D1
RESULTS AND DISCUSSION The three classes of systems that were investigated are shown in Figure 1. The first class corresponds to systems with dominant dispersion interactions, which include the sandwich configuration of the benzene dimer, the methane dimer, the hydrogen dimer, and the π-stacked adenine−thymine (Figure 1A−D). The second class corresponds to hydrogen-bonded systems, which include the water dimer, the ammonia dimer, the methanol dimer, and the adenine−thymine Watson−Crick DNA base pair (Figure 1E−H). The last class of complexes investigated is the mixed systems: the ethene−ethyne complex, the T-shape benzene dimer, and the benzene−water and benzene−ammonia complexes (Figure 1I−L). This classification was taken from ref 9. The potential energy surfaces of the first class of systems are displayed in Figure 2. The intermolecular equilibrium distances R0 and equilibrium binding energies ΔE0 of all of these systems are shown in Table 1. The benzene dimer is first investigated in a sandwich configuration. Like the semiempirical DFT-Dn methods (n = 1, 2, 3), DFT-D(EFP) and HF-D(EFP) both tend to underestimate the binding energy of the benzene sandwich in comparison to CCSD(T). MP2 overestimates the binding energy for this system. DFT-D(EFP) gives distances and binding energies that are similar to those of both DFT-D2 and DFT-D3, with an intermolecular equilibrium distance R0 of 3.82 Å and a binding energy ΔE0 of −1.11 kcal/mol. HF-D(EFP) overestimates the intermolecular equilibrium distance and underestimates the binding energy compared to CCSD(T), DFT-D2, DFT-D3, and DFT-D(EFP). In fact, the HF-D(EFP) results are very similar to those predicted by DFT-D1. For the methane dimer, all of the DFT-Dn methods underestimate the intermolecular equilibrium distance, while the HF-D(EFP) intermolecular distance is in excellent agreement with the CCSD(T) value. The binding energy at R0 for all of the -D methods is within 0.2 kcal/mol of the CCSD(T) value. The CCSD(T) interaction between two hydrogen molecules is very weak (smaller than 0.1 kcal/mol); it is therefore a challenge to model the PES accurately. As may be seen in Table 1, two basis sets were used for the H2 dimer: 6-311++G(d,p) and 6-311++G(3df,3p). The larger basis set increases the MP2 and CCSD(T) binding energies by ∼0.03 kcal/mol. The effect on DFT-Dn is small because the dependence on the basis set of the empirically determined -D methods is minimal. This is not the case for the -D(EFP) method, since it arises from first principles and therefore has an explicit basis set dependence. For the smaller basis set, DFT-D(EFP) does not have a minimum on the H2 dimer potential energy surface, while the HF-D(EFP) binding energy is smaller than that predicted by CCSD(T) by ∼0.03 kcal/mol. The binding energies predicted by both -D(EFP) methods are similar to those predicted by MP2 and CCSD(T). The π-stacking interaction between the DNA bases adenine and thymine (A−T) for the methods considered here are
dimer
Table 1. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Systems with Dominant Dispersion Contributionsa
■
CCSD(T)
CCSD(T)/CBSb
optimized at the MP2/cc-pVTZ level of theory with counterpoise correction (CP) were taken from ref 55. The equilibrium geometries of the benzene−water and benzene−ammonia complexes optimized at the MP2/cc-pVTZ level of theory with counterpoise correction were also taken from ref 55. The ethene−ethyne complex optimized at the CCSD(T)/cc-pVQZ level of theory was adapted from ref 55. The CCSD(T)/CBS energies were taken from refs 56, 54, and 57.
3.9c (−1.70) 3.71 (−0.53) 3.36d (−0.11) 3.36d (−0.11) 3.17 (−11.66)
The Journal of Physical Chemistry A
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 3. Potential energy curves of a set of hydrogen-bonded complexes. (A) Water dimer; (B) ammonia dimer; (C) methanol dimer; (D) adenine−thymine Watson−Crick pair. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations.
Table 2. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Hydrogen-Bonded Systemsa dimer (H2O)2 (NH3)2 (CH3OH)2 A···T (WC)
DFT-D1 2.90 3.28 3.58 5.96
(−6.25) (−3.89) (−6.44) (−16.55)
DFT-D2 2.88 3.28 3.54 5.94
(−6.48) (−4.19) (−6.88) (−18.17)
DFT-D3 2.90 3.28 3.56 5.96
(−6.41) (−4.06) (−6.82) (−17.90)
DFT-D(EFP) 2.84 3.20 3.50 5.84
(−6.52) (−4.24) (−6.94) (−20.11)
HF-D(EFP) 2.94 3.36 3.60 5.92
(−5.50) (−3.13) (−5.90) (−17.07)
MP2 2.92 3.32 3.58 5.98
(−5.93) (−3.75) (−6.47) (−16.03)
CCSD(T)
CCSD(T)/CBSb
2.92 (−5.85) 3.34 (−3.67) 3.58 (−6.46) −
2.91 (−5.02) 3.21 (−3.17) − 5.97 (−16.74)
a The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized CCSD(T)/cc-pVQZ geometry, except for the adenine−thymine DNA pair which was optimized with MP2/cc-pVTZ CP.56
only slightly smaller than the CCSD(T) value of 2.92 Å. The intermolecular equilibrium distance of 2.94 Å obtained with HF-D(EFP) is very close to the CCSD(T) and MP2 values. Similar trends are observed for the ammonia and methanol dimers, for which the DFT-D(EFP) method underestimates the intermolecular equilibrium distance and overestimates the binding energy. HF-D(EFP) overestimates R0 by only 0.02 Å for these two systems in comparison to CCSD(T) but underestimates the binding energy by about 0.55 kcal/mol. The adenine−thymine Watson−Crick pair contains two hydrogen bonds. The PESs obtained with all DFT-D methods for this complex are similar to the ones obtained with MP2. The values of R0 obtained with DFT-Dn (n = 1, 2, 3), and with DFT-D(EFP) and HF-D(EFP) are in good agreement with the MP2/cc-pVTZ CP optimized value. The DFT-Dn and HF-D(EFP) binding energies are all within 2 kcal/mol of the CCSD(T)/CBS value, while the DFT-D(EFP) method overestimates the binding energy by nearly 3.4 kcal/mol. The last class of systems investigated is the mixed systems. The PESs are shown in Figure 4. The values of R0 and ΔE0 are displayed in Table 3.
summarized in Table 1. CCSD(T) calculations were not performed for this system, but the value of R0 optimized at the MP2/cc-pVTZ CP level of theory55 and the value of ΔE0 obtained at the CCSD(T)/CBS level56 are reported in Table 1. MP2 greatly overestimates the binding energy of A−T: 18.2 kcal/mol vs 11.66 kcal/mol for CCSD(T)/CBS. In contrast, all of the DFT-D methods predict binding energies that are within 2 kcal/mol of the CCSD(T)/CBS value. Specifically, DFT-D(EFP) predicts a binding energy that is 0.61 kcal/mol higher than the CCSD(T)/CBS value. The HF-D(EFP) energy is just 0.68 kcal/mol smaller than the CCSD(T)/CBS value. All of the -D methods predict intermolecular distances that are in good agreement with the MP2/cc-pVTZ CP optimized value. The second class of systems is the hydrogen-bonded complexes. The PESs of these systems are shown in Figure 3, and the values of R0 and ΔE0 are reported in Table 2. For the water dimer, all DFT-D methods overestimate the binding energy in comparison to CCSD(T) by 0.5−0.7 kcal/mol, while the HF-D(EFP) binding energy is only 0.35 kcal/mol too low. The intermolecular equilibrium distances obtained with the DFT-D methods are between 2.84 and 2.90 Å, which is 2165
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
Figure 4. Potential energy curves of a set of mixed complexes. (A) Ethene−ethyne complex; (B) T-shape benzene dimer; (C) water−benzene complex; (D) ammonia−benzene complex. The functional B3LYP is used for all DFT calculations with a 6-311++G(d,p) basis set. A 6-311++G(d,p) basis set is also used for HF, MP2 and CCSD(T) calculations.
Table 3. Intermolecular Equilibrium Distances R0 in Å and Binding Energies ΔE0 in kcal/mol of Mixed Systemsa dimer ethene−ethyne benzene (T-shape) benzene−H2O benzene−NH3
DFT-D1 4.44 5.06 3.42 3.62
(−1.45) (−2.10) (−3.31) (−1.98)
DFT-D2 4.30 4.86 3.28 3.42
(−1.85) (−3.25) (−4.43) (−2.96)
DFT-D3 4.40 4.98 3.36 3.54
(−1.81) (−3.14) (−4.21) (−2.81)
DFT-D(EFP)
HF-D(EFP)
4.36 4.88 3.32 3.50
4.60 5.10 3.52 3.74
(−1.45) (−2.51) (−3.26) (−2.04)
(−1.23) (−2.16) (−2.93) (−1.80)
MP2 4.44 4.80 3.32 3.48
(−1.68) (−5.35) (−4.34) (−3.42)
CCSD(T) 4.48 4.88 3.36 3.52
(−1.53) (−4.38) (−4.04) (−3.09)
CCSD(T)/CBSb 4.42 (−1.51) 5.0c (−2.61) 3.38 (−3.29) 3.56 (−2.32)
a
The binding energy is the value indicated in parentheses. bEnergy values were obtained at the CCSD(T)/CBS(Δa(DT)Z) level of theory with the optimized MP2/cc-pVTZ geometry, except for the ethene−ethyne complex which was optimized at the CCSD(T)/cc-pVQZ level of theory.56 c These values are obtained at the CCSD(T)/aug-cc-pVQZ* level of theory from ref 54.
by ∼2 kcal/mol. The DFT-D2 and DFT-D3 methods underestimate the binding energy by ∼1 kcal/mol, while the DFT-D1 method is very similar to HF-D(EFP). For the water−benzene and ammonia−benzene complexes, DFT-D(EFP) slightly underestimates the intermolecular equilibrium distances and underestimates the equilibrium binding energies. On the other hand, HF-D(EFP) overestimates R0 for both the benzene−ammonia complex and the water− benzene complex. Overall, DFT-D(EFP) and HF-D(EFP) yield PESs that are similar to those of DFT-D1, while the DFT-D3 method is in better agreement with CCSD(T). In general, for the mixed species, the -D(EFP) method tends to underestimate binding energies slightly more than the DFT-D3 method, but the differences are typically 1 kcal/mol or less.
For the ethene−ethyne complex, the R0 values predicted by all of the DFT-D methods are smaller than those predicted by CCSD(T). DFT-D(EFP) and HF-D(EFP) slightly underestimate the binding energy (by 0.08 and 0.3 kcal/mol respectively), while the DFT-D2 and DFT-D3 methods overestimate the binding energy by 0.3 kcal/mol. For this system, the DFT-D1 method is in very good agreement with CCSD(T). For the T-shape benzene dimer, all DFT-Dn methods as well as HF-D(EFP) and DFT-D(EFP) underestimate the binding energy in comparison to CCSD(T). MP2 on the other hand overestimates the interaction between the monomers. DFT-D(EFP) predicts an equilibrium intermolecular distance that is identical to that of CCSD(T) and a binding energy that is nearly 2 kcal/mol too small. HF-D(EFP) overestimates R0 and also underestimates the equilibrium binding energy 2166
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
■
Models: The Benzene Dimer Case. J. Phys. Chem. A 2014, 118, 9561− 9567. (9) Antony, J.; Grimme, S. Density Functional Theory Including Dispersion Corrections for Intermolecular Interactions in a Large Benchmark Set of Biologically Relevant Molecules. Phys. Chem. Chem. Phys. 2006, 8, 5287−5293. (10) Grimme, S. Density Functional Theory with London Dispersion Corrections. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 211−228. (11) Klimeš, J.; Michaelides, A. Perspective: Advances and Challenges in Treating Van Der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (12) Levine, I. N. Quantum Chemistry; Pearson Prentice Hall: Upper Saddle River, NJ, 2009; Vol. 6. (13) Dreizler, R. M.; Engel, E. Density Functional Theory; Springer: Berlin, 2011. (14) Zhao, Y.; Truhlar, D. G. Hybrid Meta Density Functional Theory Methods for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions: The MPW1B95 and MPWB1K Models and Comparative Assessments for Hydrogen Bonding and Van Der Waals Interactions. J. Phys. Chem. A 2004, 108, 6908−6918. (15) Zhao, Y.; Truhlar, D. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (16) Zhao, Y.; Truhlar, D. G. Exploring the Limit of Accuracy of the Global Hybrid Meta Density Functional for Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2008, 4, 1849−1868. (17) Peverati, R.; Truhlar, D. G. Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810−2817. (18) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van Der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (19) Vydrov, O. A.; Van Voorhis, T. Nonlocal Van Der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103, 063004. (20) Vydrov, O. A.; Van Voorhis, T. Improving the Accuracy of the Nonlocal Van Der Waals Density Functional with Minimal Empiricism. J. Chem. Phys. 2009, 130, 104105. (21) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. Van Der Waals Density Functional: Self-Consistent Potential and the Nature of the Van Der Waals Bond. Phys. Rev. B 2007, 76, 125112. (22) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy Van Der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (23) Grimme, S. Accurate Description of Van Der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (24) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (25) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (26) London, F. The General Theory of Molecular Forces. Trans. Faraday Soc. 1937, 33, 8b−26. (27) Xu, P.; Zahariev, F.; Gordon, M. S. The R−7 Dispersion Interaction in the General Effective Fragment Potential Method. J. Chem. Theory Comput. 2014, 10, 1576−1587. (28) Buckingham, A. D. Theory of Long-Range Dispersion Forces. Discuss. Faraday Soc. 1965, 40, 232−238. (29) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and FreeAtom Reference Data. Phys. Rev. Lett. 2009, 102, 073005.
CONCLUSIONS The main conclusion of the present work is that the DFT-D(EFP) and HF-D(EFP) methods predict potential energy surfaces for a variety of types of intermolecular complexes with an accuracy that is comparable to that of the -D methods of Grimme and co-workers. The cost of computing the dispersion correction with this new method is slightly higher than the empirical -D methods but still small in comparison to the DFT part of the calculation. The importance of the -D(EFP) approach is that the dispersion correction is derived from first principles, without the need for empirically fitted parameters. In principle, this means that the EFP-based -D method is more broadly applicable. In the -D(EFP) methods, the intermolecular dispersion energy is calculated using the approach of the effective fragment potential method. The DFT-D(EFP) and HF-D(EFP) methods were applied to multiple test sets: complexes with dominant dispersion interactions, hydrogen-bonded complexes, and mixed complexes. For all of these systems, DFT-D(EFP) performs similarly to the semiempirical DFT-Dn (n = 1, 2, 3) corrections by Grimme. The HF-D(EFP) method also performs surprisingly well. Unlike the DFT methods, HF contains no other source of electron correlation (no correlation functional). Therefore, one cannot expect the HF-D method to perform consistently as well as DFT-D. Nonetheless, the HF-D method is appealing as a low-cost approximate alternative to MP2. The dispersion correction in DFT-D(EFP) and HF-D(EFP) only contains the R−6 term. Results could potentially be improved by including higher order terms.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Tel.: 515-294-0452. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by a U.S. National Science Foundation Software Infrastructure (SI2) grant, ACI-1047772.
■
REFERENCES
(1) Stone, A. The Theory of Intermolecular Forces, 2nd ed.; Oxford University Press: Oxford, 2013. (2) Vondrasek, J.; Kubar, T.; Jenney, F. E.; Adams, M. W. W.; Kozisek, M.; Cerny, J.; Sklenar, V.; Hobza, P. Dispersion Interactions Govern the Strong Thermal Stability of a Protein. Chem.Eur. J. 2007, 13, 9022−9027. (3) Serrano, A. L.; Waegele, M. M.; Gai, F. Spectroscopic Studies of Protein Folding: Linear and Nonlinear Methods. Protein Sci. 2012, 21, 157−170. (4) Ruiz-Blanco, Y. B.; Marrero-Ponce, Y.; Paz, W.; García, Y.; Salgado, J. Global Stability of Protein Folding from an Empirical Free Energy Function. J. Theor. Biol. 2013, 321, 44−53. (5) Pace, C. N.; Scholtz, J. M.; Grimsley, G. R. Forces Stabilizing Proteins. FEBS Lett. 2014, 588, 2177−2184. (6) Tewari, A. K.; Srivastava, P.; Singh, V. P.; Singh, P.; Khanna, R. S. Molecular Recognition Phenomenon in Aromatic Compounds. Res. Chem. Intermed. 2013, 39, 2925−2944. (7) Riley, K. E.; Hobza, P. On the Importance and Origin of Aromatic Interactions in Chemistry and Biodisciplines. Acc. Chem. Res. 2013, 46, 927−936. (8) Strutyński, K.; Gomes, J. A. N. F.; Melle-Franco, M. Accuracy of Dispersion Interactions in Semiempirical and Molecular Mechanics 2167
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
Article
The Journal of Physical Chemistry A
(53) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (54) Sinnokrot, M. O.; Sherrill, C. D. Highly Accurate Coupled Cluster Potential Energy Curves for the Benzene Dimer: Sandwich, TShaped, and Parallel-Displaced Configurations. J. Phys. Chem. A 2004, 108, 10200−10207. (55) Jurečka, P.; Šponer, J.; Č erny, J.; Hobza, P. Benchmark Database of Accurate (MP2 and CCSD(T) Complete Basis Set Limit) Interaction Energies of Small Model Complexes, DNA Base Pairs, and Amino Acid Pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985−1993. (56) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. Basis Set Consistent Revision of the S22 Test Set of Noncovalent Interaction Energies. J. Chem. Phys. 2010, 132, 144104. (57) Diep, P.; Johnson, J. K. An Accurate H2-H2 Interaction Potential from First Principles. J. Chem. Phys. 2000, 112, 4465−4473.
(30) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body Van Der Waals Interactions. Phys. Rev. Lett. 2012, 108, 236402. (31) Kronik, L.; Tkatchenko, A. Understanding Molecular Crystals with Dispersion-Inclusive Density Functional Theory: Pairwise Corrections and Beyond. Acc. Chem. Res. 2014, 47, 3208−3216. (32) Becke, A. D.; Johnson, E. R. Exchange-Hole Dipole Moment and the Dispersion Interaction. J. Chem. Phys. 2005, 122, 154104. (33) Becke, A. D.; Johnson, E. R. A Density-Functional Model of the Dispersion Interaction. J. Chem. Phys. 2005, 123, 154101. (34) Johnson, E. R.; Becke, A. D. A Post-Hartree−Fock Model of Intermolecular Interactions. J. Chem. Phys. 2005, 123, 024101. (35) Steinmann, S. N.; Corminboeuf, C. A Generalized-Gradient Approximation Exchange Hole Model for Dispersion Coefficients. J. Chem. Phys. 2011, 134, 044117. (36) Steinmann, S. N.; Corminboeuf, C. Comprehensive Benchmarking of a Density-Dependent Dispersion Correction. J. Chem. Theory Comput. 2011, 7, 3567−3577. (37) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Courier Dover Publications: Mineola, NY, 2012. (38) Conrad, J. A.; Gordon, M. S. Modeling Systems with π-π Interactions Using the Hartree−Fock Method with an Empirical Dispersion Correction. J. Phys. Chem. A 2014, DOI: 10.1021/ jp510288k. (39) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968−1986. (40) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646−9663. (41) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103, 379−387. (42) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632−672. (43) Pruitt, S. R.; Bertoni, C.; Brorsen, K. R.; Gordon, M. S. Efficient and Accurate Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2786−2794. (44) Amos, R. D.; Handy, N. C.; Knowles, P. J.; Rice, J. E.; Stone, A. J. Ab-Initio Prediction of Properties of Carbon Dioxide, Ammonia, and Carbon Dioxide···Ammonia. J. Phys. Chem. 1985, 89, 2186−2192. (45) Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. Rev. Mod. Phys. 1963, 35, 457−464. (46) Yamaguchi, Y.; Osamura, Y.; Goddard, J.; Schaefer, H., III. A New Dimension to Quantum Chemistry; Oxford University Press: New York, 1994. (47) Press, W.; Flannery, B.; Teukolsky, S.; Vetterling, W. Numerical Recipes: The Art of Scientific Computing (Fortran Version); Cambridge University Press: Cambridge, 1989. (48) Gross, E. K. U.; Ullrich, C. A.; Grossman, U. J. Density Functional Theory of Time-Dependent Systems. NATO ASI Ser., Ser. B: Phys. 1995, 337, 149−171. (49) Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107, 999− 1016. (50) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (51) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: GAMESS a Decade Later; Elsevier: Amsterdam, 2005. (52) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. 2168
DOI: 10.1021/acs.jpca.5b00379 J. Phys. Chem. A 2015, 119, 2161−2168
