Average local ionization energy generalized to correlated wavefunctions Ilya G. Ryabinkin and Viktor N. Staroverov Citation: The Journal of Chemical Physics 141, 084107 (2014); doi: 10.1063/1.4893424 View online: http://dx.doi.org/10.1063/1.4893424 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/8?ver=pdfcov Published by the AIP Publishing Articles you may be interested in A self-interaction-free local hybrid functional: Accurate binding energies vis-à-vis accurate ionization potentials from Kohn-Sham eigenvalues J. Chem. Phys. 140, 18A510 (2014); 10.1063/1.4865942 Potential energy curves and electronic structure of 3 d transition metal hydrides and their cations J. Chem. Phys. 129, 214302 (2008); 10.1063/1.2996347 First-principles T -matrix calculations of double-ionization energy spectra of atoms and molecules J. Chem. Phys. 123, 144112 (2005); 10.1063/1.2069907 Orbital-dependent correlation energy in density-functional theory based on a second-order perturbation approach: Success and failure J. Chem. Phys. 123, 062204 (2005); 10.1063/1.1904584 Open-shell localized Hartree–Fock approach for an efficient effective exact-exchange Kohn–Sham treatment of open-shell atoms and molecules J. Chem. Phys. 118, 10439 (2003); 10.1063/1.1560132

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

THE JOURNAL OF CHEMICAL PHYSICS 141, 084107 (2014)

Average local ionization energy generalized to correlated wavefunctions Ilya G. Ryabinkin1 and Viktor N. Staroverov2,a) 1

Department of Physical and Environmental Sciences, University of Toronto Scarborough, Toronto, Ontario M1C 1A4, Canada 2 Department of Chemistry, The University of Western Ontario, London, Ontario N6A 5B7, Canada

(Received 5 May 2014; accepted 7 August 2014; published online 26 August 2014) The average local ionization energy function introduced by Politzer and co-workers [Can. J. Chem. 68, 1440 (1990)] as a descriptor of chemical reactivity has a limited utility because it is defined only for one-determinantal self-consistent-field methods such as the Hartree–Fock theory and the Kohn–Sham density-functional scheme. We reinterpret the negative of the average local ionization energy as the average total energy of an electron at a given point and, by rewriting this quantity in terms of reduced density matrices, arrive at its natural generalization to correlated wavefunctions. The generalized average local electron energy turns out to be the diagonal part of the coordinate representation of the generalized Fock operator divided by the electron density; it reduces to the original definition in terms of canonical orbitals and their eigenvalues for one-determinantal wavefunctions. The discussion is illustrated with calculations on selected atoms and molecules at various levels of theory. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4893424] I. INTRODUCTION

In the context of one-determinantal self-consistent field (SCF) methods it is always possible to define the average local orbital energy function 1   |φ (r)|2 , ρ(r) i=1 i i N

¯ (r) =

(1)

where φi (r) is the spatial part of the ith canonical spinorbital,  i is the corresponding eigenvalue, N is the number  2 of electrons in the system, and ρ(r) = N i=1 |φi (r)| is the 1–5 were the total electron density. Politzer and co-workers first to recognize the significance of this quantity and studied it extensively in the framework of approximate Hartree– Fock (HF) and Kohn–Sham (KS) methods (see Refs. 6–8 for reviews). They interpreted the function I¯(r) = −¯ (r) as the average local ionization energy, that is, the average energy required to remove an electron at point r, and found that it provides useful information about molecular reactivity,6–11 electronegativity,12, 13 polarizability,3, 14 and other properties.6–8 In a parallel development, Nagy and co-workers15–17 studied the concept of local temperature of electronic systems and showed that it is closely related to ¯ KS (r), the average local KS orbital energy. More recently, Bulat et al.18 compared the average local orbital energies constructed from HF and KS atomic orbitals and observed that the difference ¯ KS (r) − ¯ HF (r) mimics the correction term that must be added19–21 to the Slater potential22 to reproduce the exact effective potential for exchange. This observation was explained in our latest work,23, 24 where we found that the quantity defined by Eq. (1) naturally arises in the derivation of a) Electronic mail: [email protected]

0021-9606/2014/141(8)/084107/8/$30.00

accurate analytic approximations to the exact exchange potential. The average local orbital energy emerges, therefore, not only as an interpretive tool of computational chemistry but also as a fundamental ingredient of density-functional approximations. The original definition1 of the average local orbital energy given by Eq. (1) has two limitations. First, it seems to lack invariance with respect to unitary transformations of the occupied orbitals. Bulat et al.18 demonstrated, however, that the lack of invariance is only apparent by deriving an equivalent expression for ¯ (r) which is invariant to unitary transformations of occupied orbitals. The second limitation is that all known expressions for ¯ (r) apply only to one-determinantal SCF methods. Given the increasing prominence of the average local ionization energy in electronic structure theory, it is desirable to extend this concept to correlated-wavefunction methods. In this work, we solve both problems by devising a generalized definition of the function ¯ (r) that is properly invariant with respect to orbital transformations and applies to correlated wavefunctions. For one-determinantal SCF wavefunctions our definition reduces to that of Bulat et al.18 We also show how one can efficiently compute the generalized ¯ (r) for atoms and molecules from a given wavefunction. Because the term “orbital energy” loses its meaning for correlated wavefunctions and because the interpretation of −¯ (r) as an average local electron removal energy is justified only in an approximate sense, we suggest a new name for the function (r): ¯ average local electron energy (ALEE), the rationale for which will be made in Sec. III A. The negative of the generalized ALEE, I¯(r) = −¯ (r), may be thought of as a particular generalization of the average local ionization energy of Politzer and co-workers.1–8

