Influence of Atoms-in-Molecules Methods on Shared-Electron

Apr 26, 2010 - The effect of using different methods to obtain atoms in molecules (AIMs) on the shared-electron distribution indices (SEDIs) and domai...
0 downloads 12 Views 3MB Size
Download PDF
8754

J. Phys. Chem. A 2010, 114, 8754–8763

Influence of Atoms-in-Molecules Methods on Shared-Electron Distribution Indices and Domain-Averaged Fermi Holes† Patrick Bultinck,‡ David L. Cooper,§ and Robert Ponec*,| Department of Inorganic and Physical Chemistry, Ghent UniVersity, Krijgslaan 281 (S3), 9000 Gent, Belgium, Department of Chemistry, UniVersity of LiVerpool, LiVerpool, L69 7ZD, United Kingdom, and Institute of Chemical Process Fundamentals, Academy of Sciences of the Czech Republic V.V.i., Prague 6, Suchdol 2, 165 02, Czech Republic ReceiVed: February 25, 2010; ReVised Manuscript ReceiVed: April 9, 2010

The effect of using different methods to obtain atoms in molecules (AIMs) on the shared-electron distribution indices (SEDIs) and domain-averaged Fermi holes (DAFHs) is examined using a test set of diatomic molecules. Use of Bader’s binary AIM model gives significantly different SEDIs as a function of internuclear distances than do self-consistent Hirshfeld-based AIM models. DAFH eigenvectors remain very similar for all AIM among the different methods. The corresponding eigenvalues are found to differ significantly, although the sums of complementary eigenvalues show only relatively small changes. The choice of a binary division of molecular space into AIM domains is found to lead to extra structure in SEDI curves, but this probably has limited chemical significance. Introduction More than 1500 years after the ancient Greek philosopher Demokritos suggested the existence of atoms as the smallest building blocks of matter, the same idea was revived and fruitfully explored by John Dalton, whose atomic hypothesis reshaped mysteriously magical alchemy into the quantitative science of chemistry.1 Because of its immense impact on modern chemistry, the idea of a molecule consisting of interacting atoms is still used up to the present day, and much, if not most, of our understanding of molecular structure and reactivity is based on the view that molecules consist of atoms held together by chemical bonds. Because of these successes, an entire jargon of chemical concepts has grown that apparently also stood the test of time despite often being short of a proper quantum mechanical derivation from first principles. One of these widely used concepts is the atom in the molecule (AIM). The first attempts to put this concept on a firmer theoretical basis are due to Coulson2 and Moffitt3,4 and their idea of the valence state of the AIM. Since then, many different approaches have been introduced to describe the AIM. In this article, we focus on methods that rely on the molecular electron density function and its 3D partitioning into atomic parts. Without doubt, one of the most popular methods is Bader’s quantum theory of the atom in the molecule (QTAIM).5-7 Despite its wide acceptance and use, the QTAIM approach to the partitioning of molecule into atomic fragments is not the only one possible, and in the past few years several other alternative partitioning schemes have been proposed. Consequently, there are now multiple ways to partition the molecule into atomic fragments, each of which can be more useful for various specific purposes. The existence of several partitioning schemes has recently been addressed in the study by Parr, Ayers, and Nalewajski, who concluded the AIM to be best considered †

Part of the “Klaus Ruedenberg Festschrift”. * Corresponding author. E-mail: [email protected]. ‡ Ghent University. § University of Liverpool. | Academy of Sciences of the Czech Republic v.v.i.

a so-called noumenon, an ideal conceptual construct ultimately unknowable by observation or by unique definition but conceivable by mind and reason.8 Our aim in this study is to perform a systematic comparison of the impact of various existing definitions of the AIM on the description of chemical bonding as provided by the numerical values of the shared-electron distribution indices (SEDIs) and by visual analysis of domain-averaged Fermi holes (DAFHs). Given that each of these approaches relies on the integration over the domains of individual atoms of certain quantities derived from the (correlated) pair density, it is thus of interest to see whether, and to what extent, the results of such an analysis depend on the particular partitioning scheme used for the definition of AIM. The motivation for the present study is the recent observation by Ponec and Cooper that depending on the AIM method used (in their case Bader’s method and the Mulliken9 approach) potentially interesting chemical information was or was not found in the SEDI versus internuclear distance curves.10 At that time, the rather structure-poor Mulliken plots were discarded because of the well-known shortcomings of the Mulliken method.11 In the present article, we wish to examine the differences in SEDI plots between Bader’s method and the recent self-consistent Hirshfeld related methods. Moreover, we wish to examine the effect of the AIM method on DAFH eigenvectors and eigenvalues to see whether the chemical conclusions drawn from inspection of DAFH pictures remain the same. To accomplish this goal, the above approaches have been applied to monitoring the process of electron reorganization accompanying the splitting of bonding electron pairs in a small set of simple diatomic molecules consisting of the isoelectronic series N2, CO, and NO+, as well as the highly polar LiF molecule. Theoretical Methods The AIM in QTAIM is derived from basic quantum mechanics using the fact that the topological condition of zero-flux serves as the boundary condition for the application of Schwinger’s principle of stationary action in the definition of an open

