A (Nearly) Universally Applicable Method for Modeling Noncovalent

Jun 2, 2012 - DCPs are utilized by simple modifications to input files and can be used in any computational package that can read effective-core poten...
0 downloads 5 Views 968KB Size
Download PDF
Letter pubs.acs.org/JPCL

A (Nearly) Universally Applicable Method for Modeling Noncovalent Interactions Using B3LYP Edmanuel Torres† and Gino A. DiLabio*,†,‡ †

National Institute for Nanotechnology, 11421 Saskatchewan Drive, Edmonton, Alberta, Canada T6G 2M9 Department of Physics, University of Alberta, Edmonton, Alberta, Canada

‡

S Supporting Information *

ABSTRACT: B3LYP is the most widely used density-functional theory (DFT) approach because it is capable of accurately predicting molecular structures and other properties. However, B3LYP is not able to reliably model systems in which noncovalent interactions are important. Here we present a method that corrects this deficiency in B3LYP by using dispersion-correcting potentials (DCPs). DCPs are utilized by simple modifications to input files and can be used in any computational package that can read effective-core potentials. Therefore, the technique requires no programming. DCPs (developed for H, C, N, and O) produce the best results when used in conjunction with 6-31+G(2d,2p) basis sets. The B3LYP-DCP approach was tested on the S66, S22, and HSG-A benchmark sets of noncovalently interacting dimers and trimers and was found to, on average, significantly outperform almost all other DFT-based methods that were designed to treat van der Waals interactions. Users of B3LYP who wish to model systems in which noncovalent interactions (viz., steric repulsion, hydrogen bonding, π-stacking) are present, should consider B3LYP-DCP. SECTION: Molecular Structure, Quantum Chemistry, and General Theory

I

noncovalent interactions. Our early (moderately successful) attempts to solve this problem were based on the fact that B3LYP (and other DFT methods) do not properly model the electron density far away (i.e., 3 − 10 Å) from the nuclei. The solution we proposed was to correct the long-range electron density through the use of atom-centered potentials named dispersioncorrecting potentials (DCPs)13−15 (see eq 1). This proposal is similar to a plane wave approach introduced by von Lilienfeld et al. for the treatment of systems using periodic DFT.16−18 The DCP approach is an attractive platform for the DFT dispersion problem because it can be utilized by any computational chemistry package that can use effective core potentials. Effective core potentials are normally used in simulations involving metals and other heavy atoms, and many computational packages have options for reading-in customized potentials. Therefore, a DCP-based solution to the dispersion problem in B3LYP can be implemented by modifying input files rather than by modifying the programs that run them, enabling many previous versions of popular computational packages (e.g., Gaussian 03, rev. D.01 and older) to perform dispersion-corrected B3LYP calculations with DCPs. Earlier generations of DCPs were developed for the carbon atom13−15,19 and corrected long-range electron density through the use of a local potential. In this report, we describe the development of a new generation of B3LYP-specific DCPs for

t is well-known among the computational methods development community that conventional density-functional theory (DFT) approaches are not capable of accurately and reliably modeling noncovalent interaction energies.1−4 This failing of DFT, which was described by Becke to arise from the fact that DFT does not contain the right physics to model these weak interactions, has inspired a great deal of research focused on the development of “dispersion-corrected” DFT methods.5−9 The importance of this effort is underscored by the fact that noncovalent interactions play a central role in many areas of science, ranging from biochemistry to condensed matter physics, and that many systems of interest are too large to be modeled using computationally expensive, highly accurate ab initio wave function techniques. Despite the fact that there has been significant progress in the development of dispersion-corrected DFTs, conventional DFTs are still widely used. For example, the B3LYP approach,10,11 which predicts repulsive long-range interactions between molecules that are attracted to each other by dispersion forces,4 continues to see widespread use, as judged by the fact that the original paper describing the approach was cited over 4000 times in 2011.12 The popularity of B3LYP lies in the fact that the method performs reasonably well for the prediction of many important molecular and reaction properties, e.g., geometries, bond energies, and barrier heights, and is especially well-validated for organic systems. No doubt there is a certain degree of name recognition and a good understanding of the pitfalls of B3LYP that also contribute to its popularity. The continuing widespread use of B3LYP motivated us to find a generally applicable solution to its inability to reliably model Published 2012 by the American Chemical Society

Received: May 2, 2012 Accepted: June 1, 2012 Published: June 2, 2012 1738

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