141, 084107-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-2

I. G. Ryabinkin and V. N. Staroverov

J. Chem. Phys. 141, 084107 (2014)

Eq. (10) we have

II. AVERAGE LOCAL ELECTRON ENERGY FOR ONE-DETERMINANTAL SCF METHODS

Our first objective is to cast Eq. (1) in a form that is properly invariant under unitary transformations of the occupied orbitals and is more general than the expressions given by Bulat et al.18 For this purpose it is convenient to adopt the Dirac bra-ket notation (see Ref. 25 for an overview). We begin by noting that canonical orbitals are solutions of the eigenvalue problem Fˆ |φi  = i |φi ,

(2)

where Fˆ is some effective one-electron Hamiltonian. In the HF theory, Fˆ is the Fock operator; in the KS method, it is the KS Hamiltonian. This Hamiltonian is diagonal in the canonical orbital basis, φi |Fˆ |φj  = δij i .

1  r|φi φi |Fˆ |φj φj |r, ρ(r) i=1 j =1 N

N

(4)

where r|φi  = φi (r) and φj |r = φj∗ (r). We recognize that γˆ =

N 

|φi φi |

(5)

i=1

is the one-electron reduced density matrix (1-RDM) operator whose coordinate representation is r|γˆ |r  ≡ γ (r, r ) =

N 

φi (r)φi∗ (r ).

(6)

i=1

This allows us to rewrite Eq. (4) as ¯ (r) =

1 r|γˆ Fˆ γˆ |r. ρ(r)

(7)

As is well known,26 in one-determinantal SCF methods the density matrix at convergence commutes with the effective Hamiltonian, [γˆ , Fˆ ] = 0,

(8)

γˆ 2 = γˆ .

(9)

and is idempotent,

¯ (r) =

1 r|Fˆ γˆ |r. ρ(r)

(10)

Equation (10) is the sought-for definition of the ALEE for all one-determinantal SCF methods made without explicit reference to orbitals. Now let us obtain an explicit formula for the orbitalinvariant form of ¯ (r). Inserting the closure relation into



r|Fˆ |r r |γˆ |r dr .

(11)

is the Hartree (electrostatic) potential of ρ(r), and 

ˆ   = − 1 γ (r, r ) r|K|r 2 |r − r |

(14)

is the coordinate representation of the integral Fock exchange operator. Substituting Eqs. (6) and (12) into Eq. (11) and using the fact that γ (r, r) = ρ(r) we obtain   1 1 ¯ (r) = − ∇r2 γ (r, r ) ρ(r) 2 r =r  1 |γ (r, r )|2  dr . + v(r) + vH (r) − (15) 2ρ(r) |r − r | The first term on the right-hand side of Eq. (15) involves the Laplacian form of the kinetic-energy density, τL (r) = −

 1 2 ∇r γ (r, r ) r =r , 2

(16)

and the last is the Slater averaged exchange-charge potential,22  1 |γ (r, r )|2  dr . vXS (r) = − (17) 2ρ(r) |r − r | Thus, we can write the ALEE of the HF theory as ¯ HF (r) =

τL (r) + v(r) + vH (r) + vXS (r), ρ(r)

(18)

where all orbital-dependent quantities on the right-hand side are defined in terms of occupied HF orbitals and are invariant with respect to their unitary transformations. The analogous expression for the KS ALEE is the same except that vXS (r) is replaced with the multiplicative exchange-correlation potential, ¯ KS (r) =

Using these two properties of γˆ we can cast Eq. (7) as

1 ρ(r)

In the HF method, Fˆ is the Fock operator whose coordinate representation is25   1 ˆ  , r|Fˆ |r  = δ(r − r ) − ∇r2 + v(r) + vH (r) + r|K|r 2 (12) where v(r) is the external potential,  ρ(r2 ) (13) dr vH (r) = |r − r2 | 2

(3)

Using Eq. (3) we can write Eq. (1) as ¯ (r) =

¯ (r) =

τL (r) + v(r) + vH (r) + vXC (r), ρ(r)

(19)

where all density-dependent quantities on the right-hand side are defined in terms of KS orbitals. The orbital-invariant definitions of the ALEE by Eqs. (18) and (19) were originally proposed by Bulat et al.18 who arrived at them by algebraic manipulation on the HF and KS equations (as did we in Ref. 24). Here, we went further and revealed that Eqs. (18) and (19) are two special cases of one definition given by Eq. (10).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-3

I. G. Ryabinkin and V. N. Staroverov

J. Chem. Phys. 141, 084107 (2014)

The pair function can be always written as25

III. AVERAGE LOCAL ELECTRON ENERGY FOR CORRELATED WAVEFUNCTIONS

