Performance of the Density Matrix Functional Theory in the Quantum

ARTICLE pubs.acs.org/JPCA

Performance of the Density Matrix Functional Theory in the Quantum Theory of Atoms in Molecules Marco García-Revilla,* E. Francisco, A. Costales, and A. Martín Pendas* Departamento de Química Física y Analítica, Facultad de Química, Universidad de Oviedo, 33006 Oviedo, Spain ABSTRACT: The generalization to arbitrary molecular geometries of the energetic partitioning provided by the atomic virial theorem of the quantum theory of atoms in molecules (QTAIM) leads to an exact and chemically intuitive energy partitioning scheme, the interacting quantum atoms (IQA) approach, that depends on the availability of second-order reduced density matrices (2-RDMs). This work explores the performance of this approach in particular and of the QTAIM in general with approximate 2-RDMs obtained from the density matrix functional theory (DMFT), which rests on the natural expansion (natural orbitals and their corresponding occupation numbers) of the first-order reduced density matrix (1-RDM). A number of these functionals have been implemented in the promolden code and used to perform QTAIM and IQA analyses on several representative molecules and model chemical reactions. Total energies, covalent intra- and interbasin exchange-correlation interactions, as well as localization and delocalization indices have been determined with these functionals from 1-RDMs obtained at different levels of theory. Results are compared to the values computed from the exact 2-RDMs, whenever possible.

1. INTRODUCTION More than 20 years have passed since the publication of R. F. W. Bader’s Atoms in Molecules,1 a text that marked a clear inflection point in the evolution of modern theories of the chemical bond. The quantum theory of atoms in molecules (QTAIM), as the theory is now widely known, has provided the core of real space analyses of wave functions, contributing to revitalize the study of reduced density matrices, and to join traditionally separated communities of researchers. A realm where the QTAIM has been particularly successful is the partitioning of physical observables into atomic or group contributions. By providing well-defined quantum mechanical subsystems in real space, the theory allows any expectation value to be written as a sum of one-, two-, or in general, many-basin contributions. Traditionally, only one-electron operators, such as the kinetic energy or the electronnucleus attraction, have been considered, so most standard practitioners of the QTAIM stick to one-basin additive decompositions. Since total electronic energies obviously contain two-electron contributions via interelectron repulsions, atomic or group energies have been defined by resort to the atomic virial theorem,2 which relates kinetic and potential terms at stationary points on potential energy surfaces, in exactly the same way as the ordinary virial theorem does. In this way, EA, the energy of quantum atom or group A, is introduced as its minus kinetic energy, EA = TA. These standard energies have been successfully used3 to theoretically justify the nature of the additivity rules taught, for instance, in organic chemistry.4 In spite of the fact that any process in which an atom’s kinetic energy increases is interpreted in the virial view as accompanied by r 2011 American Chemical Society

a concomitant energetic stabilization, there are well-known situations (steric crashes, for instance) where kinetic energy rises have been rationalized in terms of energetic destabilizations, not the contrary. This has generated an enduring controversy in the specialized literature.58 On the other hand, chemists are used to thinking in terms of interacting entities and interaction energies rather than resorting to the use of the atomic or group energies based on the atomic virial theorem. However, it is still possible to reconcile the more standard use of the QTAIM with this interacting picture by freeing the former from the one-basin decomposition constraint, directly obtaining the intra- and interbasin partitioning of the electron repulsion energy. In this way, a flexible, chemically intuitive, and exact energy decomposition is obtained that goes beyond (and is independent of) the fulfillment of the atomic virial theorem. All these ideas have been successfully implemented in the QTAIM within the so-called interacting quantum atoms (IQA) approach.9,10 All the QTAIM quantities are borrowed by IQA, which has been able to deal with the practical computation of the most difficult ones, namely, the above-mentioned one- and two-basin electron repulsion interaction terms. This has allowed the QTAIM both the possibility to be applied at arbitrary molecular geometries (not necessarily the equilibrium ones) and, more important, to restore the pair potential language that has been so useful in atomistic simulations in the theory of chemical bonding.

Received: April 29, 2011 Revised: September 12, 2011 Published: September 27, 2011 1237

dx.doi.org/10.1021/jp204001n | J. Phys. Chem. A 2012, 116, 1237–1250

The Journal of Physical Chemistry A As is clear from the last paragraph, an essential ingredient of the improved QTAIM energy partitioning scheme implemented in IQA is the real space partitioning of the electronelectron repulsion terms. A proper implementation of their computation needs from the (diagonal) second-order reduced density matrix, 2-RDM, F2(r1,r2). This limits its applicability to wave function approaches, either in the single- (e.g., HartreeFock) or multideterminant (configuration interaction, CI, complete active space, CAS) case. Neither perturbation nor standard density functional theory (DFT) methods provide a well-defined F2, and IQA cannot be applied on them. Notice that the most accurate methods available at this moment for the calculation of the electronic structure of medium size molecules, such as CCSD(T) or CASPT2, fall into this category. Moreover, 2-RDMs are not routinely available in standard computational packages, making it a painful task to extract them, if at all possible. A way to examine chemical bonding issues from highly correlated electronic structure methods relies on approximately reconstructing the 2-RDM. From the several possibilities available at this moment, the use of density matrix functional theory1127 (DMFT) is appealing. DMFT is a very active research field that, for our present purposes, writes approximate expressions for the 2-RDM in terms of functionals (DMFs) of the first-order reduced density matrix (1-RDM). In particular, from the occupation numbers, ni, and natural orbitals, χi, of its natural expansion. The aim of this work is to examine the performance of several DMFs within QTAIM. We will, on the one hand, study how well several DMFs are able to reconstruct each of the QTAIM intra- and interbasin electronelectron repulsions, as well as their associated intrabasin electron localization (LI) and interbasin electron delocalization (DI) indices, for a set of test systems for which the exact 2-RDM is also known. This may shed some light on the nature of the errors in DMFs and, perhaps, suggest ways to improve them. On the basis of what we will learn with these examples, we will present approximate DMFT QTAIM decompositions on MP2 and CC pseudo-1-RDMs, commenting on their plausible accuracies. The layout of the paper is as follows. First, we will define in section II the QTAIM quantities to be analyzed in the rest of the article. We will then discuss in section III some of the DMFs that will be used to reconstruct 2-RDMs from the natural expansion of the first order density. In section IV we will present the test systems on which we have made calculations. Then we will show in section V the performance of DMFs with respect to exact calculations made with CAS wave functions. In light of those results, section VI will be devoted to analyze the results on MP2, CCSD, and CISD calculations, for which either no exact results exist (the first two cases), or the exact calculation is not presently implemented in our code (CISD). We will end with a summary and some conclusions.

II. QTAIM DEFINITIONS AND THE IQA APPROACH Within QTAIM, the following one- and two-basin partitioning of the molecular energy is used:9,10 " # Z Z B ^ ∑ E¼∑ dr1 T F1 ðr1 ; r10 Þ r A ΩA B 1B Z Z 1 F ðr1 , r2 Þ ZA ZB dr1 dr2 2 þ ∑ ð1Þ þ ∑ 2 A, B ΩA r 12 ΩB A > B RAB where ΩA represents the atomic basin of atom A, and F1(r1;r10 ) and F2(r1;r2) are the 1-RDM and 2-RDM, respectively. Taking over the

ARTICLE

above partition, the IQA approach of QTAIM joins energetic terms such that chemically meaningful objects appear, in the light of McWeeny’s theory of electronic separabilitiy.28 Thus eq 1 becomes E ¼ ¼

∑A ðTA ∑A EAself

AA þ Ven þ VeeAA Þ þ

þ

∑ ðVnnAB þ VenAB þ VenBA

A>B

∑ EABint

þ VeeAB Þ

ð2Þ

A>B