10.1021/jp101707w  2010 American Chemical Society Published on Web 04/26/2010

Influence of AIM on SEDI and DAFH

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8755

system,5 although questions have been raised about its uniqueness.12 The resulting AIM is one that attaches to each AIM an atomic domain in which all electron density is uniquely attributed to that AIM and there can be no electron density outside this domain that is ascribed to that atom. In general, an AIM density FA(r) for atom A may be written as the product of the molecular density function F(r) with a positive weight function 0 e wA(r) e 1 as in

FA(r) ) wA(r)F(r) ) wA(r)

∑ ηiFi(r)

(1)

i

in which ηi is the occupation number of natural orbital i and Fi(r) is its density. QTAIM considers that only two values are possible for wA(r), namely, zero or one. Other techniques used in the present article allow wA(r) to vary freely between zero and one and allow every AIM to extend to infinity in every direction, although obviously the larger values of wA(r) are obtained in regions dominated by the AIM A. Several such socalled fuzzy AIM have been described in the literature, and we focus here on those recently derived from the Hirshfeld method.13 A promolecular density was constructed in these Hirshfeldbased approaches by assembling isolated atomic densities at the same geometry as that in the molecule, and a weight function is introduced as follows

wA(r) )

FA0 (NA0 , r)

∑ FB0 (NB0 , r)

(2)

B

It is then assumed that this weight function can also be used to distribute the molecular density through eq 1. In eq 2, F0A(N0A,r) is the isolated atomic density for atom A with NA0 electrons. The Hirshfeld method has been criticized on several occasions because, among other reasons, different AIM density functions can be obtained depending on the isolated atom states used in eq 2.14-18 More precisely, as shown explicitly in eq 2, the weight function depends on the isolated atomic populations used for computing the weight function. Traditionally, neutral atoms are used {N0A ) ZA}, leading to debatable atomic charges in the case of highly polar species such as LiF.15-17 Moreover, the resulting charges sometimes depend dramatically on the values used for the set {N0A}. The so-called Hirshfeld-I method, introduced recently,17,18 solves this problem by adapting these populations until self-consistency is reached. The essence is that starting from an initial arbitrary choice of {N0A}, new atomic populations for the AIM are computed via

NA )

∫ FA(r) dr ) ∫ wA(r)F(r) dr

(3)

In the next step, isolated atomic densities are computed for exactly this number of electrons and are used to obtain new weight functions according to eq 2, meaning that except at the first iteration the {N0A} equals the set {NA} obtained from the previous iteration. This entire procedure is repeated until selfconsistency is reached. The outcome of this Hirshfeld-I procedure has been shown to be unique and to exhibit marginal basis set dependence19 while leading to atomic charges that reproduce molecular electrostatic potentials20 very well.

More recently, an alternative Hirshfeld-related technique was introduced under the name iterated stockholder atoms (ISAs).21 This scheme reaches self-consistency without recourse to isolated state atomic densities, and it also leads to a unique solution.22,23 Instead of using true atomic densities, this method composes density functions attached to every atom in the following way. Starting from, for example, density functions of neutral atoms, the Hirshfeld recipe is carried out and the AIM density is computed at each point in space. Then, for a chosen set of radii around the nuclei, the average AIM density on that shell is computed. The weight function on a specific radius is then given by

w (r) )

〈FA〉|r-RA|

A

∑ 〈FB〉|r-R | B

(4)

B

The notation 〈FA〉|r-RA| denotes the average AIM density of atom A on a shell of radius |r - RA|. As in Hirshfeld-I, this process is repeated, except that there is no recourse to isolated atomic densities. The requirement in the case of ISA is that the average density for each AIM at each radius no longer changes.21 The main difference between the ISA method and Hirshfeld-I is thus that in the latter we insist on using an isolated atomic density, whereas in the ISA method, every shell around the atom is treated as a separate variable. As a consequence, the densities for the atoms used to obtain the ISA AIM need not be monotonically decreasing.22 The original all-electron ISA method can easily be extended to more orbital specific approaches. The partitioning of the molecular density can be done not only atom-by-atom, as in QTAIM or the Hirshfeld and Hirshfeld-I methods, but can also be done in the following way

F(r) )

∑ ∑ ηiwiA(r)Fi(r) A

(5)

i

where an orbital specific weight function has been introduced (cf. eq 4)