1 ρ(r)[ρ(r2 ) + ρXC (r, r2 )], (27) 2 a relation which defines ρXC (r, r2 ), the exchange-correlation hole density. If we cast Eq. (25) in terms of ρXC , we obtain the following expression  τL (r) ρXC (r, r2 ) + v(r) + vH (r) + dr2 . (28) (r) ¯ = ρ(r) |r − r2 | P (r, r2 ) =

A. Definition

Our second objective is to generalize the definition of the ALEE of Eq. (10) to correlated wavefunctions. In Sec. II, we saw that the central ingredient in devising a properly invariant ALEE expression for one-determinantal SCF methods was the operator Fˆ γˆ . Let us write out the coordinate representation of this operator explicitly. Using Eqs. (6) and (12) we obtain   ˆ r|F γˆ |r  = r|Fˆ |r2 r2 |γˆ |r  dr2 1 = − ∇r2 γ (r, r ) + v(r)γ (r, r ) 2  γ (r, r )γ (r2 , r2 ) + dr2 |r − r2 |  1 γ (r, r2 )γ (r2 , r ) dr2 . − 2 |r − r2 |

a generalization of vXS (r) to include correlation effects. Recall that in the HF theory, HF ρXC (r, r2 ) = ρX (r, r2 ) = −

(20)

The last two terms of Eq. (20) can be combined into a single term to give 1 r|Fˆ γˆ |r  = − ∇r2 γ (r, r ) + v(r)γ (r, r ) 2  HF (r, r2 ; r , r2 ) dr2 , +2 |r − r2 |

(21)

where HF is the spinless two-electron reduced density matrix (2-RDM) of the HF theory,25   1 1 γ (r, r )γ (r2 , r2 )− γ (r, r2 )γ (r2 , r ) . HF (r, r2 ; r , r2 ) = 2 2 (22) We propose to generalize the one-determinantal ALEE by replacing the HF one- and two-electron RDMs in Eq. (21) with the general one- and two-electron RDMs arising in the context of any wavefunction method. Specifically, we define the generalized ALEE as ¯ (r) =

1 ˆ r|G|r, ρ(r)

(23)

ˆ is an integral operator with the kernel where G ˆ   = − 1 ∇r2 γ (r, r ) + v(r)γ (r, r ) r|G|r 2  (r, r2 ; r , r2 ) +2 dr2 , |r − r2 |

(24)

and γ and  are, respectively, the 1-RDM and 2-RDM corresponding to a given wavefunction. Explicitly, the proposed definition is  2 P (r, r2 ) τ (r) + v(r) + dr , (25) ¯ (r) = L ρ(r) ρ(r) |r − r2 | 2 where τL (r) is the Laplacian form of the kinetic-energy density extracted from the wavefunction and P (r, r2 ) is the pair function defined by P (r, r2 ) = (r, r2 ; r, r2 ).

Here the last term is the Slater averaged exchange-correlationcharge potential,27  ρXC (r, r2 ) S vXC (r) = (29) dr2 , |r − r2 |

(26)

|γ (r, r2 )|2 , 2ρ(r)

HF (r, r2 ) at r is just so the Coulomb potential of ρXC  HF ρXC (r, r2 ) dr2 = vXS (r). |r − r2 |

(30)

(31)

This means that Eq. (28) correctly reduces to Eq. (18) for onedeterminantal wavefunctions. The exchange-correlation energy corresponding to a given wavefunction can be always written as27  1 S EXC = (r) dr. (32) ρ(r)vXC 2 Therefore, the generalized Slater potential may be formally interpreted as the potential energy of a reference electron at r averaged over all other electrons. The use of the qualifier “formally” is essential here because the quantum force field corresponding to the Slater potential is not conservative.28 Now we are in position to justify the name ALEE for the function ¯ (r). Compare Eqs. (18), (19) and (28) giving the expressions for ¯ (r) in the HF, KS, and correlated-wavefunction theories, respectively. All these equations define ¯ (r) as the sum of the average local kinetic energy per electron and three potential-energy contributions. The only difference between these expressions is in the last term representing exchangecorrelation effects. The potential-energy part of (r) ¯ has been even used to define the “characteristic atomic radius” by the classical turning-point equation.29 Observe also that the total electronic energy of the system excluding the constant nuclear repulsion, Eelec , is related to ¯ (r) by  Eelec = ρ(r)¯ (r) dr − Vee , (33) where Vee =



P (r, r2 ) 1 dr2 = |r − r2 | 2



  S (r) dr ρ(r) vH (r) + vXC

(34) is the electron-electron interaction energy. The fact that there is a double-counting of electron-electron interaction in the generalized ¯ (r) reinforces our interpretation of this quantity as the total energy of an electron at r averaged over positions

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-4

I. G. Ryabinkin and V. N. Staroverov

J. Chem. Phys. 141, 084107 (2014)

of all other electrons in the system. It is significant that this interpretation is strictly correct for any wavefunction method in any basis set, whereas the association of −¯ (r) with electron removal energies is valid only in an approximate sense. There exists yet another way of writing the generalized ALEE which resembles the original definition of ¯ (r) in terms of canonical orbitals and their eigenvalues. This expression is discussed in Sec. III B. B. Calculation of the ALEE in terms of the generalized Fock operator