(l = p for H and l = f for C, N, O) or are semilocal (l = s for H and l = s, p, or d for C, N, O), and operate on electron density in specific angular momentum channels. The DCPs presented in this work were developed by optimizing cli and ξli values such that the error in the predicted interactions chosen from a set of noncovalently bonded dimers (the “fitting” set) was minimized relative to those obtained by CCSD(T)/CBS-quality methods.21−23 Table 1 lists the dimers belonging to the fitting set. The binding energies were computed at four to six points along potential energy surfaces (PESs) that describe dimer dissociation along a conveniently defined path. Representative points on the PES were chosen partway up the repulsive wall, near the PES minima, midway to dissociation, and near complete dimer dissociation. The members of the fitting set used for DCP optimization were somewhat dynamic and depended upon the stage of the optimization process. Some final hand-tuning was also done in certain cases. A random walk-based search method was used for the optimization. A guess at the form of the initial DCPs was provided from previous generations of DCPs13 or simple optimizations performed “by hand”, with some conjugate gradient preoptimizations performed in certain cases. The quality of a generated DCP was evaluated by a cost function computed as the weighted average of the absolute error of the DCP computed binding energies (BEDCP) relative to the highlevel reference binding energy (BEHL) data of the fitting set.

the hydrogen, carbon, oxygen, and nitrogen atoms that make use of both local and semilocal potentials. The performance of the DCPs is demonstrated on various benchmark sets of noncovalently bonded dimers and compared to other dispersion-corrected DFT methods. DCPs are atom-centered potentials (U(r)) composed of Gaussian-type functions, as described in ref 13, with the following general form: Nl

Ul(r ) = r −2 ∑ clir nli e−ζlir

2

(1)

i=1 20

In eq 1, Nl is the number of Gaussian functions, nli is the power of r (set to 2 throughout this work), cli is the coefficient of the Gaussian, and ξli is its exponent. Functions are local Table 1. Dimers Comprising the DCP Fitting Seta atom(s) carbon and hydrogen nitrogen oxygen mixedc

dimers (HCCH)2, (CH4)2b, (C2H6)2, stacked (C6H6)2,c T-shaped (C6H6)2c CH4−NH3, H3N−HNH2,d H3N−NH3,e (H3CCN)2 (CO2)2, (H2CO)2, H2O−HOH,d H2O−OH2e c,f c c uracil−uracil,

acetic acid dimer, formamide dimer

a

Structural data and PESs are provided in the Supporting Information. b Two conformers: D3h and C3v. cFitting data taken from the S66-8 data available at www.bedgb.com and described in ref 27. dHydrogenbonded. eLone-pair−lone-pair repulsion surface. fπ-stacked.

Table 2. Optimized B3LYP DCPs for the H, C, N, and O Atoms for Use with and without Counterpoise Correctionsa CP

non-CP

a

function type

ζi

ci

ζi

ci

H

P and higher

0.112790145 0.043193200 0.005733838 0.084641561

0.000226456 −0.000068835 −0.000000390 −0.000052055

0.120883601 0.044528578 0.005658790 0.174740501

0.000231333 −0.000070677 −0.000000451 −0.000049845

F and higher

0.154151674 0.107153962 0.011268728

0.000107696 0.001239979 −0.000000928

S−F

0.041251699 0.113028386 0.052663045 0.121760255 0.054474481

0.000021135 0.048953906 −0.000013041 −0.024832974 −0.006696057

0.204353696 0.106552301 0.103629038 0.010449572 0.043390720 0.117045600 0.064526518 0.123222363 0.066381605

0.000103405 −0.000068461 0.001279009 −0.000001033 0.000020668 0.056519710 −0.000013344 −0.021488594 −0.010148271

0.138865910 0.029813210 0.005040192 0.067091103 0.053208698 0.013353957

0.000423807 −0.000032790 −0.000000020 −0.000960396 −0.000545012 −0.000126251

0.209342166 0.164560971 0.014906716 0.042686780 0.122439426 0.066959033

0.000283597 0.002088445 −0.000000109 0.000020569 −0.000111442 −0.005065957

0.259303647 0.091449067 0.008021004 0.199348486 0.100227232 0.056155683

0.002157399 −0.000532364 −0.000000082 −0.000169550 −0.000220388 −0.000101932

0.192168931 0.166560549 0.016734867 0.039337457 0.123780982 0.061418151

0.000242019 0.003000000 −0.000000131 0.000016041 −0.000129441 −0.003520223

S−P C

P−F D-F N

F and higher

S−F P−F D-F O

F and higher

S−F P−F D-F

a The entries are given as angular momentum function type and the ζi and ci values according to eq 1. Additional formatting of the listed potentials is required for use in computational chemistry packages. The Supporting Information contains examples of input files demonstrating the use of DCPs in the Gaussian and NWChem packages.

1739

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