wiA(r)

)

〈FiA〉|r-RA|

∑ 〈FiB〉|r-R | B

(6)

B

in which the average is taken for each molecular orbital separately. This is not possible in the Hirshfeld-I method because it would require the availability of isolated atomic densities for quite small numbers of electrons, which leads to all sorts of conceptual and computational problems. The generalized ISA method described above can obviously also be extended to any sort of partitioning per orbital or group of orbitals and will henceforth be denoted ISA-G. Because the main goal of the present article is to examine the eventual influence of the choice of a particular AIM partitioning scheme on the description of the splitting of the chemical bonds, as monitored by the numerical values of SEDI indices and by DAFH analysis, the main ideas of both of these methodologies are described here to the extent necessary for the present study. Both the DAFH and SEDI analyses are based on integration of the exchange-correlation density matrix given here as

8756

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Bultinck et al.

Fexc(r1, r′1 ;r2, r′2) ) F(r1, r′1)F(r2, r′2) - 2F(2)(r1, r′1 ;r2, r′2) (7) This exchange-correlation density measures the departure of two electrons from independent electron behavior to truly correlated behavior via the one electron density matrices F(r1,r1′) and F(r2,r2′) and the two-electron density matrix F(2)(r1,r1′;r2,r2′). In the case that Fexc(r1,r1′;r2,r2′) is integrated over the domain of a single AIM A, the resulting quantity gA(r1,r1′)

∫r )r′ wA(r2)Fexc(r1, r′1;r2, r′2) dr2 ) F(r1, r′1) ∫r )r′ wA(r2)F(r2, r′2) dr2 2 ∫r )r′ wA(r2)F(2)(r1, r′1 ;r2, r′2) dr2 ) NAF(r1, r′1) - 2 ∑ Γijklψ*i (r1)ψk(r′1)SjlA

gA(r1, r′1) )

2

2

2

2

2

2

ijkl

)

∑ ψ*i (r1)(N ηiδik - 2 ∑ ΓijklSjlA)ψk(r′1) A

ik

)

jl

∑ ψ*i (r1)Gikψk(r′1)

(8)

ik

represents the DAFH associated with the AIM domain of that atom.24-26 In previous studies, some of us have shown that the most useful feature of the DAFH method is that it straightforwardly provides useful information not only about the core and lone electron pairs associated with the atomic domain of A but also about the free or broken valences created by the formal splitting of the bonds required to isolate the corresponding domain from the rest of the molecule.10,24-29 The analysis consists, in the first step, of the diagonalization of the matrix G ){Gik} representing the hole 8 in an appropriate basis. The resulting set of eigenvectors (as well as the associated eigenvalues) is then, in a second step, subjected to the isopycnic localization procedure,27-31 whose goal it is to transform the usually delocalized eigenvectors into more localized functions that are reminiscent of chemical bonds, lone pairs, and so on in terms of which chemists traditionally think about molecules and their structures. The structural information is extracted by means of visual inspection of the resulting localized eigenvectors as well as by examining the numerical values of the associated (transformed) eigenvalues. The lone and core electron pairs present in the AIM domain of a given atom have eigenvalues close to two. Broken valences for nonpolar bonds characteristic of ideal sharing of the bonding electron pair are characterized by eigenvalues close to unity. The formal splitting of polar bonds results in free valences whose associated eigenvalues can often deviate substantially from unity. Such deviations provide information about the polarity of the formally broken bonds. Further integration of 8 over the domain of a second atom, B, yields the quantity

kAB )

∫r )r′ wA(r1)gA(r1, r′1) dr1 1

(9)

1

The SEDI ∆AB between atoms A and B is then obtained as

∆AB ) kAB + kBA

(10)

The so-called localization indices correspond to diagonal elements ∆AA ) kAA. The SEDI is the generalization of the Wiberg-Giambiagi-Mayer index32-34 to the case of 3D splitting of the molecular density function in AIM, as first suggested by Bader and Stephens.35 The following normalization condition holds, irrespective of the AIM method employed

∫r )r′ gA(r1, r′1) dr1 ) NA 1

1

(11)

However, except for calculations at the single-determinant level, the matrix G need not be positive definite. Typically, in addition to the dominant positive eigenvalues, diagonalization of G can yield a small number of low magnitude negative ones. Because the isopycnic transformation can be applied only for the positive eigenvalues, it follows that there can be (usually small) deviations from the normalization condition after application of the localization procedure, so that the hole is not entirely invariant in DAFH analysis. This is an issue to which we intend to return in future work. Computational Methods To address the issues mentioned above, full-valence CASSCF calculations were performed with the Cartesian cc-pvdz basis set for a set of small diatomic molecules, namely, LiF and the isoelectronic molecules CO, N, and NO+. All calculations were carried out in MOLPRO.36 For all molecules, the different AIM techniques described above were used to obtain atom condensed overlap matrices (AOMs). For the Bader analysis, the proaim program,37 was used and for all other methods, programs were developed in house. For the Bader, Hirshfeld-I, and ISA techniques AOM elements SAij are obtained as