The quantity defined by Eq. (24) is the coordinate repreˆ known as the generalized sentation of an integral operator G (or extended) Fock operator or orbital Lagrangian. This operator arises naturally in multiconfiguration SCF methods,30–33 in the derivation of the extended Koopmans theorem,34–42 and in other contexts.43–45 In those methods where the total energy is stationary with respect to orbital variations, the operaˆ is Hermitian.34, 41 These methods include HF SCF, full tor G configuration interaction (FCI), and multiconfiguration SCF techniques with a restricted active space (RAS) or a complete ˆ is used in active space (CAS). In fact, the Hermiticity of G multiconfiguration SCF methods as a convergence criterion for orbital optimization.32 When is an orbital-optimized wavefunction, it is possible to set up the Hermitian eigenvalue problem ˆ i (r) = λi fi (r), Gf

(35)

where, according to Eq. (24),    1 2 ˆ − ∇r + v(r) γ (r, r )fi (r ) dr Gfi (r) ≡ 2   (r, r2 ; r , r2 ) + 2 dr2 fi (r ) dr . (36) |r − r2 | Assuming that the eigenfunctions fi (r), sometimes called the canonical extended HF orbitals,45 form a complete set and using Eq. (35), we have the following diagonal representation:   ˆ ˆ i fi |r = r|G|r = r|G|f λi |fi (r)|2 , (37) i

i

ˆ and where the summation is over all eigenfunctions of G each fi (r) is assumed to be normalized. Then the generalized ALEE can be computed by the equation 1  λ |f (r)|2 , (38) ¯ (r) = ρ(r) i i i which formally resembles Eq. (1). In finite-basis-set calculations, Eq. (35) can be solved in ˆ in an matrix form starting from the matrix representation of G 34 orthonormal basis of some spin-orbitals   Gpq = hpr γrq + 2 qrst prst, (39) r

rst

where hpq are elements of the one-electron core Hamiltonian matrix, prst are antisymmetrized two-electron integrals, while γ pq and  qrst are elements of the basis-set representations of the 1-RDM and 2-RDM. The number of

TABLE I. Canonical HF orbital eigenvalues ( i ) and eigenvalues of the generalized Fock operator (λi ) computed for a Be atom at various levels of theory using the 6-31G* basis set. λi , Eh Orbital 1s 2s 1p 3s 2p 1d 4s a

 i , Eh

HF

CAS SCFa

FCI

− 4.709455 − 0.301538 0.082241 0.438941 0.464361 1.068860 2.850434

− 4.709455 − 0.301538 0.000000 0.000000 0.000000 0.000000 0.000000

− 4.695461 − 0.312265 − 0.015563 − 0.000931 − 0.000001 − 0.000041 − 0.000020

− 4.698524 − 0.312642 − 0.015681 − 0.002488 − 0.000138 − 0.000050 − 0.000568

The active space includes 2 electrons and all orbitals except 1s.

eigenfunction-eigenvalue pairs of the G matrix is finite and equal to the number of linearly independent one-electron basis functions. ˆ = Fˆ γˆ , the eigenfuncAt the HF level of theory, where G ˆ tions of G are the canonical HF orbitals (fi = φ i ), and the eigenvalues are λi = ni  i , where ni is the occupation number (1 or 0) of the canonical HF orbital φ i with the eigenvalue  i . Thus, Eq. (38) correctly reduces to Eq. (1) for the HF method. For a correlated wavefunction represented by a linear combination of determinants, configuration mixing populates virtual orbitals and depopulates the occupied ones. As a result, all eigenvalues λi in post-HF methods are generally nonzero (Table I). Note that although Eq. (38) is similar to Eq. (1), its right-hand side is not a weighted average of λi because  ρ(r) = i |fi (r)|2 even for HF wavefunctions. It has been repeatedly pointed out34, 39, 45 that the eigenvalue problem for the generalized Fock operator is closely related to the generalized eigenvalue equation expressing the extended Koopmans theorem,34–42  ˆ i (r) = ωi Gu

γ (r, r )ui (r ) dr ,

(40)

where −ωi are estimates of the ionization energies of the system and ui (r) are the corresponding amplitudes. Because the right-hand sides of Eqs. (35) and (40) are different, the eigenfunctions fi (r) and ui (r) and their eigenvalues have different properties.34, 39, 45 In particular, individual λi values do not have a physical meaning.45 But if we substitute Eq. (38) into Eq. (33) and use the fact that each fi (r) is normalized, we obtain an exact HF-like energy formula Eelec =



λi − Vee ,

(41)

i

which reveals the physical significance of the sum of λi . In those methods where the orbitals used to build the wavefunction are not variationally optimal (e.g., truncated configuration interaction, truncated coupled-cluster expanˆ is not sions, and standard perturbation theory), the operator G ˆ Hermitian. In such cases, one can compute r|G|r directly from its nondiagonal matrix representation G and values of the basis functions. Another possibility is to use Eq. (28).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-5

I. G. Ryabinkin and V. N. Staroverov

J. Chem. Phys. 141, 084107 (2014)