Table 3. Signed Deviations in Binding Energies (kcal/mol) Calculated for the Noncovalently Bonded Dimers of the S66 Benchmark Set Using B3LYP/6-31+G(2d,2p) with DCPs and with (CP) and without (non-CP) Counterpoise Correctionsa dimer

accepted high-level valueb

CP

Hydrogen Bonded Dimersc water−water 4.92 0.26 water−methanol 5.59 0.19 water−methylamine 6.91 0.32 water−peptided 8.10 0.18 methanol−methanol 5.76 0.17 methanol−methylamine 7.55 0.22 methanol−peptided 8.23 0.21 methanol−water 5.01 0.20 methylamine−methanol 3.06 −0.05 methylamine−methylamine 4.16 −0.29 methylamine−peptided 5.42 −0.23 methylamine−water 7.27 0.15 peptided−methanol 6.19 −0.20 peptided−methylamine 7.45 −0.05 peptided−peptided 8.63 −0.09 peptided−water 5.12 −0.15 uracil−uracil (base pair) 17.18 −0.10 water−pyridine 6.86 0.04 methanol−pyridine 7.41 −0.03 formic acid−formic acid 19.09 0.69 formamide−formamide 16.27 −0.13 formic acid−uracil 19.49 0.23 formamide−uracil 19.19 −0.03 mean absolute/signed error 0.18/0.07 mean absolute/signed % error 2.5/0.6 Dispersion-Dominated Dimers benzene−benzenee (π-stacked) 2.82 0.22 pyridine−pyridine (π-stacked) 3.90 0.18 uracil−uracil (π-stacked) 9.83 −0.44 benzene−pyridine (π-stacked) 3.44 0.21 benzene−uracil (π-stacked) 5.71 −0.30 pyridine−uracil (π-stacked) 6.82 −0.50 benzene−ethane 1.43 0.17 uracil−ethane 3.38 −0.11 uracil−ethyne 3.74 −0.17 pyridine−ethane 1.87 0.07 pentane−pentane 3.78 0.14 neopentane−pentane 2.61 −0.08 neopentane−neopentane 1.78 −0.15 cyclopentane−neopentane 2.40 −0.02 cyclopentane−cyclopentane 3.00 0.15

non-CP

dimer

accepted high-level valueb

CP

non-CP

Dispersion-Dominated Dimers benzene−cyclopentane 3.58 0.19 0.03 benzene−neopentane 2.90 0.00 −0.08 uracil−pentane 4.85 −0.26 0.00 uracil−cyclopentane 4.85 −0.14 0.00 uracil−neopentane 4.14 −0.20 0.04 ethene−pentane 3.71 −0.04 −0.13 ethyne−pentane 2.01 0.11 0.03 peptided−pentane 1.75 −0.08 −0.05 mean absolute/signed error 0.17/−0.05 0.09/−0.04 mean absolute/signed % error 4.7/−0.02 3.3/−1.3 Mixed Interactions benzene−benzene (T-shaped) 2.88 −0.10 0.07 pyridine−pyridine (T-shaped) 3.54 −0.29 −0.15 benzene−pyridine (T-shaped) 3.33 −0.06 0.02 benzene−ethyne (CH−π) 2.87 −0.04 0.15 ethyne−ethyne (T-shaped) 1.52 −0.04 0.09 benzene−acetic acid (OH−π) 4.71 −0.37 −0.07 benzene−formamide (NH−π) 4.36 −0.29 −0.10 benzene−water (OH−π) 3.28 −0.11 0.28 benzene−methanol (OH−π) 4.19 −0.06 0.18 benzene−methylamine (OH−π) 3.23 −0.13 −0.02 benzene−peptide (NH−π) 5.28 −0.18 −0.06 pyridine−pyridine (CH−N)f 4.15 −0.63 −0.47 ethyne−water (CH−O)f 2.85 −0.18 0.07 ethyne−acetic acid (OH−π)f 4.87 −0.07 0.13 pentane−acetic acid 2.91 0.04 0.13 pentane−formamide 3.53 −0.10 −0.06 benzene−acetic acid (π-stacked) 3.80 0.06 0.18 peptide−ethane 3.00 −0.09 −0.04 methylamine−pyridineg 3.97 −0.39 −0.11 pyridine−ethyne 3.99 −0.26 −0.06 mean absolute/signed error 0.17/-0.05 0.12/0.01 mean absolute/signed % error 4.7/-0.2 3.5/0.7 Performance Statistics for the Entire S66 Set mean absolute error 0.17 0.14 mean signed error −0.06 0.03 mean absolute percent error 3.9 3.3 mean signed percent error −1.4 0.4 error range −0.63 to −0.47 to +0.69 +0.77 percent error range −15.3 to −14.3 to +11.6 +11.2