SijA )

A (r)ψj(r) dr ∫ ψ*(r)w i

(12)

In the ISA-G technique, they are obtained as

SijA )

(

∫ ψ*(r) i

)

wiA(r) + wjA(r) ψj(r) dr 2

(13)

The SEDI and DAFH analyses are performed using the valence space only first- and second-order density matrices from the full-valence CASSCF calculations. For the QTAIM and Hirshfeld-I procedures, however, the total one-electron density is partitioned into AIMs. Results and Discussion To gain insight into the effect of the choice of the AIM method on the resulting SEDI values and on the DAFH analysis, we consider first the molecule N2, for which SEDI and DAFH analysis has been previously studied using QTAIM.29 Because the QTAIM-based results for this molecule were completely in line with chemical expectations, it serves as a good test to see whether the same is true for other AIM methods. Next, CO is studied as an example of a polar bond. The case of LiF is described as an example of a highly polar system. Lastly, the molecule NO+ is studied as an example of a molecule for which the internal charge distribution near dissociation is completely different from that near equilibrium.

Influence of AIM on SEDI and DAFH

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8757

Figure 1. N2 dissociation: QTAIM DAFH eigenvalues and eigenvectors over the domain of atom N1 as a function of internuclear distance, R.

QTAIM-based DAFH analysis for N2 yields four dominant eigenvectors whose values vary only slightly with increasing internuclear distance, R. One of the eigenvalues is close to two, and inspection of the associated eigenvector in Figure 1 shows that it corresponds to a lone pair on N. The remaining three eigenvalues remain close to one. Combined with visual inspection of the shape of the corresponding eigenvectors, these can be associated with the free valences of one σ and two π components of the formally broken NN triple bond. Whereas the eigenvalues of these free valences vary only marginally with R, increasing internuclear distance has a profound effect on the shape of the associated eigenvectors. In keeping with intuitive expectations, the observed changes can straightforwardly be correlated with the picture of monotonic splitting of the bonding electron pairs. Therefore, whereas for shorter R the eigenvectors resemble the occupied σ and π orbitals of the NN triple bond, elongation of the NN bond length results in increasing localization of the eigenvectors onto only one of the atoms, and at the limit of complete dissociation, these localized functions finally take the form expected for dissociation into two N(4S) atoms. Because the main goal of this study is to investigate the impact of various alternative definitions of AIM, a similar systematic DAFH analysis of the dissociation process was performed using the Hirshfeld-I, ISA, and ISA-G partitioning schemes. The resulting picture for the Hirshfeld-I method is

summarized in Figure 2. Comparison of Figures 1 and 2 shows that the choice of the particular partitioning scheme has practically no effect on the shape of the eigenvectors and that the description of the dissociation process emerging from the analysis is essentially the same. In addition to this striking visual resemblance, close coincidence is observed also in the range of internuclear distances where the splitting apparently takes place. Consistent with intuitive expectation, both methods agree in predicting that the splitting of the bonding electron pair corresponding to the π bond precedes the dissociation of the σ bond. The same is true for the remaining two partitioning schemes, ISA and ISA-G. (See Figures S2 and S3 in the Supporting Information.) In a sense, this result is not too surprising because in the case of ideally nonpolar bonds, like those in N2, it is clear that any physically reasonable definition of AIM must necessarily correspond to a 1:1 partitioning of molecular electron density into atomic fragments. This inherent simplicity is also reflected in the description of the dissociation process in terms of SEDI versus internuclear distance plots. As Figure 3 shows, the shapes of all of the plots are very similar, and the only very noticeable difference is the systematic shift of the QTAIM SEDI plot toward lower values. The existence of such a shift could, however, be anticipated because of the reports of previous studies,10,38 according to which the inclusion of electron correlation reduces the SEDI for the N2 molecule at

8758

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Bultinck et al.

Figure 2. N2 dissociation: Hirshfeld-I DAFH eigenvalues and eigenvectors over the domain of atom N1 as a function of internuclear distance, R.

Figure 3. Comparison of different AIM methods for the SEDI versus internuclear distance (R, in angstroms) curves for N2.

equilibrium geometry to a value close to two, which strongly deviates from the Hartree-Fock value close to three that is more consistent with the existence of the NN triple bond. As suggested by the plot, the above extreme reduction of SEDI is partially absent in the alternative partitioning schemes that consistently predict for N2 at equilibrium geometry a value around 2.4. All of the AIM techniques converge to the same limiting value of zero at large internuclear distances.