C. Asymptotic behavior of the ALEE

−0.90

lim ¯ (r) = HOMO ,

r→∞

(42)

except on the nodal surface of the HOMO. The property of the KS ALEE expressed by Eq. (42) was stated earlier by Ayers et al.17 Furthermore, the magnitude of the HOMO eigenvalue of the exact KS scheme is equal to the exact first vertical ionization energy,50–52 Imin , and for the exact functional we would have KS HOMO = −Imin .

(43)

Equation (43) is of course violated in approximate densityfunctional theory, but can be satisfied through fractional HOMO depopulation53, 54 or by other methods. The asymptotic limit of the ALEE from an exact manyelectron wavefunction can be established by analyzing the r → ∞ limit of each term on the right-hand side of Eq. (28). The external potential and the Coulomb potentials of ρ(r) and the exchange-correlation hole density all vanish at infinity, so it remains to examine the kinetic-energy contribution, τ L /ρ. It is known37, 55, 56 that all natural spin-orbitals have the same long-range exponential fall-off factor, e−cr with c = 2Imin . Using this fact and the expressions for ρ and τ L in terms of natural spin-orbitals we obtain through straightforward differentiation τL (r) = −Imin . r→∞ ρ(r)

lim ¯ (r) = lim

r→∞

(44)

We should add that in practical calculations with finite Gaussian basis sets, the apparent r → ∞ limit of ¯ (r) may depend on how this quantity is computed. If the ALEE is computed by Eq. (28), the limit will not exist because the ratio τ L /ρ for any orbital with the leading asymptotic term 2 φ(r) ∼ r β e−αr diverges at large r as −2α 2 r2 . At the same time, (38), where the leading asymptotic terms  if one uses Eq. 2 λ |f (r)| and ρ(r) are determined by the most difof M i=1 i i fuse basis function, then ¯ (r) will remain finite and approach a constant.

− ε(r) [Eh]

−0.95 1

He ( S) KS (exact)

−1.00

HF (basis−set limit) FCI / cc−pVQZ FCI / cc−pV6Z

−1.05 0

1

2

3 r [a0]

4

5

6

FIG. 1. ALEE in a He atom computed at various levels of theory. The KS KS and HF ALEEs are constants equal to HOMO = −Imin = −0.9037 Eh and HF HOMO = −0.9180 Eh , respectively.

IV. NUMERICAL EXAMPLES

We implemented Eq. (38) in the GAMESS(US) code57–59 and used it to compute the generalized ALEE for a number of atoms and molecules. In all calculations, we used Dunning’s correlation-consistent basis sets60 taken from the Environmental Molecular Sciences Laboratory (EMSL) Basis Set Library.61, 62 Correlated wavefunctions were obtained by the FCI, CAS SCF, and RAS SCF methods. The generalized Fock matrices constructed from these wavefunctions had a residual antisymmetry with max|Gpq − Gqp | ≤ 10−7 (10−6 for the CO molecule) and were symmetrized before diagonalization. The exact atomic ionization energies were taken from Refs. 63 and 64. Consider first the ground-state He atom. At the HF and KS levels of theory, this system has only one occupied orbital (1s), so its ALEE is a constant, ¯ (r) = HOMO = 1s .

(45)

The ALEEs calculated from two finite-basis-set FCI wavefunctions, however, are position-dependent (Fig. 1). Each FCI 0.0 −Imin −1.0

− ε(r) [Eh]

One of the characteristic properties of the ALEE is that it generally does not vanish at infinity but approaches a systemand method-dependent constant. The value of this constant can be deduced from the definitions. In the KS scheme, the radial dependence of each atomic or molecular orbital at large r is φi (r) ∼ r βi e−αi r , where αi = −2i and β i is some constant.46 In the HF theory, canonical orbitals decay in a similar way for atoms with only stype orbitals. In general, the asymptotically leading term of an HF orbital is φi (r) ∼ r βi e−αr , where α = −2HOMO is the same for all orbitals, and β i is largest for the highest occupied molecular orbital (HOMO).47–49 Although the large-r behavior of KS and HF orbitals is governed by different laws, the total electron density in both methods is asymptotically dominated by |φHOMO (r)|2 , and so from Eq. (1) we see that in either case

−2.0 Be (1S)

−3.0

HF / cc−pCVQZ −4.0

FCI / cc−pCVQZ

−5.0 0

1

2

3

r [a0]

FIG. 2. ALEE in a Be atom computed by the HF and FCI methods using the cc-pCVQZ basis set. The HF curve approaches the limit  2s = −0.3093 Eh ; the FCI curve approaches a limit estimated at −0.3482 Eh . The exact limit, −Imin = −0.3426 Eh , is shown by the dotted line.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-6

I. G. Ryabinkin and V. N. Staroverov

J. Chem. Phys. 141, 084107 (2014)

−Imin

−5

− ε(z) [Eh]

− ε(r) [Eh]



−0.50

−1 1

Ne (1S)

−10 10

1 +

−0.65

H2 ( Σg) FCI / cc−pV5Z

1

2

3



1.5Re

−30 0



2Re



−0.60

HF / cc−pVQZ RAS SCF / cc−pVQZ

−20



−0.55



3Re

4

