J. Phys. Chem. 1993,97, 8617-8624
8617
Crystal Orbital Hamilton Populations (COHP). Energy-Resolved Visualization of Chemical Bonding in Solids Based on Density-Functional Calculations Richard Dronskowski’ Max-Planck-Institut f i r Festkikperforschung, Heisenbergstr. I, 0-70569 Stuttgart, Germany
Peter E. Blochl IBM Research Division, Ziirich Research Laboratory, CH-8803 Riischlikon, Switzerland Received: January 20, 1993
After giving a concise overview of the current knowledge in the field of quantum mechanical bonding indicators for molecules and solids, we show how to obtain energy-resolved visualization of chemical bonding in solids by means of density-functional electronic structure calculations. On the basis of a band structure energy partitioning scheme, i.e., rewriting the band structure energy as a sum of orbital pair contributions, we derive what is to be defined as crystal orbital Hamilton populations (COHP). In particular, a COHP(e) diagram indicates bonding, nonbonding, and antibonding energy regions within a specified energy range while an energy integral of a COHP gives access to the contribution of an atom or a chemical bond to the distribution of one-particle energies. A further decomposition into specific atomic orbitals or symmetry-adapted linear combinations of atomic orbitals (hybrids) can easily be performed by using a projector technique involving unitary transformations. Because of its structural simplicity and the availability of reliable thermodynamic data, we investigate the bonding within crystalline silicon (diamond phase) first. As a basis set, both body-centered-cubic screened and diamond screened atomic-centered tight-binding linear muffin-tin orbitals (TB-LMTOs) are used throughout. The shape of COHP versus energy diagrams and the significance of COHP subcontributions (s-s, sp3-sp3) are analyzed. Specifically, the difference between the COHP energy integral (one-particle bond energy) and the experimental bond energy is critically examined. While absolute values for the one-particle bond energy show a high basis set dependence due to changing on-site (crystal field) COHP terms, the shape of off-site (bonding) COHP functions, elucidating the local bonding principle within an extended structure, remains nearly unaffected.
1. Introduction Numerical quantum chemistryfaces a challenge when it comes to the extraction of chemical information out of highly complex computations. The complete set of classic (and highly useful) bonding descriptorscannot be unambigously defined by quantum mechanics since they are not measurable quantities of a microscopic system. Nevertheless, there are some interpretational schemes for molecular wave functions that were developed in the past 40 years. 1.1. Bonding Indicators for Molecules. In his pioneering work Mulliken showed how to assign electrons both to bonds and to atomic centers.’ If there is a molecular orbital $I built up from two atomic orbitals cpr and cp,, the bonding interaction is carried through the overlap integral S,,,whereas , the overlap population C*,,C,.S,,~ may be interpreted as a kind of chemical bond order. Traditionally, overlap populations are divided up equally between atoms for the calculation of gross values. There is strong ambiguity and basis set dependency inherent in Mulliken’s scheme. There have been further developments and generalizations to polyatomics by Davidson2 and Roby,3 involving the molecular one-particle density operator and a projection operator technique, giving access to occupation numbers N and shared electron numbers u. By constructing minimal sets of modified atomic orbitals (slightly deformed Hartree-Fock atomic orbitals), Ahlrichs et al. suggested another technique to make these population analyses less basis-set depe~~dent.~.S It has been empirically verified that shared electron numbers are closely related to bond energies of small molecules, Le., they may serve as a measure of the bond strength. In the end, this result
strengthensthevalidity of the famous approximation of Wolfsberg and Helmholz.6 1.2. Bonding Indicators for Extended Structures. In contrast to the field of molecular quantum mechanics, less attention has been given to the concept of local bonds in extended systems, probably because of historic reasons: Solid-state theory has mainly been concerned with the simplest elemental structures of metals in which bonding is assumed to be metallic, i.e., widely isotropic as in Drude’s classic idea. Therefore, a local way of looking at chemical bonds in solids might seem to be ill-fitting. However, since solid-state chemistry is a burgeoning field and solid-state compounds can be looked upon as being nothing but large molecules,7 an understanding of solid-state structures and their local bonding characteristics is indispensible. In addition to the nonquantum mechanicalschemes, e.g.,variousbond length-bond strength re1ationships,”l0 we know at least two related, semiempirical formulations of bonding indicators for the solid. A very useful and popular tool has been introduced by Hughbanksand Hoffmannwithin the framework of semiempirical extended Hiickel theory,lIJ2 a special tight-binding method with overlap. Once the electronic band structure has been calculated, Mulliken’s overlap population technique is applied to the crystal, Le., all levels in a certain energy interval are investigated as to their bonding proclivities, measured by R{c*,c,S,,) (positive = bonding, zero = nonbonding, negative = antibonding).13 In this way, overlap population-weighted densities of states are defined, evaluating all pairwise orbital interactions between two atoms by the product of overlap matrixand density-of-states (DOS)matrix elements. While the energy increases through the bands, bonding characteristics appear as a function of the energy. This method
0022-365419312097-8617$04.00/0 0 1993 American Chemical Society
Dronskowski and Bliichl
8618 The Journal of Physical Chemistry, Vol. 97, No. 33, 1993
was dubbed crystal orbital overlap population (COOP), serving as the solid-state analogue to the molecular bond order. The COOP method has found a large number of interesting applications in bonding studies of extended materials,14 and it is now a standard tool for tight-binding extended Hiickel (EH) calculations. However, there are limitations. Since the COOP method is essentially based on an electron number partitioning scheme by means of Mulliken's population analysis, its extension to first-principles calculations seems to be problematic. First, the basis set dependence will generate problems if one moves away from the Slater orbitals of EH theory. Second, even within EH theory, there are problems if one tries to compare overlap populations from d-d and p p bonding interactions, for instance, simply because of the consequences of different spatial extent of atomic functions.15 The other bonding descriptor for solids goes back to the work of Sutton et al., who also used a semi-empirical theory, the tightbinding bond (TBB) model.l6 They partitioned the cohesive energy of a solid into atomic centered, bond centered, and additional electrostatic and exchange-correlation contributions, approximated by a sum of pair potentials in order to obtain interatomic forces. These energy-partitioning schemes have a long history in quantum ~hemistry,l~-'~ and Coulson was probably the first to set up such an expression.2o Within the density-functional scheme of Harris and Foulkes21-22 the total energy of a solid can be expressed as a functional of an approximate charge density and the one-electron eigenvalues obtained by solving SchrMinger's equation for the potential created by this charge density once. Local charge neutrality to eliminate long-range Coulomb terms is advocated as the simplest approximation. The covalent bond energy between two atoms is given as the product between off-diagonal density matrix elements nWvrtypically calculated by a recursion method,23-25 and the corresponding elements of the Hamiltonian matrix H,,". Plots of u-, r-,and &bond energies for canonical d bands in closed packed metals are available.16 If one wants to set up a bonding indicator within a first principles density-functional theory, it seems that an energypartitioning scheme (like in TBB) is favorable with respect to an electron-partitioning scheme (COOP). By doing so, problems arising from the principal basis-set dependence of the overlap integral may be kept to a minimum. However, to facilitate the interpretation of such a tool, it is worthwhile to make it similar to the COOP method by first generating an energy-resolved bonding descriptor which may be energy-integrated in the sequel. Thus, the graphical identification of levels with different bonding characteristics is much easier. 2. Theory
We start by constructing one-electron wave functions in the solid by a linear combination of atomic centered orbitals X R L with mixing coefficients (eigenvectors) U R L as ~ in
where j stands for the band index whereas R specifies the atomic site and L is a shortened notation for the angular momentum quantum numbers L lm. Note that stands for a double sum running first over all atoms R within the primitive unit cell and second over all atomic orbitals L centered on R. Because of the orthogonality of the wave function we have
with the abbreviation for the overlap matrix
The one-electron problem may be written as in
with the Hamiltonian matrix H R L X L ~5 ( X R L P I X R L ' )
=
( x ~ -I v2+ UO)IXRL')
(5)
Here and in the following we use Rydberg atomic units. The so-called band structure energy Eb"&is defined as the sum of all occupied one-electron eigenvalues e/, i.e.
with occupation numbersfi, 0 S f i I2 , assuming the non-spinpolarized case here and in the sequel. It can also be expressed as the energy integral of the distribution of one-electron eigenvalues:
Now if we multiply (4) from the left by the conjugate complex eigenvectors, we have
which can be shortened by applying (2) to read
Finally, by insertion of (9) into (7) we arrive at
R L RZ.'
Thus, the distribution of all one-particle energies has been rewritten into a sum of pair contributions, which we call crystal orbital Hamilton populations (COHP). While using a shortrange (tight-binding) orbital basis set, the sum of all pairwise interactions rapidly converges in real space, similar to the formalism based on approximated and improved density operators within the TBB model. For the band structure energy we have
where the density matrix is obtained by energy integration of the DOS matrix: ~ m p , L= ,
J"dS
NRL.,R,L,(e)
(13)
There is an advantage in using N ( r ) instead of n because N(t) is independent from the Fermi level e ~ Furthermore, . a graphical representation using N(e) helps to visually identify bonding, nonbonding, and antibonding states, as will be seen later. This is useful when discussing the effects on the stability when electrons are added or removed by interstitial or substitutional dopants.
Crystal Orbital Hamilton Populations The differencebetween the overlap population-weighted DOS and the energy partitioning presented here lies in the use of H instead of S.It is clear that
but on the other hand, such a relation is not valid for any single term since there will be differencesbetween the two partitionings if single terms (bonds) are investigated. The on-site COHP ( R R’) terms, corresponding to atomic contributions,consist of n X n matrices if the number of orbitals per site is n. There are also off-site COHP terms (R # R’), emerging from bonding interactions between the atoms, so that the band structure energy can be expressed as
UjeFCej - e> = C i
H ~ ~ L N R L , + ~ ( ~ ) RL
The Journal of Physical Chemistry, Vol. 97, NO. 33, 1993 8619 with an effective one-particle potential of the form
v(,) = Uext(i)+ U H ( 3 + uxcIn(i)l
(19)
which consists of the external potential, a Hartree potential UH = JdY 2n(Y)/li-YIand thepotential for exchangeandcorrelation vxc{n(i))= dE,,/dn(i). To obtain the cohesive energy we need to subtract the total energies of the atoms:
Ea‘(n(i)j=
1 Diei- s d i n(i)[-u,(i) 2
+ uxc{n(i)j] +
i
from the total energy of the crystal ECr
EE’{n(Z))= JcFdeDje,d(ej - e)
- Jdi
n(i)
i
EC
R L Z R L , w H R L ~ # L J v ~ , R L ~ ( ~ N ( 15 )
The second part represents the covalent contributionsof the bonds. If there are bondingcontributions, the system undergoes a lowering of its energy, indicated by negative off-diagonal COHPs. If a structure is destabilized by antibonding contributions,there will be positive off-diagonal COHPs. Note that according to this convention COHP and COOP have in principle different plus/ minus signs if there is a fixed relationship between H and S as expressed by the Wolfsberg-Helmholz approximation6 within EH theory. Can off-diagonal COHPs be interpreted as a measure of real bond energies? Generally, there is a complication arising from the fact that the band structure energy is not identical with the total energy, the latter including many-particle as well as electrostatic effects, which are essential for chemical bonding. Of course, it has been customary to investigate the variation of one-particle eigenvalues with changing internuclear geometry (as in Walsh diagrams). However, to obtain at least semiquantitative conclusions, we will take a closer look at the foundations of total energy calculations. Within density-functional (DF)theory, the density and the total energy of the electronic ground state are given by the minimum of a total-energy functional E ( n ( i ) )of the electron density n ( i ) . The functional
can be divided into the kinetic energy of noninteracting electrons with the same density, a classic Coulomb energy of the charge density, and a functional for exchange and correlation. veXt(i)is the Coulomb potential of the fixed nuclei. The energy due to exchange and correlation is usually expressed as
The above expression is exact, but the correct energy functional for exchange and correlation is not known. Therefore one introduces the local density approximation so that t, depends only on the local value of the electron density. The function c,(n) is usually extracted from an exact many-particlecalculation of the free electron gas. In the density functional theory, the solution of the manyparticle problem has been replaced by the self-consistent solution of the KohnSham equations:26,27
where Z R is the charge of the nucleus at site R.If we express the sum of one-particle energies by the sum of COHPs, we obtain the cohesive energy
Y off-site COHP FEF
on-site COHP
atomic eigenvalues
using eqs 15,20, and 21. We have shortened the notation for the small difference charge density and difference exchange-correlation parts. It should be noted that the onsite COHPs, the atomic eigenvalues, and the electrostatic energy difference carry undetermined constants. These terms depend on the choice of the boundary conditions for the electrostatic potential at infinite distance from the system under study. However, despite this arbitrariness, the sum of all three terms together is well defined, as is the cohesive energy. For a molecule there is a natural choice for the boundary condition, namely, natural boundary conditions. For a crystal, one usually chooses the average electrostatic potential to define the zero of the energy scale. With this choice all three terms are well defined individually. Thus, in general one may not expect integrated off-site COHPs to be identical with real bond energies, which sum up to the cohesive energy. This would only be the case if the difference between ( i ) integrated on-site COHPs, Le., the sum of atomic ‘crystal field” terms, and (ii) sum of eigenvalues of the free atom (second line) would accidentally cancel the sum of contributions from difference charge density, exchange-correlation, and Madelung terms (third line). A similar observation has already been made in the DF molecular studies by Gunnarsson et al., who stated that the nonlinear dependence of the total energy on the occupation numbers renders a unique quantitative definition of individual bonds impossible.28 However, all quantum-chemical experience suggests that the eigenvalue terms are an important indicator for the formation of bonds. So, by looking at the offsite COHPs one may nevertheless expect to detect the main bonding characteristics.
8620 The Journal of Physical Chemistry, Vol. 97, No. 33, 19’93
We do not propose to describe binding exclusively by off-site COHPs. It would not lead to an appropriate description of binding for example in salts, which gain their cohesive energy mostly from electron transfer and electrostatic attraction, whereas covalent bond formation is of secondary importance. However, we found that off-site COHPs are an extremely useful tool to analyze covalent bonding. The discussion given above explains how to relate the bonding descriptors to the cohesive energy and shows which additional contributions must be discussed to obtain a complete description of bonding. The simplest way to investigate the bonding between two interacting atoms in the solid would be to look at the complete COHP between them, taking all valence orbitals into account. However, it may sometimesbe useful to focus on pair contributions of some specific orbitals. Generally, any such orbital interaction may be extracted from the total off-site COHP by using the projection operator technique that is described in the Appendix. Thus, there is accessto symmetry-adaptedorbitals (hybrids) while analyzing bonding interactions. 2.1. The LMTO Method. For a detailed description of the linear muffin-tin orbital (LMTO) method the reader is referred to the original literature29 as well as to some recent review articles.30-32 In a nutshell, the LMTO method stands for a linearized form of the Korringa-Kohn-Rostoker (KKR) method.33934 The basis sets of the LMTO method adjust dynamically to the respective potentials, in contrast to most other methods (based on simple fixed basis functions such as Slater orbitals, plane waves, or Gaussians), which are easy to handle in principle, but tough in practice owing to huge basis sets if treating d and f level materials. The method accounts for the potential from all the electrons and is applicable to materials composed of atoms from any part of the periodic table. Even though the LMTO method is not restricted to potentials of the muffin-tin (MT) type, the use of the atomic-spheres approximation (ASA), which uses muffin tin spheres blown up to overlapping and volume filling spheres, is most common. In addition, the potential within the atomic spheres, which represent spheridized WignerSeitz spheres, is spheridized. In loosely packed structures, empty spheres must be introduced to maintain space-filling. The LMTOASA method is extraordinarily efficient and yields excellent band structures and lattice constants for systems for which a homogeneous and dense sphere packing can be achieved. This method is used throughout this paper. There are also several full-potentialmodificationsof the LMTO method which allow the treatment of anisotropic atomic deformations and which take the full potential without shape approximations into account.3s-37 Within a specificrepresentation denoted by a,an LMTO basis is constructed as almost minimal set for a given MT potential. In the interstitial regions with flat potentials, the wave functions of the valence electrons are expanded into Hankel enuelope functions K& according to
(-V2
+ K’)K&(? - fi) = 0
for i #
fi, R’, ...
(23)
whereas in the core-like regions, i.e., within the MT spheres, one seeks numerical solutions of the radial SchrMinger equation
with
Concerning the kinetic energy K~ within the interstitial region, one may satisfactorily choose the simplest case ( ~ = 2 0) for the occupied levels. The energy dependent partial waves + R L ( t , i ) are expanded into a Taylor series up to linear order in the energy, which has
Dronskowski and Blhhl been proven to yield a sufficiently accurate description of the valence band region, which is typically about 1 Ry wide. This is the central idea of all linear methods in band theory.29 TO shorten the notation we use 14) instead of 4(i) and arrive at
(26) = I$RL) + I ~ J ; L ) ( ~ e/) + ~ ( -e ed12 The abbreviations 1 4 ~ = ~ )J ~ R L ( and ~ ~ )ldE) ) = (a/&)14RL(e))lr, stand for value and energy derivative of the energy dependent partial waves at tu. Finally, the envelopefunctions are replaced within MT spheres by linear combinations of the 4’s and &s so that the resulting basis orbital is continuous with value and derivative at the sphere surface. This procedure is also called augmentation. Within the ASA, Hamiltonian and overlap matrix can then be expressed by structure constants and potential parameters. Structure constants38 S w Ll express the boundary conditionsof the envelope function centered at site R with angular momentum L at the site R’ for partial waves with angular momentum L’ . They are completely potential independent and sensitive only to the atomic structure. The potential parameters, on the other hand, reflect the boundary conditions of the partial waves at the sphere radius. They are insensitive to the atomic structure but depend on theactual potentialand thesphereradius. The potential parameters describe the phase shift of the atomic potential. 2.2. Tight-Binding LMTO.The assignment of local bonding characteristics to atomic pairs is easiest if the calculation is performed with a basis set of short-range (atomic-like) basis functions. The LMTO approach offers the possibilities of using different kinds of basis sets, which are exactly of the same quality, because they span the same Hilbert space. There are, for example, conventional orbitals (a= 0 representation), derived from atomcentered Hankel functions as envelope functions, nearly orthogonal orbitals (a = y representation) and tight-binding (TB) orbitals, easily transformable into each other. For the generation of COHPs, we search for short-ranged, localized orbitals similar to those appearing in most semiempirical tight-binding methods such as EH theory. Most semiempirical TB methods work at the price of computational inaccuracy since information on the potentials and the spatial extent of the atomic wave functions to be specified requires experience. In contrast, the recently developed exact transformation of LMTO basis orbitals into a set of localized TB orbitals offers afirst principles access to a TB method, combined with the LMTO method’s internal computational accuracy.39 Therefore, TB-LMTO can be regarded as a link between empirical TB methods on the one side and the possibly more accurate first principles methods on the other. The basic idea is the following: The TB basis orbitals are constructed by superposing Hankel functions on different sites such that the long-range tails cancel each other. These superpositions then serve as the envelope functions that yield the tight binding orbitals after augmentation. This construction is analogous to the electrostatic screening of the long-range potential of a charge multipole. The localized envelope functions can be expressed as I+RL(e))
= IKO)”(l + asa) (27) where d a plays the role of a screening charge and (KO)” are the bare Hankel functions KORL(?) = Ir - R I’-’YL(i - &, with YL(P) being the spherical harmonics for angular momentum L. The so-called screening matrix a is chosen to be diagonal a = (YR+RRISLL, (28) and the screened structure constants Sa are computed from Dyson’s equation
so
sa= + Pasa (29) The optimum values for a that yield the most localized envelope
The Journal of Physical Chemistry, Vol. 97, No. 33, 1993 8621
Crystal Orbital Hamilton Populations functions and orbitals were originally found by trial and error;39 they are as = 0.3485 &, ap * 0.0530 Bp, and a d * 0.0107 = &, for s, p, and d orbitals, respectively. They turned out to be almost independent of the crystal structure if all orbitals on all atoms are screened in the same way. On the other hand, in some cases it is convenient to use sitedependent screening constants, especially in open structures in which there are both atoms and empty spheres. If one chooses not to screen on the empty sphere sites, the screening constants need to be adjusted for the reduced density of screening sites32 and then yield
-1.0
II
,
,
,
,
,
.4
.z aR,
(UR/W)*"'a/
(30) 0.
where U R is the screeningradius at the site R and w is the average WS radius, in combination with the structure constraint
-.Z
N -.4
3. Application There are three reasons for choosingelemental Si in the diamond structure as a test case: First, the bonding principle in this fundamental structure type is very simple. In fact, point group symmetry Td lets us safely assume that each atom carries sp3 hybridized orbitals, the optimum symmetry-adapted choice for tetrahedrally coordinatedmain-groupelements. Second, reliable thermodynamic data for the bonding are available. The experimental value for the cohesive energy of a-Si was reporteda to be 447 kJ/mol, which corresponds to approximately 224 kJ/mol for the S i S i bond.41 This value is in fair agreement with an average value of 196 kJ/mol, compiled from the entire variety of molecules containing S i S i single bonds, such as c a r b o s i l a n ~ . ~ ~ Third, Si in the diamond structure has always served as benchmark material for first principles computations. A chemically more interesting study dealing with metal-rich solid-state materials will be the subject of a forthcoming paper.43 A minimal basis set of short-range atom-centered TB-LMTOs was used throughout the investigation, i.e., with one s, three p, and five d orbitals on each Si atom. To maintain space filling, additional empty spheres (ES)were placed into the remaining tetrahedral structural voids, also equipped with TB-LMTOs of s-, p-, and d-like character. Thus, the primitivetrigonal cell used for the electronic structure calculation contains two Si atoms and two ES sites. If Si atoms and ES's are regarded as being equal, the packing corresponds to a body-centered-cubic(bcc) structure type. Keeping that in mind, there are two choices for the generation of tight-binding LMTOs: We speak of bcc-screening if the screeningconstants for Si and ES are both chosen to be the optimum values (see section 2.2). Applying diamond screening, the screening values for Es's ((Y$,d) are set to zero, and the screening values for Si atoms are computed via eqs 30 and 31, Le., they are a : = 0.4391,a:i = 0.1061,and : a = 0.0340. In both cases, the electronic energy was computed using the local density-functional theory, parametrized according to von Barth and Hedin.4 The integration in k space was performed with the help of an improved37 tetrahedron method45 using 47 inequivalent k points and 148 different tetrahedra. The selfconsistent potential parameters are the same as in the paper of Andersen et a1.& For taking a look at plotted bcc-screened TB orbitals we also refer the reader to the workof Andersenet al.,%inwhichextremely compact orbitals with an effective radius of about 1 .O-1.2times the nearest-neighbor distance were obtained. With diamondscreening, the orbitals on Si are a little more diffuse, having a considerable spillover into the interstitial regions. Let us focus on the S i S i bond with the help of a few COHP diagrams, shown in Figures 1 and 2,in which the (dimensionless)
I
1
I
.5
0.
-1.0
-.I
0
-.I .I
Energy ( R y l
Figure 1. COHP (dotted line) and integrated COHP (full line) of the Si-Si bond in a-Si. Depictedare the contributionsfrom (i) s-s combination (top), (ii) sp%$ combination (center), and (iii) all Si-centeredvalence orbitals (bottom), using a bcc-screening p r d u r e . Up to the Fermi energy (vertical dashed line) the band structure energy is lowered by ( i ) -0.034 Ry (-44 kJ/mol), (ii) 4 . 3 9 6 Ry (-520 kJ/mol), and (iii) -0.468 Ry (-614 W/mol).
COHP is drawn as a dotted line (axis on the left-hand side) and its energy integral, giving the sum over all bonding, nonbonding and antibonding interactions, as a solid line (axis on the righthand side, Ry h i t s ) . The Fermi level is shown by a dashed vertical line. In the discussion, a special emphasis will be put on the shape of the COHP curves, and one-electron bond energies will be given in kJ1m01.47 In Figure 1 one can see three COHP diagrams of the S i S i bond in a&, calculated using bcc-screened TB-LMTOs. At the top, only the s-s contribution has been picked out. The first bonding region is around -0.9 Ry, visible as a sharp negative COHP spike. The COHP integral would lead to a lowering of the band structure energy of roughly -150 kJ/mol if the Fermi level were lying at about -0.65Ry. But since there are antibonding levels (sharp positive COHP spike above -0.6 Ry), which are filled up to the real Fermi level, the integrated COHP gives a total s-s bonding effect of 4 4 kJ/mol for the bond. Another unoccupied region of s-s antibonding nature around 0.1 Ry can also be seen. The middle diagram shows the bonding contributions of a pair of sp3 hybrid orbitals. These contributions have been extracted from the complete s, p, d orbital basis set using an sp3 projection operator, thus retaining the full variational freedom of the basis set within the calculations. The spatial orientation of the sp3-sp3 orbital pair is such that the orbitals are centered on neighboring silicon atoms and point towards each other along the connecting [ 1 1 11 direction. One may easily verify the existence of the s bonding contributions between -0.9 and -0.6 Ry and of the p bonding contributions between -0.4 and -0.1 Ry. Antibonding levels appear beyond 0.1 Ry in the unoccupied band region. At this point we must stress a central observation: in the present case, the integrated COHP of the symmetry-adapted hybrid
8622 The Journal of Physical Chemistry, Vol. 97, No. 33, 1993
I
! .I
0.
-.1 -.a
t
I
Dronskowski and Bliichl local-density approximation. Its size depends somewhat on the parametrization used.5’ A similar observation has been made by Gunnarsson et al. while performing DF calculations on diatomic molecules.28 They reported an empirical formula relating the real bond energy to one third of the sum of the eigenvalues. In addition there is a non-DF (TBB) recursion method calculation of a-Si by Paxton,5* using semiempirical potential parameters that had been fitted to pseudopotential calculations.53 Paxton reports the one-electron bond energy to be -785 kJ/mol, overestimating it by a factor of 3.5. The clue, we find, lies in the setup of the energy-partitioning scheme,which forms the basis of the COHP technique and related methods. By looking at eq 22, we again stress the importance of the last two lines, which contain the atomic eigenvalue sum and on-site COHPs as well as the change of the energy terms resulting from electrostatic interaction and exchange and correlation. These must be taken into account to yield reliable values for the cohesive (and bond) energy. Equation 15 shows how the on-site COHP terms enter the calculation of the band structure energy. Using bcc-screened TB orbitals, one obtains the following COHP energy integrals at the Fermi level: on-site Si, -0.232 Ry (-305 kJ/mol); on-site ES, 0.122 Ry (160kJ/mol); off-site Si-Si, -0.468 Ry (-614 kJ/mol), depicted in Figure 1, bottom panel; off-site Si-ES, -0.107 Ry (-140 kJ/mol); off-site ES-ES, -0.023 Ry (-30kJ/mol). Since a primitive trigonal cell contains two Si atoms, two ES sites,four S i S i contacts, eight Si-ES contacts, and four ES-ES contacts, we arrive at
om
-1
-4
I
-1.0
-.S
0
.5
Energy ( R y l
Figure 2. Same as in Figure 1, but using diamond-screenedTB orbitals. Up to the Fermi energy the band structure energy is lowered by (i) -0.095 Ry (-125 kJ/mol), (if) -1.002 Ry (-1315 W/mol), and (iii) -1.318 Ry (-1731 kJ/mol).
orbital is minimized for the real electron count, in contrast to the s-s contribution shown before. So, our general observation is that the best choice for the orbitals is the one that minimizes the one-electron bond energy. Also, the projector’s arbitrary fixed weighting between s and p orbitals, which one may assume to be considerably greater than 1-3 in the valence bands48 turns out to be surprisingly well chosen for describing the bonding. The COHP diagram at the bottom represents the contribution of the complete orbital set (s,p, and d orbitals) centered on the Si atoms. Regarding the shape, there is no strong difference between this diagram and the one shown before. d orbital participation leads to no qualitative difference, and the main contribution indeed comes from s and p functions. The energy integral of the COHP at the Fermi level is -614 kJ/mol, only 15% higher than without d orbital participation. Because of our previous use of the sp3projection operator in the COHP diagram shown before, this 15% difference is proven to result only from the additional d orbitals, and not from a difference in the setup of the basis set. The Si-Si one-particle bond energy, obtained as sum of the off-site COHPs (including only Si-centered orbitals) is with -0.468 Ry (-614 kJ/mol) substantially larger (by a factor of 2.7!) than the experimental bond energy, which is -0.171 Ry (-224 kJ/mol). This discrepancy can be only partly attributed to a failure of the density functional approximation, which overestimates the cohesive energy by only 13%. The theoretical cohesive energy of a-Si is obtained by subtracting ( i ) the valence total energy49 of atomic Si (-7.735 Ry) plus the spin polarization energy50 of atomic Si (-0.048 Ry) from (ii) the calculated total valence energy of crystalline a-Si (-8.166 Ry), using the same scheme for the local-density approximation. This leads to a cohesive energy of -0.383 Ry (-503 kJ/mol), corresponding to a theoretical S i s i bond energy of about -251 kJ/mol, larger by 13% than the experimental one. This slight overbonding effect is indeed due to the failure of the
‘(2 X -0.232
2
+ 2 X 0.122 + 4 X -0.468 + 8 X -0.107 + 4 X -0.023) Ry = -1.520 Ry
per atom, which is 0.262 Ry (15%) too small (in absolute value) compared to the calculated band structure energy (-1.782 Ry). This discrepancy is due to the neglect of the small second-nearestneighbor COHPs. Thus, we can trace the discrepancy between the bond energy derived from COHPs between Si-centered orbitals and the bandstructure energy per bond to the following three points. First, the on-site “crystal field” terms contribute strongly to the band structure energy. They describe the energy required to contract the atomic wave functions into localized tight-binding orbitals, the promotion of electrons from the Si-centereds orbitals to Si-centeredp orbitals consistent with sp3bond orbitals and the partial electron transfer from Si-centered orbitals to orbitals on empty spheres, which gives rise to an artificial ionic contribution to the total energy. Second, since the same TB screening constants for Si and ES sites are used, there is considerable Si-ES “bonding”, which is of course unphysical, and reflects that this basis set is not very well chosen for a chemical discussion. Third, the neglect of COHPs between second-nearestneighbors results in a moderate overestimateof the (negative) band structure energy. To address thesepoints,we recalculatedthe electronic structure, now with diamond-screened TB-LMTOs for Si and ES, fixing all ES-a’s to zero and enlarging the Si-a’s, as described above. The total energy (-8.166 Ry/atom), the band structure energy (-1.779 Ry/atom instead of -1.782 Ry/atom), and the Fermi energy (-0.043Ry instead of -0.048 Ry) are practically unaffected by the change in the basis set-as expected. Now we turn to the influence of the new basis set on the energy partitioning. A first look at corresponding COHP diagrams in Figure 2 reveals that there are only small differences in shape between diamond and bcc screening. The s-s contribution (top) to the bonding interaction is essentially the same as the one with bcc-screened orbitals, with a smaller antibonding peak in the
Crystal Orbital Hamilton Populations virtual bands. While the combination of sp3 hybrids (middle) also shows the antibonding peak at 0.1 Ry to be slightly smaller in magnitude, the bonding interaction around -0.55 Ry is more strongly emphasized when using diamond-screened orbitals. In comparing the total S i S i interaction (bottom) with the one from bcc screening, we again note that bonding levels around -0.55 Ry are stronger and that the first antibonding levels above 0.1 Ry now appear as bonding levels. The latter effect is counterintuitive, as it forces the COHP energy integral to become even more negative beyond the Fermi level. The explanation must again be sought in the interaction of on-siteand off-site COHPs. Thed orbitals havevery high energy, reflected in a strongly positive on-site contribution to the corresponding COHP if d orbitals are admixed to the wave functions. The mixing with d orbitals must be explained in two steps. First, electrons are promoted into d orbitals, which is energeticallyunfavorable. Second, this contribution is more than offset by the "bonding" hybridization, because inclusion of d orbitals must lower the energy. For a discussion of binding, however, it would be natural only to include s and p orbitals, which again emphasizes the need of a truly minimal basis set for the use of COHP de~criptors.5~ Concerning the absolute size of the one-electron S i s i bond energies, however, we find even larger values of -125, -1 3 15, and -173 1 kJ/mol for s-s contribution, sp3-sp3 contribution, and total orbital interaction, respectively, which are about 2.7 times larger in magnitude than the one-electron bond energies computed with bcc-screened TB-LMTOs and more than 7 times larger than the experimental bond energy. Since the band structure energy has remained almost constant, the reason for the absolute size of the off-site COHPs is again found in the interaction of the on-site atomic and off-site terms. For diamond screening, the COHP energy integrals at the Fermi level are as follows: on-site Si, 0.841 Ry (1 104 kJ/mol); on-site ES, 0.0 16Ry (21 kJ/mol); off-site SiSi,-l.3 18 Ry (-1 73 1 kJ/mol), see Figure 2, bottom panel; off-site Si-ES, 0.003 Ry (4 kJ/mol); off-site ES-ES, 0.002 Ry (3 kJ/mol). For the band structure energy (in real space) we have 1
-(2 2
X 0.841
+ 2 X 0.016 + 4 X -1.318 + 8 X 0.003 +
4 X 0.002) Ry = -1.763 Ry per atom, only 0.016 Ry (0.9%) too small compared to the calculated band structure energy (-1.779 Ry). There are three observations: First, energy participation from empty spheres, both on-site and off-site, has practically been eliminated. This implies that for diamond screening Si-centered orbitals alone will give a good description of the band structure, and empty sphere-centered orbitals could be excluded from the entire calculation. The best choice would be a truly minimal basis set that is still localized. Orbitals that are even closer to this ideal than the TB orbitals used here can be constructed using downfolding theory.55 These more sophisticated techniques for generating LMTO orbitals exploit theknowledgeof the potential when setting upthescreening matrix. A different method, which also arrives at a minimal tight-binding basis set, has been developed by Galli and Parrinello56 in the context of the Car-Parrinello method. Here the orbitals are completely variational and localized by a gauge potential that breaks thesymmetry ofthedensity functional theory against unitary transformation of the wave functions. Second, the restriction to nearest-neighbor interactions is justified by the negligible error in the band structure energy. Third, the "crystal field" on-siteCOHPs of Si have dramatically increased and have become positive, counterbalanced by off-site S i S i bonding COHPs in order to give the same band structure energy. This is the price paid for the greater diffuseness of the orbitals located on Si.
The Journal of Physical Chemistry, Vol. 97, No. 33, 1993 8623
To summarize,we have demonstrated how to visualizechemical bonding in solidsfrom density-functionalTB-LMTO calculations. The bonding principle can easily be extracted from the shapes of crystal orbital Hamilton populations with respect to orbital participation over the energy range specified. Thus, a COOPlike method can be set up within the framework of ajhtprinciples technique. On the one hand, the absolute values of the COHPs depend strongly on the choice of the tight-binding orbitals. This contribution is counterbalanced by "crystal field" terms, represented by on-site COHP contributions. On the other hand, the shapes of the COHPs are fairly insensitivewith respect to changing basis sets. Tocompare both shapes and absolutevalues of COHPs between different materials, the necessity to use identicalscreening constants is obvious. For loosely packed structures involving empty spheres, one has to seek for a compromise between ( i ) an appropriate magnitude of the TB orbitals and (ii) a reasonable description of the interaction between all sites in the structure. Appendix. COHP Projections via Euler Rotations While performing a COHP bond analysis, it may sometimes be interesting to know which particular atomic orbital has the strongest impact within the complete set of orbital interactions. In other cases, it will already be obvious what kind of symmetryadapted linear combination of orbitals needs to be constructed in order to describe a chemical bond satisfactorily. Perhaps the best-known examplefor the latter case is given by Pauling's valence bond treatment of chemical bonding in hydrocarbons where sp3, sp2, and sp hybrids represent the optimum combination~.~.~7 A general formalism, based on a maximum overlap criterion, was worked out by Murrell.58 In any case, the extraction of a particular orbital contribution can be performed by using a projection operator technique. Let us assume that there is an original basis set IXRL) (like s, p x ,p,,, ,..) which shall be transformed into a set of symmetry-adapted basis functions I ~ R M ) (like sp3,d2sp3,... ), according to
and T L Ma transformationvector ~~ whose entries consist of mixing coefficients from the original to the symmetry-adapted basis set. If I ~ R M ) is a complete set, the wave function reads
The COHP between two neighboring symmetry-adapted orbitals is written as coHP,~p,,+,*= ( ' ~ R MI ~ I [ P R M ' ) ~ ~ U * R M , ~ U R , M , -, ~E,)~ E= I
64.5)
8624 The Journal of Physical Chemistry, Vol. 97, No. 33, 1993
Note that Nand N(e) undergo separate transformations that are reciprocal to each other. Furthermore, there may be. the need to align hybrid orbitals to artificial directions in real space. To do so, one can easily rotate these hybrids via Euler rotations, principally equivalent with artificial axes transformations. A complete description involves the Eulerian angles a,j3,and y and the corresponding rotations P(a,B,y), equivalent to unitary transformations of functions 4~ 1 4lm. The transformation matrices do not allow combinations between functions with different ]Is, but only combinations from functions of the same ]Is with different m's: P(a,B,T)d,m = C4/mo'(a,B>T)dm
('4.7)
m'
The general matrix59 reads D'(a,B,v)mrm =
For a minimal basis set, consisting of one s, three p, and five d functions on each site, the transformation matrix would have 9 X 9 elements, and it would be blocked into one 1 X 1, one 3 X 3, and one 5 X 5 submatrix. A directed valence orbital, built by an LCAO procedure, is rotated by changing the elements of the column vector of the expansion coefficients: m
(A.11) If one uses stepwiseEuler transformationsto generate a continuous set of slightly varying projectors, one could easily scan through a bond in real space and thereby investigate its spatial extent with COHP diagrams.
Acknowledgment. We are grateful to Prof. Ole Krogh Andersen for his continuous support as well as for initiating the project. Also, we would like to express our thanks to Dr. Ove Jepsen for his help with orbital transformations. It is a pleasure to thank Dr. Uwe SchiSnberger for interesting discussions as well as Prof. Roald Hoffmann for providing us with valuable remarks on overlap population based bonding desriptors.
References and Notes (1) Mulliken, R. S. J. Chem. Phys. 1955,23, 1833. (2) Davidson, E. R. J. Chem. Phys. 1%7,46, 3320. (3) Roby, K. R. Mol. Phys. 1974,27, 81. (4) Heinzmann, R.; Ahlrichs, R. Theor. Chim. Acta 1976, 42, 33. ( 5 ) Ehrhardt, C.; Ahlrichs. R. Theor. Chim. Acta 1985,68,231.
Dronskowski and Bliichl (6) Wolfsberg, M.; Helmholz, L. J. Chem. Phys. 1952,20, 837. (7) Hoffmann, R. Angew. Chem. 1987,99,811; Angew. Chem., Int. Ed. Engl. 1987,26,846. (8) Pauling, L. The Nature of the Chemical Bond, 3rd ed.; Cornell University Press: Ithaca, NY, 1960. (9) Brown, I. D. In O'Keeffe, M., Navrotsky, A., Eds.; Structure and Bonding in Crystals, Academic Press: New York, 1981; Vol. 11. (10) Brown, I. D.; Altermatt, D. Acta Crystallogr. 1985,B41, 244. (11) Hoffmann, R.; Lipscomb, W. N. J. Chem. Phys. 1962,36, 2179. (12) Hoffmann, R. J. Chem. Phys. 1963,39,1397. (13) Hughbanks, T.; Hoffmann, R. J . Am. Chem. Soc. 1983,105,3528. (14) Hoffmann, R. Solids and Surfaces: A Chemist's View ofBonding in Extended Structures; VCH: Weinheim, New York, 1988. (15) Hoffmann, R. Personal communication, 1991. (16) Sutton, A. P.; Finnis, M. W.; Pettifor, D. G.; Ohta, Y . J . Phys. C 1988,21, 35. (17) Hall, G. G. Proc. R. Soc. London 1952,A213, 113. (18) Manning, P. P. Proc. R. SOC.London 1955,A230, 424. (19) Ruedenberg, K. Rev. Mod. Phys. 1962,34, 326. (20) Coulson, C. A. Proc. R . Soc. London 1939,A169, 413. (21) Harris, J. Phys. Rev. E 1985,31, 1770. (22) Foulkes, M. Ph.D. Thesis University of Cambridge, UK, 1987. (23) Heine, V. Solid State Phys. 1980,35, 1. (24) Jones, R.; Lewis, M. W. Philos. Mag. B 1984,49,95. (25) Ohta, Y.; Finnis, M. W.; Pettifor, D. G.; Sutton, A. P. J. Phys. F 1987,I?, L273. (26) Hohenberg, P.; Kohn, W. Phys. Rev. 1964,136,B 864. (27) Kohn, W.; Sham, L. J. Phys. Reu. 1965,140, A 1133. (28) Gunnarsson, 0.;Harris, J.; Jones, R.0. J. Chem. Phys. 1977,67, 3970. (29) Andersen, 0. K. Phys. Rev. B 1975,12,3060. (30) Andersen, 0. K.; Jepsen, 0.;GIBtzel, D. In Bassani, F., et al., Eds.; Highlights of Condensed-Matter Theory; North-Holland New York, 1985. (31) Skriver, H. L. The LMTO Method; Springer: Berlin, 1984. 132) Andersen. 0.K.: Jemen. O.:Sob. M. In Yussouff. M.. Ed.:Electronic B a d Structure and its Ap&tiohs; Springer: Berlin,'1986. (33) Korringa, J. Physica 1947,13,392. (34) Kohn, W.; Rostoker, N. Phys. Rev. 1954,94, 1111. (35) Weyrich, K. H. Phys. Rev.B 1988,37, 10269. (36) Methfessel, M. Phys. Rev. B 1988,38, 1537. (37) BlBchl, P. Ph.D. Thesis University of Stuttgart, FRG, 1989. , . (38) They should not be confused with the overlap integral S (39) Andersen, 0. K.; Jepsen, 0. Phys. Rev. Lett. 1984,53, 2571. (40) McMahan, A. K. Phys. Rev. B 1984,30,5385. (41) In the diamond structure each Si atom has four nearest Si neighbors, and each Si atom is involved in four Si-Si bonding interactions. However, since every Si atom contributes exactly 50% to each bond, there are two bonds per Si atom. (42) Fritz, G.; Matern, E. Carbosilanes. Syntheses and Reactions; Springer: Heidelberg, 1984. (43) Dronskowski, R.; Andersen, 0. K., manuscript in preparation. (44) von Barth, U.; Hedin, L. J. Phys. C 1972,5, 1629. (45) Jepsen, 0.;Andersen, 0. K. Solid Stare Commun. 1971,9, 1763. (46) Andersen, 0. K.; Pawlowska, Z.; Jepsen, 0. Phys. Rev. B 1986,34, 5253. (47) One Rydberg unit equals approximately 1313 kJ/mol. (48) We would like to thank a perceptive reviewer for spotting that the s contribution to the occupied valence levels may be considerably greater than 25%. (49) Desclaux, J. P. Hartree-Fock- Dirac-Slater Atomic Program;CEA: Paris, 1969, modified version, 1970. (50) Jepsen, 0. Personal communication, 1990. (51) Jones, R. 0.;Gunnarsson, 0. Rev. Mod. Phys. 1989,61,689. (52) Paxton, A. T. Philos. Mag. B 1988,58,603. (53) Gou-Xin, Q.;Chadi, D. J. Phys. Rev. B 1987,35, 1288. (54) Consequently, a further corrective step (not performed here) would be to downf0ld5~the d orbitals and leave only s andp orbitals for performing the self-consistent calculations. ( 5 5 ) Lambrecht, W. R. L.; Andersen, 0.K. Phys. Rev. B 1986,34,2439. (56) Galli, G.; Parrinello, M. Phys. Rev. Lett. 1992,69,3547. (57) Pauling, L. J. Am. Chem. SOC.1931,53, 1367. (58) Murrell, J. N. J. Chem. Phys. 1960,32,767. (59) Tinkham, M. Group Theory and Quantum Mechanics; McGrawHill: New York, 1964.