Having established that all methods perform well for a relatively simple case, a more challenging comparison can be made using the isoelectronic CO molecule as a test case. It is worth mentioning that Ponec and Cooper previously found for LiH that depending on the AIM technique chosen extra structure could be introduced to the SEDI plot.10 They suggested that this additional structure might be chemically significant. More precisely, when using QTAIM, a maximum in SEDI appeared for a value of R in the vicinity of significant changes in the wave function. When using the Mulliken AIM technique,9 no such maximum was retrieved.10 In the case of CO, the shape of the SEDI versus internuclear distance curve also depends significantly on the AIM method used, as is shown in Figure 4. Again, some extra structure can indeed appear when using QTAIM. The maximum in the QTAIM plot appears near a separation of 1.5 Å, which is somewhat larger than the position of the minimum in the potential energy surface (ca. 1.15 Å). In this case, no major changes in the wave function appear in the vicinity of the observed maximum of the SEDI curve. Inspection of Figure 4 confirms that the appearance of a maximum depends critically on the choice of AIM partitioning. The maximum is only found when using the QTAIM scheme, which naturally corresponds to the most extreme possible partitioning of a molecular density function into atomic parts.

Influence of AIM on SEDI and DAFH

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8759

Figure 4. Comparison of different AIM methods for the SEDI versus internuclear distance (R, in angstroms) curves for CO. Figure 6. Variation with internuclear distance (R, in angstroms) of DAFH eigenvalues associated with the broken valence of the σ component of CO.

Figure 5. Variation with internuclear distance (R, in angstroms) of DAFH eigenvalues associated with the broken valence of the π component of CO.

We find that the eigenvalues associated with the broken valences of the π and σ components of the CO bond also display such drastically different behavior between different AIM methods. To demonstrate these differences, the results from the Hirshfeld-I and QTAIM methods for the broken valence of the CO π bond are depicted in Figure 5: additional structure is again observed for QTAIM. The absence of the maximum when using other AIM techniques may indicate that the maximum in the QTAIM SEDI curve is not due in this case to underlying physical or chemical effects but rather to the very nature of QTAIM. To investigate whether this is related to the binary character of QTAIM, a new binary procedure, coined Hirshfeld-IDom, was derived from Hirshfeld-I in the following way. Once selfconsistency is reached in the Hirshfeld-I scheme, an extra step is introduced, for which at every point in space r, the weight wA(r) in eq 1 for the atom with the highest weight is changed to unity, and the weight for all other atoms is put to zero. The “Dom” of Hirshfeld-IDom is short for “dominant”. Figure 5 shows very clearly that this Hirshfeld-IDom procedure leads to extra structure that is reminiscent of that observed with QTAIM. This leads to the conclusion that the extra structure is indeed related to the use of disjoint atomic domains. It is also worth noting that the SEDI plot from the binary procedure derived from Hirshfeld-I shows similar behavior to the QTAIM case, that is, showing some extra structure. This is also related to recent results by Heyndricx et al.,39 where it was found that SEDI involving multiple centers, so-called multicenter indices,40,41 only show chemically appealing structure when using disjoint atomic domains. It was also

found in that case that gradually making the AIM method more binary introduces the chemically appealing structure. Even so, the fact that it only occurs for binary AIM methods makes attaching chemical or physical importance to the structure in the SEDI curves debatable. An important conclusion that can be deduced from Figure 5 is that although the choice of the particular partitioning scheme affects the actual form of the plots, the differences diminish with increasing internuclear distance, and, for sufficiently large R, the values for both C and O converge to the value of unity anticipated for the covalent dissociation of a π bond of CO. The deviations from unity of the π eigenvalues for shorter R indicate that for distances close to equilibrium the electron pair of the π bond is shared unevenly. Indeed, the actual values of the eigenvalues associated with the DAFH eigenvectors for the C and O domains allow an estimation of the contributions of these atoms to the unevenly shared electron pair of the π bond. As can be seen, the sum of the complementary eigenvalues from the two domains remains remarkably close to two. Such an interpretation is also straightforwardly corroborated by the striking complementarity of the pair of eigenvalues corresponding to the broken valences of the σ bonds, consistent with the shared electron pair nature of both the σ and π components of the bonding in CO. As expected, the most uneven sharing of the π-bonding pair results from the binary QTAIM partitioning, whereas the remaining schemes more or less consistently predict less polarity of the CO π bond. A similar difference in the dependence on R is found for the DAFH eigenvalues corresponding to the broken valences of the σ bond. There is, however, an important caveat. Because DAFH analysis shows one single lone pair in the valence space on each atom and the π bonds are split evenly over the two atoms, both electrons from the σ bond should become attached to the oxygen atom at large R, even though the molecule dissociates to neutral atoms. This is indeed what is found in Figure 6, and, despite the differences between the different AIM methods for the σ bond eigenvalue near equilibrium, all methods agree that at large R both electrons should be attached to oxygen. In contrast with the important impact of the particular partitioning schemes on the R dependence both of SEDI and of DAFH eigenvalues, the shapes of the associated DAFH eigenvectors exhibit remarkable insensitivity to the actual choice of the AIM definition, and so the description of the dissociation provided by the different approaches remains essentially un-