−10

5

−5

r [a0]

Re

• • 0 z [a0]

5

10

FIG. 4. ALEE for H2 computed from FCI/cc-pV5Z wavefunctions for various internuclear distances. The equilibrium H–H distance is Re = 1.4011a0 . Each curve shows ¯ (r) along the internuclear axis.

FIG. 3. ALEE for a Ne atom computed by the HF and RAS SCF methods using the cc-pVQZ basis set. Note the logarithmic scale of ¯ (r). The RAS SCF wavefunction includes all configurations generated by single and double excitations from the HF reference determinant. The HF curve approaches the HF limit HOMO = −0.8490 Eh ; the RAS curve approaches −0.8093 Eh . The exact limit, −Imin = −0.7945 Eh , is shown by the dotted line.

be tuned by changing the internuclear distance, R. As a result, the ALEE functions computed for H2 from FCI wavefunctions also change with R. At R = Re , the FCI ALEE in H2 near each nucleus resembles the FCI ALEE in the He atom (Fig. 1). As the bond is stretched, the plots of ¯ (r) approaches a straight horizontal line at −0.5 Eh . This happens because in the R → ∞ limit the H2 molecule is a pair of non-interacting hydrogen atoms, and the ALEE for an H atom is precisely (r) ¯ = −0.5Eh . Consistent with Eq. (44), the values of −¯ (z) at |z| = 10a0 in Fig. 4 closely match the vertical ionization energies of H2 computed at the FCI/cc-pV5Z level of theory for the corresponding R. The plots of ¯ (r) for the HF wavefunction are not shown in Fig. 4 to avoid clutter. For R = Re , 1.5Re , 2Re , and 3Re they would be just straight horHF = −0.5945, −0.5029, −0.4430, and izontal lines at HOMO −0.3722 Eh , respectively. As R increases, the HF ALEE deviates from the exact value (−0.5 Eh ) more and more, reflecting the increasing inability of a single Slater determinant to describe static correlation. Our last example involves a larger molecule, CO, and illustrates the changes in the ALEE accompanying qualitative improvement of the wavefunction (Fig. 5). The three

ALEE curve starts at about −1.05 Eh at the nucleus and increases monotonically, approaching a limiting value which is very close to the exact −Imin . The ALEE obtained from multideterminantal wavefunctions is not constant because the electrons in He are correlated: when the reference electron is near the nucleus, the second electron is pushed away, so the nuclear attraction felt by the reference electron becomes less shielded than in the HF picture. This deshielding lowers the electron’s local energy. For a Be atom, the ALEEs computed from the HF and FCI wavefunctions are both position-dependent (Fig. 2) and exhibit the characteristic atomic shell structure reported by Politzer et al.3 Correlation effects are small on the scale of this figure but still visible: electron correlation lowers ¯ (r) near HF the nucleus and at large r. Note that for the Be atom HOMO > −Imin . HF and correlated-wavefunction ALEEs for a Ne atom (Fig. 3) are qualitatively similar to those for Be, except HF < −Imin . that for Ne we have HOMO The H2 molecule is another interesting example (Fig. 4). Here, the amount of static electron correlation can

−0.5 CAS SCF

HF

(CAS+SD)SCF −0.6 −0.7

• C

• O

• C

• O

• C

• O

−0.8 −0.9 −1.0 −1.1 −1.2

FIG. 5. Color-coded maps of ¯ (r) computed for the ground-state CO molecule (Re = 1.1283 Å) from various wavefunctions using the cc-pVTZ basis set. The maps are for a planar elliptic coordinate grid with foci at the C and O nuclei. The complete specification of the CAS SCF and (CAS+SD)SCF wavefunctions is given in text. The dimensions of each box are 4 Å × 4 Å. The color-coding spans only the narrow, chemically relevant range of ALEE values, shown in units of Eh next to the color bar.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-7

I. G. Ryabinkin and V. N. Staroverov

wavefunctions used here are, in the order of increasing accuracy: HF, full valence CAS SCF obtained by distributing 10 active electrons over 8 valence orbitals, and (CAS+SD)SCF, which includes all full-valence CAS configurations plus all single and double excitations from this CAS to recover dynamic electron correlation. All orbitals in these three wavefunctions are variationally optimized. Figure 5 shows that the transition from HF to CAS SCF (inclusion of static correlation) lowers the ALEE mostly around the O atom, while the transition from CAS SCF to (CAS+SD)SCF (inclusion of dynamic correlation) lowers the ALEE mostly in the bonding region. This trend is in line with our earlier observation that electron correlation affects the ALEE primarily in the highdensity regions. It also reflects a reorganization of the electronic structure of the CO molecule. As is known,65 the actual dipole moment of the CO molecule is consistent with the formula Cδ − Oδ + , whereas the HF method predicts the reverse polarity, Cδ + Oδ − . The ALEE plots computed at the CAS SCF and (CAS+SD)SCF levels indicate that the O-end of the CO molecule is less prone to electrophilic attack than predicted by the HF method, in agreement with the actual polarity of the molecule. V. CONCLUSION