where we have used an easy to decrypt nomenclature, in which A,B represent atoms, i.e., atomic basins plus their corresponding nuclei.10,29,30 Intrabasin contributions define an atomic self-energy, and all interbasin ones define the pairwise-additive interaction energy between pairs of atoms. Interactions are read in the chemical scale by decomposing31 F2 into Coulombic and exchange-correlation contributions: F2 ðr1 ; r2 Þ ¼ Fðr1 ÞFðr2 Þ Fxc 2 ðr1 ; r2 Þ

ð3Þ

AA (F(r1) = F1(r1;r1)) so that VAA ee can be decomposed as Vee = AA AA + V , where V is a purely Coulombic term: VAA C xc C Z Z 1 Fðr1 ÞFðr2 Þ VCAA ¼ dr1 dr2 ð4Þ 2 ΩA r12 ΩA

and VAA xc is the exchange correlation contribution to the self-energy of atom A: Z Z 1 Fxc AA 2 ðr1 ; r2 Þ dr1 dr2 ð5Þ Vxc ¼ 2 ΩA ΩA r12 In the same way, each interatomic interaction EAB int may be divided into a classical AB AB BA þ Ven þ Ven þ VCAB Þ VclAB ¼ ðVnn

ð6Þ

and a purely quantummechanical exchangecorrelation term, VAB xc : Z Z Fxc ðr1 ; r2 Þ AB Vxc ¼ dr1 dr2 2 ð7Þ r12 ΩA ΩB AB VAA C and VC collect all the electrostatic intraatomic and interatomic electronelectron interactions, respectively. On its turn, VAB cl includes a longrange term that contains the pointcharge and multipolar interactions

Vcl,ABlong-range ¼

∞

B Q Ω1 A 1 Q Ω l2 2

∞

∑ ∑ Cl m l m ðR^ Þ Rll mþ l þm1 l m l m 1

1

1 2

1 2

2

1

2

ð8Þ

2