0.43 0.24 0.77 0.12 0.14 0.64 0.04 0.36 0.09 0.11 −0.08 0.62 −0.08 0.36 −0.10 0.20 −0.20 0.45 0.31 −0.09 −0.07 −0.19 −0.05 0.25/0.18 3.7/3.1 0.14 0.00 −0.02 0.05 −0.04 −0.21 0.09 0.08 0.06 0.04 −0.13 −0.26 −0.25 −0.19 −0.09

a The data are ordered according to interaction type, with performance statistics provided for data in each section. Overall performance statistics for the entire S66 set are provided at the bottom of the Table. bReference values are generally of CCSD(T)/CBS quality and are taken from ref 27. c All indicated dimers that have a hydrogen-bonded interaction are named in the form donor−acceptor. d“Peptide” is the nomenclature used in ref 27 and refers to methylamide. eSlipped-parallel configuration. fHydrogen-bonding interaction. gNH−N and CH(methylamine)−π interactions.

The optimization routine used the Gaussian 09 program24 to compute the BEs with the trial DCPs using the B3LYP functional25 with 6-31+G(2d,p) basis sets. Optimized DCPs were obtained with and without the inclusion of counterpoise (CP) corrections26 for basis set superposition error and are presented in Table 2. The accuracy of the BEs predicted by B3LYP increases as a function of the strength and nature of the noncovalent interactions (see, for example, Figure 2 of ref 4). The main challenge associated with the development of DCPs for B3LYP

was to improve their performance for dispersion and other weak interaction without degrading their performance for hydrogen (H)-bonding interactions. The use of fully local potentials (i.e., those that use high (f) angular momentum functions for heavy atoms) tend to affect BEs indiscriminately, and so it was necessary to introduce semilocal (i.e., s, p, and d) functions into the DCP in order to achieve the necessary selectivity for reproducing particular interactions. Similar considerations were made, for example, in developing DCPs that can accurately describe the behavior of sp, sp2 and sp3 carbon in various 1740

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

Although six members of our fitting set overlap with the validation set, the latter is quite large, making it essentially independent of the fitting set. This is an important consideration from the perspective of the generality of the optimized DCPs. The DCPs were evaluated by calculating absolute and signed errors (defined such that overbinding yields a positive number), absolute and signed percent errors, and statistics were generated by calculating mean absolute error (MAE = (1/k)Σki | k DCP − BEHL − BEDCP i i |) and mean sign error (MSE = (1/k)Σi (BEi )), along with the analogous mean percent errors. The BEHL i geometries used for the calculations associated with the S66 set were not altered from their database values.28 The data presented in Table 3 show that, in general, there is overbinding of H-bonded dimers and under-binding in dispersion-bound dimers. Nevertheless, the performance of B3LYP-DCP/6-31+G(2d,2p) is excellent, and the form of the DCPs presented in Table 2 is clearly very promising for treating noncovalent interactions in a variety of systems. We note, however, that DCPs cannot correct the underlying deficiencies in B3LYP that are associated with erroneous charge transfer within some charge-transfer complexes. The significant improvement that DCPs offer to the B3LYP approach in predicting noncovalent interactions is clearly illustrated in Figure 1, which shows the frequencies associated with the signed percent errors in predicted binding energies of the members of the S66 set. Figure 1a shows that non-CP corrected B3LYP/6-31+G(2d,2p) without dispersion corrections perform best on dimers interacting by H-bonding and worst for dimers interacting principally by dispersion. Overall, the signed percent errors range from ca. +10% to −220%. Figure 1b,c illustrate how DCPs improve the performance of B3LYP/ 6-31+G(2d,2p): All of the signed percent errors in the binding energies of the S66 set are now binned between ca. ± 14%. Table 4a demonstrates the performance of B3LYP-DCP with various basis sets against the S66 benchmark set. Additional performance assessments were made using the S22 set of dimers.29−31 The DCPs that were designed for use with CP corrections perform quite consistently as a function of basis set size. The non-CP DCPs show more variability with basis set, with performance degrading by a factor of ca. 3 when the best and worst results are compared. We did some testing of B3LYP-DCP using very small basis sets (i.e., 3-21G(d), etc) and found that the performance was unacceptably poor. We therefore recommend that basis sets no smaller than 6-31+G(d,p) be used in conjunction with B3LYP-DCP. The best performance is achieved with non-CP corrected B3LYPDCP/6-31+G(2d,2p), which gives MAEs of 0.14/0.23 kcal/mol and MAPEs of 3.3/5.4 for the S66/S22 benchmark sets. Table 4b provides some context for the quality of the predictions made by B3LYP-DCP/6-31+G(2d,2p) by comparing its performance statistics to those of other methods using 6-31+G(2d,2p) basis sets. In connection to Figure 1, B3LYP without DCPs is shown to have MAEs ranging from ca. 3.0 to 3.8 kcal/mol and MAPEs up to 93%. Commonly available dispersion-corrected DFT methods such as M06-2X,8 B97-D32 and ωB97X-D9 show good performance for both the S66 and S22 benchmark sets, and, of the three methods, M06-2X performs best. However, all three approaches are outperformed by B3LYP-DCP/6-31+G(2d,2p). Table 4c summarizes some of the recent results reported by Goerigk et al.33 using very large basis sets with their DFT-D3 approach and with the M06-2X method. The use of very large basis sets with M06-2X results in slightly worse performance as

