On the Influence of Density Functional Approximations on Some Local

ARTICLE pubs.acs.org/JPCA

On the Influence of Density Functional Approximations on Some Local Bader’s Atoms-in-Molecules Properties Vincent Tognetti* and Laurent Joubert IRCOF, CNRS UMR 6014, and FR 3038, Universit
e de Rouen et INSA de Rouen, 76821 Mont St. Aignan Cedex, France

bS Supporting Information ABSTRACT: In this article, we assess the ability of various density functionals to predict accurate values for some basic properties of the bond critical points of about 50 small molecules, including the recently proposed reduced gradient variation rates and involving typical ionic and covalent bonds, agostic interactions, and van der Waals complexes. The relation between the computed deviations and the geometric variations are discussed, as well as the topology variations. The possible correlation of these descriptors to atomization energies is considered, and the relevance of an accurate QTAIM analysis for correct descriptions of potential energy surfaces is addressed. Finally, we provide typical margins of error for the evaluation of these quantities and discuss their consequences for computational applications.

1. INTRODUCTION Computational chemistry has become an unavoidable tool for the rationalization of chemical processes and the design of new specific molecules. Notably, methods based on quantum mechanics have now reached an accuracy and a time scaling that enable them to correctly and routinely describe systems of thousands of atoms. As a consequence, a great deal of information can be extracted from these calculations, but not all of them provide clear insights into the studied phenomenon. It is thus the role of chemical models to build a small collection of meaningful descriptors from which one can easily (if possible) and unambiguously extract the most relevant information. From this viewpoint, the Quantum Theory of Atoms in Molecules (QTAIM) developed by Bader and co-workers1,2 constitutes an appealing framework, as it is deeply grounded in fundamental physics3,4 while proposing a unified and comprehensive approach of chemical bonding. Nevertheless, if the concepts and the formalism have now fully matured through useful controversies (see, for instance, refs 5
7), some practical points might still remain, from our point of view, underestimated. Among them, to the best of our knowledge, the influence of the computational protocol on the QTAIM descriptor values has not been the subject of a comprehensive study. Our goal here, however, is not to compensate for such a lack, but only to focus on one of the important computational parameters, namely, the choice of the quantum method and, more specifically, of the exchange-correlation functional in the framework of Kohn
Sham density functional theory (DFT).8 This choice is motivated by the fact that DFT is now used by a very large community of researchers, ranging from pure theoreticians to applied experimentalists. Obviously, the choice of the method is only one of the input parameters of a calculation. For instance, the influence of the basis set on QTAIM properties was studied by Volkov et al.,9 that r 2011 American Chemical Society

of electron correlation in post-Hartree
Fock methods by Gatti et al.,10 the impact of basis set superposition errors by Duran and co-workers,11 and the role of relativistic corrections by Reiher and co-workers.12 Concerning the role of DFT approximations itself, one has to mention the recent study by Rykounov et al. on substituted hydropyrimidines,13 those by Jablonski and Palusiak on hydrogen-bonded systems14 and on covalent bonds,15 and that on energetic molecules by Yau and co-workers.16 The quality of the densities produced by DFT was also assessed by Bochevarov and Friesner but from a global (integrated) perspective and not from a QTAIM approach.17 We intend here to extend the scope and the number of molecules considered (almost always less than a dozen in the above-mentioned studies), as well as the range of the considered DFT methods within a more systematic benchmark procedure. Such a work is somewhat complementary to those focusing on experimental QTAIM analysis procedures and on the comparison between experimental and calculated QTAIM properties.18,19

