FULL PAPER
WWW.C-CHEM.ORG
Halogen Bonds with Benzene: An Assessment of DFT Functionals Alessandra Forni,*[a] Stefano Pieraccini,*[a,b] Stefano Rendine,[b] and Maurizio Sironi*[a,b] The performance of an extensive set of density functional theory functionals has been tested against CCSD(T) and MP2 results, extrapolated to the complete basis set (CBS) limit, for the interaction of either DCl or DBr (D 5 H, HCC, F, and NC) with the aromatic system of benzene. It was found that double hybrid functionals explicitly including dispersion, that is, B2PLYPD and mPW2PLYPD, provide the better agreement with the CCSD(T)/CBS results on both energies and equilibrium geometry, indicating the importance of dispersive contribu-
tions in determining this interaction. Among the less expensive functionals, the better performance is provided by the xB97X and M062X functionals, while the xB97XD and B97D functionals are shown to work very well for bromine comC 2013 Wiley plexes but not so well for chlorine complexes. V Periodicals, Inc.
Introduction
is compensated by an increase of electronic charge density in a belt around the bond axis. The site of electronic depletion can hence interact with a source of electrons represented by the Lewis base A located along the direction of the DAX bond. This pictorial representation has been confirmed by quantum mechanical calculations using Natural Bond Orbitals analysis and by looking at the details of the electrostatic potential (ESP) around the halogen atom. The ESP shows in fact the presence of a positive region along the extension of the DAX bond, denoted as “r-hole,”[15] and the formation of a negative belt perpendicular to the DAX bond. The essential features of such r-hole, that is, the maximum value of the ESP and the size of the positive region, strongly depend, of course, on the moiety D and on the nature of the halogen atom. In particular, the r-hole becomes more pronounced as the electron-withdrawing capability of D and the polarizability of the halogen atom increase. As a result, for a same D group, the iodine atom shows a well-defined r-hole, which decreases on going to bromine and then to chlorine. The r-hole is not observed on fluorine, unless very strong electronegative moieties are present in the molecule.[16–18] A good correlation has been observed between the values of the ESP on the r-hole and the strength of the halogen bond interaction.[14,15,19,20] The presence of the r-hole explains also the directionality of halogen bonding. It is expected,
In recent years, halogen bonding[1,2] has received a growing attention as an important type of interaction involved in many different fields ranging from material science to biochemistry. Halogen bonding has been shown to play a key role in the development of new solid-state materials with attractive properties, like for example organic semiconductors and superconductors,[3] and materials for nonlinear optical applications.[4–7] In biochemical applications, halogen bonding has offered interesting perspectives in protein-ligand complexation with important consequences in medicinal chemistry.[8–11] Halogen bond is a highly directional interaction between a covalently bonded halogen atom X and an acceptor group A with nucleophilic character, according to the scheme DX/A. The moiety D bonded to the halogen atom has a large variability, ranging from inorganic to organic species, whereas the nucleophilic site A is usually represented by a lone pair of a heteroatom such as oxygen, nitrogen, sulfur, or by a p system as, for example, a phenyl group. Of course, the strength of halogen bonding depends on the halogen atom and on the nature of both the moieties D and A. A broad range of halogen bond interaction strengths have been reported, from about 1.6 to 43.0 kcal/mol for the ClCl interaction in Cl2[12] and the I2I2 interaction in I32,[13] respectively. The physical origin of such quite unexpected phenomenon, where an electronegative species X interacts with a nucleophilic site A, has been explained by some authors[14,15] on the basis of a simple intuitive model. The electronic structure of a halogen atom covalently bonded to another atom can be roughly described as s2px2py2pz1 assuming z as the direction of the bond. This electronic arrangement is associated with a depletion of the electron charge distribution along the bond axis with respect to a 5/3 mean population of p electrons in each direction, as found in a (spherical) isolated halogen atom. More precisely, such reduction of charge distribution is predicted in the region outwards the covalent bond DAX, and it 386
Journal of Computational Chemistry 2014, 35, 386–394
DOI: 10.1002/jcc.23507
[a] A. Forni, S. Pieraccini, M. Sironi Istituto di Scienze e Tecnologie Molecolari (ISTM), CNR, and INSTM UdR, Via Golgi 19, Milano 20133, Italy E-mail: [email protected] (A. Forni) E-mail: [email protected] (S. Pieraccini) E-mail: [email protected] (M. Sironi) [b] S. Pieraccini, S. Rendine, M. Sironi Dipartimento di Chimica and INSTM UdR, Universit a degli Studi di Milano, Via Golgi 19, Milano 20133, Italy Contract/grant sponsor: CINECA Award (N. HP10BSXIH2); Contract/grant sponsor: LISA 2012 Grant; Contract/grant sponsor: The Fondazione Banca del Monte di Lombardia C 2013 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
however, that the electrostatic contribution is only one of the energetic terms giving rise to halogen bonding. Depending on the nature of halogen, moiety D and Lewis base A, other forces, in particular dispersion interactions, can play a significant role, as evidenced by previous studies on the nature of halogen bonding.[21–23] Halogen bonding has been the subject of numerous computational studies on both molecular and macromolecular systems. At molecular level, highly correlated methods (e.g., the Mïller– Plesset approach) and density functional theory (DFT) have been largely used to study this interaction in gas-phase[19,21,24– 27] and in solution using the so-called continuum models of the solvent.[25,26,28,29] Quantum mechanics/molecular mechanics (QM/MM) methods have been used to describe the interaction of halogenated molecules with proteins.[30,31] In this context, it should be observed that a purely MM approach generally fails to correctly describe halogen bonding interactions. In fact, classical force fields are usually based on a set of single point charges centered on each atom of the system and cannot then describe the anisotropy of the ESP around the halogen atom. As we have recently shown,[31] such anisotropy can be recovered through introduction of suitable “pseudo” atoms located in the proximity of the halogen atoms, allowing to correctly reproduce the X-ray structures of complexes of proteins with halogenated ligands. Similar approaches have also been recently proposed by other authors.[32–35] The large body of theoretical calculations carried out up to now have described halogen bonding interactions in a large variety of chemical situations, but, to our knowledge, very little attention has been devoted to the interaction of a halogen atom with a p electron system.[21,36–38] Such interaction has been shown to play a key role in biochemical processes. For example, a systematic study of the interaction of Cl and Br atoms with a p electron system has allowed to identify a potent inhibitor of the serine protease factor Xa (fXa),[39] an important enzyme in the blood coagulation cascade. Such study has permitted to develop a potent bioavailable molecule (AVE-3247) which entered the clinical development.[40] Several surveys of the Cambridge Crystallographic and the PDB databases have evidenced an even increasing number of halogen bonding interactions.[39,41–43] In particular, Matter et al.,[39] analyzing the geometries of Cl(Br)/p interactions, have evidenced a strong preference toward the T-shaped conformation, and Lu et al.[43] have found that 1/3 of the halogen bonding interactions in biomolecules involve an interaction with a p system. In spite of their ubiquity in biological systems, halogen bonds with p systems are generally weaker than those with lone pairs of heteroatoms. Previous theoretical studies on the Cl(Br)/p bonding,[21] performed at MP2 level, demonstrated the importance of the dispersive term in describing this interaction, unlike charge transfer forces which play only a very minor role. For these reasons, an accurate and reliable description of the DX/p interaction would require the use of highly correlated methods. Extensive evidence has been collected over the past years indicating that coupled cluster theory represents the best choice to this aim. In particular, the CCSD(T) implementation, where single and double excitations are taken into
account in an iteratively fashion, whereas triple excitations are perturbatively computed, could be considered as the “gold standard” approach,[44] especially if evaluated with a complete basis set (CBS). However, the CCSD(T) method formally scales as M7, where M is the number of basis functions, resulting in a very CPU consuming approach which strongly limits the dimensions of the system that could be studied. The MP2 approach, which scales more favorably (M5), is known to give an overestimation of the binding energies when treating noncovalent interactions involving p systems.[45–47] Finally, the DFT approach, which scales better than (for nonhybrid functionals) or as (for hybrid functionals) the Hartree–Fock theory, faces with the difficulty to choice, among a plethora of existing functionals, the most correct one for the problem to be tackled. We have recently reported about the DX/p noncovalent bonding between NCX or PhX (X 5 F, Cl, Br, and I) and the aromatic system of benzene.[42] In that study, CCSD(T) binding energies extrapolated at the CBS limit have been reported, together with MP2 and DFT results obtained with a limited number of functionals. Such analysis allowed to evidence the failure of some of the most commonly used functionals (e.g., B3LYP and PBE) in describing the DX/p interaction, though they were demonstrated to be adequate to treat halogen bonds with lone pairs, for example, in QM/MM studies on the interaction of halogenated ligands with proteins.[30,31] Other more recently developed functionals (i.e., B97D, the xB97 family and M062X), conversely, have been found to reproduce the CCSD(T) energy profiles with good accuracy. The purpose of this work is to extend our previous analysis by examining both a much more extensive set of functionals, belonging to different categories as described in the next section, and a greater number of halogenated molecules DX, where D 5 H, F, HCC, and NC and X 5 Cl, Br, giving rise to weak to moderate halogen bonds with the aromatic system of benzene. Such different D groups have been chosen in view to investigate the electronic influence of substituents, characterized by variable electron-withdrawing strengths, on the extension of the r-hole and then on the relative importance of the dispersive versus electrostatic contributions to the interaction. Owing to the large amount of functionals considered in this study, we decided to focus our attention on only chlorine and bromine derivatives, which are known to give rise to intermediate interactions between the stronger ones involving iodine and the weaker or even negligible ones involving fluorine. According to a PDB survey,[43] over 90% of the CX/p interactions in protein-ligand complexes involve chlorine and bromine atoms. Benchmark-quality binding energies and equilibrium distances of chlorine and bromine from the p cloud of benzene, evaluated with the CCSD(T) method at the CBS limit according to different extrapolation schemes, will be also presented, together with the results obtained with MP2 and the recently proposed MP2.X method.[48]
Theoretical Methods Accurate binding energies values for the dimers of DCl or DBr (D 5 H, F, HCC, and NC) with benzene in the T-shaped Journal of Computational Chemistry 2014, 35, 386–394
387
FULL PAPER
WWW.C-CHEM.ORG
geometry have been obtained at the CCSD(T) level of theory at the CBS limit, as described below. Due to the high computational costs of the CCSD(T) method, we have used the monomers at fixed geometry as obtained at the MP2/aug-ccpVTZ level. Such rigid monomers approximation relies on the consideration that the geometry of the monomers in the dimer is only slightly affected by their reciprocal intermolecular interaction. An estimate of such approximation was carried out on all systems at the DFT level (see below) and permitted us to quantify it in 0.01 kcal/mol on average. The CCSD(T)/CBS equilibrium distances were obtained by performing singlepoint CCSD(T) energy calculations extrapolated to the CBS limit on a grid of different distances along the dissociation coordinate. Near the minimum, additional grid points were considered to localize the equilibrium geometry with a preci˚. sion of 0.05 A The extrapolation at the CBS limit of CCSD(T) results can be carried out following different protocols as reported in the literature. In this work, we have adopted three different extrapolation schemes, here briefly described. In the first one, denoted in the following as CCSD(T)/CBS1, correlation energies at the CBS limit were first obtained at the MP2 level, using the two-point extrapolation scheme of Halkier et al.[49] with the aug-cc-pVTZ and aug-cc-pVQZ basis sets. Next, such energies were summed to the interaction energies as obtained at the Hartree–Fock limit evaluated with the aug-cc-pVQZ basis set. The final CCSD(T)/CBS1 values were then obtained by adding to the previous terms a DCCSD(T) correction, estimated as the difference between the extrapolated correlation energies obtained at CCSD(T) and MP2 levels according to the Halkier scheme with the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The other two estimates of the CCSD(T) interaction energies at the CBS limit, denoted as CCSD(T)/CBS2 and CCSD(T)/CBS3, were obtained by extrapolating at the CBS limit, the interaction energies at the MP2 level, using either the two-point extrapolation of Halkier with the aug-cc-pVTZ and aug-ccpVQZ basis sets (CCSD(T)/CBS2 scheme) or the three-point extrapolation of Feller[50] with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets (CCSD(T)/CBS3 scheme). Then, we added the DCCSD(T) correction term evaluated, in both CBS2 and CBS3 schemes, as the difference of the interaction energies obtained at CCSD(T) and MP2 levels using the aug-ccpVTZ basis set. Single-point energy calculations along the dissociation coordinate within the rigid monomer approximation were also performed at MP2 and MP2.X levels of theory. The MP2.X approach has been recently proposed by Hobza and coworkers[51] to overcome the well-known problem of the tendency of MP2 to overestimate the binding energy in complexes involving aromatic moieties, in opposition to an underbinding prediction given by MP3.[52] Hobza proposed to correct the MP2 binding energies extrapolated to the CBS limit with a term given by Cx ½DEðMP3=xÞ 2DEðMP2=xÞ, where DE(MP3/x) and DE(MP2/x) are the binding energies evaluated, respectively, at the MP3 and MP2 levels with basis set x, and Cx is a coefficient depending on the chosen basis set[48] and tending to 0.5 for a CBS To obtain the MP2.X binding energies, 388
Journal of Computational Chemistry 2014, 35, 386–394
we have used the extrapolated MP2 values as obtained with the three-point Feller scheme,[50] and we have computed the correction term using the aug-cc-pVDZ basis set, corresponding to a Cx coefficient equal to 0.52. DFT calculations on the DCl/benzene and DBr/benzene dimers were performed by fully optimizing their geometry in the T-shaped approach with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets. Up to 34 functionals have been used to test their performance in describing the DX/p interaction. To evaluate the effect of the monomers relaxation on the binding energy values, single-point DFT calculations were also performed using the optimized minimum distances and the same fixed geometry of the monomers as adopted in the CCSD(T) scans. As said above, relaxation of the monomers geometry determined an energy gain of only 0.01 kcal/mol on average. In all cases, DFT calculations were performed by using a large pruned integration grid (99 radial shells and 590 angular points per shell) to avoid artifacts associated with numerical integration procedures. The need of using such ultrafine grid was in particular evidenced for meta-generalized gradient approximation (GGA) functionals.[53] The 34 XC functionals considered in the present investigation are reported in the following. They are classified according to the category to which they belong (see Supporting Information for further details), that is, (i) the functionals based on the GGA, BLYP,[54,55] PW91PW91 (or PW91),[56] HCTH407,[57] and B97D[58,59]; (ii) The hybrid GGA functionals (H-GGA), B3LYP,[55,60,61] B3P86,[60,62] O3LYP,[55,61,63] BH and HLYP,[54,55] PBE1PBE (or PBE0),[64,65] mPW1PW91,[56,66] mPW3LYP,[55,66–68] B97-2,[59,69] and B98[59,70]; (iii) The range-separated or long-range corrected GGA functionals (LC-GGA), LC-BLYP,[71] xB97X,[72] xB97XD,[58,73] and CAM-B3LYP[74]; (iv) The meta-GGA functionals (M-GGA), M06L,[75] VSXC,[76] TPSSTPSS (or TPSS),[77] and sHCTH[78]; (v) The hybrid meta-GGA functionals (HM-GGA), BMK,[79] TPSSh,[80] sHCTHhyb,[78] M05,[75] M052X,[81] M06,[82] M062X,[82] and M06HF[83,84]; (vi) The double hybrid GGA functionals (DH-GGA), B2PLYP,[54,55,85] mPW2PLYP,[66,86] B2PLYPD,[87] and mPW2PLYPD[87]; (vii) One among the generalized gradient exchange functionals, TPSSVWN5.[61,77] All CCSD(T), MP2, MP2.X, and DFT binding energies were corrected for basis set superposition error using the Boys–Bernardi[88] approach. All calculations were performed using the Gaussian 09 Rev. C.01[89] suite of programs.
Results and Discussion Weak to moderate halogen bonding interactions, such as those presently investigated, are dominated by electrostatic and dispersive contributions, unlike strong halogen bonds, where charge transfer can significantly contribute to the interaction.[90] The halogen atoms in the molecules here reported, DCl and DBr with D 5 H, HCC, F, and NC, are characterized by different r-holes as induced by substituents having different electron-withdrawing capabilities: the stronger is the electronwithdrawing moiety (in the order, H < HCC < F < NC), the larger and deeper is the r-hole. The interactions they can give rise to will be then characterized by a different balance of dispersion and electrostatic contributions. The features of the WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
r-hole can be pictorially evidenced by looking at the ESP plotted on an electron density isosurface, as already reported by Politzer et al.[14] In Figures 1a and 1b, the two cases of NCCl and NCBr, respectively, are exemplified. It can be clearly observed that the r-hole is more pronounced, as expected, on the more polarizable bromine atom with respect to chlorine atom. Recently, the origin of the halogen bond with “conventional” halogen acceptors (i.e., nitrogen and oxygen atoms) has also been rationalized using the Fukui function, f(r).[27,91] This function describes the variation of the electron [91,92] density q(r) upon changing the number of 1electrons N. @qðrÞ 1 The right derivative, namely f ðrÞ5 @N 5qN11 ðrÞ2qN ðrÞ, evidences the regions where a negative charge added to the system will tend to be distributed, that is, the regions more prone to act as electron acceptors. It has been demonstrated that it can be simply approximated by the square of the density of the LUMO orbital.[93] In Figures 1c and 1d, the Fukui function f 1 ðrÞ is plotted on the electron density isosurface for the NCCl and NCBr molecules, respectively, evidencing, similar to the ESP plot, the greater tendency of bromine with respect to chlorine atom to accept electrons along the CAX direction. The geometrical and energetical features of the interaction of DCl and DBr with the aromatic system of benzene have been first investigated by the CCSD(T) approach. We have determined the equilibrium distance R and the binding energy DE at the CBS limit by considering the dimers in the T-shaped configuration. Three different extrapolation procedures have
FULL PAPER
been used, denoted by CBS1, CBS2, and CBS3, as described in the previous section. The results reported in Tables 1 and 2 for the DCl/p and DBr/p interactions, respectively, clearly indicate a good convergence of both R and DE values as obtained by the three extrapolation procedures. The equilibrium distance of the halogen from the aromatic ring is in fact the same for the three approaches, indicating that a possible difference among them is smaller than the width of the grid used to ˚ ). Also, the binding energies locate the minimum (i.e., 0.05 A obtained with the three extrapolations are very close, with maximum differences of only 0.05 and 0.08 kcal/mol for DCl/p and DBr/p interactions, respectively, both obtained when the D moiety is the fluorine atom. For this reason, from now on we will refer to only one CCSD(T) set of extrapolated results to be used as reference for the other methods, namely CCSD(T)/ CBS1. By comparing the CCSD(T) results obtained for the different dimers, we observe a correlation between the DE values and the magnitude of the r-hole (see Supporting Information Fig. S1 for a plot of the ESP for the eight DX molecules), as already pointed out by other authors.[14,15,19] It is, however, to be evidenced a minor deviation from this trend for FX and NCX molecules, the latter showing a slightly lower maximum in the ESP than the former, albeit the corresponding interaction with benzene (at both CCSD(T) and MP2 levels) is slightly stronger for NCX. Such result could be ascribed to a greater polarization induced on the halogen by the cyano group with respect to
Figure 1. MP2/aug-cc-pVDZ ESP (first line) and Fukui function f1(r) (second line) mapped onto the isosurface of electron density (0.001 electrons/bohr-3) for NCCl (left) and NCBr (right) molecules.
Journal of Computational Chemistry 2014, 35, 386–394
389
FULL PAPER
WWW.C-CHEM.ORG
˚ ), and counterpoise-corrected interaction energies, DE Table 1. Minimum intermolecular distances between Cl and centroid of the benzene ring, R (A (kcal/mol), for dimers of DCl (D 5 H, HCC, F, NC) with benzene in the T-shaped approach.[a] HCl/p
CCSD(T)/aug-cc-pVDZ CCSD(T)/aug-cc-pVTZ CCSD(T)/CBS1 CCSD(T)/CBS2 CCSD(T)/CBS3 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ MP2/aug-cc-pVQZ MP2.X
HCCCl/p
FCl/p
NCCl/p
R
DE
R
DE
R
DE
R
DE
3.60 3.55 3.45 3.45 3.45 3.465 3.385 3.351 3.45
21.05 21.27 21.43 21.45 21.43 21.62 21.93 22.07 21.58
3.50 3.40 3.35 3.35 3.35 3.363 3.273 3.245 3.30
21.95 22.19 22.39 22.41 22.40 22.73 23.05 23.20 22.59
3.40 3.35 3.25 3.25 3.25 3.266 3.204 3.176 3.25
22.47 22.76 23.03 23.05 23.08 23.39 23.60 23.77 23.26
3.40 3.35 3.25 3.25 3.25 3.311 3.208 3.180 3.25
22.80 23.10 23.35 23.39 23.38 23.65 24.02 24.21 23.61
[a] CCSD(T) and MP2.X interaction energies evaluated on a 0.05-A˚ spaced grid.
the fluorine atom when interacting with benzene, determining a larger stabilizing contribution to DE for the NCX molecules. The electrostatic term represents in fact only the first-order approximation to the total interaction energy, but high-order terms may also play an important role. As for the CCSD(T) equilibrium distances R of the halogen from the centroid of the aromatic ring, we note that they tend to become smaller with increasing the magnitude of the r-hole, with the only exception of the FX/p system, which shows the same or even shorter R values than the stronger NCX/p system. Anyway, the equilibrium distances are always shorter than the sum of the van der Waal radius of the halogen atom[94] and the halfthickness of the aromatic ring,[95,96] corresponding to 3.5 and ˚ for the chlorine and the bromine atom, respectively. 3.6 A The MP2 results, also reported in Tables 1 and 2, show a general overestimate of the binding energies and an underestimate of the equilibrium distances. Both overestimate and underestimate tend to increase with the basis set dimensions, giving the unpleasant result that the smaller basis set, aug-ccpVDZ, happens to provide the better agreement with the extrapolated CCSD(T) results. This finding, already evidenced in our previous investigation on the DX/p interaction,[42] is a well-known problem of the MP2 methods in describing noncovalent interactions involving p systems.[45–47] An opposite tendency is usually shown by the MP3 method, and for this reason, Hobza and coworkers have recently proposed the
MP2.X method,[51] in which the results of the two methods are opportunely mixed through a weight, Cx, depending on the basis set quality, as reported in the Theoretical Methods section. The weight Cx had been determined to better reproduce the binding energies of the S66 set,[97] which does not contain any halogen bonded system. The MP2.X results, reported in Tables 1 and 2, show a significant improvement with respect to the MP2 method. The deviations from the CCSD(T)/CBS1 results are quite small, giving percentage RMS errors of 8.6 and 10.0% for the chlorine and the bromine halogen bonded complexes, respectively. Such errors are, however, slightly greater than the corresponding values reported by Hobza for dispersionbound, mixed dispersion-electrostatic, and electrostatic-bound complexes (about 7, 3, and 1%, respectively),[48] denoting a special character for halogen bonding and in particular for the DX/p interaction. DFT binding energies and equilibrium distances between the halogen atom and the centroid of the aromatic ring are reported as Supporting Information Tables S1–S4. Comparison of the DFT binding energies and equilibrium distances with the CCSD(T) reference results are graphically reported in Figures 2 and 3, respectively, where the ratios of the DFT/aug-ccpVQZ over the CCSD(T)/CBS1 values, averaged separately on the chlorine and bromine complexes, are plotted. Furthermore, the mean signed error, the mean absolute error (MAE), the mean signed percentage error, and the mean absolute
˚ ), and counterpoise-corrected interaction energies, DE Table 2. Minimum intermolecular distances between Br and centroid of the benzene ring, R (A (kcal/mol), for dimers of DBr (D 5 H, HCC, F, NC) with benzene in the T-shaped approach.[a] HBr/p
CCSD(T)/aug-cc-pVDZ CCSD(T)/aug-cc-pVTZ CCSD(T)/CBS1 CCSD(T)/CBS2 CCSD(T)/CBS3 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ MP2/aug-cc-pVQZ MP2.X
HCCBr/p
FBr/p
NCBr/p
R
DE
R
DE
R
DE
R
DE
3.65 3.60 3.50 3.50 3.50 3.524 3.438 3.402 3.50
21.56 21.83 22.03 22.06 22.04 22.30 22.70 22.87 22.27
3.55 3.50 3.40 3.40 3.40 3.424 3.358 3.314 3.40
22.45 22.71 22.96 22.99 23.00 23.39 23.73 23.92 23.26
3.45 3.40 3.30 3.30 3.30 3.334 3.260 3.260 3.30
23.16 23.51 23.89 23.94 23.97 24.23 24.57 24.82 24.24
3.45 3.40 3.35 3.35 3.35 3.369 3.280 3.251 3.35
23.36 23.69 24.00 24.04 24.05 24.37 24.78 25.01 24.35
[a] CCSD(T) and MP2.X interaction energies evaluated on a 0.05-A˚ spaced grid.
390
Journal of Computational Chemistry 2014, 35, 386–394
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 2. Ratios of DFT/aug-cc-pVQZ binding energies over the CCSD(T)/CBS1 values, averaged on the DCl/p, blue, and DBr/p, red (D 5 H, F, HCC, and NC) systems.
percentage error (MAPE) computed separately for the four chlorine and the four bromine complexes with the three considered basis sets, are reported as Supporting Information Tables S5–S8. The binding energies reported in Supporting Information Tables S1 and S3 for chlorine and bromine complexes, respectively, indicate that all functionals generally reproduce the trend provided by the CCSD(T) calculations, with increasing DE values on increasing the strength of the electron-withdrawing group and, for a same D substituent, on going from chlorine to bromine. There are, however, some remarkable exceptions, as denoted by those functionals (BLYP, B3LYP, B3P86, BH and HLYP, CAM-B3LYP, and BMK) which are not able to demonstrate the existence of the weaker HCl/p or even HBr/p halogen bonded complexes, giving repulsive curves at all intermolecular distances (such functionals are obviously excluded when calculating the average ratios plotted in Figures 2 and 3 and
the mean errors reported in Supporting Information Tables S5–S8). It also stands out the too large DE overestimates given by the VSXC functional. Finally, it is noteworthy the case of a discrete number of functionals, such as for example B97D, xB97X, and almost all the functionals of the Minnesota collection, which show a mismatch in the stability order on going from the FBr/p to the NCBr/p interaction, though they provide the correct order of stability for the corresponding chlorine complexes. The trend of the CCSD(T) equilibrium distances is instead correctly reproduced by all the functionals with only a few exceptions (see Supporting Information Tables S2 and S4 for chlorine and bromine complexes, respectively). About the dependence of the DFT results on the basis set, it is evident from Supporting Information Tables S1–S4 that binding energies usually decrease with increasing the dimensions of the basis set for almost all the functionals and that converged results are obtained with the aug-cc-pVTZ basis set.
Figure 3. Ratios of DFT/aug-cc-pVQZ equilibrium distances over the CCSD(T)/CBS1 values, averaged on the D-Cl/p, blue, and D-Br/p, red (D 5 H, F, HCC, and NC) systems. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Journal of Computational Chemistry 2014, 35, 386–394
391
FULL PAPER
WWW.C-CHEM.ORG
Equilibrium distances show a more variable behavior with increasing the dimensions of the basis set, as both a contraction and an elongation of the distance is observed depending not only on the functional but also on the system studied. Variations are more significant for the weaker HX/p interactions, but for most functionals it can be observed that “converged” results are obtained using the aug-cc-PVTZ basis set. The performance of the different examined functionals is summarized by the ratios plotted in Figures 2 and 3, where the closer is the ratio to 1, the better is the functional. From Figure 2, it is immediate to see that all the considered functionals, except VSXC, underestimate the binding energy. The large overestimate given by VSXC, which was parametrized without including noncovalent interactions,[76] was already evidenced in the investigation by Johnson et al. on purely dispersion-bound complexes,[53] in the study of the hydrogen bonded benzene-H2O system by Neves et al.[98] and in the systematic analysis by Zhao and Truhlar[82] including extensive sets of intermolecular interactions, in particular those dominated by dispersion contributions. It is also evident the better performance of almost all functionals in describing the DBr/p with respect to the weaker DCl/p interaction. The closest agreement with the CCSD(T)/CBS1 binding energy results, provided by functionals giving DE ratio > 0.80 for both chlorine and bromine complexes, is given by the functionals xB97X, M062X, B2PLYPD, and mPW2PLYP2D (MAPE values below or just above 10%, see Supporting Information Tables S5 and S7). If limiting to bromine complexes, B97D and xB97XD can be also included within the better performing functionals. It is interesting to note that all the considered functionals explicitly including dispersion, that is, B97D, xB97XD, B2PLYPD, and mPW2PLYP2D, fall within the subset of the better performing functionals, confirming the importance of the dispersion forces, besides electrostatic contributions, in determining the DX/p interaction. In particular, comparison of B2PLYPD and mPW2PLYP2D with the corresponding forms without dispersion, B2PLYP and mPW2PLYP2, respectively, indicates the substantial improving implied by the inclusion of the dispersion term (compare, e.g., the MAPEs for mPW2PLYP2D and mPW2PLYP2, below 10% and above 38%, respectively, for all basis sets). The good performance of the other two functionals, xB97X and M062X, can be ascribed to the inclusion of dispersion-bound complexes in the training sets used for the functional parametrization. It is worth noting the slightly better performance of xB97X with respect to xB97XD, a result which generalizes what previously reported in the study of the interaction of NCX and PhX with benzene,[42] where also iodine, besides chlorine and bromine, had been included in the analysis. It can be probably ascribed to the choice of the database used for parametrizing these functionals, where compounds exhibiting DX/p halogen bonding were not included. By looking at the performance of different functionals categories tested in this work, we cannot note a clear trend on going from the less accurate pure GGA functionals to the most accurate double hybrid GGA functionals. When climbing the Jacob’s ladder of DFT accuracy,[99] we find a large variability within each rung. The main origin of such variability is ascribable to 392
Journal of Computational Chemistry 2014, 35, 386–394
inclusion of the additive dispersion term into the functionals, as discussed above. Within such class of dispersion-corrected functionals, we note only a slight improvement for the double hybrid B2PLYPD and mPW2PLYP2D functionals with respect to the pure GGA (B97D) and the hybrid LC-GGA (xB97XD) functionals, according to the Jacob’s ladder. Moreover, almost uniform accuracy of the double hybrid functionals for chlorine and bromine complexes makes them more suitable when more halogen bonds are present in the investigated system. When excluding dispersion-corrected functionals, we note that pure GGA functionals, including those developed within the meta-approximation, cannot be used to treat the DX/p noncovalent interaction, though HCTH407 and M06 show a definite better performance with respect to the other functionals on this rung. Introduction of exact exchange, giving rise to H-, LC-, and HM-GGA functionals, does not involve, by itself, an improvement with respect to pure GGA functionals, however, it is evident a strong dependence on the family of functionals under consideration. For example, all the H-GGA functionals, including the popular B3LYP, are confirmed to be totally unable to describe long-range dispersive interactions (MAPEs above 60% for all basis sets). Conversely, xB97X and the Minnesota collection Myz provide in general a good description of the DX/p interaction, with M062X on the top of the hybrid functionals category. It is interesting to observe, in the Myz family, a systematic improvement of the performance when increasing the amount of exact exchange (M052X and M062X are better than M05 and M06, respectively). However, too much exact exchange appears to be detrimental, as can be deduced by the performance of M06HF, which is comparable with that of M06. Such trend matches that obtained for halogen bonding involving heteroatoms,[27] though in that case it had been demonstrated that the role of dispersion contribution is in general less important than what deduced for the DX/p halogen bonding. Finally, the double hybrids without dispersion, B2PLYP and mPW2PLYP, do not outperform the xB97X and the Myz collection, indicating that introduction of MP2 correlation, besides HF exchange, does not warrant, by itself, a better performance than that of lower rungs on Jacob’s ladder. Looking at Figure 3 and Supporting Information Tables S6 and S8, it appears that geometry is well described by a greater number of functionals, with respect to binding energy, for both chlorine and bromine complexes. In particular, besides the functionals explicitly including dispersion, all the Minnesota collection provides an excellent description of the minimum geometry. For these functionals, the MAE is less or only slightly greater than 0.05 A˚, which corresponds to the grid spacing used for the determination of the minimum distance at CCSD(T) level. The better agreement is provided once again by the B2PLYPD and mPW2PLYPD functionals, and the importance of including the dispersion contribution is further assessed.
Conclusions We have studied the DX/p halogen bonding between different halogenated molecules (D 5 H, HCC, F, and NC; X 5 Cl, Br) with WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
the aromatic system of benzene by CCSD(T), MP2, MP2.X, and DFT approaches, evaluating the minimum geometry and the binding energy with monomers in the T-shaped configuration. CCSD(T) calculations at the CBS limit have been performed to provide a reference for a plethora of DFT functionals and evaluate their performance in describing the present interaction. The results here obtained can be used for addressing the correct choice of the functional when studying more extended systems, for example, in the context of QM/MM calculations where more sophisticated post-HF methods are prohibitive. Results obtained at MP2 level confirm the previously reported tendency to overestimate the strength of the interaction when delocalized p systems are involved.[45–47] Such failure is adequately recovered by applying the recently proposed MP2.X approach,[51] though for the presently investigated DX/ p halogen bonding the improvement with respect to MP2 was not as good as that reported in the original paper.[51] The performance of up to 34 functionals, belonging to several categories of approximations, has been assessed by comparison with CCSD(T) extrapolated results. Such investigation allows to conclude that, owing to the mixed dispersive and electrostatic nature of the DX/p halogen bonding, a functional which correctly describes the dispersion contribution is mandatory to obtain good results. The double hybrid functionals B2PLYPD and mPW2PLYPD are here shown to provide the better results on both energies and equilibrium geometry. The corresponding forms without dispersion, B2PLYP and mPW2PLYP, fail to describe this interaction. Among the less expensive functionals, the better performance is provided by the xB97X and M062X functionals, while the xB97XD and B97D functionals are shown to work very well for bromine complexes but not so well for chlorine complexes. The worse performance is instead shown by the hybrid GGA functionals, including in particular the popular B3LYP functional whose use is then strongly discouraged for studying DX/p halogen bonded complexes. Keywords: halogen bond density functional theory benchmark calculations
How to cite this article: A. Forni, S. Pieraccini, S. Rendine, M. Sironi. J. Comput. Chem. 2014, 35, 386–394. DOI: 10.1002/ jcc.23507
]
Additional Supporting Information may be found in the online version of this article.
[1] P. Metrangolo, G. Resnati, Science 2008, 321, 918. [2] P. Metrangolo, G. Resnati, Halogen Bonding: Fundamentals and Applications; Springer: Berlin, 2008. [3] R. Kato, T. Imakubo, H. Yamamoto, R. Maeda, M. Fujiwara, H. Sawa, Mol. Cryst. Liq. Cryst. 2002, 380, 61. [4] E. Cariati, A. Forni, S. Biella, P. Metrangolo, F. Meyer, G. Resnati, S. Righetto, E. Tordin, R. Ugo, Chem. Commun. 2007, 2590. [5] J. Hulliger, P. J. Langley, Chem. Commun. 1998, 2557. [6] J. A. R. P. Sarma, F. H. Allen, V. J. Hoy, J. A. K. Howard, R. Thaimattam, K. Biradha, G. R. Desiraju, Chem. Commun. 1997, 101.
FULL PAPER
[7] P. K. Thallapally, G. R. Desiraju, M. Bagieu-Beucher, R. Masse, C. Bourgogne, J. F. Nicoud, Chem. Commun. 2002, 1052. [8] M. Adler, M. J. Kochanny, B. Ye, G. Rumennik, D. R. Light, S. Biancalana, M. Whitlow, Biochemistry 2002, 41, 15514. [9] E. Howard, R. Sanishvili, R. E. Cachau, B. Mitschler, P. Bart, V. Lamour, M. Van Zandt, E. Sibley, C. Bon, D. Moras, T. R. Schneider, A. Joachimiak, A. Podjarny, Proteins: Struct. Funct. Genet. 2004, 55, 792. [10] Y. Jiang, A. A. Alcaraz, J. M. Chen, H. Kobayashi, Y. J. Lu, J. P. Snyder, J. Med. Chem. 2006, 49, 1891. [11] P. Auffinger, F. A. Hays, E. Westhof, P. S. Ho, Proc. Natl. Acad. Sci. USA 2004, 101, 16789. [12] V. G. Tsirelson, P. F. Zou, T.-H. Tang, R. F. W. Bader, Acta Cryst. 1995, A51, 143. [13] P. Metrangolo, H. Neukirch, T. Pilati, G. Resnati, Acc. Chem. Res. 2005, 38, 386. [14] P. Politzer, J. S. Murray, T. Clark, Phys. Chem. Chem. Phys. 2010, 12, 7748. [15] P. Politzer, P. Lane, M. C. Concha, Y. Ma, S. Murray, J. Mol. Model. 2007, 13, 305. [16] D. Chopra, Cryst. Growth Des. 2012, 12, 541. [17] P. Metrangolo, J. S. Murray, T. Pilati, P. Politzer, G. Resnati, G. Terraneo, CrystEngComm. 2011, 13, 6593. [18] P. Metrangolo, J. S. Murray, T. Pilati, P. Politzer, G. Resnati, G. Terraneo, Cryst. Growth Des. 2011, 11, 4238. [19] K. E. Riley, J. S. Murray, P. Politzer, M. C. Concha, P. Hobza, J. Chem. Theory Comput. 2009, 5, 155. acˇ, R. J. Sola, M. C. Concha, [20] K. E. Riley, J. S. Murray, J. Fanfrlik, J. Rez F. M. Ramos, P. Politzer, J. Mol. Model. 2011, 17, 3309. [21] Y. X. Lu, J. W. Zou, Y. H. Wang, Q. S. Yu, Chem. Phys. 2007, 334, 1. [22] Y. X. Lu, J. W. Zou, Y. H. Wang, Y. J. Jiang, Q. S. Yu, J. Phys. Chem. A 2007, 111, 10781. acˇ, R. J. Sola, M. C. Concha, [23] K. E. Riley, J. S. Murray, J. Fanfrlik, J. Rez F. M. Ramos, P. Politzer, J. Mol. Model. 2013, 19, 4651. [24] Y. X. Lu, J. W. Zou, J. C. Fan, W. N. Zhao, Y. J. Jiang, Q. S. Yu, J. Comput. Chem. 2009, 30, 725. [25] A. Forni, S. Rendine, S. Pieraccini, M. Sironi, J. Mol. Graph. Model. 2012, 38, 31. [26] K. E. Riley, K. M. Merz, J. Phys. Chem. A 2007, 111, 1688. [27] S. Kozuch, J. M. L. Martin, J. Chem. Theory Comput. 2013, 9, 1918. [28] Y. Lu, H. Li, X. Zhu, W. Zhu, H. Liu, J. Phys. Chem. A 2011, 115, 4467. [29] Y. Lu, H. Li, X. Zhu, H. Liu, W. Zhu, Int. J. Quantum Chem. 2012, 112, 1421. [30] Y. Lu, T. Shi, Y. Wang, H. Yang, X. Yan, X. Luo, H. Jiang, W. Zhu, J. Med. Chem. 2009, 52, 2854. [31] S. Rendine, S. Pieraccini, A. Forni, M. Sironi, Phys. Chem. Chem. Phys. 2011, 13, 19508. [32] M. A. A. Ibrahim, J. Comput. Chem. 2011, 32, 2564. [33] M. Kolar, P. Hobza, J. Chem. Theory Comput. 2012, 8, 1325. [34] M. Carter, A. K. Rapp e, P. S. Ho, J. Chem. Theory Comput. 2012, 8, 2461. [35] W. L. Jorgensen, P. Schyman, J. Chem. Theory Comput. 2012, 8, 3895. [36] O. V. Shishkin, Chem. Phys. Lett. 2008, 458, 96. [37] D. J. R. Duarte, M. M. De Las Vallejos, N. M. Peruchena, J. Mol. Model. 2010, 16, 737. [38] H. G. Wallnoefer, T. Fox, K. R. Liedl, C. S. Tautermann, Phys. Chem. Chem. Phys. 2010, 12, 14941. [39] H. Matter, M. Nazar, S. Gussregen, D. W. Will, H. Schreuder, A. Bauer, Angew. Chem. Int. Ed. 2009, 48, 2911. [40] M. Nazare, M. Essrich, D. W. Will, H. Matter, K. Ritter, V. Wehner, PCT Int. Appl. EP 1314733, 2003. [41] C. Bissantz, B. Kuhn, M. Stahl, J. Med. Chem. 2010, 53, 5061. [42] A. Forni, S. Pieraccini, S. Rendine, F. Gabas, M. Sironi, ChemPhysChem 2012, 13, 4224. [43] Y. Lu, Y. Wang, W. Zhu, Phys. Chem. Chem. Phys. 2010, 12, 4543. [44] K. E. Riley, M. Pitonak, P. Jurecka, P. Hobza, Chem. Rev. 2010, 110, 5023. [45] C. D. Sherrill, Rev. Comput. Chem. 2009, 26, 1. [46] G. S. Tschumper, Rev. Comput. Chem. 2009, 26, 39. [47] M. O. Sinnokrot, E. F. Valeev, C. D. Sherrill, J. Am. Chem. Soc. 2002, 124, 10887. [48] K. E. Riley, J. Rezac, P. Hobza, Phys. Chem. Chem. Phys. 2011, 13, 21121. [49] A. Halkier, W. Klopper, T. Helgaker, P. Jorgensen, P. R. Taylor, J. Chem. Phys. 1999, 111, 9157. [50] D. Feller, J. Chem. Phys. 1993, 98, 7059.
Journal of Computational Chemistry 2014, 35, 386–394
393
FULL PAPER
WWW.C-CHEM.ORG
acˇ, M. Pito [51] R. Sedlak, K. E. Riley, J. Rez nak, P. Hobza, ChemPhysChem 2013, 14, 698. [52] L. Grafova, M. Pitonak, J. Rezac, P. Hobza, J. Chem. Theory Comput. 2010, 6, 2365. [53] E. R. Johnson, R. A. Wolkow, G. A. DiLabio, Chem. Phys. Lett. 2004, 394, 334. [54] A. D. Becke, Phys. Rev. A 1988, 38, 3098. [55] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. [56] J. P. Perdew, In Electronic Structure of Solids; P. Zieche, H. Eschrig, Eds.; Akademie: Berlin, Germany, 1991. [57] A. D. Boese, N. C. Handy, J. Chem. Phys. 2001, 114, 5497. [58] S. Grimme, J. Comput. Chem. 2006, 27, 1787. [59] A. D. Becke, J. Chem. Phys. 1997, 107, 8554. [60] A. D. Becke, J. Chem. Phys. 1993, 98, 5648. [61] S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. [62] J. P. Perdew, Phys. Rev. B 1986, 33, 8822. [63] A. J. Cohen, N. C. Handy, Mol. Phys. 2001, 99, 607. [64] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. [65] C. Adamo, V. Barone, J. Chem. Phys. 1999, 110, 6158. [66] C. Adamo, V. Barone, J. Chem. Phys. 1998, 108, 664. [67] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2004, 108, 6908. [68] X. Xu, W. A. Goddard, III, Proc. Natl. Acad. Sci. USA 2004, 101, 2673. [69] P. J. Wilson, T. J. Bradley, D. J. Tozer, J. Chem. Phys. 2001, 115, 9233. [70] H. L. Schmider, A. D. Becke, J. Chem. Phys. 1998, 108, 9624. [71] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540. [72] J. D. Chai, M. Head-Gordon, J. Chem. Phys. 2008, 128, 084106. [73] J. D. Chai, M. Head-Gordon, Phys. Chem. Chem. Phys. 2008, 10, 6615. [74] T. Yanai, D. Tew, N. C. Handy, Chem. Phys. Lett. 2004, 393, 51. [75] Y. Zhao, D. G. Truhlar, J. Chem. Phys. 2005, 123, 194101. [76] T. Van Voorhis, G. E. Scuseria, J. Chem. Phys. 1998, 109, 400. [77] J. M. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Phys. Rev. Lett. 2003, 91, 146401. [78] A. D. Boese, N. C. Handy, J. Chem. Phys. 2002, 116, 9559. [79] A. D. Boese, J. M. L. Martin, J. Chem. Phys. 2004, 121, 3405. [80] V. N. Staroverov, G. E. Scuseria, J. Tao, J. P. Perdew, J. Chem. Phys. 2003, 119, 12129. [81] Y. Zhao, D. G. Truhlar, J. Chem. Theory Comput. 2006, 2, 364. [82] Y. Zhao, D. G. Truhlar, Theor. Chem. Acc. 2008, 120, 215. [83] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 5121.
394
Journal of Computational Chemistry 2014, 35, 386–394
[84] [85] [86] [87] [88] [89]
[90] [91] [92] [93] [94] [95] [96] [97] [98] [99]
Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 13126. S. Grimme, J. Chem. Phys. 2006, 124, 034108. T. Schwabe, S. Grimme, Phys. Chem. Chem. Phys. 2006, 8, 4398. T. Schwabe, S. Grimme, Phys. Chem. Chem. Phys. 2007, 9, 3397. S. F. Boys, F. Bernardi, Mol. Phys. 1970, 19, 553. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision C.01; Gaussian, Inc., Wallingford CT, 2009. O. Hassel, Science 1970, 170, 497. B. Pinter, N. Nagels, W. A. Herrebout, F. D. Proft, Chem. Eur. J. 2013, 19, 519. R. G. Parr, W. Yang, J. Am. Chem. Soc. 1984, 106, 4049. W. Yang, R. G. Parr, R. Pucci, J. Chem. Phys. 1984, 81, 2862. A. Bondi, J. Phys. Chem 1964, 68, 441. A. W. Hanson, Acta Crystallogr. 1964, 17, 559. I. Saraogi, V. G. Vijay, S. Das, K. Sekar, T. N. G. Row, Cryst. Eng. 2003, 6, 69. acˇ, K. E. Riley, P. Hobza, J. Chem. Theory Comput. 2011, 7, 2427. J. Rez A. R. Neves, P. A. Fernandes, M. J. Ramos, J. Chem. Theory Comput. 2011, 7, 2059. J. P. Perdew, K. Schmidt, AIP Conf. Proc. 2001, 577, 1.
Received: 17 August 2013 Revised: 29 October 2013 Accepted: 21 November 2013 Published online on 13 December 2013
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Halogen Bonds with Benzene: An Assessment of DFT Functionals Alessandra Forni,*[a] Stefano Pieraccini,*[a,b] Stefano Rendine,[b] and Maurizio Sironi*[a,b] The performance of an extensive set of density functional theory functionals has been tested against CCSD(T) and MP2 results, extrapolated to the complete basis set (CBS) limit, for the interaction of either DCl or DBr (D 5 H, HCC, F, and NC) with the aromatic system of benzene. It was found that double hybrid functionals explicitly including dispersion, that is, B2PLYPD and mPW2PLYPD, provide the better agreement with the CCSD(T)/CBS results on both energies and equilibrium geometry, indicating the importance of dispersive contribu-
tions in determining this interaction. Among the less expensive functionals, the better performance is provided by the xB97X and M062X functionals, while the xB97XD and B97D functionals are shown to work very well for bromine comC 2013 Wiley plexes but not so well for chlorine complexes. V Periodicals, Inc.
Introduction
is compensated by an increase of electronic charge density in a belt around the bond axis. The site of electronic depletion can hence interact with a source of electrons represented by the Lewis base A located along the direction of the DAX bond. This pictorial representation has been confirmed by quantum mechanical calculations using Natural Bond Orbitals analysis and by looking at the details of the electrostatic potential (ESP) around the halogen atom. The ESP shows in fact the presence of a positive region along the extension of the DAX bond, denoted as “r-hole,”[15] and the formation of a negative belt perpendicular to the DAX bond. The essential features of such r-hole, that is, the maximum value of the ESP and the size of the positive region, strongly depend, of course, on the moiety D and on the nature of the halogen atom. In particular, the r-hole becomes more pronounced as the electron-withdrawing capability of D and the polarizability of the halogen atom increase. As a result, for a same D group, the iodine atom shows a well-defined r-hole, which decreases on going to bromine and then to chlorine. The r-hole is not observed on fluorine, unless very strong electronegative moieties are present in the molecule.[16–18] A good correlation has been observed between the values of the ESP on the r-hole and the strength of the halogen bond interaction.[14,15,19,20] The presence of the r-hole explains also the directionality of halogen bonding. It is expected,
In recent years, halogen bonding[1,2] has received a growing attention as an important type of interaction involved in many different fields ranging from material science to biochemistry. Halogen bonding has been shown to play a key role in the development of new solid-state materials with attractive properties, like for example organic semiconductors and superconductors,[3] and materials for nonlinear optical applications.[4–7] In biochemical applications, halogen bonding has offered interesting perspectives in protein-ligand complexation with important consequences in medicinal chemistry.[8–11] Halogen bond is a highly directional interaction between a covalently bonded halogen atom X and an acceptor group A with nucleophilic character, according to the scheme DX/A. The moiety D bonded to the halogen atom has a large variability, ranging from inorganic to organic species, whereas the nucleophilic site A is usually represented by a lone pair of a heteroatom such as oxygen, nitrogen, sulfur, or by a p system as, for example, a phenyl group. Of course, the strength of halogen bonding depends on the halogen atom and on the nature of both the moieties D and A. A broad range of halogen bond interaction strengths have been reported, from about 1.6 to 43.0 kcal/mol for the ClCl interaction in Cl2[12] and the I2I2 interaction in I32,[13] respectively. The physical origin of such quite unexpected phenomenon, where an electronegative species X interacts with a nucleophilic site A, has been explained by some authors[14,15] on the basis of a simple intuitive model. The electronic structure of a halogen atom covalently bonded to another atom can be roughly described as s2px2py2pz1 assuming z as the direction of the bond. This electronic arrangement is associated with a depletion of the electron charge distribution along the bond axis with respect to a 5/3 mean population of p electrons in each direction, as found in a (spherical) isolated halogen atom. More precisely, such reduction of charge distribution is predicted in the region outwards the covalent bond DAX, and it 386
Journal of Computational Chemistry 2014, 35, 386–394
DOI: 10.1002/jcc.23507
[a] A. Forni, S. Pieraccini, M. Sironi Istituto di Scienze e Tecnologie Molecolari (ISTM), CNR, and INSTM UdR, Via Golgi 19, Milano 20133, Italy E-mail: [email protected] (A. Forni) E-mail: [email protected] (S. Pieraccini) E-mail: [email protected] (M. Sironi) [b] S. Pieraccini, S. Rendine, M. Sironi Dipartimento di Chimica and INSTM UdR, Universit a degli Studi di Milano, Via Golgi 19, Milano 20133, Italy Contract/grant sponsor: CINECA Award (N. HP10BSXIH2); Contract/grant sponsor: LISA 2012 Grant; Contract/grant sponsor: The Fondazione Banca del Monte di Lombardia C 2013 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
however, that the electrostatic contribution is only one of the energetic terms giving rise to halogen bonding. Depending on the nature of halogen, moiety D and Lewis base A, other forces, in particular dispersion interactions, can play a significant role, as evidenced by previous studies on the nature of halogen bonding.[21–23] Halogen bonding has been the subject of numerous computational studies on both molecular and macromolecular systems. At molecular level, highly correlated methods (e.g., the Mïller– Plesset approach) and density functional theory (DFT) have been largely used to study this interaction in gas-phase[19,21,24– 27] and in solution using the so-called continuum models of the solvent.[25,26,28,29] Quantum mechanics/molecular mechanics (QM/MM) methods have been used to describe the interaction of halogenated molecules with proteins.[30,31] In this context, it should be observed that a purely MM approach generally fails to correctly describe halogen bonding interactions. In fact, classical force fields are usually based on a set of single point charges centered on each atom of the system and cannot then describe the anisotropy of the ESP around the halogen atom. As we have recently shown,[31] such anisotropy can be recovered through introduction of suitable “pseudo” atoms located in the proximity of the halogen atoms, allowing to correctly reproduce the X-ray structures of complexes of proteins with halogenated ligands. Similar approaches have also been recently proposed by other authors.[32–35] The large body of theoretical calculations carried out up to now have described halogen bonding interactions in a large variety of chemical situations, but, to our knowledge, very little attention has been devoted to the interaction of a halogen atom with a p electron system.[21,36–38] Such interaction has been shown to play a key role in biochemical processes. For example, a systematic study of the interaction of Cl and Br atoms with a p electron system has allowed to identify a potent inhibitor of the serine protease factor Xa (fXa),[39] an important enzyme in the blood coagulation cascade. Such study has permitted to develop a potent bioavailable molecule (AVE-3247) which entered the clinical development.[40] Several surveys of the Cambridge Crystallographic and the PDB databases have evidenced an even increasing number of halogen bonding interactions.[39,41–43] In particular, Matter et al.,[39] analyzing the geometries of Cl(Br)/p interactions, have evidenced a strong preference toward the T-shaped conformation, and Lu et al.[43] have found that 1/3 of the halogen bonding interactions in biomolecules involve an interaction with a p system. In spite of their ubiquity in biological systems, halogen bonds with p systems are generally weaker than those with lone pairs of heteroatoms. Previous theoretical studies on the Cl(Br)/p bonding,[21] performed at MP2 level, demonstrated the importance of the dispersive term in describing this interaction, unlike charge transfer forces which play only a very minor role. For these reasons, an accurate and reliable description of the DX/p interaction would require the use of highly correlated methods. Extensive evidence has been collected over the past years indicating that coupled cluster theory represents the best choice to this aim. In particular, the CCSD(T) implementation, where single and double excitations are taken into
account in an iteratively fashion, whereas triple excitations are perturbatively computed, could be considered as the “gold standard” approach,[44] especially if evaluated with a complete basis set (CBS). However, the CCSD(T) method formally scales as M7, where M is the number of basis functions, resulting in a very CPU consuming approach which strongly limits the dimensions of the system that could be studied. The MP2 approach, which scales more favorably (M5), is known to give an overestimation of the binding energies when treating noncovalent interactions involving p systems.[45–47] Finally, the DFT approach, which scales better than (for nonhybrid functionals) or as (for hybrid functionals) the Hartree–Fock theory, faces with the difficulty to choice, among a plethora of existing functionals, the most correct one for the problem to be tackled. We have recently reported about the DX/p noncovalent bonding between NCX or PhX (X 5 F, Cl, Br, and I) and the aromatic system of benzene.[42] In that study, CCSD(T) binding energies extrapolated at the CBS limit have been reported, together with MP2 and DFT results obtained with a limited number of functionals. Such analysis allowed to evidence the failure of some of the most commonly used functionals (e.g., B3LYP and PBE) in describing the DX/p interaction, though they were demonstrated to be adequate to treat halogen bonds with lone pairs, for example, in QM/MM studies on the interaction of halogenated ligands with proteins.[30,31] Other more recently developed functionals (i.e., B97D, the xB97 family and M062X), conversely, have been found to reproduce the CCSD(T) energy profiles with good accuracy. The purpose of this work is to extend our previous analysis by examining both a much more extensive set of functionals, belonging to different categories as described in the next section, and a greater number of halogenated molecules DX, where D 5 H, F, HCC, and NC and X 5 Cl, Br, giving rise to weak to moderate halogen bonds with the aromatic system of benzene. Such different D groups have been chosen in view to investigate the electronic influence of substituents, characterized by variable electron-withdrawing strengths, on the extension of the r-hole and then on the relative importance of the dispersive versus electrostatic contributions to the interaction. Owing to the large amount of functionals considered in this study, we decided to focus our attention on only chlorine and bromine derivatives, which are known to give rise to intermediate interactions between the stronger ones involving iodine and the weaker or even negligible ones involving fluorine. According to a PDB survey,[43] over 90% of the CX/p interactions in protein-ligand complexes involve chlorine and bromine atoms. Benchmark-quality binding energies and equilibrium distances of chlorine and bromine from the p cloud of benzene, evaluated with the CCSD(T) method at the CBS limit according to different extrapolation schemes, will be also presented, together with the results obtained with MP2 and the recently proposed MP2.X method.[48]
Theoretical Methods Accurate binding energies values for the dimers of DCl or DBr (D 5 H, F, HCC, and NC) with benzene in the T-shaped Journal of Computational Chemistry 2014, 35, 386–394
387
FULL PAPER
WWW.C-CHEM.ORG
geometry have been obtained at the CCSD(T) level of theory at the CBS limit, as described below. Due to the high computational costs of the CCSD(T) method, we have used the monomers at fixed geometry as obtained at the MP2/aug-ccpVTZ level. Such rigid monomers approximation relies on the consideration that the geometry of the monomers in the dimer is only slightly affected by their reciprocal intermolecular interaction. An estimate of such approximation was carried out on all systems at the DFT level (see below) and permitted us to quantify it in 0.01 kcal/mol on average. The CCSD(T)/CBS equilibrium distances were obtained by performing singlepoint CCSD(T) energy calculations extrapolated to the CBS limit on a grid of different distances along the dissociation coordinate. Near the minimum, additional grid points were considered to localize the equilibrium geometry with a preci˚. sion of 0.05 A The extrapolation at the CBS limit of CCSD(T) results can be carried out following different protocols as reported in the literature. In this work, we have adopted three different extrapolation schemes, here briefly described. In the first one, denoted in the following as CCSD(T)/CBS1, correlation energies at the CBS limit were first obtained at the MP2 level, using the two-point extrapolation scheme of Halkier et al.[49] with the aug-cc-pVTZ and aug-cc-pVQZ basis sets. Next, such energies were summed to the interaction energies as obtained at the Hartree–Fock limit evaluated with the aug-cc-pVQZ basis set. The final CCSD(T)/CBS1 values were then obtained by adding to the previous terms a DCCSD(T) correction, estimated as the difference between the extrapolated correlation energies obtained at CCSD(T) and MP2 levels according to the Halkier scheme with the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The other two estimates of the CCSD(T) interaction energies at the CBS limit, denoted as CCSD(T)/CBS2 and CCSD(T)/CBS3, were obtained by extrapolating at the CBS limit, the interaction energies at the MP2 level, using either the two-point extrapolation of Halkier with the aug-cc-pVTZ and aug-ccpVQZ basis sets (CCSD(T)/CBS2 scheme) or the three-point extrapolation of Feller[50] with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets (CCSD(T)/CBS3 scheme). Then, we added the DCCSD(T) correction term evaluated, in both CBS2 and CBS3 schemes, as the difference of the interaction energies obtained at CCSD(T) and MP2 levels using the aug-ccpVTZ basis set. Single-point energy calculations along the dissociation coordinate within the rigid monomer approximation were also performed at MP2 and MP2.X levels of theory. The MP2.X approach has been recently proposed by Hobza and coworkers[51] to overcome the well-known problem of the tendency of MP2 to overestimate the binding energy in complexes involving aromatic moieties, in opposition to an underbinding prediction given by MP3.[52] Hobza proposed to correct the MP2 binding energies extrapolated to the CBS limit with a term given by Cx ½DEðMP3=xÞ 2DEðMP2=xÞ, where DE(MP3/x) and DE(MP2/x) are the binding energies evaluated, respectively, at the MP3 and MP2 levels with basis set x, and Cx is a coefficient depending on the chosen basis set[48] and tending to 0.5 for a CBS To obtain the MP2.X binding energies, 388
Journal of Computational Chemistry 2014, 35, 386–394
we have used the extrapolated MP2 values as obtained with the three-point Feller scheme,[50] and we have computed the correction term using the aug-cc-pVDZ basis set, corresponding to a Cx coefficient equal to 0.52. DFT calculations on the DCl/benzene and DBr/benzene dimers were performed by fully optimizing their geometry in the T-shaped approach with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets. Up to 34 functionals have been used to test their performance in describing the DX/p interaction. To evaluate the effect of the monomers relaxation on the binding energy values, single-point DFT calculations were also performed using the optimized minimum distances and the same fixed geometry of the monomers as adopted in the CCSD(T) scans. As said above, relaxation of the monomers geometry determined an energy gain of only 0.01 kcal/mol on average. In all cases, DFT calculations were performed by using a large pruned integration grid (99 radial shells and 590 angular points per shell) to avoid artifacts associated with numerical integration procedures. The need of using such ultrafine grid was in particular evidenced for meta-generalized gradient approximation (GGA) functionals.[53] The 34 XC functionals considered in the present investigation are reported in the following. They are classified according to the category to which they belong (see Supporting Information for further details), that is, (i) the functionals based on the GGA, BLYP,[54,55] PW91PW91 (or PW91),[56] HCTH407,[57] and B97D[58,59]; (ii) The hybrid GGA functionals (H-GGA), B3LYP,[55,60,61] B3P86,[60,62] O3LYP,[55,61,63] BH and HLYP,[54,55] PBE1PBE (or PBE0),[64,65] mPW1PW91,[56,66] mPW3LYP,[55,66–68] B97-2,[59,69] and B98[59,70]; (iii) The range-separated or long-range corrected GGA functionals (LC-GGA), LC-BLYP,[71] xB97X,[72] xB97XD,[58,73] and CAM-B3LYP[74]; (iv) The meta-GGA functionals (M-GGA), M06L,[75] VSXC,[76] TPSSTPSS (or TPSS),[77] and sHCTH[78]; (v) The hybrid meta-GGA functionals (HM-GGA), BMK,[79] TPSSh,[80] sHCTHhyb,[78] M05,[75] M052X,[81] M06,[82] M062X,[82] and M06HF[83,84]; (vi) The double hybrid GGA functionals (DH-GGA), B2PLYP,[54,55,85] mPW2PLYP,[66,86] B2PLYPD,[87] and mPW2PLYPD[87]; (vii) One among the generalized gradient exchange functionals, TPSSVWN5.[61,77] All CCSD(T), MP2, MP2.X, and DFT binding energies were corrected for basis set superposition error using the Boys–Bernardi[88] approach. All calculations were performed using the Gaussian 09 Rev. C.01[89] suite of programs.
Results and Discussion Weak to moderate halogen bonding interactions, such as those presently investigated, are dominated by electrostatic and dispersive contributions, unlike strong halogen bonds, where charge transfer can significantly contribute to the interaction.[90] The halogen atoms in the molecules here reported, DCl and DBr with D 5 H, HCC, F, and NC, are characterized by different r-holes as induced by substituents having different electron-withdrawing capabilities: the stronger is the electronwithdrawing moiety (in the order, H < HCC < F < NC), the larger and deeper is the r-hole. The interactions they can give rise to will be then characterized by a different balance of dispersion and electrostatic contributions. The features of the WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
r-hole can be pictorially evidenced by looking at the ESP plotted on an electron density isosurface, as already reported by Politzer et al.[14] In Figures 1a and 1b, the two cases of NCCl and NCBr, respectively, are exemplified. It can be clearly observed that the r-hole is more pronounced, as expected, on the more polarizable bromine atom with respect to chlorine atom. Recently, the origin of the halogen bond with “conventional” halogen acceptors (i.e., nitrogen and oxygen atoms) has also been rationalized using the Fukui function, f(r).[27,91] This function describes the variation of the electron [91,92] density q(r) upon changing the number of 1electrons N. @qðrÞ 1 The right derivative, namely f ðrÞ5 @N 5qN11 ðrÞ2qN ðrÞ, evidences the regions where a negative charge added to the system will tend to be distributed, that is, the regions more prone to act as electron acceptors. It has been demonstrated that it can be simply approximated by the square of the density of the LUMO orbital.[93] In Figures 1c and 1d, the Fukui function f 1 ðrÞ is plotted on the electron density isosurface for the NCCl and NCBr molecules, respectively, evidencing, similar to the ESP plot, the greater tendency of bromine with respect to chlorine atom to accept electrons along the CAX direction. The geometrical and energetical features of the interaction of DCl and DBr with the aromatic system of benzene have been first investigated by the CCSD(T) approach. We have determined the equilibrium distance R and the binding energy DE at the CBS limit by considering the dimers in the T-shaped configuration. Three different extrapolation procedures have
FULL PAPER
been used, denoted by CBS1, CBS2, and CBS3, as described in the previous section. The results reported in Tables 1 and 2 for the DCl/p and DBr/p interactions, respectively, clearly indicate a good convergence of both R and DE values as obtained by the three extrapolation procedures. The equilibrium distance of the halogen from the aromatic ring is in fact the same for the three approaches, indicating that a possible difference among them is smaller than the width of the grid used to ˚ ). Also, the binding energies locate the minimum (i.e., 0.05 A obtained with the three extrapolations are very close, with maximum differences of only 0.05 and 0.08 kcal/mol for DCl/p and DBr/p interactions, respectively, both obtained when the D moiety is the fluorine atom. For this reason, from now on we will refer to only one CCSD(T) set of extrapolated results to be used as reference for the other methods, namely CCSD(T)/ CBS1. By comparing the CCSD(T) results obtained for the different dimers, we observe a correlation between the DE values and the magnitude of the r-hole (see Supporting Information Fig. S1 for a plot of the ESP for the eight DX molecules), as already pointed out by other authors.[14,15,19] It is, however, to be evidenced a minor deviation from this trend for FX and NCX molecules, the latter showing a slightly lower maximum in the ESP than the former, albeit the corresponding interaction with benzene (at both CCSD(T) and MP2 levels) is slightly stronger for NCX. Such result could be ascribed to a greater polarization induced on the halogen by the cyano group with respect to
Figure 1. MP2/aug-cc-pVDZ ESP (first line) and Fukui function f1(r) (second line) mapped onto the isosurface of electron density (0.001 electrons/bohr-3) for NCCl (left) and NCBr (right) molecules.
Journal of Computational Chemistry 2014, 35, 386–394
389
FULL PAPER
WWW.C-CHEM.ORG
˚ ), and counterpoise-corrected interaction energies, DE Table 1. Minimum intermolecular distances between Cl and centroid of the benzene ring, R (A (kcal/mol), for dimers of DCl (D 5 H, HCC, F, NC) with benzene in the T-shaped approach.[a] HCl/p
CCSD(T)/aug-cc-pVDZ CCSD(T)/aug-cc-pVTZ CCSD(T)/CBS1 CCSD(T)/CBS2 CCSD(T)/CBS3 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ MP2/aug-cc-pVQZ MP2.X
HCCCl/p
FCl/p
NCCl/p
R
DE
R
DE
R
DE
R
DE
3.60 3.55 3.45 3.45 3.45 3.465 3.385 3.351 3.45
21.05 21.27 21.43 21.45 21.43 21.62 21.93 22.07 21.58
3.50 3.40 3.35 3.35 3.35 3.363 3.273 3.245 3.30
21.95 22.19 22.39 22.41 22.40 22.73 23.05 23.20 22.59
3.40 3.35 3.25 3.25 3.25 3.266 3.204 3.176 3.25
22.47 22.76 23.03 23.05 23.08 23.39 23.60 23.77 23.26
3.40 3.35 3.25 3.25 3.25 3.311 3.208 3.180 3.25
22.80 23.10 23.35 23.39 23.38 23.65 24.02 24.21 23.61
[a] CCSD(T) and MP2.X interaction energies evaluated on a 0.05-A˚ spaced grid.
the fluorine atom when interacting with benzene, determining a larger stabilizing contribution to DE for the NCX molecules. The electrostatic term represents in fact only the first-order approximation to the total interaction energy, but high-order terms may also play an important role. As for the CCSD(T) equilibrium distances R of the halogen from the centroid of the aromatic ring, we note that they tend to become smaller with increasing the magnitude of the r-hole, with the only exception of the FX/p system, which shows the same or even shorter R values than the stronger NCX/p system. Anyway, the equilibrium distances are always shorter than the sum of the van der Waal radius of the halogen atom[94] and the halfthickness of the aromatic ring,[95,96] corresponding to 3.5 and ˚ for the chlorine and the bromine atom, respectively. 3.6 A The MP2 results, also reported in Tables 1 and 2, show a general overestimate of the binding energies and an underestimate of the equilibrium distances. Both overestimate and underestimate tend to increase with the basis set dimensions, giving the unpleasant result that the smaller basis set, aug-ccpVDZ, happens to provide the better agreement with the extrapolated CCSD(T) results. This finding, already evidenced in our previous investigation on the DX/p interaction,[42] is a well-known problem of the MP2 methods in describing noncovalent interactions involving p systems.[45–47] An opposite tendency is usually shown by the MP3 method, and for this reason, Hobza and coworkers have recently proposed the
MP2.X method,[51] in which the results of the two methods are opportunely mixed through a weight, Cx, depending on the basis set quality, as reported in the Theoretical Methods section. The weight Cx had been determined to better reproduce the binding energies of the S66 set,[97] which does not contain any halogen bonded system. The MP2.X results, reported in Tables 1 and 2, show a significant improvement with respect to the MP2 method. The deviations from the CCSD(T)/CBS1 results are quite small, giving percentage RMS errors of 8.6 and 10.0% for the chlorine and the bromine halogen bonded complexes, respectively. Such errors are, however, slightly greater than the corresponding values reported by Hobza for dispersionbound, mixed dispersion-electrostatic, and electrostatic-bound complexes (about 7, 3, and 1%, respectively),[48] denoting a special character for halogen bonding and in particular for the DX/p interaction. DFT binding energies and equilibrium distances between the halogen atom and the centroid of the aromatic ring are reported as Supporting Information Tables S1–S4. Comparison of the DFT binding energies and equilibrium distances with the CCSD(T) reference results are graphically reported in Figures 2 and 3, respectively, where the ratios of the DFT/aug-ccpVQZ over the CCSD(T)/CBS1 values, averaged separately on the chlorine and bromine complexes, are plotted. Furthermore, the mean signed error, the mean absolute error (MAE), the mean signed percentage error, and the mean absolute
˚ ), and counterpoise-corrected interaction energies, DE Table 2. Minimum intermolecular distances between Br and centroid of the benzene ring, R (A (kcal/mol), for dimers of DBr (D 5 H, HCC, F, NC) with benzene in the T-shaped approach.[a] HBr/p
CCSD(T)/aug-cc-pVDZ CCSD(T)/aug-cc-pVTZ CCSD(T)/CBS1 CCSD(T)/CBS2 CCSD(T)/CBS3 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ MP2/aug-cc-pVQZ MP2.X
HCCBr/p
FBr/p
NCBr/p
R
DE
R
DE
R
DE
R
DE
3.65 3.60 3.50 3.50 3.50 3.524 3.438 3.402 3.50
21.56 21.83 22.03 22.06 22.04 22.30 22.70 22.87 22.27
3.55 3.50 3.40 3.40 3.40 3.424 3.358 3.314 3.40
22.45 22.71 22.96 22.99 23.00 23.39 23.73 23.92 23.26
3.45 3.40 3.30 3.30 3.30 3.334 3.260 3.260 3.30
23.16 23.51 23.89 23.94 23.97 24.23 24.57 24.82 24.24
3.45 3.40 3.35 3.35 3.35 3.369 3.280 3.251 3.35
23.36 23.69 24.00 24.04 24.05 24.37 24.78 25.01 24.35
[a] CCSD(T) and MP2.X interaction energies evaluated on a 0.05-A˚ spaced grid.
390
Journal of Computational Chemistry 2014, 35, 386–394
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 2. Ratios of DFT/aug-cc-pVQZ binding energies over the CCSD(T)/CBS1 values, averaged on the DCl/p, blue, and DBr/p, red (D 5 H, F, HCC, and NC) systems.
percentage error (MAPE) computed separately for the four chlorine and the four bromine complexes with the three considered basis sets, are reported as Supporting Information Tables S5–S8. The binding energies reported in Supporting Information Tables S1 and S3 for chlorine and bromine complexes, respectively, indicate that all functionals generally reproduce the trend provided by the CCSD(T) calculations, with increasing DE values on increasing the strength of the electron-withdrawing group and, for a same D substituent, on going from chlorine to bromine. There are, however, some remarkable exceptions, as denoted by those functionals (BLYP, B3LYP, B3P86, BH and HLYP, CAM-B3LYP, and BMK) which are not able to demonstrate the existence of the weaker HCl/p or even HBr/p halogen bonded complexes, giving repulsive curves at all intermolecular distances (such functionals are obviously excluded when calculating the average ratios plotted in Figures 2 and 3 and
the mean errors reported in Supporting Information Tables S5–S8). It also stands out the too large DE overestimates given by the VSXC functional. Finally, it is noteworthy the case of a discrete number of functionals, such as for example B97D, xB97X, and almost all the functionals of the Minnesota collection, which show a mismatch in the stability order on going from the FBr/p to the NCBr/p interaction, though they provide the correct order of stability for the corresponding chlorine complexes. The trend of the CCSD(T) equilibrium distances is instead correctly reproduced by all the functionals with only a few exceptions (see Supporting Information Tables S2 and S4 for chlorine and bromine complexes, respectively). About the dependence of the DFT results on the basis set, it is evident from Supporting Information Tables S1–S4 that binding energies usually decrease with increasing the dimensions of the basis set for almost all the functionals and that converged results are obtained with the aug-cc-pVTZ basis set.
Figure 3. Ratios of DFT/aug-cc-pVQZ equilibrium distances over the CCSD(T)/CBS1 values, averaged on the D-Cl/p, blue, and D-Br/p, red (D 5 H, F, HCC, and NC) systems. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Journal of Computational Chemistry 2014, 35, 386–394
391
FULL PAPER
WWW.C-CHEM.ORG
Equilibrium distances show a more variable behavior with increasing the dimensions of the basis set, as both a contraction and an elongation of the distance is observed depending not only on the functional but also on the system studied. Variations are more significant for the weaker HX/p interactions, but for most functionals it can be observed that “converged” results are obtained using the aug-cc-PVTZ basis set. The performance of the different examined functionals is summarized by the ratios plotted in Figures 2 and 3, where the closer is the ratio to 1, the better is the functional. From Figure 2, it is immediate to see that all the considered functionals, except VSXC, underestimate the binding energy. The large overestimate given by VSXC, which was parametrized without including noncovalent interactions,[76] was already evidenced in the investigation by Johnson et al. on purely dispersion-bound complexes,[53] in the study of the hydrogen bonded benzene-H2O system by Neves et al.[98] and in the systematic analysis by Zhao and Truhlar[82] including extensive sets of intermolecular interactions, in particular those dominated by dispersion contributions. It is also evident the better performance of almost all functionals in describing the DBr/p with respect to the weaker DCl/p interaction. The closest agreement with the CCSD(T)/CBS1 binding energy results, provided by functionals giving DE ratio > 0.80 for both chlorine and bromine complexes, is given by the functionals xB97X, M062X, B2PLYPD, and mPW2PLYP2D (MAPE values below or just above 10%, see Supporting Information Tables S5 and S7). If limiting to bromine complexes, B97D and xB97XD can be also included within the better performing functionals. It is interesting to note that all the considered functionals explicitly including dispersion, that is, B97D, xB97XD, B2PLYPD, and mPW2PLYP2D, fall within the subset of the better performing functionals, confirming the importance of the dispersion forces, besides electrostatic contributions, in determining the DX/p interaction. In particular, comparison of B2PLYPD and mPW2PLYP2D with the corresponding forms without dispersion, B2PLYP and mPW2PLYP2, respectively, indicates the substantial improving implied by the inclusion of the dispersion term (compare, e.g., the MAPEs for mPW2PLYP2D and mPW2PLYP2, below 10% and above 38%, respectively, for all basis sets). The good performance of the other two functionals, xB97X and M062X, can be ascribed to the inclusion of dispersion-bound complexes in the training sets used for the functional parametrization. It is worth noting the slightly better performance of xB97X with respect to xB97XD, a result which generalizes what previously reported in the study of the interaction of NCX and PhX with benzene,[42] where also iodine, besides chlorine and bromine, had been included in the analysis. It can be probably ascribed to the choice of the database used for parametrizing these functionals, where compounds exhibiting DX/p halogen bonding were not included. By looking at the performance of different functionals categories tested in this work, we cannot note a clear trend on going from the less accurate pure GGA functionals to the most accurate double hybrid GGA functionals. When climbing the Jacob’s ladder of DFT accuracy,[99] we find a large variability within each rung. The main origin of such variability is ascribable to 392
Journal of Computational Chemistry 2014, 35, 386–394
inclusion of the additive dispersion term into the functionals, as discussed above. Within such class of dispersion-corrected functionals, we note only a slight improvement for the double hybrid B2PLYPD and mPW2PLYP2D functionals with respect to the pure GGA (B97D) and the hybrid LC-GGA (xB97XD) functionals, according to the Jacob’s ladder. Moreover, almost uniform accuracy of the double hybrid functionals for chlorine and bromine complexes makes them more suitable when more halogen bonds are present in the investigated system. When excluding dispersion-corrected functionals, we note that pure GGA functionals, including those developed within the meta-approximation, cannot be used to treat the DX/p noncovalent interaction, though HCTH407 and M06 show a definite better performance with respect to the other functionals on this rung. Introduction of exact exchange, giving rise to H-, LC-, and HM-GGA functionals, does not involve, by itself, an improvement with respect to pure GGA functionals, however, it is evident a strong dependence on the family of functionals under consideration. For example, all the H-GGA functionals, including the popular B3LYP, are confirmed to be totally unable to describe long-range dispersive interactions (MAPEs above 60% for all basis sets). Conversely, xB97X and the Minnesota collection Myz provide in general a good description of the DX/p interaction, with M062X on the top of the hybrid functionals category. It is interesting to observe, in the Myz family, a systematic improvement of the performance when increasing the amount of exact exchange (M052X and M062X are better than M05 and M06, respectively). However, too much exact exchange appears to be detrimental, as can be deduced by the performance of M06HF, which is comparable with that of M06. Such trend matches that obtained for halogen bonding involving heteroatoms,[27] though in that case it had been demonstrated that the role of dispersion contribution is in general less important than what deduced for the DX/p halogen bonding. Finally, the double hybrids without dispersion, B2PLYP and mPW2PLYP, do not outperform the xB97X and the Myz collection, indicating that introduction of MP2 correlation, besides HF exchange, does not warrant, by itself, a better performance than that of lower rungs on Jacob’s ladder. Looking at Figure 3 and Supporting Information Tables S6 and S8, it appears that geometry is well described by a greater number of functionals, with respect to binding energy, for both chlorine and bromine complexes. In particular, besides the functionals explicitly including dispersion, all the Minnesota collection provides an excellent description of the minimum geometry. For these functionals, the MAE is less or only slightly greater than 0.05 A˚, which corresponds to the grid spacing used for the determination of the minimum distance at CCSD(T) level. The better agreement is provided once again by the B2PLYPD and mPW2PLYPD functionals, and the importance of including the dispersion contribution is further assessed.
Conclusions We have studied the DX/p halogen bonding between different halogenated molecules (D 5 H, HCC, F, and NC; X 5 Cl, Br) with WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
the aromatic system of benzene by CCSD(T), MP2, MP2.X, and DFT approaches, evaluating the minimum geometry and the binding energy with monomers in the T-shaped configuration. CCSD(T) calculations at the CBS limit have been performed to provide a reference for a plethora of DFT functionals and evaluate their performance in describing the present interaction. The results here obtained can be used for addressing the correct choice of the functional when studying more extended systems, for example, in the context of QM/MM calculations where more sophisticated post-HF methods are prohibitive. Results obtained at MP2 level confirm the previously reported tendency to overestimate the strength of the interaction when delocalized p systems are involved.[45–47] Such failure is adequately recovered by applying the recently proposed MP2.X approach,[51] though for the presently investigated DX/ p halogen bonding the improvement with respect to MP2 was not as good as that reported in the original paper.[51] The performance of up to 34 functionals, belonging to several categories of approximations, has been assessed by comparison with CCSD(T) extrapolated results. Such investigation allows to conclude that, owing to the mixed dispersive and electrostatic nature of the DX/p halogen bonding, a functional which correctly describes the dispersion contribution is mandatory to obtain good results. The double hybrid functionals B2PLYPD and mPW2PLYPD are here shown to provide the better results on both energies and equilibrium geometry. The corresponding forms without dispersion, B2PLYP and mPW2PLYP, fail to describe this interaction. Among the less expensive functionals, the better performance is provided by the xB97X and M062X functionals, while the xB97XD and B97D functionals are shown to work very well for bromine complexes but not so well for chlorine complexes. The worse performance is instead shown by the hybrid GGA functionals, including in particular the popular B3LYP functional whose use is then strongly discouraged for studying DX/p halogen bonded complexes. Keywords: halogen bond density functional theory benchmark calculations
How to cite this article: A. Forni, S. Pieraccini, S. Rendine, M. Sironi. J. Comput. Chem. 2014, 35, 386–394. DOI: 10.1002/ jcc.23507
]
Additional Supporting Information may be found in the online version of this article.
[1] P. Metrangolo, G. Resnati, Science 2008, 321, 918. [2] P. Metrangolo, G. Resnati, Halogen Bonding: Fundamentals and Applications; Springer: Berlin, 2008. [3] R. Kato, T. Imakubo, H. Yamamoto, R. Maeda, M. Fujiwara, H. Sawa, Mol. Cryst. Liq. Cryst. 2002, 380, 61. [4] E. Cariati, A. Forni, S. Biella, P. Metrangolo, F. Meyer, G. Resnati, S. Righetto, E. Tordin, R. Ugo, Chem. Commun. 2007, 2590. [5] J. Hulliger, P. J. Langley, Chem. Commun. 1998, 2557. [6] J. A. R. P. Sarma, F. H. Allen, V. J. Hoy, J. A. K. Howard, R. Thaimattam, K. Biradha, G. R. Desiraju, Chem. Commun. 1997, 101.
FULL PAPER
[7] P. K. Thallapally, G. R. Desiraju, M. Bagieu-Beucher, R. Masse, C. Bourgogne, J. F. Nicoud, Chem. Commun. 2002, 1052. [8] M. Adler, M. J. Kochanny, B. Ye, G. Rumennik, D. R. Light, S. Biancalana, M. Whitlow, Biochemistry 2002, 41, 15514. [9] E. Howard, R. Sanishvili, R. E. Cachau, B. Mitschler, P. Bart, V. Lamour, M. Van Zandt, E. Sibley, C. Bon, D. Moras, T. R. Schneider, A. Joachimiak, A. Podjarny, Proteins: Struct. Funct. Genet. 2004, 55, 792. [10] Y. Jiang, A. A. Alcaraz, J. M. Chen, H. Kobayashi, Y. J. Lu, J. P. Snyder, J. Med. Chem. 2006, 49, 1891. [11] P. Auffinger, F. A. Hays, E. Westhof, P. S. Ho, Proc. Natl. Acad. Sci. USA 2004, 101, 16789. [12] V. G. Tsirelson, P. F. Zou, T.-H. Tang, R. F. W. Bader, Acta Cryst. 1995, A51, 143. [13] P. Metrangolo, H. Neukirch, T. Pilati, G. Resnati, Acc. Chem. Res. 2005, 38, 386. [14] P. Politzer, J. S. Murray, T. Clark, Phys. Chem. Chem. Phys. 2010, 12, 7748. [15] P. Politzer, P. Lane, M. C. Concha, Y. Ma, S. Murray, J. Mol. Model. 2007, 13, 305. [16] D. Chopra, Cryst. Growth Des. 2012, 12, 541. [17] P. Metrangolo, J. S. Murray, T. Pilati, P. Politzer, G. Resnati, G. Terraneo, CrystEngComm. 2011, 13, 6593. [18] P. Metrangolo, J. S. Murray, T. Pilati, P. Politzer, G. Resnati, G. Terraneo, Cryst. Growth Des. 2011, 11, 4238. [19] K. E. Riley, J. S. Murray, P. Politzer, M. C. Concha, P. Hobza, J. Chem. Theory Comput. 2009, 5, 155. acˇ, R. J. Sola, M. C. Concha, [20] K. E. Riley, J. S. Murray, J. Fanfrlik, J. Rez F. M. Ramos, P. Politzer, J. Mol. Model. 2011, 17, 3309. [21] Y. X. Lu, J. W. Zou, Y. H. Wang, Q. S. Yu, Chem. Phys. 2007, 334, 1. [22] Y. X. Lu, J. W. Zou, Y. H. Wang, Y. J. Jiang, Q. S. Yu, J. Phys. Chem. A 2007, 111, 10781. acˇ, R. J. Sola, M. C. Concha, [23] K. E. Riley, J. S. Murray, J. Fanfrlik, J. Rez F. M. Ramos, P. Politzer, J. Mol. Model. 2013, 19, 4651. [24] Y. X. Lu, J. W. Zou, J. C. Fan, W. N. Zhao, Y. J. Jiang, Q. S. Yu, J. Comput. Chem. 2009, 30, 725. [25] A. Forni, S. Rendine, S. Pieraccini, M. Sironi, J. Mol. Graph. Model. 2012, 38, 31. [26] K. E. Riley, K. M. Merz, J. Phys. Chem. A 2007, 111, 1688. [27] S. Kozuch, J. M. L. Martin, J. Chem. Theory Comput. 2013, 9, 1918. [28] Y. Lu, H. Li, X. Zhu, W. Zhu, H. Liu, J. Phys. Chem. A 2011, 115, 4467. [29] Y. Lu, H. Li, X. Zhu, H. Liu, W. Zhu, Int. J. Quantum Chem. 2012, 112, 1421. [30] Y. Lu, T. Shi, Y. Wang, H. Yang, X. Yan, X. Luo, H. Jiang, W. Zhu, J. Med. Chem. 2009, 52, 2854. [31] S. Rendine, S. Pieraccini, A. Forni, M. Sironi, Phys. Chem. Chem. Phys. 2011, 13, 19508. [32] M. A. A. Ibrahim, J. Comput. Chem. 2011, 32, 2564. [33] M. Kolar, P. Hobza, J. Chem. Theory Comput. 2012, 8, 1325. [34] M. Carter, A. K. Rapp e, P. S. Ho, J. Chem. Theory Comput. 2012, 8, 2461. [35] W. L. Jorgensen, P. Schyman, J. Chem. Theory Comput. 2012, 8, 3895. [36] O. V. Shishkin, Chem. Phys. Lett. 2008, 458, 96. [37] D. J. R. Duarte, M. M. De Las Vallejos, N. M. Peruchena, J. Mol. Model. 2010, 16, 737. [38] H. G. Wallnoefer, T. Fox, K. R. Liedl, C. S. Tautermann, Phys. Chem. Chem. Phys. 2010, 12, 14941. [39] H. Matter, M. Nazar, S. Gussregen, D. W. Will, H. Schreuder, A. Bauer, Angew. Chem. Int. Ed. 2009, 48, 2911. [40] M. Nazare, M. Essrich, D. W. Will, H. Matter, K. Ritter, V. Wehner, PCT Int. Appl. EP 1314733, 2003. [41] C. Bissantz, B. Kuhn, M. Stahl, J. Med. Chem. 2010, 53, 5061. [42] A. Forni, S. Pieraccini, S. Rendine, F. Gabas, M. Sironi, ChemPhysChem 2012, 13, 4224. [43] Y. Lu, Y. Wang, W. Zhu, Phys. Chem. Chem. Phys. 2010, 12, 4543. [44] K. E. Riley, M. Pitonak, P. Jurecka, P. Hobza, Chem. Rev. 2010, 110, 5023. [45] C. D. Sherrill, Rev. Comput. Chem. 2009, 26, 1. [46] G. S. Tschumper, Rev. Comput. Chem. 2009, 26, 39. [47] M. O. Sinnokrot, E. F. Valeev, C. D. Sherrill, J. Am. Chem. Soc. 2002, 124, 10887. [48] K. E. Riley, J. Rezac, P. Hobza, Phys. Chem. Chem. Phys. 2011, 13, 21121. [49] A. Halkier, W. Klopper, T. Helgaker, P. Jorgensen, P. R. Taylor, J. Chem. Phys. 1999, 111, 9157. [50] D. Feller, J. Chem. Phys. 1993, 98, 7059.
Journal of Computational Chemistry 2014, 35, 386–394
393
FULL PAPER
WWW.C-CHEM.ORG
acˇ, M. Pito [51] R. Sedlak, K. E. Riley, J. Rez nak, P. Hobza, ChemPhysChem 2013, 14, 698. [52] L. Grafova, M. Pitonak, J. Rezac, P. Hobza, J. Chem. Theory Comput. 2010, 6, 2365. [53] E. R. Johnson, R. A. Wolkow, G. A. DiLabio, Chem. Phys. Lett. 2004, 394, 334. [54] A. D. Becke, Phys. Rev. A 1988, 38, 3098. [55] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. [56] J. P. Perdew, In Electronic Structure of Solids; P. Zieche, H. Eschrig, Eds.; Akademie: Berlin, Germany, 1991. [57] A. D. Boese, N. C. Handy, J. Chem. Phys. 2001, 114, 5497. [58] S. Grimme, J. Comput. Chem. 2006, 27, 1787. [59] A. D. Becke, J. Chem. Phys. 1997, 107, 8554. [60] A. D. Becke, J. Chem. Phys. 1993, 98, 5648. [61] S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. [62] J. P. Perdew, Phys. Rev. B 1986, 33, 8822. [63] A. J. Cohen, N. C. Handy, Mol. Phys. 2001, 99, 607. [64] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. [65] C. Adamo, V. Barone, J. Chem. Phys. 1999, 110, 6158. [66] C. Adamo, V. Barone, J. Chem. Phys. 1998, 108, 664. [67] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2004, 108, 6908. [68] X. Xu, W. A. Goddard, III, Proc. Natl. Acad. Sci. USA 2004, 101, 2673. [69] P. J. Wilson, T. J. Bradley, D. J. Tozer, J. Chem. Phys. 2001, 115, 9233. [70] H. L. Schmider, A. D. Becke, J. Chem. Phys. 1998, 108, 9624. [71] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540. [72] J. D. Chai, M. Head-Gordon, J. Chem. Phys. 2008, 128, 084106. [73] J. D. Chai, M. Head-Gordon, Phys. Chem. Chem. Phys. 2008, 10, 6615. [74] T. Yanai, D. Tew, N. C. Handy, Chem. Phys. Lett. 2004, 393, 51. [75] Y. Zhao, D. G. Truhlar, J. Chem. Phys. 2005, 123, 194101. [76] T. Van Voorhis, G. E. Scuseria, J. Chem. Phys. 1998, 109, 400. [77] J. M. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Phys. Rev. Lett. 2003, 91, 146401. [78] A. D. Boese, N. C. Handy, J. Chem. Phys. 2002, 116, 9559. [79] A. D. Boese, J. M. L. Martin, J. Chem. Phys. 2004, 121, 3405. [80] V. N. Staroverov, G. E. Scuseria, J. Tao, J. P. Perdew, J. Chem. Phys. 2003, 119, 12129. [81] Y. Zhao, D. G. Truhlar, J. Chem. Theory Comput. 2006, 2, 364. [82] Y. Zhao, D. G. Truhlar, Theor. Chem. Acc. 2008, 120, 215. [83] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 5121.
394
Journal of Computational Chemistry 2014, 35, 386–394
[84] [85] [86] [87] [88] [89]
[90] [91] [92] [93] [94] [95] [96] [97] [98] [99]
Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 13126. S. Grimme, J. Chem. Phys. 2006, 124, 034108. T. Schwabe, S. Grimme, Phys. Chem. Chem. Phys. 2006, 8, 4398. T. Schwabe, S. Grimme, Phys. Chem. Chem. Phys. 2007, 9, 3397. S. F. Boys, F. Bernardi, Mol. Phys. 1970, 19, 553. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision C.01; Gaussian, Inc., Wallingford CT, 2009. O. Hassel, Science 1970, 170, 497. B. Pinter, N. Nagels, W. A. Herrebout, F. D. Proft, Chem. Eur. J. 2013, 19, 519. R. G. Parr, W. Yang, J. Am. Chem. Soc. 1984, 106, 4049. W. Yang, R. G. Parr, R. Pucci, J. Chem. Phys. 1984, 81, 2862. A. Bondi, J. Phys. Chem 1964, 68, 441. A. W. Hanson, Acta Crystallogr. 1964, 17, 559. I. Saraogi, V. G. Vijay, S. Das, K. Sekar, T. N. G. Row, Cryst. Eng. 2003, 6, 69. acˇ, K. E. Riley, P. Hobza, J. Chem. Theory Comput. 2011, 7, 2427. J. Rez A. R. Neves, P. A. Fernandes, M. J. Ramos, J. Chem. Theory Comput. 2011, 7, 2059. J. P. Perdew, K. Schmidt, AIP Conf. Proc. 2001, 577, 1.
Received: 17 August 2013 Revised: 29 October 2013 Accepted: 21 November 2013 Published online on 13 December 2013
WWW.CHEMISTRYVIEWS.COM