We showed that Politzer’s average local ionization energy index, I¯(r) = −¯ (r), can be extended to correlatedwavefunction methods. The proposed generalized definition of (r) ¯ is given by Eq. (28). This equation suggests a new interpretation of ¯ (r) as the ALEE—the total energy of an electron at r averaged over positions of all other electrons in the system. ALEEs extracted from one- and multideterminantal wavefunctions may differ qualitatively because electron correlation can induce chemically significant changes in ¯ (r) relative to the HF theory. In particular, the HF and KS ALEEs for He and H2 are just constants, whereas correlatedwavefunction ALEEs for the same systems have significant spatial variations. The ALEEs obtained from an exact wavefunction and from exact KS orbitals are distinct but have the same r → ∞ limit equal to the negative of the exact first ionization energy of the system. The r → ∞ limit of an ALEE computed from a HF wavefunction is the negative of the HF HOMO energy. It should be noted that the original definition of the function I¯(r) can be extended to correlated wavefunctions in more than one way. Forinstance, one could define the generalized index I¯(r) = n In |gn (r)|2 /ρ(r), where gn (r) are Dyson orbitals66 and In are the corresponding ionization energies or by a similar expression in terms of the amplitudes and eigenvalues of the extended Koopmans theorem. Nevertheless, once the HF ¯ (r) of Eq. (1) is cast as Eq. (18) and reinterpreted as an ALEE, the transition to Eq. (28) is arguably more natural and practicable than other possibilities. Our generalization of the ALEE concept paves the way for systematic studies and practical applications of this quantity. It would be especially interesting to explore the generalized ALEE as a descriptor of chemical reactivity at various post-HF levels of theory. In addition, because the ALEE

J. Chem. Phys. 141, 084107 (2014)

expressions given by Eqs. (28) and (38) represent enticingly simple reduced forms of the many-electron Schrödinger equation, we expect these expressions to be useful for studying exact effective potentials of the Kohn–Sham density-functional theory. ACKNOWLEDGMENTS

I.G.R. thanks Dr. Alex Granovsky for help with the GAMESS code and Sviataslau Kohut for supplying auxiliary data. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grants Program. 1 P.

Sjoberg, J. S. Murray, T. Brinck, and P. Politzer, Can. J. Chem. 68, 1440 (1990). 2 J. S. Murray, J. M. Seminario, P. Politzer, and P. Sjoberg, Int. J. Quantum Chem., Symp. 24, 645 (1990). 3 P. Politzer, J. S. Murray, M. E. Grice, T. Brinck, and S. Ranganathan, J. Chem. Phys. 95, 6699 (1991). 4 J. S. Murray, T. Brinck, and P. Politzer, J. Mol. Struct.: THEOCHEM 255, 271 (1992). 5 P. Politzer, F. Abu-Awwad, and J. S. Murray, Int. J. Quantum Chem. 69, 607 (1998). 6 J. S. Murray, and P. Politzer, in Theoretical Organic Chemistry, edited by C. Párkányi (Elsevier, Amsterdam, 1996), pp. 189–202. 7 P. Politzer, and J. S. Murray, in Theoretical Aspects of Chemical Reactivity, edited by A. Toro-Labbé (Elsevier, Amsterdam, 2007), pp. 119–137. 8 P. Politzer, J. S. Murray, and F. A. Bulat, J. Mol. Model. 16, 1731 (2010). 9 A. Toro-Labbé, P. Jaque, J. S. Murray, and P. Politzer, Chem. Phys. Lett. 407, 143 (2005). 10 F. A. Bulat, A. Toro-Labbé, T. Brinck, J. S. Murray, and P. Politzer, J. Mol. Model. 16, 1679 (2010). 11 J. S. Murray, Z. Peralta-Inga Shields, P. Lane, L. Macaveiu, and F. A. Bulat, J. Mol. Model. 19, 2825 (2013). 12 P. Politzer, J. S. Murray, and M. E. Grice, Collect. Czech. Chem. Commun. 70, 550 (2005). 13 P. Politzer, Z. Peralta-Inga Schields, F. A. Bulat, and J. S. Murray, J. Chem. Theory Comput. 7, 377 (2011). 14 P. Jin, J. S. Murray, and P. Politzer, Int. J. Quantum Chem. 96, 394 (2004). 15 A. Nagy, R. G. Parr, and S. Liu, Phys. Rev. A 53, 3117 (1996). 16 T. Gál and Á. Nagy, Mol. Phys. 91, 873 (1997). 17 P. W. Ayers, R. G. Parr, and A. Nagy, Int. J. Quantum Chem. 90, 309 (2002). 18 F. A. Bulat, M. Levy, and P. Politzer, J. Phys. Chem. A 113, 1384 (2009). 19 J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 45, 101 (1992). 20 O. Gritsenko, R. van Leeuwen, E. van Lenthe, and E. J. Baerends, Phys. Rev. A 51, 1944 (1995). 21 A. D. Becke and E. R. Johnson, J. Chem. Phys. 124, 221101 (2006). 22 J. C. Slater, Phys. Rev. 81, 385 (1951). 23 I. G. Ryabinkin, A. A. Kananenka, and V. N. Staroverov, Phys. Rev. Lett. 111, 013001 (2013). 24 S. V. Kohut, I. G. Ryabinkin, and V. N. Staroverov, J. Chem. Phys. 140, 18A535 (2014). 25 R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, New York, 1989), pp. 20–46. 26 R. McWeeny, Rev. Mod. Phys. 32, 335 (1960). 27 J. C. Slater, Phys. Rev. 91, 528 (1953). 28 V. Sahni, Quantal Density Functional Theory (Springer, Berlin, 2004), pp. 205–213. 29 Z.-Z. Yang and E. R. Davidson, Int. J. Quantum Chem. 62, 47 (1997). 30 J. Hinze, J. Chem. Phys. 59, 6424 (1973). 31 T. L. Gilbert, J. Chem. Phys. 60, 3835 (1974). 32 P. E. M. Siegbahn, J. Almlöf, A. Heiberg, and B. O. Roos, J. Chem. Phys. 74, 2384 (1981). 33 R. H. A. Eade and M. A. Robb, Chem. Phys. Lett. 83, 362 (1981). 34 O. W. Day, D. W. Smith, and C. Garrod, Int. J. Quantum Chem., Symp. 8, 501 (1974). 35 D. W. Smith and O. W. Day, J. Chem. Phys. 62, 113 (1975). 36 O. W. Day, D. W. Smith, and R. C. Morrison, J. Chem. Phys. 62, 115 (1975).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