Figure 1. Histograms showing the frequency of signed percent errors of calculated binding energies for the members of the S66 set of noncovalently bonded dimers using (a) B3LYP/6-31+G(2d,2p) without corrections for dispersion, and (b) B3LYP/6-31+G(2d,2p) used with DCPs. The binding energies were computed without CP corrections. (c) Close-up view of the plot in panel b.

bonding scenarios and intermolecular interactions (e.g., methane dimer vs π-stacked benzene dimer). Nevertheless, as can be seen in Table 2, the optimized DCPs are dominated by high-angular momentum functions. To a significant extent, the DCPs designed for use without CP corrections also compensate for basis set superposition errors, and this generally results in functions in the high-angular momentum channel having more repulsive character compared to the DCPs designed for use with CP corrections. We point out that the CP and non-CP DCPs are not interchangeable. That is, the application of DCPs that were designed for use with CP corrections in non-CP calculations would result in poor performance. The converse is also true (vide inf ra). The performance of the DCPs were assessed by computing the BEs of the compounds in the S66 set27 of noncovalently bonded dimers with B3LYP-DCP/6-31+G(2d,2p) and comparing those values to the BEHL values tabulated for the set. This set contains 66 dimers (see Table 3) that are bound by dispersion, hydrogen-bonding, and/or mixed interactions, and represent interaction strengths ranging over ca. 1.75 to 19.5 kcal/mol. 1741

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

Table 4. MAE (MAPE) of Binding Energies for the Noncovalently Bonded Dimers of the S66a and S22b Sets Calculated Using Selected DFT-Based Methods,c with (CP) and without (non-CP) Counterpoise Corrections (a) DCPs: B3LYP-DCP/Basis S66 basis set

CP

6-31+G(d,p) aug-cc-pVDZ 6-31+G(2d,p)d 6-31+G(2d,2p)e,f 6-311+G(2d,2p) aug-cc-pVTZ

0.20 (4.2) 0.19 (4.5) 0.16 (3.7) 0.17 (3.9)g 0.16 (3.8) 0.16 (3.9)

S22 non-CP

0.29 (6.0) 0.35 (9.7) 0.16 (3.4) 0.14 (3.3)h 0.22 (5.2) 0.23 (5.3) (b) Other DFT/6-31+G(2d,2p)

CP

non-CP

0.29 (5.9) 0.34 (6.8) 0.32 (6.0) 0.31 (6.1i) 0.30 (5.8) 0.30 (5.9)

0.31 (7.8) 0.23 (7.3) 0.28 (6.0) 0.23 (5.4)j 0.31 (5.9) 0.37 (6.9)

S66 method

CP

B3LYP B97-D M06−2X ωB97XD

3.48 (92.6) 0.39 (8.3) 0.30 (7.6) 0.64 (15.8)

S22 non-CP 3.06 (82.6) 0.39 (9.6) 0.25 (5.8) 0.30 (8.2) (c) Other DFT/def2-QZVP

CP

non-CP

3.87 (88.6) 0.46 (5.6) 0.44 (7.2) 0.18 (5.3)

3.35 (78.0) 0.45 (9.3) 0.43 (6.0) 0.58 (13.3)

S66 method

CP

BLYP-D3k,l B3LYP-D3k,l M06−2Xj

S22 non-CP 0.39m 0.40n 0.28

CP

non-CP 0.24o 0.36o 0.40

a

