THE JOURNAL OF CHEMICAL PHYSICS 142, 204311 (2015)
Real-space representation of electron correlation in π -conjugated systems Jian Wang1,a) and Evert Jan Baerends2,a) 1
School of Science, Huzhou University, Zhejiang 10083, China Afdeling Theoretische Chemie, FEW, Vrije Universiteit, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands 2
(Received 5 March 2015; accepted 6 May 2015; published online 29 May 2015) π-electron conjugation and aromaticity are commonly associated with delocalization and especially high mobility of the π electrons. We investigate if also the electron correlation (pair density) exhibits signatures of the special electronic structure of conjugated systems. To that end the shape and extent of the pair density and derived quantities (exchange-correlation hole, Coulomb hole, and conditional density) are investigated for the prototype systems ethylene, hexatriene, and benzene. The answer is that the effects of π electron conjugation are hardly discernible in the real space representations of the electron correlation. We find the xc hole to be as localized (confined to atomic or diatomic regions) in conjugated systems as in small molecules. This result is relevant for density functional theory (DFT). The potential of the electron exchange-correlation hole is the largest part of v xc , the exchange-correlation Kohn-Sham potential. So the extent of the hole directly affects the orbital energies of both occupied and unoccupied Kohn-Sham orbitals and therefore has direct relevance for the excitation spectrum as calculated with time-dependent DFT calculations. The potential of the localized xc hole is comparatively more attractive than the actual hole left behind by an electron excited from a delocalized molecular orbital of a conjugated system. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921725] I. INTRODUCTION
π-conjugated systems are characterized by alternating single (σ) and double (σ and π) bonds. Benzene is a prototype for the cyclic conjugated structures which are known for their special properties (cf. the heme group, chlorophyll, and graphene). The associated concept of aromaticity has received much interest.1,2 The π electrons are commonly conceived as particularly mobile, and for instance when exposed to magnetic fields, a large ring current can be induced around the benzene ring.3–5 Kekulé6 in the 19th century conceived the cyclic structure with alternating single bond and double bonds for benzene. Depending on whether a given bond is a single bond or double bond, there are two possibilities for the Kekulé structure. To comply with the experimental evidence that bond lengths in benzene are the same, Kekulé suggested the two alternative structures are in “oscillation” or “resonating.” This notion found a natural explanation in quantum mechanics, which allows enhanced stability by superposition of (trial) wave functions. Pauling and Wheland7,8 found that the superposition of the Kekulé electronic structures and other possible electronic structures, such as Dewar structures, results in energetic stabilization. This “resonance” effect in benzene had been previously proposed by Slater.9 It provides an energetic mechanism of enhanced stability for the symmetric structure. In the early molecular orbital based calculations,10–12 benzene was modelled by separating the electrons into orthogonal σ and π subspaces. The σ electrons occupy a closedshell configuration and make up the molecular σ core, while the π electrons move in orbitals orthogonal to this core. The a)Electronic addresses: [email protected] and [email protected]
0021-9606/2015/142(20)/204311/12/$30.00
molecular orbitals are constructed according to the molecular symmetry, which is D6h for benzene. The electronic structure from the molecular orbital picture provides a roughly correct explanation for the experimental spectra. In conjugated chains, the experimental spectra justify the delocalized picture of the π molecular orbitals.13–15 There have been conflicting views on the delocalized nature of the π electron subsystem in aromatic compounds. Ab initio self-consistent field (SCF) molecular orbital calculations16,17 and various high level configuration interaction (CI) calculations18–23 provide, as expected, detailed comparison with experiments. However, Cooper, Gerratt, and Raimondi24 have used valence bond models to argue that the π electrons can be pictured as residing in almost localized orbitals on the individual carbon atoms of benzene, defying the text-book picture of its delocalized π ring structure. In their calculation, they started from initially delocalized symmetric molecular orbitals for the π-electrons but lifted symmetry restrictions during subsequent optimization. Shaik, Hiberty, and co-workers25,26 have argued that the π-electrons of benzene tend to pair with each other which will distort the equal sided hexagon, and the molecule maintains 6-fold symmetry only because of the σsystem (cf. also the Kohn-Sham MO analysis of Pierrefixe and Bickelhaupt27). Of course, orbitals are not observables, and for instance, the Hartree-Fock wave function is invariant under arbitrary unitary transformation of the orbitals.28 Both molecular orbital (Hartree-Fock plus CI) and valence bond type representations of the wave function may lead to different pictures of the localized versus delocalized nature of the bonding in conjugated compounds. Similarly in solids, one may choose the delocalised Bloch function or the localised Wannier functions to describe the electronic structure.29–31 In contrast,
142, 204311-1
© 2015 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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-2
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
the pair density and other quantities derived from the twoparticle density matrix are defined, unambiguous, real-space quantities. It will be interesting to see whether conjugation, delocalization, and aromaticity can be traced in the description of electron correlation with these spatial probabilities. We focus on the pair density and derived quantities, such as Fermi hole, exchange-correlation (XC) hole, Coulomb hole, and conditional density. We note that such a study has immediate relevance for the much used Kohn-Sham model of density functional theory (DFT). Often the XC energy of DFT Exc (which contains a kinetic energy contribution) is expressed through the coupling constant integrated XC hole, 1 hole d 3r 1 n(1) d 3r 2 n¯ xc (2|1)/r 12. (1) Exc[n] = 2 hole The coupling constant integrated XC hole n¯ xc (2|1) in DFT contributes not only to the electron-electron Coulomb energy but also takes care of the kinetic energy difference between the exact kinetic energy and that modeled by the Kohn-Sham orbitals, T − Ts . The kinetic energy contribution to the correlahole tion is folded into the XC hole n¯ xc by the so-called coupling32–34 constant integration, 1 hole n¯ xc (2|1) = dλ nhole (2) xc,λ(2|1), 0
(2|1) is the hole for a system in which the Coulomb where nhole xc,λ interaction between electrons is given by λ/r 12, but the external
potential is adapted in such a way that the density remains equal to the exact one. The coupling constants with 0 ≤ λ ≤ 1 connect the Kohn-Sham system (λ = 0) to the physical system (λ = 1). It is very difficult to construct the XC hole with the coupling constant integration.35–37 However, a large part of the energy density for Exc and of the Kohn-Sham potential vxc(r) can be obtained directly from the full hole nhole (2|1),38–40 xc,λ=1 1 hole v + vc,kin, 2 xc hole vxc = vxc + vc,kin + vresp, hole vxc = d 3r 2 nhole xc,λ=1(2|1)/r 12, ϵ xc =
(3)
where ϵ xc is the XC energy density per particle, from which the XC energy may be obtained by Exc[n] = n(1)ϵ xc([n]; 1)d1. Here, vc,kin is the energy density that accounts for the correlation correction in kinetic energy, T − Ts = n(1)vc,kin(1)d1 and vresp is the so-called response part of vxc. The hole potential hole vxc is usually the largest part of vxc. Notably in situations of bond breaking, the components vc,kin and vresp may also become significant,41–43 and this is in particular the case when timedevelopment is studied with time-dependant DFT (TDDFT) (beyond the linear response regime).44 However, here we will exclusively deal with the static ground state and equilibrium hole geometries. The more or less attractive nature of the vxc potential, depending on the localized or delocalized nature of the hole, affects directly the orbital shapes and the orbital energies of the Kohn-Sham orbitals. If there is a reasonable correspondence between the size of the exchange-correlation hole and the size of the electron hole that is created in an excitation, the orbital energy difference between the occupied and virtual orbitals involved in the excitation will be a very
good first approximation to the energy of the excitation.45 On the other hand, if there is a large mismatch, as is the case in solids and as may be the case in π-conjugated systems if the exchange-correlation hole is less delocalized than the π orbitals involved in an excitation, the orbital energy difference may not be a good first approximation to the excitation energy, and the time-dependent DFT calculations may be expected to be problematic. The paper is organized as follows. First, we outline the computational method for the pair density and its derived quantities. Then in Secs. III A and III B, the holes in benzene are analyzed. Next, in Sec. III C, the holes are compared to those of related systems such as hexatriene and ethene, and of cyclohexane, and in Sec. III D to holes from electron gas models. Finally, the conclusions are drawn.
II. COMPUTATIONAL METHOD
The pair density is defined as46,47 Γ(1, 2) = N(N − 1) × Ψ∗(123 · · · N)Ψ(123 · · · N)d3 · · · dN,
(4)
where N is the total number of electrons. The normalization follows McWeeny’s convention,47 which is a factor 1/2 different from Löwdin’s.46,48 The quantity can be interpreted as the probability that two electrons be simultaneously found at positions r 1 and r 2 with spins s1, and s2, respectively. Each electron has three spatial coordinates and one spin coordinate. The integral implies integration over space coordinates and summation over spin coordinates. If the spin coordinates are summed out, the pair density is then a function of the spatial coordinates of the two electrons, Γ(r1, r2). To visualize Γ(r1, r2) in three dimensions, we may fix the position of one electron, say r1 (the reference electron), and then look at the pair density (or the exchange-correlation hole) as a function of the position (r2) of the second electron. The pair density has previously also been studied through the reduced quantity: intracule.49–51 The XC hole surrounding an electron at 1 is defined as hole (2|1) = nxc
Γ(1, 2) − n(1)n(2) , n(1)
where n(1) is the density defined as n(1) = N Ψ∗(12 · · · N)Ψ(12 · · · N)d2 · · · dN.
(5)
(6)
The effects of exchange and correlation on the electron distribution are described by the XC hole. For the model system of the electron gas, the hole has been extensively studied,35,52–58 and the results have been used to construct XC functionals in DFT.59 The XC holes in a few prototype systems, such as H2,41,60–63 He,64,65 atomic Si,66 and crystal Si,67–69 have been carefully studied. Here, we study the XC hole for π-conjugated systems, previously the Coulomb hole around the bond midpoint has been reported.70 The density is related to the pair density as 1 Γ(1, 2)d2. (7) n(1) = N −1
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-3
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
From this, one can derive the important sum rule for the XC hole, hole nxc (2|1)d2 = −1. (8) A closely related quantity is the conditional density defined as the probability to find an electron at 2 when one is known to be at 1, hole ncond(2|1) = n(2) + nxc (2|1).
(9)
It integrates to N − 1 electrons because the XC hole integrates to −1, and the density to N. The spherically averaged conditional density in the Hartree-Fock approximation can be expanded as71 1 Dσ s2 + · · ·, (10) 3 where the arguments (r, s) denote the spherical average on a shell of radius s around the reference point r, and σσ (r, s) = ncond
1 (∇ρσ )2 , Dσ = τσ − 4 ρσ
(11)
τσ = Σi |∇ψi,σ |2 is the kinetic energy of electrons with spin σ, 2 and 14 (∇ρρ σσ ) is the von Weizsäcker kinetic energy functional. A smaller Dσ implies smaller probability to find a second electron around the reference electron, which may be interpreted as a stronger electron localisation. On this basis, the quantity Dσ has been used by Becke to propose an “electron localization function” (ELF).71 A possible connection between ELF and the phenomenon of aromaticity has been explored by several authors.72–74 The connection is however not very direct, and we will focus on the quantities directly related to the pair density, i.e., the exchange and correlation holes. Another quantity of interest is the one-particle density matrix defined as ′ γ(1 , 1) = N Ψ∗(1′2 · · · N)Ψ(12 · · · N)d2 · · · dN. (12) Its eigenfunctions are the natural orbitals, which were believed to have certain optimal properties for the description of electron correlation (shortest expansion)46 although this does not appear to be the case for general systems,75 γ (1′, 1) = λi χ∗i (1′) χ i (1), (13) i
where λi are the occupation numbers of the natural orbitals χ i . The density is just the diagonal of γ, n(1) = γ(1, 1). In the method of CI, one usually selects a set of M spin orbitals (M > N). One can then build various determinants (or configurations), Φi , each from N of those orbitals. The variational wave function can be written as a linear combination of those configurations, Ψ= Ci Φi . (14) i
The simplest case is M = N spin orbitals (Hartree-Fock), the wave function in this case is a single determinant. The pair density in this case has been worked out by Dirac76,77 ΓHF(1, 2) = γHF(1, 1)γHF(2, 2) − γHF(1, 2)γHF(2, 1),
(15)
where γHF (1, 1′) =
N
ϕi (1)ϕ∗i (1′).
(16)
i
Here, ϕi are the Hartree-Fock spin orbitals. The XC hole in this case reflects the effect of electron exchange, it is by definition the exchange hole, also called the Fermi hole, nhole x (2|1) = −
|γHF1, 2|2 . nHF (1)
(17)
The negative sign implies the reduction of probability to find the second electron at position 2 while a reference electron is at 1. The Fermi hole also integrates to −1, i.e., obeys the same sum rule as the total XC hole (see Eq. (8)). The Coulomb correlation hole can be defined as the difference between the XC hole and Fermi hole,60 hole hole nhole c (2|1) = nxc (2|1) − n x (2|1) .
(18)
Since both the XC hole and Fermi hole integrate to −1, the Coulomb hole integrates to zero. It only represents a redistribution in space of the electron density, i.e., reduced density in one region will show up in other regions.
III. CALCULATION RESULTS
We use the GAMESS package78 to calculate the wave function. The CI wave function is obtained using the multiconfiguration self-consistent field (MCSCF) or complete active space self-consistent field (CASSCF) method. The wave function is calculated in the graphical unitary group approach;19 as a byproduct, the elements of the two-matrix can be exported. Dunning’s correlation-consistent basis set cc-pVDZ79,80 is chosen. In terms of contracted functions, the basis set is 2s1p for H and 3s2p1d for C. A. Fermi hole in benzene
The optimal benzene structure has a C–C bond length 1.397 Å and C–H bond length of 1.084 Å. The Hartree-Fock energy is −230.722 a.u., as compared to Ermler and Kern’s −230.749 a.u.16 using an extended basis set (2s1p for H and 4s2p1d for C). Fig. 1 shows the Fermi hole where the reference electron is chosen at 1 bohr directly above a carbon atom. The atomic 2pz orbital has its maximum just below 1 bohr. The upper panel is a picture of the hole in three dimensions, from which it is clear the hole is quite localized around the reference electron, with only small contributions below the plane and at other atomic sites (ortho, para). The lower panel of Fig. 1 gives more quantitative information of the Fermi hole in the plane passing through the reference position and parallel to the molecular plane. The Fermi hole in benzene and other π-conjugated systems has been previously studied by Bader et al.81 They emphasized the feature at the para site as a signature of π delocalization. But in fact the amplitude at the para-site is very small (the magnitude there is in the order of 10−3 e/(bohrs)3). The surprise is rather that the Fermi hole is actually not delocalized.
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-4
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
FIG. 1. Benzene: The Fermi hole from the Hartree-Fock wave function. The reference electron is at 1 bohr above a C atom. The upper panel shows the 3-dimensional distribution of the hole, while the lower panel gives the cross section in the plane through the reference electron and parallel to the molecular plane.
We will compare to a non-conjugated system (cyclohexane, see Sec. III C) and note that the Fermi hole is not more delocalized in the π-conjugated system than in the non-conjugated system. Physically, a localized Fermi hole is favorable because it maximizes the attractive electron-hole interaction. Compared to the spherical hole in the uniform electron gas, the Fermi hole in benzene has a remarkable structure. It will be interesting to identify how its features are related to the π and σ bondings. For the Hartree-Fock approximation, it is straightforward to decompose the π and σ contributions since there are distinct σ and π orbitals. We may write the oneparticle density matrix as a sum of π and σ orbital contributions, γHF (1′, 1) = ϕ∗i,σ (1′)ϕi,σ (1) + ϕ∗i,π (1′)ϕi,π (1) i ∈σ
i ∈π
= γσ (1′, 1) + γπ (1′, 1) .
(19)
So the density can also be decomposed into two parts, n (1) = nσ (1) + nπ (1) .
(20)
The Fermi hole (Eq. (17)) can be decomposed as X hole X hole X hole nhole x (2|1) = nσ−σ (2|1) + nσ−π (2|1) + nπ−π (2|1),
(21)
where, for example, X hole nσ−σ (2|1) = −
|γσ (1, 2)|2 n (1)
(22)
and a similar form for the π part. The cross term is γσ (1, 2)γπ (2, 1) + γσ (2, 1)γπ (1, 2) . (23) n (1) It is to be noted that the partial sums of these holes yield nσ (1) X hole nσ−σ (2|1)d2 = − , n(1) nπ (1) X hole nπ−π (2|1)d2 = − , (24) n(1) X hole nσ−π (2|1)d2 = 0. X hole nσ−π (2|1) = −
Fig. 2 shows this decomposition for the Hartree-Fock wave function along the hexagon lines at 1 bohr above the molecular plane while the reference electron is placed at 1 bohr above the atom C3. We see that curves for the σ − σ and σ − π parts quickly decrease to zero, being already very small at the nearest carbon atoms. The contribution from the π − π part has a small tail at the para site, but the magnitude at meta sites is zero. Fig. 3 is a similar decomposition, but along the vertical line perpendicular to the molecular plane, passing through the reference position at 1 bohr above C3. One observes that the greatest (negative) depth of the hole is not exactly at the referX hole ence position, but a little below 1 bohr. Both the nσ−σ (2|1) X hole and nπ−π (2|1) are definite negative (cf. Eq. (22)) and are symmetric with respect to the plane as a consequence of their X hole definitions. But nσ−π (2|1) is negative above the molecular
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-5
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
the plane, the Fermi hole takes the shape of the localized orbital above the plane. B. XC hole and Coulomb hole in benzene
FIG. 2. Benzene: The decomposition of the Fermi hole into contributions from the σ − σ, σ − π, and π − π pairs along the hexagon lines parallel the molecular plane pass through the reference electron at 1 bohr above the C3 atom.
plane and positive below the plane. This positive part almost completely annihilates the negative π − π and σ − σ contributions below the plane. Because of this compensating effect X hole from nσ−π (2|1), the Fermi hole below the molecular plane is much shallower than above the plane. This is in agreement with what is known about the Fermi hole behavior: if the density at the reference electron position is predominantly contributed by a single localized orbital, resulting from an orbital localizing transformation of the Hartree-Fock orbitals, then the Fermi hole will take approximately the shape of this localized orbital squared.82 Now orbital localization leads, in conjugated systems, to loss of the (anti)symmetrical π orbitals, instead it produces orbitals that localize either above or below the plane of the molecule, cf. the well known “banana bonds” above and below the plane in ethene. With our reference position above
We perform correlated calculations with the traditional CASSCF method, where the active space of course has to be limited. A π-only correlation, with only the 6 π electrons in the active space and the rest inactive, is perfectly feasible. The σ − π correlation can be studied by taking σ electrons in the active space. For example, we have done CASSCF calculations with inclusion of 6 π electrons + 4 σ electrons in the active space. A CASSCF wave function of 10 electrons in 14 active orbitals resulted in a CI wave function with 1 002 001 configuration state functions. The setup yields a total energy of −230.839 a.u., which is lower than the Hartree-Fock energy by 0.117 a.u. To be able to discern the different contributions from the π and σ electrons, in the following, we report results of a CASSCF calculation of the 6 π electrons in an active space of 3 occupied π orbitals and 21 virtual orbitals. All the σ electrons are in the inactive core, their interaction with the π electrons is obtained through the anti-symmetrization of the Φσ and Φπ group functions, Ψ = N AΦσ Φπ ,
(25)
where Φσ is the determinant with the Hartree-Fock σ orbitals, and Φπ is a CI wave function consisting of the determinant of the Hartree-Fock π orbitals as leading term, with many single, double, triple, etc., substitutions added (excitations from occupied π orbitals to both virtual π and σ orbitals). Clearly, the σ − π correlation interaction is approximated in this case. But the π − π correlation can be studied more extensively by inclusion of more virtual orbitals in the active space. A wave function with 1 163 800 configuration state functions constructed from the 24 active orbitals is obtained with a total energy of −230.814 a.u. The results for group functions of Ref. 47 are then immediately applicable (cf. also Ref. 83). The one-particle density matrix is just the sum of the one-particle density matrices of the two group functions. The one-particle density matrix is in diagonal form in the basis of Hartree-Fock σ orbitals for Φσ , and in the basis of natural orbitals for Φπ , γΨ (1′, 1) = γσ (1′, 1) + γπCI (1′, 1) = ϕ∗i,σ (1′)ϕi,σ (1) + λ j ϕ∗j (1′)ϕ j (1). i ∈σ
(26)
j
Real-space representation of electron correlation in π -conjugated systems Jian Wang1,a) and Evert Jan Baerends2,a) 1
School of Science, Huzhou University, Zhejiang 10083, China Afdeling Theoretische Chemie, FEW, Vrije Universiteit, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands 2
(Received 5 March 2015; accepted 6 May 2015; published online 29 May 2015) π-electron conjugation and aromaticity are commonly associated with delocalization and especially high mobility of the π electrons. We investigate if also the electron correlation (pair density) exhibits signatures of the special electronic structure of conjugated systems. To that end the shape and extent of the pair density and derived quantities (exchange-correlation hole, Coulomb hole, and conditional density) are investigated for the prototype systems ethylene, hexatriene, and benzene. The answer is that the effects of π electron conjugation are hardly discernible in the real space representations of the electron correlation. We find the xc hole to be as localized (confined to atomic or diatomic regions) in conjugated systems as in small molecules. This result is relevant for density functional theory (DFT). The potential of the electron exchange-correlation hole is the largest part of v xc , the exchange-correlation Kohn-Sham potential. So the extent of the hole directly affects the orbital energies of both occupied and unoccupied Kohn-Sham orbitals and therefore has direct relevance for the excitation spectrum as calculated with time-dependent DFT calculations. The potential of the localized xc hole is comparatively more attractive than the actual hole left behind by an electron excited from a delocalized molecular orbital of a conjugated system. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921725] I. INTRODUCTION
π-conjugated systems are characterized by alternating single (σ) and double (σ and π) bonds. Benzene is a prototype for the cyclic conjugated structures which are known for their special properties (cf. the heme group, chlorophyll, and graphene). The associated concept of aromaticity has received much interest.1,2 The π electrons are commonly conceived as particularly mobile, and for instance when exposed to magnetic fields, a large ring current can be induced around the benzene ring.3–5 Kekulé6 in the 19th century conceived the cyclic structure with alternating single bond and double bonds for benzene. Depending on whether a given bond is a single bond or double bond, there are two possibilities for the Kekulé structure. To comply with the experimental evidence that bond lengths in benzene are the same, Kekulé suggested the two alternative structures are in “oscillation” or “resonating.” This notion found a natural explanation in quantum mechanics, which allows enhanced stability by superposition of (trial) wave functions. Pauling and Wheland7,8 found that the superposition of the Kekulé electronic structures and other possible electronic structures, such as Dewar structures, results in energetic stabilization. This “resonance” effect in benzene had been previously proposed by Slater.9 It provides an energetic mechanism of enhanced stability for the symmetric structure. In the early molecular orbital based calculations,10–12 benzene was modelled by separating the electrons into orthogonal σ and π subspaces. The σ electrons occupy a closedshell configuration and make up the molecular σ core, while the π electrons move in orbitals orthogonal to this core. The a)Electronic addresses: [email protected] and [email protected]
0021-9606/2015/142(20)/204311/12/$30.00
molecular orbitals are constructed according to the molecular symmetry, which is D6h for benzene. The electronic structure from the molecular orbital picture provides a roughly correct explanation for the experimental spectra. In conjugated chains, the experimental spectra justify the delocalized picture of the π molecular orbitals.13–15 There have been conflicting views on the delocalized nature of the π electron subsystem in aromatic compounds. Ab initio self-consistent field (SCF) molecular orbital calculations16,17 and various high level configuration interaction (CI) calculations18–23 provide, as expected, detailed comparison with experiments. However, Cooper, Gerratt, and Raimondi24 have used valence bond models to argue that the π electrons can be pictured as residing in almost localized orbitals on the individual carbon atoms of benzene, defying the text-book picture of its delocalized π ring structure. In their calculation, they started from initially delocalized symmetric molecular orbitals for the π-electrons but lifted symmetry restrictions during subsequent optimization. Shaik, Hiberty, and co-workers25,26 have argued that the π-electrons of benzene tend to pair with each other which will distort the equal sided hexagon, and the molecule maintains 6-fold symmetry only because of the σsystem (cf. also the Kohn-Sham MO analysis of Pierrefixe and Bickelhaupt27). Of course, orbitals are not observables, and for instance, the Hartree-Fock wave function is invariant under arbitrary unitary transformation of the orbitals.28 Both molecular orbital (Hartree-Fock plus CI) and valence bond type representations of the wave function may lead to different pictures of the localized versus delocalized nature of the bonding in conjugated compounds. Similarly in solids, one may choose the delocalised Bloch function or the localised Wannier functions to describe the electronic structure.29–31 In contrast,
142, 204311-1
© 2015 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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-2
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
the pair density and other quantities derived from the twoparticle density matrix are defined, unambiguous, real-space quantities. It will be interesting to see whether conjugation, delocalization, and aromaticity can be traced in the description of electron correlation with these spatial probabilities. We focus on the pair density and derived quantities, such as Fermi hole, exchange-correlation (XC) hole, Coulomb hole, and conditional density. We note that such a study has immediate relevance for the much used Kohn-Sham model of density functional theory (DFT). Often the XC energy of DFT Exc (which contains a kinetic energy contribution) is expressed through the coupling constant integrated XC hole, 1 hole d 3r 1 n(1) d 3r 2 n¯ xc (2|1)/r 12. (1) Exc[n] = 2 hole The coupling constant integrated XC hole n¯ xc (2|1) in DFT contributes not only to the electron-electron Coulomb energy but also takes care of the kinetic energy difference between the exact kinetic energy and that modeled by the Kohn-Sham orbitals, T − Ts . The kinetic energy contribution to the correlahole tion is folded into the XC hole n¯ xc by the so-called coupling32–34 constant integration, 1 hole n¯ xc (2|1) = dλ nhole (2) xc,λ(2|1), 0
(2|1) is the hole for a system in which the Coulomb where nhole xc,λ interaction between electrons is given by λ/r 12, but the external
potential is adapted in such a way that the density remains equal to the exact one. The coupling constants with 0 ≤ λ ≤ 1 connect the Kohn-Sham system (λ = 0) to the physical system (λ = 1). It is very difficult to construct the XC hole with the coupling constant integration.35–37 However, a large part of the energy density for Exc and of the Kohn-Sham potential vxc(r) can be obtained directly from the full hole nhole (2|1),38–40 xc,λ=1 1 hole v + vc,kin, 2 xc hole vxc = vxc + vc,kin + vresp, hole vxc = d 3r 2 nhole xc,λ=1(2|1)/r 12, ϵ xc =
(3)
where ϵ xc is the XC energy density per particle, from which the XC energy may be obtained by Exc[n] = n(1)ϵ xc([n]; 1)d1. Here, vc,kin is the energy density that accounts for the correlation correction in kinetic energy, T − Ts = n(1)vc,kin(1)d1 and vresp is the so-called response part of vxc. The hole potential hole vxc is usually the largest part of vxc. Notably in situations of bond breaking, the components vc,kin and vresp may also become significant,41–43 and this is in particular the case when timedevelopment is studied with time-dependant DFT (TDDFT) (beyond the linear response regime).44 However, here we will exclusively deal with the static ground state and equilibrium hole geometries. The more or less attractive nature of the vxc potential, depending on the localized or delocalized nature of the hole, affects directly the orbital shapes and the orbital energies of the Kohn-Sham orbitals. If there is a reasonable correspondence between the size of the exchange-correlation hole and the size of the electron hole that is created in an excitation, the orbital energy difference between the occupied and virtual orbitals involved in the excitation will be a very
good first approximation to the energy of the excitation.45 On the other hand, if there is a large mismatch, as is the case in solids and as may be the case in π-conjugated systems if the exchange-correlation hole is less delocalized than the π orbitals involved in an excitation, the orbital energy difference may not be a good first approximation to the excitation energy, and the time-dependent DFT calculations may be expected to be problematic. The paper is organized as follows. First, we outline the computational method for the pair density and its derived quantities. Then in Secs. III A and III B, the holes in benzene are analyzed. Next, in Sec. III C, the holes are compared to those of related systems such as hexatriene and ethene, and of cyclohexane, and in Sec. III D to holes from electron gas models. Finally, the conclusions are drawn.
II. COMPUTATIONAL METHOD
The pair density is defined as46,47 Γ(1, 2) = N(N − 1) × Ψ∗(123 · · · N)Ψ(123 · · · N)d3 · · · dN,
(4)
where N is the total number of electrons. The normalization follows McWeeny’s convention,47 which is a factor 1/2 different from Löwdin’s.46,48 The quantity can be interpreted as the probability that two electrons be simultaneously found at positions r 1 and r 2 with spins s1, and s2, respectively. Each electron has three spatial coordinates and one spin coordinate. The integral implies integration over space coordinates and summation over spin coordinates. If the spin coordinates are summed out, the pair density is then a function of the spatial coordinates of the two electrons, Γ(r1, r2). To visualize Γ(r1, r2) in three dimensions, we may fix the position of one electron, say r1 (the reference electron), and then look at the pair density (or the exchange-correlation hole) as a function of the position (r2) of the second electron. The pair density has previously also been studied through the reduced quantity: intracule.49–51 The XC hole surrounding an electron at 1 is defined as hole (2|1) = nxc
Γ(1, 2) − n(1)n(2) , n(1)
where n(1) is the density defined as n(1) = N Ψ∗(12 · · · N)Ψ(12 · · · N)d2 · · · dN.
(5)
(6)
The effects of exchange and correlation on the electron distribution are described by the XC hole. For the model system of the electron gas, the hole has been extensively studied,35,52–58 and the results have been used to construct XC functionals in DFT.59 The XC holes in a few prototype systems, such as H2,41,60–63 He,64,65 atomic Si,66 and crystal Si,67–69 have been carefully studied. Here, we study the XC hole for π-conjugated systems, previously the Coulomb hole around the bond midpoint has been reported.70 The density is related to the pair density as 1 Γ(1, 2)d2. (7) n(1) = N −1
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-3
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
From this, one can derive the important sum rule for the XC hole, hole nxc (2|1)d2 = −1. (8) A closely related quantity is the conditional density defined as the probability to find an electron at 2 when one is known to be at 1, hole ncond(2|1) = n(2) + nxc (2|1).
(9)
It integrates to N − 1 electrons because the XC hole integrates to −1, and the density to N. The spherically averaged conditional density in the Hartree-Fock approximation can be expanded as71 1 Dσ s2 + · · ·, (10) 3 where the arguments (r, s) denote the spherical average on a shell of radius s around the reference point r, and σσ (r, s) = ncond
1 (∇ρσ )2 , Dσ = τσ − 4 ρσ
(11)
τσ = Σi |∇ψi,σ |2 is the kinetic energy of electrons with spin σ, 2 and 14 (∇ρρ σσ ) is the von Weizsäcker kinetic energy functional. A smaller Dσ implies smaller probability to find a second electron around the reference electron, which may be interpreted as a stronger electron localisation. On this basis, the quantity Dσ has been used by Becke to propose an “electron localization function” (ELF).71 A possible connection between ELF and the phenomenon of aromaticity has been explored by several authors.72–74 The connection is however not very direct, and we will focus on the quantities directly related to the pair density, i.e., the exchange and correlation holes. Another quantity of interest is the one-particle density matrix defined as ′ γ(1 , 1) = N Ψ∗(1′2 · · · N)Ψ(12 · · · N)d2 · · · dN. (12) Its eigenfunctions are the natural orbitals, which were believed to have certain optimal properties for the description of electron correlation (shortest expansion)46 although this does not appear to be the case for general systems,75 γ (1′, 1) = λi χ∗i (1′) χ i (1), (13) i
where λi are the occupation numbers of the natural orbitals χ i . The density is just the diagonal of γ, n(1) = γ(1, 1). In the method of CI, one usually selects a set of M spin orbitals (M > N). One can then build various determinants (or configurations), Φi , each from N of those orbitals. The variational wave function can be written as a linear combination of those configurations, Ψ= Ci Φi . (14) i
The simplest case is M = N spin orbitals (Hartree-Fock), the wave function in this case is a single determinant. The pair density in this case has been worked out by Dirac76,77 ΓHF(1, 2) = γHF(1, 1)γHF(2, 2) − γHF(1, 2)γHF(2, 1),
(15)
where γHF (1, 1′) =
N
ϕi (1)ϕ∗i (1′).
(16)
i
Here, ϕi are the Hartree-Fock spin orbitals. The XC hole in this case reflects the effect of electron exchange, it is by definition the exchange hole, also called the Fermi hole, nhole x (2|1) = −
|γHF1, 2|2 . nHF (1)
(17)
The negative sign implies the reduction of probability to find the second electron at position 2 while a reference electron is at 1. The Fermi hole also integrates to −1, i.e., obeys the same sum rule as the total XC hole (see Eq. (8)). The Coulomb correlation hole can be defined as the difference between the XC hole and Fermi hole,60 hole hole nhole c (2|1) = nxc (2|1) − n x (2|1) .
(18)
Since both the XC hole and Fermi hole integrate to −1, the Coulomb hole integrates to zero. It only represents a redistribution in space of the electron density, i.e., reduced density in one region will show up in other regions.
III. CALCULATION RESULTS
We use the GAMESS package78 to calculate the wave function. The CI wave function is obtained using the multiconfiguration self-consistent field (MCSCF) or complete active space self-consistent field (CASSCF) method. The wave function is calculated in the graphical unitary group approach;19 as a byproduct, the elements of the two-matrix can be exported. Dunning’s correlation-consistent basis set cc-pVDZ79,80 is chosen. In terms of contracted functions, the basis set is 2s1p for H and 3s2p1d for C. A. Fermi hole in benzene
The optimal benzene structure has a C–C bond length 1.397 Å and C–H bond length of 1.084 Å. The Hartree-Fock energy is −230.722 a.u., as compared to Ermler and Kern’s −230.749 a.u.16 using an extended basis set (2s1p for H and 4s2p1d for C). Fig. 1 shows the Fermi hole where the reference electron is chosen at 1 bohr directly above a carbon atom. The atomic 2pz orbital has its maximum just below 1 bohr. The upper panel is a picture of the hole in three dimensions, from which it is clear the hole is quite localized around the reference electron, with only small contributions below the plane and at other atomic sites (ortho, para). The lower panel of Fig. 1 gives more quantitative information of the Fermi hole in the plane passing through the reference position and parallel to the molecular plane. The Fermi hole in benzene and other π-conjugated systems has been previously studied by Bader et al.81 They emphasized the feature at the para site as a signature of π delocalization. But in fact the amplitude at the para-site is very small (the magnitude there is in the order of 10−3 e/(bohrs)3). The surprise is rather that the Fermi hole is actually not delocalized.
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-4
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
FIG. 1. Benzene: The Fermi hole from the Hartree-Fock wave function. The reference electron is at 1 bohr above a C atom. The upper panel shows the 3-dimensional distribution of the hole, while the lower panel gives the cross section in the plane through the reference electron and parallel to the molecular plane.
We will compare to a non-conjugated system (cyclohexane, see Sec. III C) and note that the Fermi hole is not more delocalized in the π-conjugated system than in the non-conjugated system. Physically, a localized Fermi hole is favorable because it maximizes the attractive electron-hole interaction. Compared to the spherical hole in the uniform electron gas, the Fermi hole in benzene has a remarkable structure. It will be interesting to identify how its features are related to the π and σ bondings. For the Hartree-Fock approximation, it is straightforward to decompose the π and σ contributions since there are distinct σ and π orbitals. We may write the oneparticle density matrix as a sum of π and σ orbital contributions, γHF (1′, 1) = ϕ∗i,σ (1′)ϕi,σ (1) + ϕ∗i,π (1′)ϕi,π (1) i ∈σ
i ∈π
= γσ (1′, 1) + γπ (1′, 1) .
(19)
So the density can also be decomposed into two parts, n (1) = nσ (1) + nπ (1) .
(20)
The Fermi hole (Eq. (17)) can be decomposed as X hole X hole X hole nhole x (2|1) = nσ−σ (2|1) + nσ−π (2|1) + nπ−π (2|1),
(21)
where, for example, X hole nσ−σ (2|1) = −
|γσ (1, 2)|2 n (1)
(22)
and a similar form for the π part. The cross term is γσ (1, 2)γπ (2, 1) + γσ (2, 1)γπ (1, 2) . (23) n (1) It is to be noted that the partial sums of these holes yield nσ (1) X hole nσ−σ (2|1)d2 = − , n(1) nπ (1) X hole nπ−π (2|1)d2 = − , (24) n(1) X hole nσ−π (2|1)d2 = 0. X hole nσ−π (2|1) = −
Fig. 2 shows this decomposition for the Hartree-Fock wave function along the hexagon lines at 1 bohr above the molecular plane while the reference electron is placed at 1 bohr above the atom C3. We see that curves for the σ − σ and σ − π parts quickly decrease to zero, being already very small at the nearest carbon atoms. The contribution from the π − π part has a small tail at the para site, but the magnitude at meta sites is zero. Fig. 3 is a similar decomposition, but along the vertical line perpendicular to the molecular plane, passing through the reference position at 1 bohr above C3. One observes that the greatest (negative) depth of the hole is not exactly at the referX hole ence position, but a little below 1 bohr. Both the nσ−σ (2|1) X hole and nπ−π (2|1) are definite negative (cf. Eq. (22)) and are symmetric with respect to the plane as a consequence of their X hole definitions. But nσ−π (2|1) is negative above the molecular
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: 61.16.135.116 On: Sat, 30 May 2015 06:23:46
204311-5
J. Wang and E. J. Baerends
J. Chem. Phys. 142, 204311 (2015)
the plane, the Fermi hole takes the shape of the localized orbital above the plane. B. XC hole and Coulomb hole in benzene
FIG. 2. Benzene: The decomposition of the Fermi hole into contributions from the σ − σ, σ − π, and π − π pairs along the hexagon lines parallel the molecular plane pass through the reference electron at 1 bohr above the C3 atom.
plane and positive below the plane. This positive part almost completely annihilates the negative π − π and σ − σ contributions below the plane. Because of this compensating effect X hole from nσ−π (2|1), the Fermi hole below the molecular plane is much shallower than above the plane. This is in agreement with what is known about the Fermi hole behavior: if the density at the reference electron position is predominantly contributed by a single localized orbital, resulting from an orbital localizing transformation of the Hartree-Fock orbitals, then the Fermi hole will take approximately the shape of this localized orbital squared.82 Now orbital localization leads, in conjugated systems, to loss of the (anti)symmetrical π orbitals, instead it produces orbitals that localize either above or below the plane of the molecule, cf. the well known “banana bonds” above and below the plane in ethene. With our reference position above
We perform correlated calculations with the traditional CASSCF method, where the active space of course has to be limited. A π-only correlation, with only the 6 π electrons in the active space and the rest inactive, is perfectly feasible. The σ − π correlation can be studied by taking σ electrons in the active space. For example, we have done CASSCF calculations with inclusion of 6 π electrons + 4 σ electrons in the active space. A CASSCF wave function of 10 electrons in 14 active orbitals resulted in a CI wave function with 1 002 001 configuration state functions. The setup yields a total energy of −230.839 a.u., which is lower than the Hartree-Fock energy by 0.117 a.u. To be able to discern the different contributions from the π and σ electrons, in the following, we report results of a CASSCF calculation of the 6 π electrons in an active space of 3 occupied π orbitals and 21 virtual orbitals. All the σ electrons are in the inactive core, their interaction with the π electrons is obtained through the anti-symmetrization of the Φσ and Φπ group functions, Ψ = N AΦσ Φπ ,
(25)
where Φσ is the determinant with the Hartree-Fock σ orbitals, and Φπ is a CI wave function consisting of the determinant of the Hartree-Fock π orbitals as leading term, with many single, double, triple, etc., substitutions added (excitations from occupied π orbitals to both virtual π and σ orbitals). Clearly, the σ − π correlation interaction is approximated in this case. But the π − π correlation can be studied more extensively by inclusion of more virtual orbitals in the active space. A wave function with 1 163 800 configuration state functions constructed from the 24 active orbitals is obtained with a total energy of −230.814 a.u. The results for group functions of Ref. 47 are then immediately applicable (cf. also Ref. 83). The one-particle density matrix is just the sum of the one-particle density matrices of the two group functions. The one-particle density matrix is in diagonal form in the basis of Hartree-Fock σ orbitals for Φσ , and in the basis of natural orbitals for Φπ , γΨ (1′, 1) = γσ (1′, 1) + γπCI (1′, 1) = ϕ∗i,σ (1′)ϕi,σ (1) + λ j ϕ∗j (1′)ϕ j (1). i ∈σ
(26)
j