8760

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Bultinck et al.

Figure 7. CO dissociation: Hirshfeld-I DAFH eigenvalues and eigenvectors over the domains of C and O as a function of internuclear distance, R. CO dissociation: Hirshfeld-I DAFH eigenvalues and eigenvectors over the domains of C and O as a function of internuclear distance, R.

Influence of AIM on SEDI and DAFH

Figure 8. Nitrogen atomic charge q(N) as a function of internuclear distance, R (in angstroms), in NO+.

changed. To demonstrate the effect of the systematic variation of the internuclear distance on the shape of the DAFH eigenvectors and to connect the picture of electron reorganization emerging from the numerical data in Figures 4-6 with visual insights from DAFH eigenvectors, the progress of dissociation of CO, as monitored by the DAFH analysis using the Hirshfeld-I partitioning scheme, is displayed in Figure 7 for different values of R. Practically the same description results from all of the partitioning procedures. (See Figures S6-S8 in the Supporting Information.) As is clear from Figure 7, the general trends in the individual DAFH eigenvectors are analogous to those observed earlier for N2. For short internuclear distances close to equilibrium, the eigenvectors corresponding to broken valences of the σ and π bonds are reminiscent of the corresponding occupied molecular orbitals, and the corresponding eigenvalues provide information about the contributions of C and O to the unevenly shared electron pairs of the polar bonds. This picture of the bonding changes, however, with the elongation of the internuclear distance, and the increasing localization of the eigenvectors onto individual atoms is again consistent with the picture of the progressive splitting of the bonding electron pairs of σ and π bonds. Unlike the case of N2, the progress of the dissociation appears to be reflected to different extents by the DAFHs averaged over C and O atoms, but this could be partially due to the fact that the corresponding complementary pictures are displayed using the same common isovalue so that because of different orbital exponents on C and O the shape of the eigenvectors in the vicinity of individual nuclei can be slightly biased. Completely consistent with the qualitative expectation is also the picture of the bonding for internuclear distances close to complete dissociation where the populations and the shapes of individual eigenvectors completely correspond to the anticipated covalent dissociation into neutral C and O atoms. Another polar system that is worthy of study here is the NO+ molecule: it is isoelectronic with N2 and CO but also shows interesting charge shifts during dissociation. The overall positive charge on this ion is predominantly located on the nitrogen atom near equilibrium geometry, whereas the molecule dissociates into N and O+.42,43 As is shown in Figure 8, all AIM methods succeed in yielding the qualitatively correct charge on the N atom near equilibrium and at dissociation. The charges derived from individual AIM methods start to deviate significantly for short R, with the most dramatic difference corresponding to QTAIM. Indeed, QTAIM exhibits remarkable behavior at smaller internuclear distances, with nitrogen charges reaching highly positive values. Much the same was observed previously when using variationally optimized density matrices, and so

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8761

Figure 9. SEDI versus internuclear distance (R, in Å) for NO+.

Figure 10. SEDI versus internuclear distance (R, in angstroms) for LiF.

these high QTAIM charges are certainly not due to the use of the present CASSCF calculations.43 The behavior of the SEDI as a function of internuclear distance is very reminiscent of CO, as is shown in Figure 9. Again, it is only with QTAIM that extra structure appears. Interestingly, the maximum in the QTAIM SEDI for NO+ appears near a separation of 1.10 Å, which is somewhat shorter than the equilibrium distance of 1.55 Å, and also not in a region that drastic changes in the wave function are expected. We conclude that any supposedly chemically meaningful reasons for the position of the SEDI maxima for LiH,10 and for other molecules, are at best debatable. The close resemblance in the SEDI plots between the isoelectronic CO and NO+ molecules is also reflected in a close similarity of the DAFH descriptions of the two dissociation processes. The DAFH eigenvectors and eigenvalues for NO+ are available as Supporting Information (Figures S11-S14). The first qualitative indication of the resemblance is that the DAFH analysis yields the same the number of dominant eigenvalues (and associated eigenvectors) for the two systems. In the case of the hole averaged over the N atom of the NO+ ion, the analysis yields four dominant eigenvalues, as in the case of the hole averaged over the C atom in CO. The same is also true for the holes averaged over the O atoms in NO+ and CO, for which four dominant eigenvalues are detected in each case. Consistent with what was observed above for N2 and for CO, the choice of the particular AIM method has only a marginal effect on the shapes of the corresponding DAFH eigenvectors. Instead, the most noteworthy changes occur in the DAFH eigenvalues. Similarly to what was found for CO and reported in Figures 5 and 6, the QTAIM-based DAFH eigenvalues for a σ and π bond (Figures S15 and S16 in the Supporting Information) show a