See ref 27. bSee ref 29, with binding energies taken from ref 31. cThe units for MAE and MAPE are kcal/mol and %, respectively. dThis basis set was used to optimize the DCPs. eBasis set recommended for use with the B3LYP DCPs reported in this work. fWhen CP corrections are not applied to the binding energy data obtained from DCPs that were developed for use with CP corrections, the MAE and MAPE, respectively, increase to the following: g0.38 kcal/mol and 9.0% and i0.41 kcal/mol and 10.2%. When CP corrections are applied to the binding energy data obtained from DCPs that were developed for use without CP corrections, the MAE and MAPE, respectively, increase to the following: h0.41 kcal/mol and 10.1% and j0.61 kcal/mol and 11.5%. kTaken from ref 33. lUsing the original damping function implementation of D3, see ref 33 and references contained therein. MAEs decrease to the following: m0.19 and n0.28 when the short-range damping of Becke and Johnson34 is used. oData from ref 36.

approaches produce MAEs of 0.15 kcal/mol. The largest error in the BEs predicted by the non-CP approach for the acetate ion−water−acetic acid trimer is large (−1.0 kcal/mol, 6.1%). The extensive review of dispersion-corrected DFT methods provided by Burns et al.35 allows us to place the performance of B3LYP-DCP/6-31+G(2d,2p) very nearly equivalent to the B3LYP/aug-cc-pVTZ approach corrected for dispersion using the exchange-hole dipole moment (XDM) technique of Johnson and Becke (viz., MAE = 0.14 kcal/mol).5 For comparison, Burns et al. report MAEs for the M06-2X, ωB97X-D, and B3LYP-D3 methods with aug-cc-pVTZ of 0.47, 0.32, and 0.48 kcal/mol, respectively. In summary, we have developed new DCPs for use with the B3LYP functional. The DCPs can be used with any basis set. However, performance degrades significantly if bases smaller than 6-31+G(d,p) are employed. Optimal performance is achieved with 6-31+G(2d,2p) basis sets, which are much small than those recommend for use with DFT-D3 methods. The B3LYP-DCP/6-31+G(2d,2p) method is one of the best performing dispersion-corrected DFT methods currently available. This approach should be considered by users studying systems in which noncovalent interactions and/or steric repulsions are present. DCPs provide a nearly universal solution to the dispersion problem in B3LYP because they can be used with any computational chemistry packages (e.g., NWChem and Gaussian09/03/98, etc.) that can read effective core potential input. The DCPs are specified by modifying the input files used in these packages and therefore do not require reprogramming.

compared to the results obtained using 6-31+G(2d,2p) basis sets (see Table 4b). The BLYP-D3 and B3LYP-D3 methods display fair performance on the S66 and S22. However, DFT-D3 methods generally improve when they incorporate the shortrange damping function proposed by Becke and Johnson,34 rather than the one used by Grimme.6 It is not obvious from the results presented in ref 33 how much the performance of the DFT-D3 methods would degrade when smaller basis sets are utilized. However, degradation in performance is expected because the D3 corrections were fitted for use with def2-QZVP basis sets. For example, Burns et al., in their extensive review of dispersion-corrected DFT methods, report that B3LYP-D3/ aug-cc-pVDZ produces an MAE of 0.79 kcal/mol for the S66 set,35 which is more than a factor of 3 higher than the comparable value obtained by B3LYP-DCP/aug-cc-pVDZ (see Table 4a). A final assessment of the performance of B3LYP-DCP/ 6-31+G(2d,2p) was made by applying the method to the HSG-A benchmark set (Table 5). This set of noncovalently bound dimers and trimers is comprised of 21 neutral and charged fragments obtained from the HIV−II protease crystal structure having a bound inhibitor molecule, indinavir. Highlevel binding energies were reported by Faver et al.37 and recently refined by Marshall et al.31 For our calculations, we used structures graciously provided to us by Professor David Sherrill. In general, we find that B3LYP-DCP/6-31+G(2d,2p) performs well for the HSG-A set. Both the CP and non-CP DCPs produce similar BEs, as reflected by the fact that both 1742

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

Table 5. Signed Deviations in Binding Energies (in kcal/mol) Calculated for the Noncovalently Bonded Dimers of the HSG-A Set Using B3LYP/6-31+G(2d,2p) with DCPs and with (CP) and without (non-CP) Counterpoise Correctionsa complex

accepted high-level valueb

CP

non-CP

0.518 2.283 2.478 16.526 19.076 5.998 3.308 0.581 5.066 7.509 6.274 −0.302 2.103 1.378 0.856 1.100 1.534 0.472 1.598 −0.378 9.538

0.07 −0.11 0.27 0.43 0.12 0.35 0.13 0.26 −0.01 0.27 0.13 −0.06 0.00 0.02 0.03 0.16 −0.12 0.04 −0.14 0.13 0.26 0.15/0.11 9.2/−3.8

0.11 −0.07 −0.02 1.01 0.28 0.19 0.07 0.22 0.02 0.29 0.15 0.01 0.06 0.05 0.10 0.14 −0.03 0.08 −0.15 0.04 0.02 0.15/0.12 7.5/−4.7