084107-8 37 M.

I. G. Ryabinkin and V. N. Staroverov

M. Morrell, R. G. Parr, and M. Levy, J. Chem. Phys. 62, 549 (1975). C. Ellenbogen, O. W. Day, D. W. Smith, and R. C. Morrison, J. Chem. Phys. 66, 4795 (1977). 39 J. M. O. Matos and O. W. Day, Int. J. Quantum Chem. 31, 871 (1987). 40 R. C. Morrison and G. Liu, J. Comput. Chem. 13, 1004 (1992). 41 J. Cioslowski, P. Piskorz, and G. Liu, J. Chem. Phys. 107, 6804 (1997). 42 D. Vanfleteren, D. V. Neck, P. W. Ayers, R. C. Morrison, and P. Bultinck, J. Chem. Phys. 130, 194104 (2009). 43 P.-O. Löwdin, Phys. Rev. 97, 1474 (1955). 44 D. H. Kobe, Phys. Rev. 188, 1583 (1969). 45 L. J. Holleboom, J. G. Snijders, and E. J. Baerends, Int. J. Quantum Chem. 34, 289 (1988). 46 E. Engel and R. M. Dreizler, Density Functional Theory: An Advanced Course (Springer, Berlin, 2011), pp. 84–88. 47 N. C. Handy, M. T. Marron, and H. J. Silverstone, Phys. Rev. 180, 45 (1969). 48 G. S. Handler, D. W. Smith, and H. J. Silverstone, J. Chem. Phys. 73, 3936 (1980). 49 T. Ishida and K. Ohno, Theor. Chim. Acta 81, 355 (1992). 50 C.-O. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985). 51 J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Jr., Phys. Rev. Lett. 49, 1691 (1982). 52 M. Levy, J. P. Perdew, and V. Sahni, Phys. Rev. A 30, 2745 (1984). 53 A. P. Gaiduk, D. S. Firaha, and V. N. Staroverov, Phys. Rev. Lett. 108, 253005 (2012). 38 J.

J. Chem. Phys. 141, 084107 (2014) 54 A.

P. Gaiduk, D. Mizzi, and V. N. Staroverov, Phys. Rev. A 86, 052518 (2012). 55 M. Levy and R. G. Parr, J. Chem. Phys. 64, 2707 (1976). 56 J. Katriel and E. R. Davidson, Proc. Natl. Acad. Sci. U.S.A. 77, 4403 (1980). 57 M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, Jr., J. Comput. Chem. 14, 1347 (1993). 58 M. S. Gordon, and M. W. Schmidt, in Theory and Applications of Computational Chemistry: The First Forty Years, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Elsevier, Amsterdam, 2005), pp. 1167–1189. 59 The General Atomic and Molecular Electronic Structure System (GAMESS), Version of May 2013, see http://www.msg.chem.iastate.edu/ GAMESS/GAMESS.html. 60 T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 61 D. Feller, J. Comput. Chem. 17, 1571 (1996). 62 K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, and T. L. Windus, J. Chem. Inf. Model. 47, 1045 (2007). 63 E. R. Davidson, S. A. Hagstrom, S. J. Chakravorty, V. Meiser Umar, and C. Froese Fischer, Phys. Rev. A 44, 7071 (1991). 64 S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A. Parpia, and C. Froese Fischer, Phys. Rev. A 47, 3649 (1993). 65 F. Grimaldi, A. Lecourt, and C. Moser, Int. J. Quantum Chem. 1S, 153 (1967). 66 B. T. Pickup, Chem. Phys. 19, 193 (1977).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.26.31.3 On: Mon, 22 Dec 2014 17:59:41

Average local ionization energy generalized to correlated wavefunctions.

The average local ionization energy function introduced by Politzer and co-workers [Can. J. Chem. 68, 1440 (1990)] as a descriptor of chemical reactiv...
2MB Sizes 4 Downloads 10 Views