^ ) is a coupling coefficient,32 Q Ω where Cl1m1l2m2 (R lmA are the spherical atomic multipoles defined as 1=2Z 4π A ¼ r l Slm ð^r ÞFðrÞdr ð9Þ QΩ lm 2l þ 1 ΩA and Slm (^r ) is a real spherical harmonic. Also, a generally small short-range interaction VAB cl,short-range is present that measures how strongly the atoms A and B overlap in the sense defined in ref 32. On the other hand, the quantummechanical term VAB xc is associated with a covalent interaction between atoms A and B. In fact, the DI between atoms A and B, δA,B, which is given by an expression very similar to eq 7, Z Z A, B δ ¼2 dr1 dr2 Fxc ð10Þ 2 ðr1 ; r2 Þ ΩA

1238

ΩB

dx.doi.org/10.1021/jp204001n |J. Phys. Chem. A 2012, 116, 1237–1250

The Journal of Physical Chemistry A

ARTICLE

is a standard descriptor of the covalent bond order between both atoms in real space techniques,33,34 while the LI, λA, given as Z Z A λ ¼ dr1 dr2 Fxc ð11Þ 2 ðr1 ; r2 Þ ΩA

ΩA

determines the number of electrons that atom A does not share with other atoms, i.e., that do not appreciably delocalize to regions of space beyond ΩA. Further evidence of the close connection between covalency and the exchange-correlation interaction energy has recently been found.35 As different examples have shown, the bond paths of the QTAIM signal preferred quantum-mechanical exchange channels or, in other words, that the existence of a bond path between an atom A and the atoms B in its environment appears to be largely determined by the relative values of the different AB exchange-correlation interaction energies. To evaluate all the energetic contributions except VAA xc and AB Vxc , we only need the 1-RDM, usually given in terms of the natural molecular orbitals (NMOs). However, the exchangecorrelation terms, as well as the LI and DI also require Fxc 2 (r1;r2). In the case of HartreeFock mono-determinantal wave functions, Fxc 2 (r1;r2) is given by the FockDirac exchange function Fxc 2 ðr1 ; r2 Þ ¼ F1 ðr1 ; r2 ÞF1 ðr2 ; r1 Þ

ð12Þ

and the evaluation of the exchange energies, LIs, and DIs, may be carried out using the algorithms described in ref 32. These algorithms may also be used with multi-determinantal wave functions provided that Fxc 2 (r1;r2) in terms of the canonical orbitals is known.31 Most times, however, this six-dimensional function is either unavailable from the output of standard packages or not well-defined in the quantum-mechanical method that is being used. As a consequence, the (exact) computation of A,B AB , and λA for correlated multi-determinantal wave VAA xc , Vxc , δ functions is far from trivial, and has been restricted up to now to systems and methods for which Fxc 2 (r1;r2) is easily available. An interesting possibility not yet explored, which is the main objective of this work, consists of using approximate 2-RDMs derived from 1-RDMs. Actually, some approximate formulas for δA,B in terms of the NMOs and their respective occupation numbers have been used in the literature. For instance, for a closed-shell molecule, Fulton36 proposed δA, B ¼

M=2

∑i, j 4ðni njÞ1=2 SΩij SΩij A

B

ð13Þ

where M is the number of spinorbitals (χi), ni (constrained to A lie 0 and 1) is the occupation number of χi, and SΩ ij = R between * ΩAχi (r)χj(r) dr is an element of the atomic overlap matrix between natural spin-orbitals. The above expression for δA,B has also recently used by Wang, Werstiuk, and Matta as an approximate DI for correlated wave functions.37,38 Even though no reference to the 2-RDM was made by Fulton in his original work, eq 13 for δA,B may be derived by substituting in eq 10 the following expression for Fxc 2 (r1;r2): Fxc 2 ðr1 ; r2 Þ

Baerends17,18 rederived the same approximation from modeling exchange and correlation holes. Their derivation corresponds to a setting of the M€uller parameter equal to 1/2. These and other approximate 2-RDMs are considered in the following section.

=

M=2

∑i, j 2ðni njÞ

1=2

χij ðr1 Þχij ðr2 Þ

ð14Þ

where χij(r) is an abbreviature for χ*i (r1)χj(r1). A generalized form of the above expression that contained one free parameter was formerly suggested by M€uller.11 More recently, Buijse and

III. 2-RDMS BASED ON THE DMFT The exchange-correlation part (Fxc 2 ) of the approximate 2-RDMs based on DMFT are built in terms of the NMOs and their occupation numbers of the molecule using the following general formula: Fxc 2 ðr1 ; r2 Þ ¼

∑i, j f ðni , nj Þχij ðrÞχij ðrÞ

ð15Þ

All of these functionals have been implemented in our promolden code, which has also been used to perform all the calculations of this work. The NMOs (or pseudo-NMOs) and ni’s may have different origins. They may arise, for example, from second-order MollerPlesset (MP2), coupled cluster, or CAS calculations. This last possibility is particularly interesting in this work, since the exact 2-RDM of CI expansions (e.g., CAS wave functions) is also implemented in promolden, allowing us to gauge the quality of the different approximate 2-RDMs by comparing their results with those of the exact 2-RDM corresponding to the same wave function. We will restrict this work to N-electron closed shell molecules with M NMOs (M/2 NMOs with spin α and M/2 NMOs with spin β), and use the convention that the ni’s are real numbers between 0 and 1. The indices i and j in eq 15 run from 1 to M/2, and the case i = j must be explicity included in the summation. Different choices of f(ni,nj) yield different approximate 2-RDMs. For instance, f BB ðni , nj Þ ¼ 2ðni nj Þ1=2

ð16Þ

corresponds to the M€uller 2-RDM (named by this author the corrected Hartree (CH)) or BuijseBaerends functional.17,18 A more recent DMF was proposed by Goedecker and Umrigar14 (GU) as a correction of the M€uller functional, in which f GU ðni , nj Þ ¼ 2ðni nj Þ1=2 2ni ð1 nj Þδi, j

ð17Þ

Another correction to the M€uller functional, named the corrected HartreeFock (CHF) functional, was given by Csanyi and Arias (CA)15 f CA ðni , nj Þ ¼ 2½ni nj þ ½ni ð1 ni Þnj ð1 nj Þ1=2

ð18Þ

In order to deal with the overestimation of the correlation energy of the CHF functional, Csanyi, Goedecker, and Arias16 (CGA) introduce the functional f CGA ðni , nj Þ ¼ ½ni nj þ ½ni ð2 ni Þnj ð2 nj Þ1=2

ð19Þ

A hybrid functional (HYB) was proposed by Staroverov and Scuseria39 that empirically combines the CH and CHF functionals with the following formula: f HYB ðni , nj Þ ¼ c0 f GU ðni , nj Þ þ ð1 c0 Þf CA ðni , nj Þ

ð20Þ

were c0 is a coefficient taken equal to 0.6. In an attempt to correct the overcorrelation of the M€uller functional, while preserving the proper description in the dissociation limit, Gritchenko et al.19 proposed a hierarchy of three levels of repulsive corrections (only the first two will be 1239



ARTICLE

considered here). The first repulsive correction to the BB functional (BBC1) restores the contribution of the cross products of weakly occupied orbitals (i,j > N/2): 8 < þ 2ðni nj Þ1=2 , i, j e N=2 f BBC1 ðni , nj Þ ¼ : 2ðni nj Þ1=2 , i, j > N=2 ði 6¼ jÞ ð21Þ BBC1

The second repulsive correction (BBC2) changes f (ni,nj) for the different strongly occupied natural orbitals χi and χj, where i,j e N/2: 8 1=2 > > < þ 2ðni nj Þ , i, j e N=2 i ¼ j i, j e N=2 ði 6¼ jÞ þ2ni nj , f BBC2 ðni , nj Þ ¼ > > : 2ðni nj Þ1=2 , i, j > N=2 ði 6¼ jÞ ð22Þ The last functional that we will use in this work is the third Piris natural orbital functional (PNOF3),40 very close to the GU approximation f PNOF3 ðni , nj Þ ¼ 2ðni nj Þ1=2 2ni ð1 nj Þδi, j þ 2ðhi hj Þ1=2 λi, j where hi = 1 ni, and ( 1, i, j e N=2 and i 6¼ j λi, j ¼ 0, otherwise

ð23Þ

ð24Þ

We should mention that an improved Piris et al. functional not yet included in our promolden code has recently been published.41 Using eqs 10, 11, and 15, the LIs and DIs are given by λA ¼

∑i, j f model ðni , njÞSΩij SΩij

δA, B ¼ 2

A

ð25Þ

A

∑i, j f model ðni , nj ÞSΩij SΩij A

B

ð26Þ

where model = BB, BBC1, BBC2, CA, GU, CGA, HYB, or PNOF3. Finally, we must note from the outset that all the functionals except GU and PNOF3 satisfy the sum rule Z Z ð27Þ dr1 dr2 Fxc 2 ðr1 ; r2 Þ ¼ N where the integrals are extended over the whole space. In the case of the GU and PNOF3 functionals, the above integral is ∑i2n2i e ∑i2ni = N. As a direct consequence of this behavior, these two functionals do not satisfy δA,B f 0 and VAB xc f 0 when the distance between the nuclei A and B tends to infinity.

IV. TEST SYSTEMS AND COMPUTATIONAL DETAILS To test the performance of the DMFT-reconstructed 2-RDMs shown above, we have examined several molecules at different levels of theory, as well as a pair of model chemical reactions along the reaction coordinate. The systems studied in more detail are a prototype of σ-bonded molecule (H2), and a very polar molecule (LiH). Both of them have been analyzed in a wide range of interatomic distances using approximate 2-RDMs derived

from 1-RDMs that, in turn, were obtained from HF, CAS, MP2, CCSD, or CISD calculations. The analyses corresponding to the HF and CAS wave functions were also carried out with the exact 2-RDMs. In this way, these uncorrelated and correlated wave functions, respectively, may be used to gauge the goodness of the approximate 2-RDMs, while the comparison betwen the results derived from the exact CAS 2-RDMs and those coming from the approximate MP2, CCSD, and CISD 2-RDMs may be useful to measure the quality of the latter. Other diatomic molecules with 12 (BeO, BN, C2, LiF), 14 (BF, CO), and 20 (AlN, MgO, SiC) electrons have also been explored using exact HF and CAS 2-RDMs and approximate CAS, MP2, CCSD, CISD 2-RDMs. In these systems, the internuclear distances were fixed to their theoretical CAS values. In order to explore mainly how the quality of the exchangecorrelation (covalent) interactions computed with the approximate 2-RDMs change with the polarity of the bond, we have finally studied the second-row hydrides (LiH, BeH2, BH3, CH4, NH3, H2O, and HF). As in the above isoelectronic series, the analysis in these molecules (except LiH) were restricted to their theoretical CAS equilibrium geometries. On the other hand, the more relevant energy components and topological indices related to chemical bonding of the reactions (i) H2S h H2 + S and (ii) 12 proton migration of [C6H7]+42 have been analyzed and compared with each other using both exact and approximate 2-RDMs, employing 6-311G++(d,p) basis sets. All the wave functions and/or 1-RDMs used in this work have been generated by means of the GAMESS43 or Gaussian0344 codes, and all the calculations using either exact or approximate 2-RDMs have been performed with our promolden program with 6-311G(d,p) or comparable/higher quality basis sets. A relatively large exploratory work has been done in each case to find the numerical quadrature and the number of radial and angular grid points that are needed to compute the energy components to a precision of about 1.03.0 kcal/mol. In some cases, but not always, a β sphere around each atom, with a radius smaller that the shortest distance to its nearest bonded atom, was used to increase this precision.

V. APPROXIMATE 2-RDM DMFT RESULTS VERSUS EXACT 2-RDM RESULTS USING CAS WAVE FUNCTIONS A. H2. Figure 1 shows the intraatomic, interatomic, and total H H HH0 tot ,V ,V exchange-correlation energies (V xc xc xc ), and the LI (λ ) H,H0 0 and DI (δ ) for the H2 (HH ) molecule calculated with the exact full-CI/6-311G(d,p) 2-RDM, the approximate DMFT 2-RDMs derived from the full-CI/6-311G(d,p) 1-RDM, and the HartreeFock 2-RDM. It is observed that all the approximations except GUPNOF3 display a qualitatively right behavior. tends to 0.3123 VH xc exhibits a maximum at about 2.00 bohr and 0 increases continuhartree in the limit RHH0 f ∞, while VHH xc ously and tends smoothly to zero at long RHH0 0 distances. According to this asymptotic behavior, λH and δH,H tend to 1 and 0 in the dissociation limit, corresponding to a single electron located on each hydrogen atom and no electron pair shared by the two atoms. Regarding the different DMFs, we can see that the BBC1 and BBC2 approximations almost fit the exact profiles of the four properties represented in the Figure, 0 while the BB 0 functional H H,H , but gives works as well as the above two in VHH xc , λ , and δ H Vxc values slightly more negative than those of BBC1 and BBC2. It is also evident that the CA model considerably underestimates 1240



0

0

ARTICLE

0

HH tot H HH H H,H Figure 1. VH values for H2 (HH0 ) molecule according to a full-CI/6-311G(d,p) calculation using the exact xc, Vxc , Vxc = 2Vxc + Vxc , λ , and δ 2-RDM and different approximate 2-RDMs derived from the exact full-CI/6-311G(d,p) 1-RDM.

the intraatomic exchangecorrelation energy and the LI. For instance, at the equilibrium geometry (RHH0 = 1.4 bohr), VH xc and λH in the CA calculation are 14% and 10% smaller than the exact values, respectively. On the contrary, (the absolute values 0 H,H0 in this functional are higher than the exact of) VHH xc and δ interatomic exchange-correlation energy and DI, respectively. In summary, the CA model predicts an overbinded and overbonded H2 molecule. We can also see that, except in VH xc, the CGA and HYB profiles are very similar. The values of the four properties in these two functionals are between those corresponding to the CA and BBC1 or BBC2 approximations. It is clear from eqs 17, 23, and 24 that the GU and PNOF3 models are equivalent for the H2 molecule. At RHH0 distances shorter

than 2.0 bohr, these two functionals predict reasonable values H for the intraatomic properties0 VH xc and 0λ but too small values HH H,H for the interatomic ones Vxc and δ ; i.e., the H2 molecule is poorly bounded and binded in these two models. Moreover, as they do not satisfy eq 27, the molecule dissociates incorrectly, and the four properties show a wrong behavior in the limit RHH0 f ∞. Finally, we observe in Figure 1 that the results calculated using all the approximate 2-RDMs are consistently better than those using the HartreeFock 2-RDM (exactly expressible in this case in terms of the HartreeFock 1-RDM). As the GU and PNOF3 calculations, all the HartreeFock exchange-correlation energies exhibit a wrong behavior in the RHH0 f ∞ limit. Similarly, the 1241



ARTICLE

Li,H H LiH tot Li H LiH Li H Figure 2. VLi values for LiH according to a CAS[2,2]//6-311G(d,p) calculation using the exact xc, Vxc, Vxc , Vxc = Vxc + Vxc + Vxc , λ , λ , and δ 2-RDM and different approximate 2-RDMs derived from the exact CAS[2,2]//6-311G(d,p) 1-RDM. In the VLi xc graph, the average number of electrons in Li atom has also been plotted and should be read on the right scale. 0

HartreeFock λH and δH,H values (not plotted in the figure) are permanently 0.5 and 1.0, respectively. All these facts are due to

the well-known shortcomings of the HartreeFock wave function at long distances for this molecule. 1242



ARTICLE

A,H Figure 3. VAH relative errors for the AHn hydrides (A = Li, Be, B, C, N, O, F) according to different approximate 2-RDMs derived from the xc and δ exact CAS[n,m]//6-311G(d,p) (n active electrons and m active orbitals) 1-RDMs. [n,m] = [4,6] for LiH; [n,m] = [6,7] for BeH2; [n,m] = [8,8] for BH3; and [n,m] = [10,9] for CH4, NH3, H2O, and HF. Errors in the HartreeFock data are given with respect to the CAS calculation using the exact 2-RDM.

H LiH B. LiH. The intraatomic (VLi xc,Vxc), interatomic (Vxc ), and

(Vtot xc

VLi xc

VH xc

VLiH xc )

= + + exchange-correlation energies, the total LIs of Li and H atoms, and the DI for the LiH molecule calculated with the exact CAS[2,2]//6-311G(d,p) 2-RDM, the approximate DMFT 2-RDMs derived from the CAS[2,2]//6-311G(d,p) 1-RDM, and the HartreeFock 2-RDM are shown in Figure 2. A first relevant fact is that, as expected, the magnitude of both the exact and the approximate VLi xc is proportional to the average number of electrons of the Li atom. We observe also that, at the equilibrium geometry (RLiH = 3.0 bohr), all the approximate Li LiH Li,H very close to functionals predict values for VLi xc, λ , Vxc , and δ those derived from the exact 2-RDM. The same happens with λH for all except the GU and PNOF3 models. However, there is a noticeable difference between the different VH xc’s. As in H2, the BB, BBC1, and BBC2 functionals give the best results, followed by CGA and HYB, and finally by CA and GUPNOF3. The Li H Li,H LiH at the limit of the values of VH xc, Vxc , λ , λ , and δ dissociation are correctly predicted by all except the GU and PNOF3 models. In particular, VH xc tends in the limit RLiH f ∞ to a value almost equal to that found in H2 when RHH0 f ∞. The small difference is due to the different level of calculation in both cases. Concerning VLi xc, we can see that the CA model has the same asymptotic limit as the exact calculation when RLiH f ∞, while all the other functionals, again except GU and PNOF3, predict a VLi xc value slightly below the exact one. It is interesting to notice that λLi, the number of electrons exclusive of the Li atom, increases linearly with RLiH from 1.0 up to 4.0 bohr, taking a value almost equal to 2.0 at RLiH = 3.0 bohr, both in the exact and all the approximate functionals. This fact is not surprising and simply highlights that the inner 1s core is well represented by all the functionals. Below RLiH = 3.0 bohr, λLi is smaller than 2.0, showing that even the 1s electrons of Li start to participate in the bonding at very small internuclear separations. Contrarily, the marginally greater than 2.0 value for RLiH > 3.0 bohr means that a very small fraction of the number of Li valence electrons are not shared with the H atom. An analogous discussion can be made with λH. At RLiH = 3.0 bohr, the DI reaches a minimum. For RLiH > 3.0 bohr, λH decreases continuously, taking the (right) asymptotic value of 1.0 at the limit of dissociation in all but the GU and PNOF3 functionals. It is very satisfying to see that all the approximate DFMT models

are able to reproduce the well-known30 maximum (minimum) of δLi,H (VLiH xc ) at about 5.5 bohr. Even though GU and PNOF3 are not an exception in this case, the behavior of δLi,H and VLiH xc with these two functionals is very different from that found in all the other cases. As a consequence of the failure to satisfy eq 27, VLiH xc and δLi,H in the GU and PNOF3 models behave radically wrong at the limit of dissociation, tending to a value above and below zero at the limit RLiH f ∞, respectively. The HartreeFock results shown in Figure 2 are clearly worse than all those obtained using approximate 2-RDMs based on the DMFT. They are reasonably correct below the equilibrium geometry (=3.0 bohr), but soon deteriorate when RLiH increases. As in the H2 molecule, this is a consequence of the wrong dissociation behavior of the HartreeFock wave function. C. The AHn Hydrides (A = Li, Be, B, C, N, O, F). Besides LiH, we have also explored the remaining second-row hydrides AHn (A = Li, Be, B, C, N, O, F) at their respective CAS[n,m]//6311G(d,p) (n active electrons and m active orbitals, [n,m] = [4,6], [6,7], [8,8], and [10,9] for LiH, for BeH2, BH3, and (CH4, NH3, H2O, and HF), respectively) equilibrium geometries. The A,H values are shown in Figure 3. The HartreeFock VAH xc and δ results at the correlated geometries are also plotted for comparison. Different conclusions arise when observing this figure. First, the more polar the molecule is, the better the DMFs work. For instance, the error for VAH xc in LiH is smaller than 1% in all the functionals except GU and PNOF3 (error =2%), whereas this error increases up to 16% in CH4 when the CA functional is used. Similar results can be applied to δA,H. If we focus exclusively on the BB, BBC1, and BBC2 functionals, their quality decreases in the order LiH > BH3 > HF > BeH2 = H2O > NH3 > CH4 in the case of VAH xc , while the order is LiH = HF > BeH2 > BH3 > H2O > NH3 > CH4 in the case of δA,H. It is also clear from the figure that the best performance is precisely achieved with the BB, BBC1, and BBC2 family of functionals. The three approximations are of similar quality in the prediction of δA,H, while BBC2 is slightly better than BBC1, and BBC1 slightly better than AH BB in the case of VAH xc . Roughly speaking, the HartreeFock Vxc results do not reproduce too accurately the correlated numbers using the exact 2-RDM. The good prediction for VOH xc in H2O is probably fortuitous. However, the HartreeFock δA,H’s for all except CH4 and NH3 are relatively good approximations to the 1243



A,H Figure 4. Linear correlation between VAH for the AHn xc and δ hydrides (A = Li, Be, B, C, N, O, F) according to different 2-RDMs derived from the exact CAS[n,m]//6-311G(d,p) (n active electrons and m active orbitals) 1-RDMs. n and m defined as in Figure 3.

corresponding CAS results using the exact 2-RDM, with relative errors smaller than 5%. Finally, it is interesting to observe that, regardless of the functional used and the polarity of the molecule, the interatomic exchangecorrelation energy and the AH DI are linearly correlated. This behavior can be seen in Figure 4 and is already known. D. Diatomic Molecules with 12, 14, and 20 Electrons. The differences with respect to the exact results in the intraatomic and interatomic properties for several AB molecules with 12 (BeO, BN, C2), 14 (BF, CO), and 20 (AlN, NaF, SiC, MgO) electrons are plotted in Figure 5. Since VAxc (λA) and VBxc (λB) expand the same range of values, it is clear from the figure that VAxc and λA, corresponding to the less electronegative atom in the pair (A), are slightly better predicted than VBxc and λB. The CA functional clearly provides too small |VAxc| and |VBxc| values, while the opposite is true for the PNOF3 approximation. As in the AHn hydrides, the BB, BBC1, and BBC2 functionals are the best, this last approximation working particularly well in the intraatomic contributions VAxc and VBxc. A detailed analysis of the data shows that, in general, the best performance of all the DMFs occurs when the AB bond is very polar. This fact seems to be related to the greater proximity of the occupation numbers of the natural MOs, at least at the equilibrium geometry, to the limiting values 0.0 and 1.0 when atoms A and B have very different electronegativities. As these limiting values occur when the system is accurately described by a single Slater determinant, the conclusion is that, as a rule, the best described the molecules are at the HF level the more accurate the DMFT results are. In line with the above statements, the HartreeFock results for NaF have, except in the case of VBxc, negligible differences with respect to the CAS results using the exact 2-RDM. Interestingly, the HartreeFock A,B are relatively values for the two-center properties VAB xc and δ good approximations to the corresponding CAS numbers using the exact 2-RDM in all except the C2 molecule. This system is also the one for which the HartreeFock intrabasin properties VAxc, VBxc, λA, and λB present more differences with respect to the CAS value with the exact 2-RDM. E. Chemical Reactions. We have also analyzed with the DMFs a couple of chemical reactions. The first one consists of the approach of a S atom in its 1D state to a ground state H2 molecule to form an H2S molecule, an emblematic example of the change in connectivity along a chemical process in the QTAIM. The exact 2-RDMs were obtained at the CAS[6,8]//6-311G++(d,p)

ARTICLE

level of theory and the approximate DMFT 2-RDMs derived from the CAS[6,8]//6-311G++(d,p) 1-RDM. The C2v symmetry has been maintained at all points, such that the S atom evolves along a line perpendicular to the H2 axis that passes through its center. Calling d the perpendicular distance between S and the HH axis, d = 1.710 bohr at the equilibrium geometry of the H2S molecule. The second reaction studied is the 12 intramolecular hydrogen migration of the benzenium cation, C6H7+.42 Exact 2-RDMs were calculated at the CAS[6,6]//6-311G++(d,p) level, and the approximate DMFT 2-RDMs derived from the CAS[6,6]//6-311G++(d,p) 1-RDM in points along the intrinsic reaction coordinate (IRC) of the migration process. To shorten the discussion, we will only discuss the intra- and the interatomic results of relevant atoms (H and S atoms of the H2S molecule, and C1 and migrating H atoms of C6H+7 ). The DMFT study in these two model reactions complements in a sense our recent analysis on the evolution of the DIs along the IRC.45 Figure 6 shows the intra- and interatomic exchange-correlaS H SH tion energies (VSxc, VH xc, Vxc ), together with the LIs (λ , λ ) and S,H DI (δ ) for the H2S dissociation process. The first general observation is that profiles of the DMFT intraatomic energies display big deviations from the exact profile at long distances, at variance with the DMFT interatomic energies, that show small deviations with respect to the exact 2-RDM curve for all H2S distances. In particular, the VSxc profile is reproduced qualitatively right only by the GU functional, while the biggest deviation, with ΔVSxc 0.20 hartree, is presented by the BB, BBC1 and CGA functionals, followed by the HYB, BBC2, and CA approximations. The behavior of VH xc is totally different: the biggest difference is displayed by CA (ΔVH xc = 0.030 au), H = 0.010 au), CGA (ΔV followed by HYB (ΔVH xc xc = 0.007 au) H and GU (ΔVxc = 0.003 au), while BB, BBC1 and BBC2 functionals fit the exact profile quite well. In the case of the λS LI, all the approximations show a good behavior except the GU functional, displaying a difference of c.a. one unit at long H2S distances. Regarding the λH LI, all the profiles present deviation patterns of the same order shown by VH xc. Finally, the δS,H DI for the different DMFs agrees fairly well with the exact profile for all the H2S distances. This example illustrates one of the difficulties of 2-RDM reconstructions. At large d values, the H2S system is a singlet combination of a free 1Σ+g dihydrogen molecule and a 1D openshell singlet sulfur atom. As the VS and λS plots show, different DMFs may reproduce the appropiate sum rule that makes λS equal to the number of electrons of a free S, while failing to describe the electronelectron couplings of the local openshell singlet. Thus, at dissociation, most of the functionals predict a spurious stabilization of the S atom, which plausibly approaches the 3P ground-state energy. This remains to be further investigated. The behavior of the intra- and interatomic exchangeC H CH correlation energies (V Cxc , V H xc , V xc ), the LIs (λ , λ ), and the δC,H DI of the 12 intramolecular hydrogen migration of the benzenium cation are shown in Figure 7. It is observed that the differences of the intraatomic exchange-correlation energies of the carbon atom with respect to the 2-RDM reference are smaller than in the case of the sulfur of H 2 S molecule, a fact related to the smaller size of the former. The smallest ΔV Cxc is shown by the CA functional, whereas all the other approximations give negative ΔV Cxc ’s and show deviation peaks at the transition state. The (positive) peak of ΔV Cxc 1244



ARTICLE

A B A,B Figure 5. VAxc, VBxc, VAB differences with respect to the exact values for several AB molecules (AB = BeO, BN, C2, BF, CO, AlN, NaF, SiC, xc , λ , λ , and δ MgO) according to several approximate 2-RDMs derived from the exact CAS[n,m]//TZV(d) 1-RDMs. [n,m] = [8,8] for BeO, BN, C2, AlN, NaF, SiC, MgO; [n,m] = [10,8] for BF; and [n,m] = [10,10] for CO. Relative errors in the HartreeFock data as computed with respect to the CAS calculation that uses the exact 2-RDM are also plotted for comparison.

in the CA calculation is due to the maximum shown by V Cxc at the transition state. With respect to VH xc , we can say that, as expected, the differences are smaller due to the smaller size of the hydrogen atom. Again, the CA approximation shows the best fitting to the exact profile, followed by the HYB and CGA approximations. The GU, BB, BBC1, and BBC2 functionals display the biggest errors, and their behavior is quantitatively similar. The interatomic exchange-correlation energies, VCH xc , are well reproduced by all the DMFs. In particular, the best fitting is shown by the CA functional followed by the CGA, GU, and HYB functionals, whereas the

BB-class functionals show the biggest differences with respect to the exact profile. Regarding the LIs, the best fitting for λC near the transition state is obtained with the HYB functional, whereas at negative and positive IRC’s the CA approximation gives the smallest deviation with respect to the 2-RDM calculation. In the case of λH, the CA functional is the best for all values of the IRC. Finally, the behavior of the δC,H profiles is well reproduced by all the DMFs. The best fitting is showed by CA, followed by the HYB, CGA, GU, and the BB-class functionals, that display the biggest errors with respect to the exact profile. 1245



ARTICLE

SH S H S,H Figure 6. VSxc, VH values for the H2 + S reaction discussed in the text versus the H2S distance, which is taken as the reaction xc, Vxc , λ , λ , and δ coordinate associated to the perpendicular approach/detachment of the S atom to/from the H2 molecule.

A general observation is that the interatomic exchangecorrelation energies are well described by DMFs and that the errors of the intraatomic energies are proportional to the size of the atom. In this sense, the DMFs tested for the discussed reactions are sufficient to reproduce the interactions between pairs of atoms, but are inefficient at recovering the intraatomic exchange-correlation energies: DMFs underestimate the correlation energy of the electrons inside an atom, at variance with the good description of interatomic correlations. Finally, we think the analysis of the intra- and interatomic exchangecorrelation energies is useful to test DMFs. The different approximations can be tested in a few points along the reaction

coordinate of a chemical reaction avoiding the calculation of the entire profile.

VI. PERFORMANCE OF DFMT FUNCTIONALS WITH 2-RDMS DERIVED FROM MP2, CCSD, OR CISD 1-RDMS A. H2 and AHn Hydrides. We have analyzed the quality of the DMFs with the approximate 2-RDMs derived from CCSD, CISD, or MP2 1-RDMs. To avoid repetition of arguments and an excessive length of the article, we will focus in this section only on the discussion of the BBC1 results. Overall, the same comments may be applied to the other functionals as well. Our 1246



ARTICLE

C,H CH C H Figure 7. VCxc, VH values versus the IRC associated with the migration of a H atom from a carbon atom (C) to the nearest carbon xc, Vxc , λ , λ , and δ atom (C0 ) in C6H+7 .

results for H2, LiH, and AHn (A = Be, B, C, N, O, F) hydrides are collected in Tables 1, 2, and 3, respectively. The results obtained for H2 by using the CCSD 1-RDM to construct the approximate BBC1 2-RDM (BBC1//CCSD results) are very close to both the exact full-CI numbers and the BBC1 results inspired on the fullHH0 CI 1-RDM. The intra- (VH xc) and interatomic (Vxc ) interactions in the BBC1//CCSD calculation are predicted with errors smaller than 0.0002 and 0.0014 hartree, respectively. Compaired to the exact full-CI results, this calculation increases the number of localized electrons (λH) by 0.0011, at the expense of a reduc0 tion of 0.0018 in the number delocalized electrons (δH,H ). The BBC1//MP2 model increases the DI of H2, from 0.8505 in

the exact full-CI calculation, up to 0.8989. As a consequence, λH is considerably reduced with respect to all the other correlated (non-HF) calculations. The interatomic exchange-correlation 0 energy, VHH xc , a measure of the degree of covalency of the molecule is thus exacerbated at the BBC1//MP2 level. In LiH, the CCSD//BBC1 and MP2//BBC1 calculations give an energetic binding term, VLiH xc , very similar to that obtained at the CAS[2,2] 2-RDM and CAS[2,2] BBC1 levels. However, the DI is predicted slightly smaller than in these two last calculations. All the correlated calculations offer the same image in relation to the number of localized electrons. This number (1.995) is only marginally smaller than 2.0, clearly showing that the 1s electrons 1247



ARTICLE 0

HH Table 1. Intrabasin (VH ) Exchangexc) and Interbasin (Vxc H H,H0 Correlation Energies, LI (λ ) and DI (δ ) at RHH0 = 1.4 bohr for the H2 (HH0 ) Molecule According to Exact and Approximate BBC1 2-RDMs Derived from HF, CAS[2,2], Full-CI, CCSD, and MP2 Calculations Using a 6-311G(d,p) Basis Set

DMF

calc

VH xc

λH

VHH xc

0

δH,H

0

exact 2-RDM

RHF

0.1972

0.5000

0.2556

1.0000

exact 2-RDM

CAS[2,2]

0.2367

0.5834

0.2246

0.8335

exact 2-RDM

full-CI

0.2462

0.5747

0.2383

0.8508

BBC1 BBC1

CAS[2,2] full-CI

0.2373 0.2469

0.5834 0.5747

0.2247 0.2384

0.8335 0.8508

BBC1

CCSD

0.2461

0.5758

0.2370

0.8490

BBC1

MP2

0.2338

0.5507

0.2459

0.8989

AH Table 3. Intrabasin (VAxc, VH xc) and Interbasin (Vxc ) ExA H change-Correlation Energies, LIs (λ , λ ), and DI (δA,H) for AHn Hydrides (A = Be, B, C, N, O, F) according to CAS[n,m]/ 6-311G(d,p) Calculations (n Active Electrons and m Active Orbitals) Using the Exact 2-RDM (First Entry), the Approximate BBC1 2-RDM Derived from the Exact CAS[n,m]/6311G(d,p) 1-RDM (Second Entry), and the Approximate BBC1 2-RDM Derived from the MP2/6-311G(d,p) 1-RDM (Third Entry)a

VAxc BeH2

BH3

H LiH Table 2. Intrabasin (VLi xc, Vxc) and Interbasin (Vxc ) ExLi H change-Correlation Energies, LIs (λ , λ ) and DI (δLi,H) at RLiH = 3.0 bohr for the LiH Molecule According to RHF 2-RDM (First Entry), CAS[2,2] 2-RDM (Second Entry), CAS[2,2] BBC1 (Third Entry), CCSD BBC1 (Fourth Entry), and MP2 BBC1 Calculations, Respectively, Using a 6-311G (d,p) Basis Set

VLi xc

VH xc

VLiH xc

λ

λ

1.6610

0.4466

0.0350

1.9899

1.8013

Li

H

NH3

H2O

δ

Li,H

0.1926

1.6632

0.4671

0.0383

1.9951

1.7768

0.2275

1.6649

0.4683

0.0383

1.9953

1.7770

0.2270

1.6676

0.4942

0.0392

1.9951

1.7870

0.2173

1.6665

0.4851

0.0383

1.9930

1.7974

0.2093

of this atom are very localized. Regarding λ , the CCSD//BBC1 and MP2//BBC1 values are about 0.01 and 0.02 units greater than in the CAS[2,2] BBC1 calculation, respectively. It is thus clear that δLi,H at the CCSD//BBC1 and MP2//BBC1 levels decrease at the expense of an increase in λH, leaving λLi practically unchanged. In summary, one can say that bonding and binding images resulting from the CCSD and MP2 calculations at the equilibrium geometry of this molecule are very reasonable. In the other AHn (A = Be, B, C, N, O, F) hydrides, the δA,H MP2//BBC1 values closely follow those obtained at the CAS// BBC1 level. An exception is H2O, where δO,H is 0.6205 and 0.4829 in the MP2//BBC1 and CAS//BBC1 calculations, respectively. Actually, δO,H in the MP2//BBC1 calculation is very close to the exact 2-RDM value (0.6368). We have checked that this large difference between MP2//BBC1 and CAS//BBC1 calculations in this molecule can not be attributed to the BBC1 approximation. In fact, all the other approximate DMFs give MP2 δO,H’s between 0.590527 (PNOF3) and 0.665329 (CA). As in LiH, the MP2//BBC1 and CAS//BBC1 δA,H’s in BeH2, BH3, and HF do not differ too much from those found when the exact 2-RDM is used. However, as the covalency of the molecule increases, the approximate DIs progressively separate from the exact ones (NH3, CH4). It is certainly unsatisfactory to see that δA,H in the CAS//BBC1 and MP2/BBC1 calculations is smaller than its corresponding 2-RDM value in NH3, while the opposite happens in CH4. Consequently, it is neither possible to say that the approximate DMFs predict higher values than the exact 2-RDM calculation for a given property nor the contrary. H

CH4

HF

a

VH xc

VAH xc

λA

λH

δA,H

2.3288

0.5094

0.0712

2.0252

1.6387

0.3025

2.3378

0.5109

0.0714

2.0266

1.6349

0.3015

2.3418 3.1137

0.5254 0.4588

0.0708 0.1529

2.0226 2.2857

1.6425 1.2571

0.2913 0.5229

3.1307

0.4595

0.1542

2.2663

1.2549

0.5357

3.1307

0.4797

0.1565

2.2456

1.2505

0.5382

4.8452

0.2197

0.2363

4.4631

0.5097

0.7615

4.8321

0.2233

0.2506

4.2801

0.5167

0.8526

4.8918

0.2285

0.2618

4.3039

0.4923

0.8697

6.7816

0.1185

0.2334

6.8482

0.2400

0.7563

6.7868 6.9085

0.1218 0.1272

0.2438 0.2488

6.7741 6.8062

0.2453 0.2412

0.8051 0.8004

8.5804

0.0663

0.1996

8.4508

0.1195

0.6368

8.7856

0.0357

0.1538

8.8979

0.0589

0.4829

8.8620

0.0650

0.2027

8.5276

0.1078

0.6205

10.4892

0.0295

0.1493

9.4796

0.0462

0.4739

10.5109

0.0323

0.1509

9.4834

0.0502

0.4656

10.6804

0.0375

0.1562

9.4734

0.0567

0.4695

The n and m values are as in Figure 3.

The binding exchangecorrelation interaction VAH xc is well reproduced by the CAS//BBC1 and MP2//BBC1 calculations in BeH2, BH3, and HF, while these two calculations are slightly overbinding in NH3 and CH4. As in the case of δO,H, VOH xc in the CAS//BBC1 calculation (0.1538) is much smaller (in absolute value) than its corresponding 2-RDM value (0.1996), whereas the MP2//BBC1 value (0.2027) differs from the exact 2-RDM one by less than 2.0 kcal/mol. As a general comment, we can say that the DMFs are able to reproduce the main features of the calculations using the exact 2-RDM, particularly the interatomic ones. The best agreement is obtained when the atoms involved in the pair have very different electronegativities, and worsens as they form less polarized bonds. The absolute errors in the intra-atomic exchange-correlation energies and DIs are proportional to the atomic size. Overall, the MP2 calculations, for which it is not possible to construct a wave function, are an interesting alternative to the CAS ones, for which it is possible to generate a wave function. The main difficulty in the latter lies in the construction of the 2-RDM, usually not available from the output of the most widely used computational codes. B. AB Molecules with 12, 14, and 20 Electrons. Finally, we compare in Table 4 the intra- and interatomic exchangeA B correlation energies (VAxc, VBxc, VAB xc ), the LIs (λ , λ ), and A,B the DI (δ ) at the CAS//2-RDM, CAS//BBC1, and CISD/ BBC1 levels for the diatomic molecules previously studied. It is important to recall that the CISD approximation captures a considerably larger part of the correlation energy than the present CAS expansions. 1248



ARTICLE

Table 4. Intrabasin (VAxc, VBxc) and Interbasin (VAB xc ) ExchangeCorrelation Energies, LIs (λA, λB), and DI (δA,B) for Several AB Molecules (AB = BeO, BN, C2, BF, CO, AlN, NaF, SiC, MgO) According to CAS[n,m]/TZV(d) Calculations (n Active Electrons and m Active Orbitals) Using the Exact 2-RDM (First Entry), the Approximate BBC1 2-RDM Derived from the Exact CAS[n,m]/TZV(d) 1-RDM (Second Entry), and the Approximate BBC1 2-RDM Derived from the CISD/TZV(d) 1-RDM (Third Entry)a VAxc BeO

BN

C2

BF

CO

AlN

VAB xc

λA

λB

δA,B

2.3216

8.7388

0.1534

2.0247

9.2856

0.6894

2.3295

8.7674

0.1580

2.0182

9.2791

0.7024

2.3428 3.3691

8.8727 6.8887

0.1721 0.2958

2.0357 3.2010

9.1980 7.6019

0.7683 1.1970

3.4017

6.9638

0.3118

3.1536

7.5545

1.2918

3.3501

6.9430

0.4043

2.9195

7.2159

1.8670

4.9557

4.9557

0.5160

5.1048

5.1048

1.7892

4.9872

4.9872

0.5678

4.9454

4.9454

2.1080

4.9653

4.9653

0.6173

4.7714

4.7714

2.4577

3.4390

10.4451

0.1845

3.7091

9.5029

0.7868

3.4517 3.5364

10.4740 10.6614

0.1891 0.2000

3.7084 3.7157

9.5022 9.4687

0.7881 0.8144

4.4664

8.6371

0.4292

4.1306

8.4350

1.4346

4.4996

8.6755

0.4498

4.0821

8.3864

1.5317

4.5590

8.8521

0.4519

4.0514

8.4125

1.5348

17.8563

6.7487

0.1768

11.5559

7.4249

1.0188

17.9185

6.8807

0.1779

11.5794

7.4485

0.9717

17.8847

6.8714

0.2316

11.2370

7.3162

1.4468

NaF

13.9035 13.9084

10.4956 10.5110

0.0422 0.0429

9.9641 9.9647

9.8274 9.8280

0.2082 0.2070

13.9003

10.6611

0.0470

9.9552

9.7736

0.2407

SiC

19.9735

5.3078

0.2564

12.3838

6.3325

1.2832

MgO

a

VBxc

20.0318

5.3755

0.2944

12.2411

6.1897

1.5687

20.0814

5.2785

0.3469

12.1631

5.8657

1.9715

15.8044

8.5075

0.1243

10.5390

8.6583

0.8026

15.8233

8.5959

0.1236

10.5268

8.6941

0.7783

15.8211

8.6807

0.1454

10.4194

8.6362

0.9435

The n and m values are as in Figure 3.

Overall, the intra-atomic properties are best reproduced for the less electronegative atom in each pair. As expected, the CAS//BBC1 results are generally closer than the CISD/BBC1 numbers to the reference (CAS//2-RDM), even though there are several exceptions to this rule. It can also be seen that the BBC1 functional works fairly well in very ionic systems such as NaF. However, for covalent molecules, the interatomic binding A,B DI are not too accurately reproduced. energy VAB xc and the δ There are, for instance, big differences between the reference CAS//2-RDM numbers and those obtained at the CISD/BBC1 level. See, in particular, the results for BN, C2, AlN, and SiC. Moreover, the dispersion between the CAS//2-RDM, CAS// BBC1, and CISD//BBC1 results for these systems are significantly higher than in the others. Taking C2 as an example, 0 0 and δC,C range from 0.5160 to 0.6173, and from VC,C xc 1.7892 to 2.4577 in going from the accurate CAS//2-RDM calculation to the approximate CISD/BBC1 model. If we plausibly suppose that the behavior of BBC1 is comparable for the CAS and CISD wave functions, the CISD/BBC1 values

might be approximately corrected from the BBC1 deficiencies through the CAS/BBC1CAS/2-RDM difference. Our results highlight a general fact, at least for all the molecules listed in the table: The CISD//BBC1 model increases covalency A,B values. With a few as measured in terms of the VAB xc and δ exceptions, this also implies a decrease in the number of electrons localized exclusively in each of the atoms. A quite general conclusion to be drawn from the analysis of these and other results is that it is not possible to establish general guidelines of the behavior of DMFT approaches as applied to the cases where the 2-RDM is not available or difficult to obtain.

VII. SUMMARY AND CONCLUSIONS The calculation of state-of-the-art chemical bonding indicators for correlated wave functions depends on the availability of 2-RDMs, which properly take into account electronelectron interactions. 2-RDMs are either not well-defined, as, for instance, in perturbative schemes, or unavailable from standard electronic structure packages. This fact severly limits the applicability of modern real space theories of the chemical bond. A possible way out consists of approximately reconstructing 2-RDMs from first-order densities (1-RDMs), a very active research field within DMFT. We have investigated in this work the performance of DFMTreconstructed 2-RDMs within the QTAIM. We have examined both intra- and interbasin exchange-correlation energies, as well as LIs and DIs for a set of test molecules ranging from very ionic to basicaly covalent, and also for few model chemical processes. DMFT approximated exchange-correlation energies and (de)localization indices have been gauged to exact values obtained with explicit 2-RDMs obtained from either full-CI or CAS expansions. As a rule, DMFs work better for highly ionic systems, where correlation corrections are smaller, and for interbasin contributions as opposed to intrabasin ones, supporting the idea that they are only able to correct small deviations of the natural occupancies from either the limiting 0 or 1 values. We have found it difficult to find general trends from our data, but it seems rather clear that the BBC1 or BBC2 functionals have the best overall accuracy and stability. All of the functionals seem to be unable to correctly describe local open-shell regions in real space, as the catastrophic results for our H2S test case in the limit of dissociation shows. We think that the ability of IQA analyses to separate the long and short-range (via inter- and intrabasin contributions) performance of DMFs may be useful to DMFT developers. A note of caution has to be made related to the generalized used of BB DIs for correlated systems after the work of Wang, Werstiuk, and Matta. 37,38 As our work shows, the accuracy of these indices may not be good enough in strong covalent interactions. We have also computed QTAIM quantities from DMFT reconstructed 2-RDMs taken from CCSD and MP2 calculations, where a proper exact reference is lacking. In most cases, the results are pretty similar, and their quality resembles that of the full-CI or CAS DMFT data. A crude cross-comparison of all our calculations shows that static correlation (as introduced by CAS) decreases sharply the covalency of most interactions (except in very ionic compounds, where a slight increase is found), but that dynamic contributions partially restore the single determinant situation. In summary, provided that advices are generally wellcome, we think that BBCn DMFs may provide a reasonable alternative to expensive CI expansions in order to obtain real space chemical 1249


The Journal of Physical Chemistry A bonding indices. Their accuracy has to be taken, nevertheless, with care. This situation will certainly improve in the near future, for new functionals come out every now and again. If this is the case, 6D integrated properties within the QTAIM might be obtained at a 3D computational cost in the near future.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected] (M.G.-R.); angel@ fluor.quimica.uniovi.es (A.M.P.).

’ ACKNOWLEDGMENT The authors acknowledge the financial support from the Spanish MICINN, Project No. CTQ2009-08376, the European Union FEDER funds, the MALTA- Consolider program (CSD2007-00045), and the FICyT Project No. IB09-019. ’ ADDITIONAL NOTE Originally submitted for the “Richard F. W. Bader Festschrift”, published as the November 17, 2011, issue of J. Phys. Chem. A (Vol. 115, No. 45). ’ REFERENCES (1) Bader, R. F. W. Atoms in Molecules; Oxford University Press: Oxford, 1990. (2) Srebrenik, S.; Bader, R. F. W. J. Chem. Phys. 1974, 61, 2536. (3) Bader, R. F. W.; Carroll, M. T.; Cheeseman, J. R.; Chang, C. J. Am. Chem. Soc. 1987, 109, 7968. (4) Cohen, N.; Benson, S. W. The Chemistry of Alkanes and Cycloakanes; Patai, S., Rappoport, Z., Eds.; John Wiley & Sons Ltd.: New York, 1992. (5) Poater, J.; Sola, M.; Bickelhaupt, F. M. Chem.—Eur. J. 2006, 12, 2889. (6) Poater, J.; Sola, M.; Bickelhaupt, F. M. Chem.—Eur. J. 2006, 12, 2902. (7) Poater, J.; Visser, R.; Sola, M.; Bickelhaupt, F. M. J. Org. Chem. 2007, 72, 1134. (8) Bader, R. F. W. Chem.—Eur. J. 2006, 12, 2902. (9) Blanco, M. A.; Martín Pendas, A.; Francisco, E. J. Chem. Theory Comput. 2005, 1, 1096–1109. (10) Martín Pendas, A.; Blanco, M. A.; Francisco, E. J. Comput. Chem. 2007, 28, 161. (11) M€uller, A. M. K. Phys. Lett. 1984, 105A, 446. (12) Hollas, A. Phys. Rev. 1999, A59, 3454. (13) Herbert, M.; Harriman, J. E. J. Chem. Phys. 2003, 382, 142. (14) Goedecker, S.; Umrigar, C. Phys. Rev. Lett. 1998, 81, 886. (15) Csanyi, G.; Arias, T. A. Phys. Rev. 2000, B61, 7348. (16) Csanyi, G.; Goedecker, S.; Arias, T. A. Phys. Rev. 2002, A65, 032510. (17) Buijse, M. A.; Baerends, E. J. Mol. Phys. 2002, 100, 401. (18) Buijse, M. A. Ph.D. Thesis. Vrije Universiteit, Amsterdam, 1991. (19) Gritchenko, O.; Pernal, K.; Baerends, E. J. J. Chem. Phys. 2005, 122, 204102. (20) Baerends, E. J. Phys. Rev. Lett. 2001, 87, 133004. (21) A., M. D. Adv. Chem. Phys. 2007, 134, 21. (22) Kollmar, C.; Heβ, B. A. J. Chem. Phys. 2003, 119, 4655. (23) Lathiotakis, N. N.; Helbig, N.; Zacarias, A.; Gross, E. K. U. J. Chem. Phys. 2009, 130, 064109. (24) Leiva, P.; Piris, M. J. Chem. Phys. 2005, 123, 214102. (25) Leiva, P.; Piris, M. J. Mol. Struct.: THEOCHEM 2006, 770, 45. (26) Gritchenko, O.; Pernal, K.; Baerends, E. J. J. Chem. Phys. 2005, 122, 204102.

ARTICLE

(27) Feixas, F.; Matito, E.; Duran, M.; Sola, M.; Silvi, B. J. Chem. Theory. Comput. 2010, 6, 2736. (28) McWeeny, R. Methods of Molecular Quantum Mechanics, 2nd ed.; Academic Press: London, 1992; Chapter 14. (29) Blanco, M. A.; Martín Pendas, A.; Francisco, E. J. Chem. Theory Comput. 2005, 1, 1096–1109. (30) Martín Pendas, A.; Francisco, E.; Blanco, M. A. Faraday Discuss. 2007, 135, 423–438. (31) Martín Pendas, A.; Francisco, E.; Blanco, M. A. J. Comput. Chem. 2004, 26, 344. (32) Martín Pendas, A.; Blanco, M. A.; Francisco, E. J. Chem. Phys. 2004, 120, 4581. ngyan, J. G.; Loos, M.; Mayer, I. J. Phys. Chem. 1994, (33) A 98, 5244–5248. (34) Fradera, X.; Austen, M. A.; Bader, R. F. W. J. Phys. Chem. A 1999, 103, 304–314. (35) Martín Pendas, A.; Francisco, E.; Blanco, M. A.; Gatti, C. Chem.— Eur. J. 2007, 13, 9362. (36) Fulton, R. L. J. Phys. Chem. 1993, 97, 7516–7529. (37) Wang, Y. G.; Werstiuk, N. H. J. Comput. Chem. 2003, 24, 379–385. (38) Wang, Y. G.; Matta, C.; Werstiuk, N. H. J. Comput. Chem. 2003, 24, 1720–1729. (39) Staroverov, V. N.; Scuseria, G. E. J. Chem. Phys. 2002, 117, 2489. (40) Piris, M. Int. J. Quantum Chem. 2006, 106, 1093. (41) Piris, M.; Matxain, J. M.; Lopez, X.; Ugalde, J. M. J. Chem. Phys. 2010, 133, 111101. (42) García-Revilla, M.; Hernandez Trujillo, J. Phys. Chem. Chem. Phys. 2009, 11, 8425. (43) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (44) Frisch, M. J. et al. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (45) García-Revilla, M.; Popelier, P. L. A.; Francisco, E.; Martín Pendas, A. J. Chem. Theory. Comput. 2011, 7, 1704.

’ NOTE ADDED AFTER ASAP PUBLICATION This paper was published ASAP on September 27, 2011. Due to production errors, the wrong version of the Table of Content and Abstract graphics were published and the paper did not appear in the Richard F. W. Bader Festschrift. The corrected Article was posted on January 20, 2012, with a note indicating the special issue assignment.

1250


Performance of the Density Matrix Functional Theory in the Quantum

Recommend Documents