8762

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Bultinck et al.

Figure 11. QTAIM DAFH eigenvectors and eigenvalues at 1.20 Å for LiF.

clearly different behavior as a function of the internuclear distance relative to those from other AIM techniques, except that the complementarity conditions are maintained. Turning now to LiF as an example of a highly polar molecule, a somewhat peculiar SEDI plot is obtained. As Figure 10 shows, the SEDI first decreases with increasing internuclear distance and then increases again before dropping dramatically around roughly 4 Å. This sort of remarkable behavior is well-known for LiF: it is due to the crossing of a covalent and ionic surface.44 This drastic change in the nature of the wave function is clearly reflected in the SEDI plot. The main features are basically the same, but, strangely, all of the Hirshfeld-related procedures indicate that the dramatic change in SEDI occurs at the same point, whereas QTAIM puts it at a significantly larger internuclear distance and also makes the transition somewhat more smooth. The underlying wave function is naturally the same for all of the AIM methods used, and thus one expects the drastic changes also to occur at the same geometry, which is clearly not the case here. DAFH analysis for LiF detects the existence of a very polar σ bond (as well as three lone pairs on fluorine) that is consistent for shorter R with predominantly ionic bonding that is close to the idealized limit of the interaction between Li+ and F-. As an example, Figure 11 shows the results of DAFH analysis at a separation of 1.20 Å when using QTAIM; there is a clear polarization toward the F atom for all of the eigenvectors with significant eigenvalues. (See Figures S17-S20 in the Supporting Information for DAFH eigenvectors on both atoms for all interatomic distances and AIM methods.) Conclusions The effect on SEDIs and on the results of DAFH analysis of using different methods to describe the AIM has been examined. Using a test set representing a homonuclear diatomic (N2), a

polar neutral molecule (CO), a charged polar molecule (NO+) exhibiting significantly different charge distributions among the atoms at equilibrium distance compared with the dissociation limit, and a highly polar molecule (LiF), it was found that the different AIM methods can yield significantly different SEDI versus internuclear distance plots. Only in the case of Bader’s QTAIM method could extra structure be found in the plots, and we concluded that the positions of any maxima could not be assigned any chemical significance. The extra structure was found to be related to the use of procedures yielding disjoint atomic domains in a molecule: a binary procedure derived from Hirshfeld-I leads to results that are more in line with those from QTAIM. Only in the case of LiF, for which a sudden drop in SEDI appears for all AIM methods, could one reasonably infer physical significance: it is well-known for this molecule that the nature of the wave function changes drastically. However, the fact remains that the QTAIM SEDI values show the dramatic decrease at a significantly longer internuclear distance than when using any of the other AIM approaches. The eigenvalues of the DAFH eigenvectors can also differ significantly between different AIM methods, although they tend to the same dissociation limit values, and complementarity conditions for the broken valences of a given bond are preserved. Practically no visually observable differences appear between the eigenvectors obtained using the different AIM techniques so that information gleaned by visual inspection of DAFH functions is essentially unchanged. Acknowledgment. P.B. acknowledges the Fund for Scientific Research - Flanders (FWO-Vlaanderen) for continuous support to his group. R.P. acknowledges the support of this work from the Grant Agency of the Czech Republic (grant no: 203/118/ 2009).