e

methane−n-t-butylformamide ethane−indan-2-olf benzene−n-methylguanidine(h+)c,g acetate ion−water−acetic acidh 3-acetylamino-propioniate ion−methanoli benzene−acetate ionj dimethylaminomethanol−2-faad,k pyridine−2-formylaminoacetamidel n-methylacetamide−2-faam n-methylacetamide−n-methylacetamiden n-methylacetamide−2-faao propane−n-t-butylformamidep butane−benzeneq butane−n-t-butylformamider ethane−ethane-1,2-diamines toluene−2-methylpropanet ethane−3-methylpyridineu ethane−n-t-butylformamidev propane−benzenew propane−benzenex acetate ion−water−n-methylacetamidey mean absolute/signed error mean absolute/signed % error a

Statistics associated with the performance of the methods are also provided. bTaken from reference 35. cProtonated N-methylguanidine. 2-Formylaminoacetamide. The names of the complexes associated with the HSG-A database are as follows: eala29-big; fala128-small; garg8; hash26asp125; iasp129-big; jaso130; kgly28-big; lgly50-ring-big; mgly50-v1; ngly127; ogly148; pile48-big; qile147; rile150-big; sile184; tleu23-big; upro181; v val33-big; wval83; xval132; ywat200. d

and thank Prof. C. David Sherrill, Dr. Lori R. Burns, and Dr. Michael S. Marshall for providing data and assistance with the HSG-A benchmark set.

Examples of how to utilize the B3LYP-DCP approach are given in the Supporting Information and at www.ualberta.ca/ ∼gdilabio. A program for generating DCP input files for Gaussian XX and NWChem is also available on our Web site. Future work will focus on the expansion of the DCP library to other atoms in the periodic table and on the development of DCPs for other DFT approaches.

■

■

ABBREVIATIONS DFT, density-functional theory; DCP, dispersion-correcting potentials; CP, counterpoise corrections/counterpoise corrected

■

ASSOCIATED CONTENT

S Supporting Information *

Sample input files illustrating the use of DCPs with the Gaussian and NWChem suite of programs, tables listing the signed deviations in predicted binding energies for the S66 and S22 sets using B3LYP-DCP/basis for the basis sets listed in Table 4a, and details of the fitting set used to optimized the DCPs. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

(1) Kristyán, S.; Pulay, P. Can (Semi)Local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175−180. (2) Tsuzuki, S.; Lüthi, H. P. Interaction Energies of van der Waals and Hydrogen Bonded Systems Calculated Using Density Functional Theory: Assessing the PW91 Model. J. Chem. Phys. 2001, 114, 3949− 3957. (3) Johnson, E. R.; Wolkow, R. A.; DiLabio, G. A. Application of 25 Density Functionals to Dispersion-Bound Homomolecular Dimers. Chem. Phys. Lett. 2004, 394, 334−338. (4) Johnson, E. R.; DiLabio, G. A. Structure and Binding Energies in van der Waals Dimers: Comparison Between Density Functional Theory and Correlated Ab Initio Methods. Chem. Phys. Lett. 2006, 419, 333−339. (5) Becke, A. D.; Johnson, E. R. Exchange-Hole Dipole Moment and the Dispersion Interaction Revisited. J. Chem. Phys. 2007, 127, 154108. (6) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (7) Grimme, S.; Antony, J.; Ehrlich, S.; Kreig, 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.

AUTHOR INFORMATION

Corresponding Author

*Phone: +1-780-641-1729; FAX: +1-780-641-1601; E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are extremely grateful to the Centre for Oil Sands Innovation, based at the University of Alberta, for funding, and to Compute/Calcul Canada and WestGrid for a generous allocation of computing resources. We benefitted from many enlightening discussions with Prof. Erin R. Johnson (UC Merced) 1743

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744

The Journal of Physical Chemistry Letters

Letter