2. COMPUTATIONAL DETAILS 2.1. Methods and Molecular Sets. We considered the performances of 10 exchange-correlation functionals, classified following Perdew’s Jacob’s scale:20 local [local density approximation (LDA)] in the SVWN parametrization,21,22 gradientcorrected [generalized gradient approximation (GGA)] functionals (PBE,23 BLYP,24,25 HCTH26), meta-GGA (MGGA) functionals (TPSS,27 VSXC28) incorporating a kinetic energy density contribution (the so-called τ variable), global hybrids [mixing Hartree
Fock (HF) exchange, where the mixing coefficient is space-independent: PBE0,29 B3LYP,30,31 BHHLYP (also Received: April 5, 2011 Published: May 06, 2011 5505

dx.doi.org/10.1021/jp2031384 | J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A known as BHandHLYP), and hyperGGA functionals (τ þ HF), here B1B9532. To have some hints about wave function theories (that have frequently been used by Bader and co-workers), we also considered HF theory and the Møller
Plesset approach, up to second perturbative order (MP2). Reference values were obtained at the couple-cluster level including single and double substitutions (CCSD), even when we were aware that this level of calculation might not fully represent a golden standard for some particular cases (notably in cases of state degeneracies for which multideterminantal methods are required). All calculations were carried out with the Gaussian 03 package,33 using the aug-cc-pVTZ triple-ζ all-electron basis set. Notice that our purpose is not to discuss the influence of the basis set. Even if larger basis sets might be required to reach full convergence for some QTAIM descriptors,9 our goal is to give error margins corresponding to protocols that could also be used for realistic modeling of larger systems. The main test set consisted of 46 molecules (see the tables in the Supporting Information for the full list). Thirty-six were taken from the frequently used G2 database of Pople and coworkers34 and were either diatomic or of the symmetric AX2 formula. Six of the Se
Y or Se
Yþ type were also included: they belong to a set, called C2, that we recently designed35 and found useful for the assessment of DFT functionals.36 Finally, the radical cations (RCs) based on the first four noble gas dimers were considered, as they are challenging systems for most DFT methods [notably due to self-interaction errors (SIEs)37]. Notice that all of these molecules exhibit only one kind of bond at a time. This choice was driven by sampling considerations: every molecule had the same weight in the total error contribution, thus avoiding biases from exaggerated statistical importance (molecules with many bonds). To broaden our scope, but without more systematically investigating this family, we also considered a vanadium(III) complex38 (in its singlet spin state) belonging to the very popular bis(imino)pyridyl family (Brokkhart
Gibson catalysts)39,40 and bearing a propyl chain so that a γ-agostic interaction could be created (cf. Figure S1 in the Supporting Information). For this purpose, the all-electron Wachtersþf basis set41 was used for the metal atom, in conjunction with the standard 6-31þþG(d,p) basis set. Moreover, a β-agostic niobium complex that is problematic for DFT was investigated, using the computational protocol given in ref 42. Thus, the G2, C2, and RC sets contain (more or less) typical covalent, iono-covalent, and ionic bonds; agostic bonds are epitomized by the vanadium and niobium complexes; and, last, a dispersion-bonded system (the methane dimer) was studied. 2.2. Local QTAIM Quantities. We now review the studied QTAIM descriptors. Only local ones were retained, all of which are related to bond critical points (BCPs). Recall that a BCP is a point in real space where the density gradient vanishes that corresponds to a minimum along the bond path (BP) and to a maximum in the perpendicular plane. Importantly, the BCP and BP concepts apply only for minimum-energy geometries, even if, at every point in the potential energy surface (PES), critical points (but not BCPs) and atomic interaction lines (AILs, but not BPs) can be obviously localized [“the existence of an electron density accumulation between nuclei along the length of the AIL is not enough to call them bonded. One external condition (i.e., not observable in F) that needs to be imposed on the nuclei is

ARTICLE

that the forces on them vanish.”2]. We show two examples (in sections 3.2 and 3.3) of the fundamental role of this requirement. Accordingly, one cannot decouple geometry optimizations and QTAIM calculations; in other words, for two given functionals, the properties for a given BCP are evaluated at two (in general) different geometries. This point was recently discussed by Matta,43 who stated that “the agreement between H
F and both DFT and MP2 seems to be much better than generally appreciated” for geometries, concluding that “we can anticipate, thus, that H
F will perform very well in recovering the molecular and atomic properties calculated at MP2 and DFT”. It is known that GGAs tend to predict bonds that are too much long whereas HF contrarily overbinds, but, on average, the relative deviations for equilibrium distances are very low ( Au ¼ Ru ¼ 1 > > > > > 1 4 > < As ¼ , Rs ¼ 3 2ð3π2 Þ1=3 > >
 > 1=6 > 1 π 7 > > > At ¼ , Rt ¼ : 4 3 6

for ρð B r Þ ¼ uð B rÞ for ρð B r Þ ¼ sð B rÞ rÞ for ρð B r Þ ¼ tð B

ð1Þ one can evaluate the corresponding variation rates at a given BCP46,47 as Aρ  4πFð B r c ÞRρ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi λ1 2 sin2 φ cos2 θ þ λ2 2 sin2 φ sin2 θ þ λ3 2 cos2 φ sin φ dθ dφ Æρ0 æð B r cÞ ¼

Z

2π

Z

π

θ¼0 φ¼0

ð2Þ where λ1 e λ2 < λ3 denote the eigenvalues of the density Hessian matrix. In cases where the ellipticity (ε = λ1/λ2
1) equals 0, eq 2 can be expressed analytically (with e = λ3/λ1)48 as Aρ jλ1 j 8πFð B r c ÞRρ 8 pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi > <
e þ arcsinh e2
1 = e2
1 for e e
1  pffiffiffiffiffiffiffiffiffiffiffiffi > :
e þ ½π=2 þ arcsinðeÞ= 1
e2 for
1 < e < 0 Æρ0 æð B r cÞ ¼

ð3Þ For every functional and each BCP property [F(rBc), r F(rBc), Æu0 æ(rBc), Æs0 æ(rBc), Æt0 æ(rBc)], the relative mean signed deviations (MSDs with respect to CCSD) and the mean absolute deviations 2

5506

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

(MADs) were evaluated (where N is the number of retained molecules) MSDðPÞ ð%Þ ¼

Pi
Pi, CCSD 100 N 1 e i e N Pi, CCSD

MADðPÞ ð%Þ ¼

jPi
Pi, CCSD j 100 N 1 e i e N jPi, CCSD j

∑

∑

ð4Þ

Finally, the overall QTAIM performances were assessed as MSDðhQTAIMiÞ ¼ MADðhQTAIMiÞ ¼

1  5 1  5

P ∈ fFð

∑

MSDðPÞ

∑

MADðPÞ

! r c Þ, r2 Fð! r c Þ, hu0 ið! r c Þ, hs0 ið! r c Þ, ht 0 ið! r c Þg

P ∈ fFð

! r c Þ, r2 Fð! r c Þ, hu0 ið! r c Þ, hs0 ið! r c Þ, ht 0 ið! r c Þg

ð5Þ Atomic units are used exclusively for all of the QTAIM descriptors values. On the contrary, distances are expressed in angstroms. These Bader analyses were carried out using the AIMAll and MORPHY softwares.49
52 Obviously, these local BCP properties constitute only one facet of QTAIM, and other quantities (more refined and sometimes more suitable, mainly based on integrations over atomic basins) could be chosen as alternatives,53
55 for instance, delocalization indices56 or the source function.57,58 All calculated values are gathered in the Supporting Information (cf. Tables S1
S7).

3. RESULTS AND DISCUSSION 3.1. “Typical” Bonds (G2, C2, and RC Sets). 3.1.1. Topologies. In some cases, non-nuclear attractor critical points

(NNACPs) have been localized. As a consequence of the Poincar
e
Hopf relation, two or more BCPs are present. Such a situation is encountered for the well-known Li2 case (and, more generally, for Lin and Nan clusters),59 for which one NNACP and two BCPs are predicted by all methods, confirming (with more tested DFT approximations) that it must not be a computational artifact.60,61 NNACPs also appears for Na2 but only with wavefunction-based methods [HF, MP2, and CCSD (one NNACP)], with all density functionals (even those including a portion of HF exchange) finding only one BCP. A more complicated topology is observed for Si2: all DFT and MP2 methods predict only one BCP, whereas HF predicts one NNACP and two BCPs and CCSD even features two NNACPs and three BCPs. Finally, one NNACP appears with all methods for the P2 molecule, except SVWN, PBE, and BLYP. As shown by these four examples, no general trend can be found on the basis of functional classes. Similarly, two different post-Hartree
Fock methods such as MP2 and CCSD can lead to different topologies. In conclusion, it can be asserted that the absence or presence of an NNACP generally cannot be easily rationalized. There is nevertheless one important rule verified by all of these sytems (although the reciprocal one is false): there is no case where a DFT approximation predicts an NNACP that CCSD does not.62 A particular case also deserves to be mentioned. For the CH, OH, SH, ClO, NO, NSe, PSe, and AsSe radicals, all methods give a BCP with a nonzero ellipticity (from 0.001 to 0.1), which is contradictory to the expected cylindrical symmetry. This remains when the tightness of the self-consistent-field procedure is increased. The stability of the single-determinant wave functions has also been checked, but they were confirmed to be true energy

Figure 1. Mean absolute deviations (%) for the evaluation of the BCP electron density (top) and Laplacian (bottom), depending on the chosen computational method.

minima. This unphysical result can be cured by the use of fractional orbital occupancies to avoid this symmetry breaking. However, this does not imply that molecular properties will be badly described. For instance, we showed in ref 35 that CCSD(T), some functionals, and complete-active-space self-consistentfield (CASSCF) approaches give almost identical geometries and atomization energies for the C2 set (including NSe, PSe, and AsSe), probably due to fortuitous error cancellations. 3.1.2. Densities and Laplacian Values at BCPs. Does a global good geometry description imply accurate BCP densities? (See Figure 1.) The linear regression coefficient r2 for MAD[F(rBc)] = f [MAD(rBeq)] seems to only partly (i.e., qualitatively) confirm this conclusion, r2 = 0.90, with: MAD½Fð B r c Þ ¼ 0:893 þ 3:071MADð B r eq Þ

ð6Þ

Some trends can actually be identified: GGAs and MGGAs tend to give BCP density values that are too small (negative MSDs) regardless of the subset. On the contrary, the inclusion of HF exchange (hybrids and HGGAs) leads to a global balanced cancellation, as corroborated also by the fact that HF gives the highest MSD values for the C2 and G2 sets. Actually, PBE0 always gives higher BCP densities than PBE, regardless of 5507

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A the molecule. Similarly, for each molecule, one has, without exception, F(rBc)(BLYP) < F(rBc)(B3LYP) < F(rBc)(BHHLYP), following an increasing function of the HF mixing. As a consequence, the MSD variations also increase with the HF contribution: MSD(B3LYP)
MSD(BLYP) = 3.6% (20% HF), MSD(PBE0)
MSD(PBE) = 3.9% (25% HF), MSD(BHHLYP)
MSD(BLYP) = 5.0% (20% HF).63 Nevertheless, there is no correlation between the rungs of Perdew’s DFT scale and the average performances: a GGA can perform better than a hybrid one [MAD(HCTH) = 4.4%, MAD(BHHLYP) = 5.1%] or worse [MAD(PBE0) = 3.5%]. Similarly, the hierarchies are also different depending on the subset (HCTH is best for C2 and G2, whereas MP2 and BHHLYP are best for RC). Concerning post-HF methods, one always has FHF(rBc) g MP2 F (rBc) (except for H2, ClO, AsSe, He2þ, and Ar2þ) and FHF(rBc) g FCCSD(rBc) (except for AsSe and Ar2þ), so that correlation (in the post-HF sense) tends to lower BCP density values. On the other hand, there are about as many cases for which the CCSD values are lower than the MP2 ones as there are where the CCSD values are higher. We now examine the density Laplacian values (bottom part of Figure 1): the range in values is considerably larger, as it goes from 14.9 (B1B95) to 66.9 (SVWN). Moreover, performances for the BCP Laplacian evaluations are not related to those for BCP densities. Indeed, the linear regression coefficient between MAD[F(rBc)] and MAD[r2F(rBc)] is small, at r2 = 0.41, confirming the independence of the global errors in the density and Laplacian values. As a corollary, the r2 value between MAD[F(rBeq)] and MAD[r2F(rBc)] is also small (0.53). From another general viewpoint, every method gives a balanced description of the whole data set in the sense that they feature about as many underestimations as overestimations; the MSD values are accordingly quite low (in absolute value), but not the MAD values (typically greater than 20). For instance, for TPSS, MSD =
2.1%, whereas MAD = 34.3%. The only clear trend is that all methods strongly underestimate the Laplacian values for RCs, contrary to wave function methods (with the notable exception of BHHLYP). On average, hybrid functionals tend to provide smaller MAD values for all molecules, as evidenced by the interval separations: [33.4, 64.5] for GGAs, [15.6, 36.4] for hybrids. This is not only because the hybrids give a better RC Laplacian description: indeed, the intervals are also separated for the G2 set ([32.3, 58.7] for GGAs, [17.2, 32.8] for hybrids). Contrary to their position in the Perdew scale, the MGGAs are not “between” the two, as VSXC provides better results than hybrids and TPSS is very close to GGA for the G2 subset, but not for the C2 set. The role of the HF exchange is also difficult to explain. Indeed, BHHLYP (50% HF exchange) gives the best hybrid result for the C2 set (12.4%), whereas HF has a very high MAD (111.8%). Thus, the performances are not clearly related to the HF exchange ratio. Even for the PBE/PBE0 couple, the situation is far from being as simple as in the density case. Indeed, at first glance, PBE0 tends to give higher absolute values for Laplacian at BCPs. However, exceptions are not rare and cannot be neglected: F2, PH2, NSe, PSeþ. Some cases are even more problematic: the sign can be different, for instance, for Cl2: þ0.0002 with PBE and
0.0426 for PBE0 (even though, clearly, the values are close to 0). This sign problem is also encountered, more consequently, for CS (PBE,
0.025; PBE0, 0.207). In addition, if other methods than

ARTICLE

PBE/PBE0 are considered, the same problems arise for PH2 (only HF predicts a positive value), Cl2 again (positive values for BLYP and TPSS), CS (negative values for PBE, BLYP, and SVWN), and PSeþ (positive values for HF and VSXC). Once more, it is not directly linked to a possible “bad” geometry, because, for instance, the difference between the equilibrium bond lengths obtained with TPSS and B3LYP equals 0.0012 Å for Cl2, and that between VSXC and B3LYP is 0.0001 Å for PSeþ; that is, the same geometry leads to opposite signs for the Laplacian values, even if small differences are not necessarily significant. Finally, deviations greater than 100% are far from being the exception, regardless of the method, which confirms the high sensitivity of Laplacian calculations. 3.1.3. Reduced Density Gradient Variation Rates. We now consider the average (reduced) density gradient variation rates. As for the Laplacian, they involve the Hessian eigenvalues, but through their squares, so that compensation between negative and positive eigenvalues cannot similarly occur. One can wonder, however, whether they are complementary or carry the same information. To this aim, the left part of Figure 2 presents the computed Æu0 æ(rBc) values with respect to the Laplacian values at the CCSD level: the two quantities, albeit built using the same eigenvalues, can be considered uncorrelated. Compared to Laplacians, the deviations for Æu0 æ(rBc) are moderate, comparable to those for the densities. However, there is no parallel evolution between these two last quantities: r2 = 0.64 for MAD[F(rBc)] with respect to MAD[Æu0 æ(rBc)]. Interestingly, with the sole exception of BHHLYP, all DFT methods tend to underestimate Æu0 æ(rBc) at the BCP (negative MSD values), and hybrid functionals clearly outperform GGAs for every subset as demonstrated by the respective almost separate intervals: [3.6, 5.7] and [5.7, 8.1] for G2, [2.0, 8.5] and [6.9, 13.3] for C2, and [8.4, 22.1] and [30.4, 36.7] for RC (this last set being sensitive to self-interaction errors, accounting for the poor behavior of GGAs). As Æu0 æ(rBc) and the BCP densities have similar mean deviations, one might expect that the deviations for Æs0 æ(rBc) and Æt0 æ(rBc) would be of the same magnitude. The situation is actually more subtle, as we notice much more positive values for the MSD, so that no definitive trend can be drawn. This fact thus constitutes indirect confirmation that the deviations for the densities and for the reduced density gradients are not correlated. Other evidence is that the regression coefficients between MAD[Æs0 æ(rBc)] and MAD[Æu0 æ(rBc)] and between MAD[Æt0 æ(rBc)] and MAD[Æu0 æ(rBc)] are 0.46 and 0.38, respectively. However, one has a nearly perfect correlation between MAD[Æs0 æ(rBc)] and MAD[Æt0 æ(rBc)], with r2 = 0.984 (right part of Figure 2): MAD½ht 0 ið B r c Þ ¼
0:3154 þ 0:9523MAD½hs0 ið B r c Þ ð7Þ Finally, the role of the HF mixing can be briefly discussed. Adding a portion of exact exchange tends to increase Æu0 æ(rBc). This is actually the case when comparing PBE and PBE0 for all molecules, excepted for H2S, Na2, SH (although the differences are very small, at 0.002, 0.0003, and 0.001, respectively), ClO, and PSe. Similarly, the Æu0 æ(rBc) values given by HF are almost always higher than those provided by PBE0 (with the only exceptions being H2S, ClF, SH, ClO, and NSeþ). The same trends are shared for BLYP and its global hybrids. As HF mixing at the same time increases densities and Æu0 æ(rBc), one can wonder whether these two effects will compensate each other when considering the quantity Æs0 æ(rBc) or Æt0 æ(rBc). In almost 50% of the 5508

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

Figure 2. (Left) Plot of Æu0 æ(rBc) with respect to the Laplacian values at the BCPs at the CCSD level. (Right) Evolution of MAD[Æt0 æ(rBc)] with respect to MAD[Æs0 æ(rBc)] considering all methods (in red, linear regression corresponding to eq 6). Atomic units are used.

Figure 3. Mean absolute deviations (%) for the evaluation of all considered QTAIM local quantities, depending on the chosen computational method.

cases, PBE0 gives higher Æs0 æ(rBc) and Æt0 æ(rBc) values than PBE, so that the two effects reciprocally balance. 3.1.4. Overall Performances. We are now able to assess the global QTAIM performances of all of these methods. The total MAD[ÆQTAIMæ] values, which show a “profile” similar to that of the Laplacian (because r2F has the highest deviations), are collected in Figure 3. It appears that HF is comparable to the less accurate density functionals whereas MP2 is close to the best ones. The differences between HF and MP2 are not negligible (a factor 2). However, MP2 cannot be considered as the most reliable method for QTAIM analyses. Moreover, because of its high computational costs, it can be advantageously replaced by DFT, taking into consideration that hybrid functionals, MGGAs, and HGGAs are more accurate than GGAs, with the best performances being obtained with the B1B95 functional. Note that some of the tested functional were parametrized on the atomization energies of the extended G2 set, for instance, VSXC and HCTH (best GGA for our tests). Are these good

Figure 4. CCSD atomization energies (in kcal/mol) with respect to BCP densities (in atomic units).

performances biased that by the fact we also used a part of this set? 3.1.5. Atomization Energies and QTAIM Descriptors. In other words, are QTAIM BCP properties linked to molecular atomization energies? To answer this question, linear regressions between these energies and each individual QTAIM descriptor were made.64 The obtained regression coefficients r2 (at the CCSD level) are equal to 0.27, 0.22, 0.24, 0.01, and 0.00 for F(rBc) (see Figure 4), |r2F(rBc)|, Æu0 æ(rBc), Æs0 æ(rBc), and Æt0 æ(rBc), respectively. Thus, none of these QTAIM descriptors is able to predict these atomization energies alone. Multilinear regressions were then tested with the five QTAIM descriptors, but the correlation coefficient remained very low (r2 = 0.32). This indirectly confirms the findings of Bochevarov and Friesner about the possible absence of a relation between good densities and good thermochemistry.17 Indeed, good performances for thermochemistry often rely on (uncontrolled) error cancellations. Recall, however, that the BCP local properties (especially the density and the Laplacian values) have been found to be 5509

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Deviations for the BCP density and Laplacian values with respect to the equilibrium bond length deviation obtained using B3LYP. Atomic units are used for QTAIM quantities; distances are in angstroms.

correlated to the strength of some bonds and, notably, to the important case of H bonding.65,66 Bader and Ess
en’s observations are also fundamental: in the homonuclear series from B2 to F2, the F(rBc) and |r2F(rBc)| variations closely parallel the variations of the corresponding binding energies.67 One thus has effective correlations between these primary local BCP properties and energetic data, when staying inside the same family of molecules. 3.1.6. Other Considerations. As shown before, there does not exist a universal relation (i.e., valid regardless of the method) between the quality of the optimized geometry and the performance for QTAIM evaluation (except maybe for densities). However, one might wonder whether, for a given and fixed functional, correlations can be established. As a test example, B3LYP was considered, because it is one of the most used functionals in computational studies. Toward this end, for every molecule, the geometrical deviation (δreq) and the QTAIM deviation for property P (δP) were plotted, as exemplified by Figure 5 for BCP density and Laplacian values. The plots can be divided in four quadrants: In the one labeled a, δreq is negative, and δP is positive. δreq and δP are both positive in quadrant b, and both are negative in quadrant c. Finally, δreq is positive, and δP is negative in quadrant d. We consider that a consistent trend exists if (almost) all points (i.e., molecules) belong to opposite quadrants: a and d (“negative” correlation) or b and c (“positive” correlation). In the first case, an overestimation of one of two quantities leads to the underestimation of the other, whereas in the second case, quantities are simultaneously overestimated or underestimated. It clearly appears that, even if quantitative relations are illusory, one qualitative trend can still be proposed: geometrical and densities deviations exhibit a negative correlation. The correlation is very roughly positive between Laplacian values and equilibrium bond lengths, but many points vary from this trend. Notably, some points fall along the dashed vertical or horizontal lines: for them, one of the two quantities is very accurately described, whereas the other is far from the reference values. For Æu0 æ(rBc), no convincing real qualitative trend can be formulated, whereas Æs0 æ(rBc) and Æt0 æ(rBc) follow a positive

correlation. In any case, no quantitative relationship can be obtained. Another important computational parameter is linked to the various thresholds involved in calculations, one of which is related to the quality of the integration grid for the energy numerical quadrature. All of the previous results were obtained using a so-called “fine” grid (75 radial shells and 302 angular points per shell). In this case, we also used an “ultrafine” (99,590) grid. Another threshold concerns the cutoffs on forces and step sizes to determine the geometrical convergence. By default, the “tight” option was used; the “very tight” one was considered for this analysis as well. Finally, the self-consistent-field procedure has its own threshold (10
8 on the density by default); a tighter one (10
11) was also investigated. Then, to assess the importance of these parameters, we extracted three representative molecules, namely, H2, CO, and NaCl, because they represent the pure covalent, iono-covalent, and ionic bonding types, respectively, for two functionals: B3LYP and TPSS. It has indeed already been shown that MGGAs can be prone to spurious oscillations or random noise due to the second-order derivative reduced variables they use.68,69 The results are gathered in the Supporting Information (Table S8). One can see that the magnitudes of the maximum deviations are 10
4 Å for the bond lengths, 10
4 au for the densities, and 10
3 au for the Laplacian. In all cases, they represent only less than 0.1% of the contribution. Numerical errors are thus definitely negligible compared to the errors associated with the quantum methods. 3.2. Agostic Bonds. Agostic bonds are essential in organometallic chemistry but have posed a challenge for theoreticians for decades.70
74 Thus they should not be excluded from this bencharmk study. Like previously, molecular graphs can depend on the used functionals, as well illustrated by the mentioned vanadium(III) complex. TPSS indeed predicts a bond path between the γ-hydrogen and the metal (Figure 6). On the other hand, PBE0 predicts a bond path between the metal and the γ-carbon. In this latter case, the bond path is also consequently curved and distorted (Figure 6). This discrepancy can be geometrically rationalized by examining the VCRCβCγ dihedral 5510

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

Figure 6. Molecular graphs for the vanadium complex. The BCPs are denoted by squares, and the BPs of interest are represented by a continuous bold line. The labels of the nuclei that lie in the projection plane are in bold, whereas those that are outside are outlined. For the left graph, corresponding to the PBE0 optimization, every point within 1.0 bohr of the projection plane was considered to belong to this plane; for the left panel, resulting from a TPSS calculation, this threshold is equal to 0.3 bohr. A view of the alkyl chain conformation is provided below each of the graphs.

angle (lower parts of Figure 6): it is smaller for PBE0 (27) than for TPSS (30), so that Cγ is better oriented to directly interact with the metal in the PBE0 geometry, whereas Hγ is better positioned to do so in the TPSS-optimized one. Nevertheless, these two different schemes can be reunified in the same framework, notably by invocating the metal
methyl agostic concept that we recently developed.75 Another example of the impact of the exchange-correlation functional on the agostic geometry was recently given.42,46,76 Pantazis et al. studied a critical agostic niobium complex and showed that some density functionals are unable to locate the agostic minimum (that has been experimentally characterized).42 For instance, BP86 succeeds, whereas BLYP fails. We carried out topological analyses of this complex. First, we used the BLYP density at the BP86 equilibrium geometry (close to the geometry obtained by X-ray diffraction): we obtained a (3,
1) critical point between the niobium and the β-hydrogen atoms. However, this critical point disappeared when the BLYP density was used at the BLYP-optimized geometry. Thus, it would be meaningless to study the properties of an agostic critical point from a singlepoint approach when no agostic bond is predicted. Also note that we showed in ref 46 that no simple rationalization can be invoked to predict the presence or absence of an agostic bond because the behaviors of the density functional at the two asymptotic boundaries of the reduced density gradient intervals have to be considered. 3.3. QTAIM and van der Waals Interactions. Finally, we want to deepen the possible links between QTAIM and PES, by tackling the following problem: Does a good description of QTAIM properties along a path in a PES ensure a good

Figure 7. Unrelaxed potential energy surface for the methane dimer depending on the used method. Interaction energies are in kcal/mol, and distances are in angstroms.

description of the studied PES? To answer this question, a (dispersion) van der Waals (ubiquitous in chemistry and biology) complex was studied. It is well-known that many exchange-correlation functionals are not able to locate a minimum.77 Also recall that, according to Slater78 and as regularly stressed by Bader,79 there is no fundamental distinction 5511

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

Figure 9. Evolution (in atomic units) of the (reduced) density gradient variation rates depending on the distance (in Å) between the two methane moieties forming the methane dimer.

Figure 8. Deviations for the C
C critical point QTAIM properties (top: density, bottom: Laplacian) with respect to MP2 reference values (in %), depending on the distance (in Å) between the two methane moieties forming the methane dimer.

between van der Waals bonding and covalent bonding: they are governed by the same physics. Now, as a particular but important case, consider a dimer composed of two protomolecules (for instance, two benzene molecules, two ethylene molecules): each one obeys the Poincar
e
Hopf (PH) equality, ni(3,
3)
ni(3,
1) þ ni(3,1)
ni(3,3) = 1 (for i = 1, 2), so that K = ∑2i=1[ni(3,
3)
ni(3,
1) þ ni(3,1)
ni(3,3)] = 2. Thus, to fulfill the PH equality, at least one more (3,
1) critical point or one more (3, 3) critical point must be created (such that K is decreased by one unit) when the interaction is switched on, regardless of the interaction distance. As a special case, we investigated the methane dimer (CH4)2, for which only one (3,
1) critical point is created by the interaction of the two monomers at every separation distance.80 Moreover, for symmetry reasons (D3d space group), this critical point is located in the middle of the C
C line. Note that this critical point can be called a BCP only when an energy minimum exists. This is the case for HCTH and PBE, but not for B3LYP, BLYP, and BPW91 (Figure 7).81,82 Now consider the densities and Laplacian values at this point with respect to the C
C distance and to the computational method. The corresponding results (Figure 8)83 are disparate: Whereas the signed relative deviation, ΔF(rBc) (%) = 100[F(rBc)
FMP2(rBc)]/FMP2(rBc), for the critical point density increases when the distance between the two moieties increases for

HCTH, PBE, BLYP, and B3LYP, it descreases for BPW91. Similarly, no trends about the sign of the deviations can be formulated: PBE, HCTH, and BLYP overestimate it, whereas there is a sign change for both HCTH and BPW91. Interestingly, B3LYP provides the best densities (at all C
C distances) but fails at the same time to find an energy minimum. In addition, the existence of a minimum is not due to the fact that the density is either too weak or too high at this van der Waals critical point. The situation is slightly different in the case of the Laplacian (also Figure 8). All DFT approximations provide values that are too low, regardless of the C
C separation. However, the method that predicts the lowest values (HCTH) is one of the few that correctly find a minimum. Once again, one cannot ascribe the absence of a minimum to a too-important charge depletion or concentration. Once again, B3LYP gives quite accurate values. Finally, Figure 9 represents the evolution of the reduced density gradients at the van der Waals critical point. For MP2, Æs0 æ(rBc) and Æt0 æ(rBc) exhibit an almost increasing linear behavior with respect to the dimer distance. In contrast, as expected, Æu0 æ(rBc) decreases, but not linearly. Once again, HCTH, which provides a qualitative good PES shape (unlike B3LYP), shows quite important deviations, whereas B3LYP is very close to MP2; for instance, Æt0 æHCTH(rBc) is not even a monotonic function of the distance, whereas Æs0 æB3LYP(rBc) and Æt0 æB3LYP(rBc) are almost linear.84 As for the previously mentioned niobium agostic complex, if one had used only a single-point approach based on the reference equilibrium geometry, one would have concluded that B3LYP (and others) is very accurate, even though it poorly describes the underlying physics. As a final remark, it is worth noticing that van der Waals corrections of the Grimme type,85 which enable accurate energy profiles to be obtained, do not modify the electron densities, whereas the PES will be completely different. However, as demonstrated by Langreth and co-workers,86 even a correct DFT description of the van der Waals interactions, based on a fully nonlocal correlation,87,88 will lead to very small density changes, being several orders of magnitude smaller than the error margins of the applied methods. In our opinion, this makes an accurate QTAIM description of large dispersion complexes 5512

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A (i.e., those that cannot be studied by reliable post-Hartree
Fock methods) not currently feasible.

4. CONCLUSIONS In this article, we have assessed the ability of DFT approximations to predict accurate local QTAIM properties for a large set of small molecules, as well as various and representative bonding types. The following conclusions can be drawn from this study: - Topology can strongly depend on the method, but the differences cannot be rationalized in terms of DFT approximation features. In addition, QTAIM analysis on singlepoint calculations can lead to erroneous conclusions when potential energy surfaces strongly differ. - GGAs and MGGAs tend to provide densities that are too small at BCPs, whereas the inclusion of a fraction of HF exchange leads to a global overestimation of these densities, although it is not directly proportional to the HF ratio. - Conversely, all methods can overestimate or underestimate Laplacian values, but hybrid functionals tend to give smaller absolute deviations. The sign of the Laplacian at a BCP can even depend on the functional used, but a wrong sign is not necessarily linked to a wrong geometry. Moreover, there exists no relation for a given method between the performance for the evaluation of BCP densities and that for the corresponding Laplacians. - The reduced density gradient variation rates are much less sensitive to the computational protocol. Furthermore, they provide additional information on the Laplacian. DFT methods tend to underestimate Æu0 æ(rBc), but hybrid functionals substantially reduce this overestimation because HF mixing increases Æu0 æ(rBc) values (but not Æs0 æ(rBc) and Æt0 æ(rBc) values). As these variation rates seem stable, their use in QTAIM analyses certainly deserves to be fostered. - Regardless of the local QTAIM property studied, the deviations cannot be quantitatively correlated to the deviations in equilibrium bond lengths, so that DFT approximations can provide, on average, excellent geometries and relatively bad QTAIM values, and vice versa. The only reported qualitative trend is a “negative” correlation between BCP densities and bond length deviations. - Similarly, accurate local QTAIM descriptions do not imply accurate energetic or thermochemical properties. More specifically, it seems that atomization energies cannot generally be simply predicted by looking at BCP properties. - Overall, no direct relation between Perdew’s Jacob’s ladder rungs and global QTAIM performances can be invoked. However, it has been shown that HF performances are close to those of the less accurate DFT approximations, whereas MP2 performances are similar to those of the most accurate functional, with B1B95 being the best for our tests. - Computational thresholds on SCF and geometrical convergences have a definitely negligible impact when compared to the impact of the method choice. - The evolution of local QTAIM properties with respect to the separation distance in dispersion complexes cannot provide information on the PES shape. The failure of some DFT approximations to locate energy minima is not due to poor-quality electron densities around the critical points. The only identified trends are that MP2 tends to give higher Laplacian values than DFT and that reference Æs0 æ(rBc) and Æt0 æ(rBc) evolutions seem linear.

ARTICLE

We now give typical margins of error for the considered local QTAIM properties (i.e., typical MAD values): - 5% for densities, - 25% for Laplacians, and - 7% for the reduced density gradient variation rates. Some of the most popular DFT approximations are determined by fitting on energetic (atomization energies, activation barriers, ionization potentials, etc.) or geometric properties.89 One can allegedly wonder whether QTAIM properties can or must be integrated in the used fitting databases. In principle, a functional must provide high-quality density properties, but an approximation can be “chemically” or computationally efficient without accurate local BCP descriptions. For these reasons, the fitting of parameters on only QTAIM properties cannot ensure that good physicochemical properties will be obtained. On the other hand, adding some BCP requirements to existing thermochemical databases is probably a nice way to improve the versatility and range of applications of a functional. We are still working on this point. Finally, as said in the Introduction, integrated properties are at least as important as local BCP properties. Recent studies have investigated their sensitivity to the chosen method.90,91 There is no doubt that such issues still deserve more analysis.

’ ASSOCIATED CONTENT

bS

Supporting Information. Tables with optimized bond lengths and calculated density, Laplacian, and (reduced) density gradient values at bond critical points, depending on the used functionals. Tables with statistical analysis and threshold dependences. View of the vanadium(III) complex. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT The Centre de Ressources Informatiques de Haute-Normandie (CRIHAN) and the Institut du D
eveloppement et des Ressources en Informatique Scientique (IDRIS) are thanked for computer time. One of the reviewers is kindly thanked for drawing our attention to ref 43 and for the noteworthy discussion about it. ’ REFERENCES (1) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, U.K., 1990. (2) Popelier, P. L. A Atoms in Molecules: An Introduction; Pearson Education: Harlow, U.K., 2000. (3) Bader, R. F. W. Phys. Rev. B 1994, 49, 13348. (4) Martín Pend
as, A.; Blanco, M. A.; Francisco, E. Chem. Phys. Lett. 2006, 417, 16. (5) Cassam-Chenai, P.; Jayatilaka, D. Theor. Chem. Acc. 2001, 105, 276. (6) Bader, R. F. W. Theor. Chem. Acc. 2002, 107, 381. (7) Delle Site, L. J. Math. Chem. 2004, 35, 289. (8) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: NewYork, 1989. (9) Volkov, A.; Koritsanszky, T.; Chodkiewicz, M.; King, H. F. J. Comput. Chem. 2008, 30, 1379. (10) Gatti, C.; MacDougall, P. J.; Bader, R. F. W. J. Chem. Phys. 1988, 88, 3792. 5513

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A (11) Salvador, P.; Fradera, X.; Duran, M. J. Chem. Phys. 2000, 112, 10106. (12) Eickerling, G.; Mastalerz, R.; Herz, V.; Scherer, W.; Himmel, H.-J.; Reiher, M. J. Chem. Theory Comput. 2007, 3, 2182. (13) Rykounov, A. A.; Tsirelson, V. G. J. Mol. Struct. (THEOCHEM) 2009, 906, 11. (14) Jablonski, M.; Palusiak, M. J. Phys. Chem. A 2010, 114, 2240. (15) Jablonski, M.; Palusiak, M. J. Phys. Chem. A 2010, 114, 12498. (16) Yau, A. D.; Byrd, E. F. C.; Rice, B. M. J. Phys. Chem. A 2009, 113, 6166. (17) Bochevarov, A. D.; Friesner, R. A. J. Chem. Phys. 2008, 128, 034102. (18) Volkov, A.; Abramov, Y.; Coppens, P.; Gatti, C. Acta Crystallogr. A 2000, 56, 332. (19) Destro, R.; Leconte, L.; Presti, L.; Roversi, P.; Soave, R. Acta Crystallogr. A 2004, 60, 365. (20) Perdew, J. P.; Ruzsinszky, A.; Tao, J.; Staroverov, V. N.; Scuseria, G. E.; Csonka, G. J. Chem. Phys. 2005, 123, 062201. (21) Slater, J. C. Quantum Theory of Molecular and Solids. Vol. 4: The Self-Consistent Field for Molecular and Solids; McGraw-Hill: New York, 1974. (22) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (23) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (24) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (25) Lee, C.; Yang, W.; R. G. Parr, R. G. Phys. Rev. B 1988, 37, 785. (26) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 1998, 109, 6264. (27) Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (28) Van Voorhis, T.; Scuseria, G. E. J. Chem. Phys. 1998, 109, 400. (29) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (30) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (31) Stephens, P. J.; Devlin, J. F.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (32) Becke, A. D. J. Chem. Phys. 1996, 104, 1040. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02, Gaussian, Inc.: Wallingford, CT, 2004. (34) Johnson, B.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1993, 98, 5612. (35) Rekhis, M.; Ouamerali, O.; Joubert, L.; Tognetti, V.; Adamo, C. J. Mol. Struct. (THEOCHEM) 2008, 863, 79. (36) Tognetti, V.; Cortona, P.; Adamo, C. Chem. Phys. Lett. 2008, 460, 536. (37) Chermette, H.; Ciofini, I.; Mariotti, F.; Daul, C. J. Chem. Phys. 2001, 114, 1447. (38) Reardon, D.; Conan, F; Gambarotta, S.; Yap, Y.; Wang, Q. J. Am. Chem. Soc. 1999, 121, 9318. (39) Small, B. L.; Brookhart, M.; Bennett, A. M. A. J. Am. Chem. Soc. 1998, 120, 4049. (40) Britovsek, G. J. P.; Bruce, M.; Gibson, V. C.; Kimberley, B. S.; Maddox, P. J.; McTavish, S. J.; Solan, G. A.; White, A. J. P.; Williams, D. J. Chem. Commun. 1998, 849.

ARTICLE

(41) Wachters, A. J. H. J. Chem. Phys. 1970, 52, 1033. Bauschlicher, C. W., Jr.; Langhoff, S. R.; Barnes, L. A. J. Chem. Phys. 1989, 91, 2399. (42) Pantazis, D. A.; McGrady, J. E.; Maseras, F.; Etienne, M. J. Chem. Theory Comput. 2007, 3, 1329. (43) Matta, C. F. J. Comput. Chem. 2010, 31, 1297. (44) Koch, U.; Popelier, P. L. A J. Phys. Chem. 1995, 99, 9747. (45) Popelier, P. L. A.; Logothetis, G. J. Organomet. Chem. 1998, 555, 101. (46) Tognetti, V.; Joubert, L.; Cortona, P.; Adamo, C. J. Phys. Chem. A 2009, 113, 12322. (47) Tognetti, V.; Joubert, L.; Adamo, C. J. Chem. Phys. 2010, 132, 211101. (48) Note that this function is continuous for e =
1 and can be written with only one of the two formulas when the arcsinh or arcsin functions are correctly extended on the imaginary axis. (49) Keith, T. A. AIMAll, version 10.12.08; TK Gristmill Software, 2010; available at aim.tkgristmill.com. (50) MORPHY98, a topological analysis program written by Popelier, P. L. A. with a contribution from Bone, R. G. A.; UMIST: Manchester, U.K.. (51) Popelier, P. L. A. Comput. Phys. Commun. 1996, 93, 212. (52) Note that analyses of Gaussian 03 post-HF calculations require the density=current keyword, as well as the pop=noab keyword, when they are performed within an unrestricted formalism. (53) Gatti, C. Z. Kristallogr. 2005, 220, 399. (54) Mebs, S.; Grabowsky, S.; F€orster, D.; Kickbusch, R.; Hartl, M.; Daemen, L. L.; Morgenroth, W.; Luger, P.; Paulus, B.; Lentz, D. J. Phys. Chem. A 2010, 114, 10185. (55) Mebs, S.; Kalinowski, R.; Grabowsky, S.; F€orster, D.; Kickbusch, R.; Justus, E.; Morgenroth, W.; Paulmann, C.; Luger, P.; Gabel, D.; Lentz, D. Inorg. Chem. 2011, 50, 90. (56) Fradera, X.; Austen, M. A.; Bader, R. F. W. J. Phys. Chem. A 1999, 103, 304. (57) Bader, R. F. W.; Gatti, C. Chem. Phys. Lett. 1998, 287, 233. (58) Lo Presti, L.; Gatti, C. Chem. Phys. Lett. 2009, 476, 308. (59) Cao, W. L.; Gatti, C.; MacDougall, P. J.; Bader, R. F. W. Chem. Phys. Lett. 1987, 141, 380. (60) Edgecombe, K. E.; Esquivel, R. O.; Smith, V. H., Jr.; M€ullerPlathe, F. J. Chem. Phys. 2002, 97, 2593. (61) Cioslowski, J. J. Phys. Chem. 1990, 94, 5496. (62) This fact has driven our choice for the number of molecules used in the calculations of averages: only the molecules that do not display NNACP with CCSD have been retained (N = 42 in eq 4) because all methods give the same topology for all them. (63) This is somewhat fallacious, because the correlation part is not the same (no VWN contribution in BLYP and BHHLYP, contrary to B3LYP) and because the scaling factor for the B88 gradient correction differs in BLYP and BHHLYP (1 and 0.5, respectively). Furthermore, if this reasoning were perfectly true, one would expect HF to predict higher BCP densities than BHHLYP, which is true in some cases (less than about 80%, e.g., H2, CH) but not in others (more than 20%, including ionic species such as LiH, BeH, and LiF, but also covalent species such as PH2). (64) As Laplacian values can be either positive or negative whereas atomization energies all have the same sign, we used the absolute value of the computed Laplacians. (65) Parthasarathi, R.; Subramanian, V.; Sathyamurthy, N. J. Phys. Chem. A 2005, 109, 843. (66) Parthasarathi, R.; Subramanian, V.; Sathyamurthy, N. J. Phys. Chem. A 2006, 110, 3349. (67) Bader, R. F. W.; Ess
en, H. J. Chem. Phys. 1984, 80, 1943. (68) Johnson, E. R.; Becke, A. D.; Sherrill, C. D.; DiLabio, G. A. J. Chem. Phys. 2009, 131, 034111. (69) Wheeler, S. E.; Houk, K. N. J. Chem. Theory Comput. 2010, 6, 395. (70) Brookhart, M.; Green, M. L. H. J. Organomet. Chem. 1983, 250, 395. 5514

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515

The Journal of Physical Chemistry A

ARTICLE

(71) Vogt, D. In Applied Homogeneous Catalysis with Organometallic Compounds; Cornils, B., Hermann, W. A., Eds.; Wiley: Weinheim, Germany, 1996; Vol. 1, pp 245
258. (72) Scherer, W.; McGrady, G. S. Angew. Chem, Int. Ed. 2004, 43, 1782. (73) Clot, E.; Eisenstein, O. Struct. Bonding (Berlin) 2004, 113, 1. (74) Lein, M. Coord. Chem. Rev. 2009, 253, 625. (75) Tognetti, V.; Joubert, L.; Raucoules, R; De Bruin, T.; Adamo, C., manuscript in preparation. (76) Pantazis, D. A.; McGrady, J. E.; Besora, M.; Maseras, F.; Etienne, M. Organometallics 2008, 27, 1128. (77) Riley, K. E.; Piton
ak, M.; Jurecka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023. (78) Slater, J. C. J. Chem. Phys. 1972, 57, 2389. (79) Bader, R. F. W. J. Phys. Chem. A 2010, 114, 7431. (80) This is not at all a general result. For example, in the case of the ethylene dimer in D2d symmetry, the topology is much richer: four ring critical points, one cage critical point, and four BCPs between the hydrogens of each C2H4 moiety. (81) Tsuzuki, S.; Uchimaru, T.; Tanabe, K. Chem. Phys. Lett. 1998, 287, 202. (82) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985. (83) Unrelaxed scans on the C
C distance, using the reference geometry given in ref 66 and the aug-cc-pVTZ basis set, MP2 for reference. (84) The (reduced) density gradient rates here are mainly driven by the positive Hessian eigenvalues (λ3) because the λ3/λ1,2 ratio is equal to 66 for C
C = 3.30 Å and to 37 for C
C = 4.30 Å (MP2). If one neglects λ1,2, eq 2 reads simply Æρ0 æ(rBc) = Aρλ3/[2F(rBc)Rρ]. (85) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (86) Thonhauser, T.; Cooper, V. R.; Li, S.; Pudzer, A.; Hyldgaard, P.; Langreth, D. C. Phys. Rev. B 2007, 76, 125112. (87) Dion, M.; Rydberg, H.; Schr€oder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (88) Vydrov, O. A; Van Voorhis, T. J. Chem. Phys. 2010, 132, 164113. (89) Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2010, 6, 107. (90) Henn, J.; Ilge, D.; Leusser, D.; Stalke, D.; Engels, B. J. Phys. Chem. A 2004, 108, 9442. (91) Phillips, J.; Hudspeth, M.; Browne, P.; Peralta, J. E. Chem. Phys. Lett. 2010, 495, 146.

5515

dx.doi.org/10.1021/jp2031384 |J. Phys. Chem. A 2011, 115, 5505–5515