Influence of AIM on SEDI and DAFH Supporting Information Available: DAFH eigenvectors and eigenvalues as a function of internuclear distance for all molecules and all AIM methods as well as plots of σ and π DAFH eigenvalues as a function of internuclear distance. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Dalton, J. A New System of Chemical Philosophy; R. Bickerstaff: London, 1808. (2) Coulson, C. A. Valence, 2nd ed.; Oxford University Press: Oxford, U.K., 1961. (3) Moffitt, W. E. Proc. R. Soc. London A 1950, 202, 534. (4) Moffitt, W. E. Proc. R. Soc. London A 1950, 202, 548. (5) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Clarendon: Oxford, U.K., 1990. (6) Bader, R. F. W. Chem. ReV. 1991, 91, 893. (7) Popelier, P. Atoms In Molecules: An Introduction; Prentice-Hall: Harlow, Great Britain, 2000. (8) Parr, R. G.; Ayers, P. W.; Nalewajski, R. F. J. Phys. Chem. A 2005, 109, 3957. (9) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (10) Ponec, R.; Cooper, D. L. THEOCHEM 2005, 727, 133. (11) Jensen, F. Introduction to computational chemistry, John Wiley & Sons: Chichester, Great Britain, 1999. (12) Cioslowski, J.; Karwowski, J. In Fundamentals of Molecular Similarity; Carbo-Dorca, R., Girones, X., Mezey, P. G., Eds.; Kluwer Academic: Dordrecht, The Netherlands, 2001, p 101. (13) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129. (14) Davidson, E. R.; Chakravorty, S. Theor. Chim. Acta 1992, 83, 319. (15) Bader, R. F. W.; Matta, C. F. J. Phys. Chem. A 2004, 108, 8385. (16) Matta, C. F.; Bader, R. F. W. J. Phys. Chem. A 2006, 110, 6365. (17) Bultinck, P.; Van Alsenoy, C.; Ayers, P. W.; Carbo´-Dorca, R. J. Chem. Phys. 2007, 126, 144111. (18) Bultinck, P. Faraday Discuss. 2007, 135, 244. (19) Bultinck, P.; Ayers, P. W.; Fias, S.; Tiels, K.; Van Alsenoy, C. Chem. Phys. Lett. 2007, 444, 205. (20) Van Damme, S.; Bultinck, P.; Fias, S. J. Chem. Theory Comput. 2009, 5, 334. (21) Lillestolen, T. C.; Wheatley, R. J. Chem. Commun. 2008, 45, 5909. (22) Bultinck, P.; Cooper, D. L.; Van Neck, D. Phys. Chem. Chem. Phys. 2009, 11, 3424.

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8763 (23) Lillestolen, T. C.; Wheatley, R. J. J. Chem. Phys. 2009, 131, 144101. (24) Ponec, R. J. Math. Chem. 1997, 21, 323. (25) Ponec, R. J. Math. Chem. 1998, 23, 85. (26) Ponec, R.; Duben, A. J. Comput. Chem. 1999, 20, 760. (27) Ponec, R.; Yuzhakov, G.; Cooper, D. L. Theor. Chem. Acc. 2004, 112, 419. (28) Ponec, R.; Cooper, D. L. J. Phys. Chem. A 2007, 111, 11294. (29) Ponec, R.; Cooper, D. L. Faraday Discuss. 2007, 135, 31. (30) Cioslowski, J. Int. J. Quant. Chem. Symp. 1990, S24, 15. (31) Cioslowski, J. J. Math. Chem. 1991, 8, 169. (32) Wiberg, K. B. Tetrahedron 1968, 24, 1083. (33) Giambiagi, M.; Giambiagi, M. S.; Grempel, D. R.; Heymann, C. D. J. Chim. Phys. 1975, 72, 15. (34) Mayer, I. Chem. Phys. Lett. 1983, 97, 270. (35) Bader, R. F. W.; Stephens, M. E. J. Am. Chem. Soc. 1975, 97, 7391. (36) Werner, H.-J.; Knowles, P. J.; Lindh, R.; Manby, F. R.; Schu¨tz, M.; Celani, P.; Korona, T.; Mitrushenkov, A.; Rauhut, G.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hetzer, G.; Hrenar, T.; Knizia, G.; Ko¨ppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pflu¨ger, K.; Pitzer, R.; Reiher, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. MOLPRO, version 2009.1; Cardiff, U.K., 2009. (37) Biegler-Ko¨nig, F. W.; Bader, R. F. W.; Tang, T. H. J. Comput. Chem. 1982, 13, 317. (38) Fradera, X; Austin, M. A.; Bader, R. F. W. J. Phys. Chem. A 1999, 103, 304. (39) Heyndrickx, W.; Matito, E.; Salvador, P.; Sola, M.; Bultinck, P. J. Comput. Chem., submitted for publication. (40) Bultinck, P.; Ponec, R.; Van Damme, S. J. Phys. Org. Chem. 2005, 18, 706. (41) Bultinck, P.; Rafat, M.; Ponec, R.; Carbo´-Dorca, R.; Van Gheluwe, B.; Popelier, P. J. Phys. Chem. A 2006, 110, 7642. (42) Van Aggelen, H.; Bultinck, P.; Verstichel, B.; Van Neck, D.; Ayers, P. W. Phys. Chem. Chem. Phys. 2009, 11, 5558. (43) Van Aggelen, H.; Verstichel, B.; Ayers, P. W.; Bultinck, P.; Cooper, D. L.; Van Neck, D. J. Chem. Phys. 2010, 132, 114112. (44) Bauschlicher, C. W.; Langhoff, S. R. J. Chem. Phys. 1988, 89, 4246.

JP101707W