(29) Jurečka, P.; Šponer, J.; Č erný, 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 See also reference 28. The “B” designation in the set name refers to revised binding energies as computed by Marshall et al., reference 31. (30) The S22 set contains 22 dimers that interact principally by dispersion, hydrogen-bonding, and/or mixed interactions. As this set is somewhat older, we relegate the details of the DCP performance to the Supporting Information and present only performance statistics in the text. (31) Marshall, M. S.; Burns, L. A.; Sherrill, C. D. Basis Set CCSD(T) : Best Convergence of the Coupled-Cluster Correction, δMP2 Practices for Benchmarking Noncovalent Interactions and the Attendant Revision of the S22, NBC10, HBC6, and HSG Databases. J. Chem. Phys. 2011, 135, 194102. (32) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (33) Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking Density Functional Methods against the S66 and S66 × 8 Datasets for Noncovalent Interactions. ChemPhysChem 2011, 12, 3421−3433. (34) Johnson, E. R.; Becke, A. D. A Post-Hartree−Fock Model of Intermolecular Interactions: Inclusion of Higher-Order Corrections. J. Chem. Phys. 2006, 124, 174104 and references therein. (35) Burns, L. A.; Vázquez-Mayagoitia, Á .; Sumpter, B. G.; Sherrill, C. D. Density-Functional Approaches to Noncovalent Interactions: A Comparison of Dispersion Corrections (DFT-D), Exchange-Hole Dipole Moment (XDM) Theory, and Specialized Functionals. J. Chem. Phys. 2011, 134, 084107. (36) Goerigk, L.; Grimme, S. A Thorough Benchmark of Density Functional Methods for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. Phys. Chem. Chem. Phys. 2011, 13, 6670−6668. (37) Faver, J. C.; Benson, M. L.; He, X. A.; Roberts, B. P.; Wang, B.; Marshall, M. S.; Kennedy, M. R.; Sherrill, C. D.; Merz, K. M. Formal Estimation of Errors in Computed Absolute Interaction Energies of Protein−Ligand Complexes. J. Chem. Theory Comput. 2011, 7, 790− 797.

(8) Zhao, Y.; Truhlar, D. G. 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. (9) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom−Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (10) Becke, A. D. Density-Functional Thermochemistry. 3. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (11) 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. (12) Web of Science citation report. www.webofknowledge.com (accessed April 2012). (13) DiLabio, G. A. Accurate Treatment of van der Waals Interactions Using Standard Density Functional Theory Methods With Effective Core-Type Potentials: Application to CarbonContaining Dimers. Chem. Phys. Lett. 2008, 455, 348−353. (14) Mackie, I. D.; DiLabio, G. A. Interactions in Large, Polyaromatic Hydrocarbon Dimers: Application of Density Functional Theory with Dispersion Corrections. J. Phys. Chem. A 2008, 112, 10968−10976. (15) Mackie, I. D.; DiLabio, G. A. Accurate Dispersion Interactions From Standard Density-Functional Theory Methods with Small Basis Sets. Phys. Chem. Chem. Phys. 2010, 12, 6092−6098. (16) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Optimization of Effective Atom Centered Potentials for London Dispersion Forces in Density Functional Theory. Phys. Rev. Lett. 2004, 93, 153004. (17) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Performance of Optimized Atom-Centered Potentials for Weakly Bonded Systems Using Density Functional Theory. Phys. Rev. B 2005, 71, 195119. (18) Lin, I.-C.; Coutinho-Neto, M. D.; Felsenheimer, C.; von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Library of DispersionCorrected Atom-Centered Potentials for Generalized Gradient Approximation Functionals: Elements H, C, N, O, He, Ne, Ar, and Kr. Phys. Rev. B 2007, 75, 205131. (19) DCPs were also developed for the silicon atom. See Johnson, E. R.; DiLabio, G. A. Theoretical Study of Dispersion Binding of Hydrocarbon Molecules to Hydrogen-Terminated Silicon(100)-2 × 1. J. Phys. Chem. C 2009, 113, 5681−5689. (20) Christiansen, P. A.; Lee, Y. S.; Pitzer, K. S. Improved Ab Initio Effective Potentials for Ar, Kr, and Xe with Applications to Their Homonuclear Dimers. J. Chem. Phys. 1979, 71, 4445. (21) Mackie, I. D.; DiLabio, G. A. Approximations to Complete Basis Set-Extrapolated, Highly Correlated Noncovalent Interaction Energies. J. Chem. Phys. 2011, 135, 134318. (22) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. An Assessment of Theoretical Methods for Nonbonded Interactions: Comparison to Complete Basis Set Limit Coupled-Cluster Potential Energy Curves for the Benzene Dimer, the Methane Dimer, Benzene−Methane, and Benzene−H2S. J. Phys. Chem. A 2009, 113, 10146−10159. (23) In some cases, an approach outlined in ref 21 was used to generate high-level binding energies. (24) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, revision C.01; Gaussian, Inc.: Wallingford, CT, 2009. (25) As implemented in reference 24. (26) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies − Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (27) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. S66: A Well-Balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures. J. Chem. Theory Comput. 2011, 7, 2427−2438. (28) Interested readers are directed to www.begdb.com. 1744

dx.doi.org/10.1021/jz300554y | J. Phys. Chem. Lett. 2012, 3, 1738−1744