Effect of on-site Coulomb repulsion term U on the band-gap states of the reduced rutile (110) TiO2 surface

Effect of on-site Coulomb repulsion term Uon the band-gap states of the reduced rutile (110) TiO2surface Carmen J. Calzado,*Norge Cruz Hernández, and Javier Fdez. Sanz Departamento de Química Física, Facultad de Química, Universidad de Sevilla, E-41012 Sevilla, Spain 共Received 26 October 2006; revised manuscript received 31 July 2007; published 15 January 2008兲 We present a study concerning the effect of the on site d-dCoulomb interaction energy Uon the band-gap states of nonstoichiometric rutile 共110兲TiO2surface. As well known, the excess electrons resulting from the formation of oxygen vacancies localize on the Ti 3dorbitals forming band-gap states. Local density approximation 共LDA兲does not give a correct description of these band-gap states, either with or without gradient corrections. The failure of LDA is often attributed to an inadequate treatment of electron correlation in systems with localized orbitals and is commonly corrected with an empirical local Coulomb repulsion term, i.e., the LDA+Umethod. This study provides a completely general strategy to estimate the Uvalue in this kind of systems, illustrated here for reduced 共110兲TiO2surface, well characterized from experiments. From ab initio embedded cluster configuration interaction calculations, combined with the effective Hamiltonian theory, a value of Uof 5.5±0.5 eV is obtained, in good agreement with those reported for this system from x-ray photoemission spectroscopy experiments 共U=4.5±0.5 eV兲. It is observed that when the ab initio estimate of U is injected into the periodic LDA+Ucalculations, a correct description of the gap states is obtained from the periodic LDA+Ucalculations. Additionally, the results indicate that the position of these states on the band gap strongly depends on the level at which lattice relaxation is taken into account, with significant differences between the density of states curves at the LDA+Ulevel obtained using the optimal generalized gradient approximation or LDA+Ugeometries. These results suggest that this combined strategy could be a useful tool for those systems where electron correlation plays a key role, and no experimental data are available for the on-site Coulomb repulsion. DOI: 10.1103/PhysRevB.77.045118 PACS number共s兲: 71.15.Mb, 31.15.V⫺, 31.15.A⫺ I. INTRODUCTION Titanium dioxide is probably one of the most studied compounds in material science, due to its versatility and large range of applications from white pigments to catalyst support 共Ref. 1and references therein兲. Together with the improvement of the experimental techniques, the knowledge of this material has taken benefit of the development of sophisticated theoretical approaches, which have helped in understanding the structure and reactivity of TiO2surfaces. Since most of the real applications deal with substoichiometric TiO2, there has been an increasing effort in characterizing both the geometrical and electronic structures of reduced TiO2, specially the rutile 共110兲surface which is the most stable.2 Reduction can be induced by removal of oxygen atoms as well as by deposition of alkali metal atoms. Both processes produce Ti+3 ions and localized band-gap states formed by Ti 3dorbitals.1Different spectroscopic techniques support this description such as ultraviolet photoemission spectra 共UPS兲,2–4resonant photoemission,5–7electron energy loss spectroscopy 共EELS兲,8,9or x-ray photoelectron spectroscopy.6Indeed all these features disappear after adsorption of molecular oxygen at room temperature. The defect states lie in the upper half of the gap, 0.75–1.18 eV below the conduction-band minimum,10 although their position depends on the vacancy concentration, moving toward the conduction band as the defect concentration increases. Also broad peaks have been observed in the region of 450–750 nm, attributed to defect levels due to oxygen vacancies, placed at about 2 eV below the conduction-band edge.11–14 However, Henderson et al.15 have recently proposed an alternative assignment of these peaks based on Franck-Condon arguments. They considered that these peaks can correspond to excitations from electron trap states located at about 1 eV below the conduction-band 共CB兲edge to Ti 3d-derived levels placed at 1 eV above the CB edge, where See and Bartynski16 have situated the first maximum in Ti 3ddensity of states. In the past years these band-gap states have renewed the interest in this material since they enhance the photocatalytic activity observed in titanium oxide on several processes such as the degradation of organic pollutants in water 共Ref. 17 and references therein兲or the solar energy conversion.18,19 The presence of these defect states reduces the energy needed for photoexcitation with respect to the stoichiometric material, making possible the use of visible light, which opens the way to low-cost applications. There have been previous theoretical studies of the reduced 共110兲surface of TiO2focusing on the band-gap states. Most of them predict that Ti 3d-derived conduction-band states become occupied after reduction, but in general they fail in the position of these new states. Among those using periodical approaches we can mention the work by Ramamoorthy et al.,20 in which using plane-wave pseudopotential techniques based on local density functional approximation 共LDA兲they found Ti 3d-like surface states but lying in the conduction band. The same result has been obtained from spin polarized calculation performed by means of the fullpotential linear muffin-tin orbital method by Paxton and Thiên-Nga.21 Lindan et al.22 using gradient corrections to the local density approximation and spin polarized calculations PHYSICAL REVIEW B 77, 045118 共2008兲 1098-0121/2008/77共4兲/045118共10兲©2008 The American Physical Society045118-1 show results which are qualitatively correct but the band-gap states are too close to the conduction-band edge. As far as we know, the best theoretical description of these states has been recently provided by Di Valentin et al.,23 achieved only when the lattice distortion induced by the extra electrons is accounted for with the Becke 3-parameter Lee-Yang-Parr 共B3LYP兲functional. The failure of pure density functional methods 关LDA or generalized gradient approximation 共GGA兲兴 could be associated with an inadequate description of the strong Coulomb interaction between 3delectrons localized on Ti atoms. This could be corrected by introducing an additional term U, which takes into account the repulsion between two electrons placed on the same 3dorbital. This approach, named LDA +U共or GGA+U兲method, has provided descriptions in a reasonable agreement with experiments for systems where LDA solutions where systematically wrong, such as the EELS spectra of nickel oxide.24 However, the main limitation of this method is that Uis conceived as an empirical term, which is varied until convergence with the experimental results. Many examples of the application of this pragmatic approach can be found in literature. A very recent one is the work by Loschen et al.25 dedicated to the study of the electronic structures of CeO2 and Ce2O3. This procedure is not comfortable, especially when experimental Uvalues are not available, and it is not possible to check whether the Uvalue giving the correct answer is or not physically meaningful. An alternative to this trial-and-error procedure is the method proposed by Cococcioni and de Gironcoli26 to determine Uin a self-consistent way. However, while quite promising for some systems, its success is not universal. For instance, the application of this method to the study of the electronic structure of CeO2provides Uvalues which are underestimated with respect to the optimal Uvalue,27 while overestimated Uvalues have been obtained for iron heme complexes.28 In this context, any independent evaluation of Uprevious to the LDA+Ucalculations would be desirable, and this is the aim of the present work. We report on a general strategy combining ab initio embedded cluster calculations and effective Hamiltonian theory which yields estimates of the on-site Coulomb repulsion term. This value is later injected on the periodic LDA+U calculations and provides, as will be discussed below, a correct description of the band-gap states of reduced 共110兲TiO2 surface as well as their dependence on the vacancy concentration. It is worth to note that TiO2is one of the simplest of the systems where electron correlation could play a crucial role on the description of the electronic structure. It has been chosen as an illustrative example since there is a vast amount of available experimental data which permit us to check the validity of our procedure. The whole strategy, however, is completely general and can be extended to systems containing more than one delectron per site. The paper is organized as follows. Section II presents the methodology to obtain Ufrom ab initio embedded cluster calculations, as well as the computational details of the periodic LDA+Ucalculations. Section III discusses the results and Sec. IV contains the main conclusions from this work. II. COMPUTATIONAL STRATEGY A. Ab initio embedded cluster calculations 1. Evaluation of U from ab initio embedded cluster calculations The on-site Coulomb repulsion Ufor reduced TiO2has been determined from ab initio embedded cluster calculations. When TiO2is reduced by creating oxygen vacancies or by the deposition of alkali metals, the electron excess localizes on Ti 3dorbitals. Let us consider a fragment of the reduced TiO2system, containing two neighbor Ti centers 共Fig. 1兲each of them with one 3delectron. Let us call aand bthe 3doccupied orbitals placed at TiAand TiBcenters, respectively. Three different situations can be conceived for these two electrons in two orbitals: 共1兲one electron per orbital, with opposite spins, represented by the determinants 兩ab ¯ 兩and 兩ba ¯ 兩; 共2兲one electron per orbital, with parallel spins, that is, the determinants 兩ab兩and 兩a ¯ b ¯ 兩; 共3兲both electrons on the same orbital, 兩aa ¯ 兩and 兩bb ¯ 兩. In the valence-bond 共VB兲framework the two first arrangements correspond to neutral determinants, the latter one to ionic ones. The combination of these determinants gives four different configurations: 共a兲a neutral singlet state of gsymmetry, Sg N, 兩Sg N典=共兩ab ¯ 兩+兩ba ¯ 兩兲/冑2; 共1兲 共b兲a neutral triplet state of usymmetry, Tu N, 兩Tu N典=共兩ab ¯ 兩−兩ba ¯ 兩兲/冑2; 共2兲 共c兲an ionic singlet state of gsymmetry, Sg I, 兩Sg I典=共兩aa ¯ 兩+兩bb ¯ 兩兲/冑2; 共3兲 共d兲an ionic singlet state of usymmetry, Su I, 兩Su I典=共兩aa ¯ 兩−兩bb ¯ 兩兲/冑2. 共4兲 Ucan be defined as the energy associated with the pairing of two electrons on the same 3dorbital, [110] [ 001 ] [110 ] FIG. 1. 共Color online兲Ti2O10 fragment used in ab initio embedded cluster calculations. Gray and large dark spheres correspond to Ti and O atoms, respectively. Ti total ion potentials 共TIPs兲are represented as small dark spheres. CALZADO, HERNÁNDEZ, AND SANZ PHYSICAL REVIEW B 77, 045118 共2008兲 045118-2 U=E兩aa ¯ 兩−E兩ab ¯ 兩共5兲 and can be written in terms of the energy difference between the ionic and neutral singlet configurations: U=ESg I−ESg N.共6兲 So, if we can evaluate the energies of the neutral and ionic states, Ucan be extracted from their energy difference. However, the energies of these two states cannot be directly obtained from standard quantum chemistry calculations. Since both states belong to the same irreducible symmetry representation, any ab initio quantum chemistry calculation will provide solutions where both states are mixed, with different cI/cNratios: 1⌽g 1=cN兩Sg N典+cI兩Sg I典, 1⌽g 2=cN ⬘兩Sg N典+cI ⬘兩Sg I典.共7兲 Then we need to design a procedure to extract this information from the eigenvalues and eigenvectors of the electronic Hamiltonian. This procedure makes use of the effective Hamiltonian theory, and it has been previously used in the study of magnetic systems, as well as the evaluation of hopping integrals in mixed-valence systems.29–41 A detailed description of the method can be found in Ref. 41, and a summary enlightening the most striking points will be provided here. Let us consider a model space Sspanned by the four VB determinants or by their four combinations 关Eqs. 共1兲–共4兲兴. Its projector is P ˆS=兺 i苸S 兩 ␾ i典具 ␾ i兩.共8兲 From ab initio embedded cluster calculations we can obtain napproximated solutions to the exact Hamiltonian, which hereafter will be considered as exact. These solutions have the largest components in the model space S,兵⌽k, k=1,n其, with energies 兵Ek其. They constitute the target space S⬘. Now we define an effective Hamiltonian in Ssuch as its neigenvalues are exact 共then equal to 兵Ek其兲 and its eigenvectors are projections of the corresponding exact eigenvectors in the model space. This is the definition of Bloch effective Hamiltonian42: H ˆeff Bloch兩P ˆS⌽k典=Ek兩P ˆS⌽k典.共9兲 This basic equation leads to the spectral definition of the Bloch effective Hamiltonian42: H ˆeff Bloch =兺 k=1,n 兩P ˆS⌽k典Ek具P ˆS⌽k ⬜兩,共10兲 where 兩P ˆS⌽k ⬜典represents the biorthogonal vector associated with 兩P ˆS⌽k典. Actually the projections of the 共orthogonal兲 eigenvectors of H ˆonto the model space have in general no reason to be orthogonal. They define an overlap matrix Sas follows: Sij =具P ˆS⌽i兩P ˆS⌽j典共11兲 and 兩P ˆS⌽k ⬜典=S−1兩p ˆS⌽k典.共12兲 The Bloch effective Hamiltonian is non-Hermitian. Its n2 matrix elements are defined from the n2conditions imposed by Eq. 共10兲. In the two-electron/two-orbital problem, nis equal to 4, and the Bloch effective Hamiltonian takes the following form, where the triplet state is the energy origin: 兩Sg N典 兩Sg I典 兩Tu N典 兩Su I典 冢 2Kab 2tab ⬘00 2tab U+2Kab 00 000 0 000 U+2共Kab −Kab ⬘兲 冣 ,共13兲 where Kab and Kab ⬘are the exchange integrals: Kab =具ab ¯ 兩Heff Bloch兩ba ¯ 典, Kab ⬘=具aa ¯ 兩Heff Bloch兩bb ¯ 典,共14兲 and tab and tab ⬘the hopping integrals between centers aand b: tab =具ab ¯ 兩Heff Bloch兩aa ¯ 典, tab ⬘=具aa ¯ 兩Heff Bloch兩ab ¯ 典.共15兲 The differences between Kab and Kab ⬘and tab and tab ⬘result from the nonhermiticity of the Bloch Hamiltonian. The hermiticity can be restored by means of the des Cloizeaux43 or Gram-Schmidt44 procedures, based on the orthogonalization of the projections of the eigenvectors on the model space. It is worth to notice that this procedure requires the knowledge of four states. The neutral ones are in most of the cases the lowest states in their symmetry, but the ionic ones are excited states. For those cases where the excited ionic states cannot be unambiguously determined or they are strongly contaminated by other configurations such as ligand-to-metal charge transfer excitations, the use of an intermediate effective Hamiltonian,45 instead of the procedure described above, is more pertinent. The intermediate Hamiltonian is built on a four dimensional model space but it is only asked to reproduce the energies of the neutral singlet and triplet states and projections of the singlet state onto the model space 共cIand cNcoefficients兲. It is possible to demonstrate that the on-site Coulomb repulsion on this scheme is equal to41 Uint =2tab cI 2−cN 2 cIcN ,共16兲 where tab keeps the valence value. This approach is reasonable since tab is only moderatly dependent on the electronic correlation as it has been shown in the past for several authors in quite different systems.33,37–39,41,46–51 EFFECT OF ON-SITE COULOMB REPULSION TERM U…PHYSICAL REVIEW B 77, 045118 共2008兲 045118-3 2. Physical model A bulk fragment of formula Ti2O10 has been chosen to evaluate U. It is composed of two Ti atoms and the nearest neighbor oxygen atoms. The remainder of the crystal is represented by means of a set of point charges, which approximately reproduce the Madelung potential of the infinite crystal. In order to avoid an artificial polarization of the electrons of the cluster, the Ti atoms in the immediate surroundings of the cluster are represented by total ion potentials 共TIPs兲共Fig. 1兲. The geometries of all the atoms are obtained from a previous periodic density functional theory 共DFT兲calculation carried on the bulk of TiO2rutile,52 which correctly reproduces the experimental lattice parameters. Notice that this fragment contains two extra electrons, resulting from the removal of an oxygen atom on the surface or any other chemical process such as the deposition of a donor atom or molecule such as an alkali metal. In order to take into account the effect on the Uvalue of the geometrical relaxations induced by these extra electrons, partial geometry optimizations have also been carried out. Details of them are provided on next section. Three sets of basis functions have been employed. Oxygen atoms are represented by means of an ANO-L-type basis functions with a contraction 共4s3p1d兲.53 For Ti atoms, the core electrons are replaced by pseudopotentials. We deal with 12 valence electrons when the pseudopotential and basis set proposed by Hay and Wadt54 are used 共basis 1兲and 10 valence electrons if Barandiarán and Seijo’s set is employed,55 with a 共3s3p4d兲contraction scheme 共basis 2兲. Finally, basis 3 consists of Hay and Wadt representation for Ti atoms, and ANO-L-type basis functions with an additional dshell for oxygen atoms 关contraction 共4s3p2d兲兴. Regarding the embedding, several choices are possible for the values of the point charges. There is a general consensus with respect to the partial covalent character of TiO2at difference with other systems extensively used in ab initio embedded cluster calculations such as MgO.56,57 For this reason a completely ionic representation for this system could be not adequate.58,59 In order to check the impact of the environment representation on the Uvalue, we have used two different models for the embedding: 共a兲Model 1: a fully ionic representation, with charges 兵+4其and 兵−2其for Ti and O atoms, respectively, 共model 1兲; 共b兲Model 2: a completely covalent representation, with charges 兵+2其and 兵−1其,共model 2兲. The bare cluster, with only TIPs in the neighborhood, gives meaningless results, as expected, due to a significant contamination of the ionic states, which impedes us a univocal characterization of the ionic states. On the other hand, as will be shown below, the effect on the Uvalue of the representation used to model the crystal is practically negligible. 3. Approximation to the exact N-electron wave functions Highly correlated N-electron wave functions are obtained from state-of-the-art quantum chemistry techniques. The onsite repulsion is evaluated from the energies and wave functions obtained from three CI spaces: 共a兲The bare valence complete active space 共CAS兲, that is a CASCI space. 共b兲The CAS+Sspace, including all the single excitations on the top of the active space. This expansion includes three types of determinants: 共1兲1h, where an electron moves from an inactive occupied orbital to an active orbital; 共2兲1p, where an electron moves from an active to an inactive virtual orbital; and 共3兲1h+1p, where an electron moves from an occupied to a virtual orbital. It is worth to note that these excitations can be coupled with simultaneous single excitations inside the active space. 共c兲The DDCI space, acronym of difference dedicated CI space,60 where all the single and double excitations on the top of the CAS are included in the CI calculations, except those which do not involve an active orbital. Thus, a double excitation from two inactive occupied orbitals to two virtual ones is not included in the calculation. While these excluded excitations, which are the most numerous, contain the most part of the electronic correlation energy, they do not contribute to the energy difference among the states involved in the evaluation of U. So, the use of the DDCI approach considerably reduces the computational cost with respect to single and double configuration interaction 共SDCI兲calculations, but it assures the introduction of the component of the correlation energy which plays a differential role in the stabilization of the states under study. This strategy has been extensively used in the recent past, especially in the evaluation of magnetic coupling constants as well as the determination of hopping integrals in numerous systems.29–41,46,48,49,61,62 The whole procedure can be summarized as follows. 共1兲A fragment of the bulk structure of TiO2has been chosen, containing two neighbor Ti centers, and the ten oxygen atoms around them, as shown in Fig. 1. This cluster is embedded in a set of point charges which model the effect of the Madelung potential of the infinite crystal on the centers of the cluster. This cluster is doped with two extra electrons, representing the effect of the reduction which takes place somewhere in the crystal surface. 共2兲Accurate ab initio embedded cluster CI calculations have been performed over a complete active space 共CAS兲 composed by two electrons and the two 3dorbitals. These calculations provide the four states 共3⌽u,1⌽u,1⌽g 1,1⌽g 2兲 with largest projections on the model space. These four states constitute the target space. The projections of the eigenvectors 3⌽uand 1⌽uonto the model space are fixed by symmetry. It is not the case for the singlet states 1⌽g 1and 1⌽g 2whose projections in the model space have a degree of freedom, namely, the ratio of the coefficients on Sg Nand Sg I共cIand cN兲. From the energies of these four states and the ratios of the coefficients cI/cNwe can univocally fix the values of the five parameters on the Bloch Hamiltonian 关Eq. 共13兲兴. The GramSchmidt procedure is less demanding, as it only requires the energies of these four states and the ratio cI/cNof the ground 1⌽g 1state. 共3兲For all the DDCI calculations, where ionic states are too high in energy, we use the intermediate effective Hamiltonian theory and Uis evaluated from the projection of the singlet state on the model space as shown in Eq. 共16兲. CALZADO, HERNÁNDEZ, AND SANZ PHYSICAL REVIEW B 77, 045118 共2008兲 045118-4 B. LDA+Ucalculations Once an estimate of the Uvalue is obtained from our ab initio embedded cluster calculations, we proceed to the study of reduced 共110兲TiO2surface by means of periodic DFT calculations. The calculations were performed using the projected augmented wave approach63 and the VASP4.6 code,64–66 with a cutoff for the plane waves of 500 eV. The electrons explicitly included in the calculations are the 3p,4s, and 3d shells of Ti and the 2sand 2pshells of O.67,68 For k-point sampling we used the lowest order Monkhorst-Pack69 set of 4⫻4⫻1kpoints, including the ⌫point. The LSDA+Uapproximation introduced by Dudarev et al.24 was used, where our Uvalue is equivalent to the Ueff parameter 共Ueff=U−J兲proposed by these authors. To describe the TiO2共110兲rutile surface a slab model is used. The slabs are obtained through replication along the three directions of a supercell that includes a portion of vacuum. The reduced surfaces are finally modeled by removing oxygen atoms from the bridge positions which have been experimentally established as the most stable vacancies.1,70–72 Two different vacancy concentrations have been considered: ␪ =0.25 and ␪ =0.5, where ␪ represents the mean number of oxygen vacancies per surface unit cell. The corresponding unit cells are composed of one primitive surface unit cell in the 关11 ¯ 0兴direction and four 共 ␪ =0.25兲or two 共 ␪ =0.50兲primitive surface unit cells in the 关001兴direction, as schematically shown in Figs. 2and 5. Regarding the slab thickness, it has been argued by one of us52 that energies are almost converged at the GGA level for those calculations using at least a five-layer slab. Assuming a similar behavior for LDA+Ucalculations, we keep the slab thickness in five layers, the analysis of the influence of the slab on the position of the band-gap state being outside the scope of this work. The notations used hereafter for the resulting super cells are 共2⫻1兲and 共4⫻1兲for vacancy concentrations ␪ =0.5 and ␪ =0.25, respectively. We perform spin polarized calculations on the triplet spin state, as suggested by previous works.22,23 The starting geometries are those obtained for the relaxed surface after oxygen removal from GGA periodic calculations by Oviedo et al.,52 and then we perform geometry optimizations in order to analyze the dependence of the band-gap states on the atomic structure, as recently suggested by Di Valentin et al.23 for the hydroxylated TiO2共110兲surface. These authors have found a close relationship between a correct structural relaxation and the occurrence of electron trapping, when hybrid exchange functionals are used. For each considered Uvalue, geometry optimizations have been carried out, maintaining the atoms in the two bottom layers fixed at their original positions, until the largest energy difference between two successive points was less than 1⫻10−3 eV. III. RESULTS A. Ab initio evaluation of U The Uvalues obtained from ab initio calculations are reported in Table I. We have analyzed 共i兲the effect of the electronic correlation, 共ii兲the dependency on the basis sets, and 共iii兲the impact of the environment model on the U value. Also the effect of the geometrical relaxation induced by the extra electrons on the Uvalue is discussed in this section. As shown in Table I, both the basis sets and the environment have a minor effect on the Uvalues, at least for the accuracy demanded for the subsequent periodic calculations, as discussed below. Only the electronic correlation plays an important role, and we focus the discussion on its effect on the Uvalue. Regarding the impact of the wave function complexity, in all cases the Uvalues obtained from the bare valence space are extremely large 共Uvalues around 13–14 eV兲, due to an artificial destabilization of the ionic states. A remarkable decrease of the Uvalue is observed when the single excitations are included in the CI expansion 共CAS+Sresults in Table I兲. This important modification of Uis due to the dynamical polarization of the ionic forms, which is introduced by the 1h+1pdeterminants. The rest of contributions included in the DDCI space produce a minor modification of U. This result suggests that for those systems where the DDCI calculation is not available, a reasonable Uvalue could be ob- [001] [110] θ=0.5 θ=0.25 VV FIG. 2. 共Color online兲Schematic representation of the two cells employed in periodic calculations: the 2⫻1共on the left兲and 4 ⫻1共on the right兲cells, corresponding to vacancy concentrations of ␪ =0.5 and ␪ =0.25, respectively. “V” represents the bridging oxygen vacancy. TABLE I. Uvalues 共eV兲obtained with different environment representations and several CI wave functions. Model Basis set CAS+S DDCI Int.Bloch Gram-Schmidt 共+4,−2兲 Model 1 1 5.80 5.84 6.16 2 5.63 5.66 6.39 3 5.43 5.47 5.54 共+2,−1兲 Model 2 1 5.37 5.45 5.95 2 5.56 5.63 6.48 3 4.93 5.02 5.24 Mean U ¯ values 5.45 5.51 5.96 EFFECT OF ON-SITE COULOMB REPULSION TERM U…PHYSICAL REVIEW B 77, 045118 共2008兲 045118-5 tained at the CAS+Slevel, at a remarkable lower computational cost. The same general trend has been found in the recent past for several magnetic binuclear Cu共II兲systems, where Uis essentially affected by the single excitations, the double excitations producing only a minor change.41 At the CAS+Slevel, both Bloch and Gram-Schmidt U values are very close, independently of the basis sets and representation of the environment. This means that the nonhermiticity inherent to the Bloch approach does not affect the Uvalue 共while it changes significantly the hopping integrals t, for instance兲. The mean Uvalue obtained at this level is around 5.5 eV. However, for the CAS+Ssolutions, the agreement with the intermediate Uvalues is not so good 共differences as large as 1.5 eV are found兲. The comparison with previous works on different systems41 revels that Uvalues provided by the intermediate Hamiltonian are always larger than those coming from Bloch or Gram-Schmidt Hamiltonians. The simplifications inherent to the intermediate Hamiltonian produce an enhancement of the Uvalue, especially for low-level wave functions, such as CAS+S ones, where the cI/cNratio of the ground 1⌽g 1state is underestimated with respect to the DDCI wave function. The differences are not so pronounced for DDCI solutions, although the same trend is observed. In spite of this limitation, it is worth to notice that the procedure based on the use of intermediate Hamiltonian is the only one available for those systems where there exists a large contamination of the ionic states 共that prevents a correct identification of the excited states兲or those with a large number of low-lying ligand-tometal charge transfer excitations. Then, the DDCI values reported in Table I, which have been obtained by means of the intermediate Hamiltonian approach, must be considered as an upper limit of the Uparameter. An additional aspect of the method used concerns whether the active space employed on the CI calculations is large enough to evaluate U. In order to analyze the suitability of the active space we have determined the singlet-triplet energy difference by means of CAS self-consistent field calculations. We have compared the results obtained when an enlarged CAS containing two active electrons and ten moleculor orbitals 共MOs兲with large Ti 3dcharacter is employed with those coming from a minimal CAS 共two electrons in two 3dMOs兲. The so-calculated singlet-triplet gaps differ by less than 1 meV. Indeed, the enlarged CAS wave functions are significantly dominated by the configurations contained in the minimal CAS 共projections larger than 98%兲. The same negligible effect is observed when the active space contains also molecular orbitals centered on the ligands 关CAS 共12MOs/6e兲兴. These results support the use of a minimal active space in our extended CI calculations. Finally, the effect of the geometry relaxation on the U value has been addressed. It is expected that the extra electrons induce a relaxation not only on the first neighbor atoms to the two Ti centers, but also to the second and perhaps third coordination spheres. However, a full geometry optimization on an extended cluster, containing several Ti centers is meaningless in the present case because, as discussed in previous section, the theoretical evaluation of Urequires a wave function distributing two electrons in two neighbor Ti sites. It has been observed that when an extended cluster is employed, the ground state corresponds to two electrons located in noncontiguous positions. This situation probably represents the most stable electronic distribution in the real system, but it is useless for our purposes, since it introduces the geometrical relaxation in quite distant Ti positions and not in two neighbor Ti centers as desired. On the other hand, a full optimization of the Ti2O10 cluster is not possible, due to the frozen positions of the total ion potentials and point charges. The nearest neighborhood of the cluster is not affected by the relaxation, which produces unrealistic final geometries. For the above reasons, we have estimated the impact of the geometrical relaxation on the Uvalue from a series of test calculations, in which a progressive breathing expansion of the Ti2O10 cluster is allowed, as well as the nearest neighbor atoms. The breathing of the cluster enlarges the Ti-Ti distance and leads to a stabilization of the neutral configuration, while the ionic ones are practically unaffected. Then the ionic-neutral energy difference, i.e., Uterm, increases. Even though significant changes in the total energy of the system are observed, the Uvalue increases by no more than 10% with respect to the original geometry. In summary, our ab initio calculations suggest a Uvalue around 5.5±0.5 eV for reduced TiO2. This value is in agreement with those extracted from XPS and EELS experiments ranging from 4 to 5 eV.57,73–75 B. Periodic LDA+Ucalculations The density-of-states 共DOS兲curves obtained from our LDA+Ucalculations are depicted in Fig. 3for the 2⫻1 unit cell. We have employed several values of Uin the LDA +Ucalculations, in order to check the sensibility of the periodic approach to the Uvalue used and to discard a fortuitous agreement between the ab initio U value and the U value which provides a correct description of the band-gap states at the LDA+Ulevel. The aim is not to suggest a U value based on empirical fitting of the band-gap states, but to test the capabilities of the strategy proposed. The left panels in Fig. 3show the DOS curves obtained at the LDA+Ulevel with GGA optimized structures. The right panels report the DOS curves obtained when also the geometry is optimized at the LDA+Ulevel. As usual for DFT calculations,76 the band gap is underestimated for all the Uvalues considered 共band gap around 2.2 eV instead of 3.1 eV兲.3As far as we know among DFT based approaches only the hybrid B3LYP functional gives a better agreement to the band-gap energy, although slightly overestimated 共3.4 eV兲.23 Examination of DOS reveals that for U=0 the Ti 3d-like surface state lays in the conduction band, as previously reported.20,21 For Uvalues ranging from 5 to 7 eV, there are localized states on the band gap. In contrast to the B3LYP calculations of Di Valentin et al.,23 these states appear both when the GGA relaxed surface is used as well as when the geometry optimization is performed at the LDA+Ulevel, probably due to the effect of the Uterm. However, the position of these states is quite dependent on the precise structural relaxation carried out, even when changes in Ti-Ti distances once the optimization has been completed never exceed 0.06 Å. As a general trend, in all cases, the electron trapping states move CALZADO, HERNÁNDEZ, AND SANZ PHYSICAL REVIEW B 77, 045118 共2008兲 045118-6 toward the valence-band edge. This is an important point to be considered in future applications of the LDA+Uprocedure, since so far most of the published LDA+Ustudies of transition metal oxides where Uvalues are fitted in an empirical way does not take into account the possible effects of the specific structural relaxation introduced for each Uvalue. For the suggested Uvalue 共5.5±0.5 eV兲two distinct peaks are observed in the density of states, 0.88 and 1.27 eV below the conduction-band edge for the LDA+Uoptimized geometry, which nicely correlate with the experimental data. For Ularger than 6.0 eV, the localized states are placed in the lower part of the band gap, close to the valence-band edge. Even when it could be possible to relate these states to the broad peaks observed in the region of 450–750 nm for irradiated rutile,11–14 the controversy regarding their assignment impedes us from doing. For U=8 eV and higher, the band gap disappears, the valence-band edge mixes with the bottom of the conduction band. Figure 4presents the DOS curves obtained for a lower vacancy concentration 共 ␪ =0.25兲when Uvalues on the range provided by the ab initio calculations are used. As for ␪ =0.5, the positions of the band-gap states are affected by the geometrical relaxation, but a nice agreement with experimental data is found for the ab initio U value. A closer inspection of Figs. 3and 4reveals that the gap states move toward the conduction-band edge as the defect concentration increases 共defect states lie at 0.54 and 0.98 eV above the valence-band edge for a ␪ =0.25 vacancy concentration, which move to 0.80 and 1.19 eV when ␪ =0.50兲in agreement with UPS measurements.77,78 The inset of Fig. 5shows this effect in detail. The curves in this figure correspond to the total density of states obtained for the 2⫻1共thin line兲and 4⫻1共solid line兲cells from LDA+Ucalculations, with U =5.5 eV on the LDA+Uoptimized geometry. In order to clarify the composition of the band-gap states and the impact of the geometry relaxation on their localized nature, we obtain the projection of the DOS 共PDOS兲on the Ti atoms for both cells, shown in Fig. 5共middle panel: 2 ⫻1 cell; bottom panel: 4⫻1 cell兲for the GGA 共left兲and LDA+U共right兲relaxed geometries 共U=5.5 eV兲. For the LDA+Ugeometry, the band-gap states are essentially centered in two specific Ti atoms in both cells, but with different nature depending on the vacancy concentration. For low concentrations, the electrons trap on subsurface Ti atoms, while for high concentrations, half of the electrons localizes on subsurface Ti atoms, the rest on surface Ti atoms, close to the oxygen defect. The relationship between vacancy concentration and Ti+3 distribution could affect the catalytic activity of the material, and even play a crucial role on subsequent apLDA+U /GG A LDA+U / LDA+U Dens i ty o f states (ar b .un i ts) -6 -4 -2 0 2 4 Energy (eV) -6 -4 -2 0 2 4 Energy (eV) U=0.0 U=5.0 U=5.5 U=6.0 U=6.5 U=7.0 FIG. 3. Density-of-states curves obtained from LDA+Ucalculations with different values of Uusing the 2⫻1 cell 共vacancy concentration ␪ =0.5兲. The zero of energy corresponds to the valence-band edge. Left panel: LDA+Uenergies with GGA optimized geometry 共LDA+U/GGA兲. Right panel: the geometry is also optimized by using the LDA+Uapproach 共LDA+U/LDA+U兲. Solid and dashed lines correspond, respectively, to the majority and minority spin components. LDA+U/GGA LDA+U/LDA+U Density of states (arb. units) -6 -4 -2 0 2 4 Energy (eV) -6 -4 -2 0 2 4 Energy (eV) U=5.5 U=6.0 U=5.0 FIG. 4. Density-of-states curves obtained from LDA+Ucalculations with three different Uvalues for a vacancy concentration of 0.25 共4⫻1 cell兲. The zero of energy corresponds to the valenceband edge. Left panel: LDA+Uenergies with GGA optimized geometry 共LDA+U/GGA兲. Right panel: the geometry is also optimized by using the LDA+Uapproach 共LDA+U/LDA+U兲. Solid and dashed lines correspond, respectively, to the majority and minority spin components. EFFECT OF ON-SITE COULOMB REPULSION TERM U…PHYSICAL REVIEW B 77, 045118 共2008兲 045118-7 plications. It is worth to mention that it could be possible to conceive several localized solutions with approximately the same energy, but in all of our calculations, independently of both the Uused and the starting spin configuration, we obtain the Ti+3 distribution described above, and we have not found an alternative way to converge over different solutions. Regarding the effect of the lattice relaxation, the comparison of the PDOS obtained from the GGA and LDA+Ugeometries indicates that the localized nature of these states is strongly related with the relaxation effects induced by the extra electrons. When the GGA geometry is used the bandgap states present a lower Ti 3dcharacter 共smaller projections on the Ti atoms兲, and larger delocalization, especially for low vacancy concentration, where only one of the bandgap states is well separated from the conduction band. IV. CONCLUSIONS The failure of LDA based methods in describing the transition metal oxide surfaces has motivated the development of alternative methods able to deal with electron correlation effects. Among those methods, LDA+Uapproach is a promising one, the main limitation is that Uis generally considered as an empirical parameter to be fixed during the calculations. In this context, a general protocol to estimate a reasonable value of Uto be injected in the LDA+Ucalculations has been presented here, and illustrated in the case of reduced 共110兲rutile surface. The procedure provides a Uvalue of 5.5±0.5 eV, which once introduced in the LDA+Ucalculations gives a correct description of the band-gap states. The 3dcharacter of these states is confirmed and the effect of the vacancy concentraTisub3 Ti sub1 Tisurf Ti sub LDA+U/GGA LDA+U/LDA+U FIG. 5. 共Color online兲Total and projected density of states for reduced TiO2共110兲, with two different vacancy concentrations ␪ =0.5 共middle panel兲and ␪ =0.25 共bottom panel兲, calculated using the LDA+Umethod on the GGA geometry 共left兲or on the relaxed LDA+U geometry 共right兲. In all cases a value of U=5.5 eV has been employed. On the top: the five-layer slabs used in the calculations 共left, vacancy concentration: ␪ =0.5; right, vacancy concentration: ␪ =0.25兲. The Ti+3 states are localized on surface and subsurface Ti ions for ␪ =0.5 and only subsurface Ti ions for ␪ =0.25. The inset shows how the increase of the vacancy concentration ␪ pushes the band-gap states toward the conduction-band edge for the LDA+Ugeometry. Solid and dashed lines correspond, respectively, to the majority and minority spin components. CALZADO, HERNÁNDEZ, AND SANZ PHYSICAL REVIEW B 77, 045118 共2008兲 045118-8 tion on the relative position is also reproduced. Moreover, our results show a strong dependence between the position of the band-gap states and the approach used to account for the lattice relaxation, with marked differences between the DOS obtained from LDA+Ucalculations using the optimal GGA geometry or the corresponding relaxed LDA+Ustructure. So, the procedure is quite promising, and can be considered as an alternative tool for those systems where electronic correlation plays a crucial role, and LDA+Udescriptions could help in understanding their properties. ACKNOWLEDGMENTS This work was funded by the Spanish DGESIC, Project No. MAT2005-1872. N.C.H. thanks the Ramón y Cajal Program from Spanish Ministerio de Ciencia y Tecnología. *Author to whom correspondence should be addressed. FAX: ⫹34954-557174. [email protected] 1U. Diebold, Surf. Sci. Rep. 48,53共2003兲. 2V. E. Henrich and P. A. Cox, The Surface Science of Metal Oxides 共Cambridge University Press, Cambridge, 1994兲. 3V. E. Henrich and R. L. Kurtz, Phys. Rev. B 23, 6280 共1981兲. 4Y. Aiura, Y. Nishihara, Y. Haruyama, T. Komeda, S. Kodaira, Y. Sakisaka, T. Maruyama, and H. Kato, Physica B 196, 1215 共1994兲. 5Z. Zhang, S. P. Jeng, and V. E. Henrich, Phys. Rev. B 43, 12004 共1991兲. 6R. Heise, R. Courths, and S. Witzel, Solid State Commun. 84, 599 共1992兲. 7U. Diebold, H. S. Tao, N. D. Shinn, and T. E. Madey, Phys. Rev. B50, 14474 共1994兲. 8W. Göpel, J. A. Anderson, D. Frankel, M. Jaehnig, K. Phillips, J. A. Schäfer, and G. Rocker, Surf. Sci. 139, 333 共1984兲. 9M. A. Henderson, Surf. Sci. 400, 203 共1998兲. 10D. C. Cronemeyer, Phys. Rev. 113, 1222 共1959兲. 11 W.-T. Kim, C.-D. Kim, and Q. W. Choi, Phys. Rev. B 30, 3625 共1984兲. 12A. K. Ghosh, F. G. Wakim, and R. R. Addiss, Phys. Rev. 184, 979 共1969兲. 13T. Asahi, A. Furube, and H. Masuhara, Chem. Phys. Lett. 275, 234 共1997兲. 14K. R. Gopidas, M. Bohorquez, and P. V. Kamat, J. Phys. Chem. 94, 6435 共1990兲. 15M. A. Henderson, W. S. Epling, C. H. F. Peden, and C. L. Perkins, J. Phys. Chem. B 107, 534 共2003兲. 16A. K. See and R. A. Bartynski, J. Vac. Sci. Technol. A 10, 2591 共1992兲. 17I. Justicia, P. Ordejón, G. Canto, J. L. Mozos, J. Fraxedas, G. A. Battiston, R. Gerbasi, and A. Figueras, Adv. Mater. 共Weinheim, Ger.兲14, 1399 共2002兲. 18B. O’Regan and M. Grätzel, Nature 共London兲353, 737 共1991兲. 19M. K. Nazeeruddin, A. Kay, I. Rodicio, R. Humphry-Baker, E. Müller, P. Liska, N. Vlachopaulos, and M. Grätzel, J. Am. Chem. Soc. 115, 6383 共1993兲. 20M. Ramamoorthy, D. Vanderbilt, and R. D. King-Smith, Phys. Rev. B 49, 16721 共1994兲. 21A. T. Paxton and L. Thiên-Nga, Phys. Rev. B 57, 1579 共1998兲. 22P. J. D. Lindan, N. M. Harrison, M. J. Gillan, and J. A. White, Phys. Rev. B 55, 15919 共1997兲. 23C. Di Valentin, G. Pacchioni, and A. Selloni, Phys. Rev. Lett. 97, 166803 共2006兲. 24S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 共1998兲. 25C. Loschen, J. Carrasco, K. M. Neyman, and F. Illas, Phys. Rev. B75, 035115 共2007兲. 26M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 共2005兲. 27S. Fabris, S. de Gironcoli, S. Baroni, G. Vicario, and G. Balducci, Phys. Rev. B 71, 041102共R兲共2005兲. 28D. A. Scherlis, M. Cococcioni, P. Sit, and N. Marzari, J. Phys. Chem. B 111, 7384 共2007兲. 29C. J. Calzado and J. P. Malrieu, Phys. Rev. B 63, 214520 共2001兲. 30C. J. Calzado and J. P. Malrieu, Eur. Phys. J. B 21, 375 共2001兲. 31J. Cabrero, C. J. Calzado, D. Maynau, R. Caballol, and J. P. Malrieu, J. Phys. Chem. A 106, 8146 共2002兲. 32I. de P. R. Moreira and F. Illas, Phys. Chem. Chem. Phys. 8, 1645 共2006兲. 33E. Bordas, C. de Graaf, R. Caballol, and C. J. Calzado, Phys. Rev. B71, 045108 共2005兲. 34E. Bordas, R. Caballol, C. de Graaf, and J. P. Malrieu, Chem. Phys. 309, 259 共2005兲. 35C. de Graaf, L. Hozoi, and R. Broer, J. Chem. Phys. 120, 961 共2004兲. 36C. J. Calzado, C. de Graaf, E. Bordas, R. Caballol, and J. P. Malrieu, Phys. Rev. B 67, 132409 共2003兲. 37N. Suaud, A. Gaita-Ariño, J. M. Clemente-Juan, and E. Coronado, Chem.-Eur. J. 10, 4041 共2004兲. 38N. Suaud, A. Gaita-Ariño, J. M. Clemente-Juan, J. SánchezMarín, and E. Coronado, Polyhedron 22, 2331 共2003兲. 39N. Suaud, A. Gaita-Ariño, J. M. Clemente-Juan, J. SánchezMarín, and E. Coronado, J. Am. Chem. Soc. 124, 15134 共2002兲. 40C. J. Calzado, J. Cabrero, J. P. Malrieu, and R. Caballol, J. Chem. Phys. 116, 2728 共2002兲. 41C. J. Calzado, J. Cabrero, J. P. Malrieu, and R. Caballol, J. Chem. Phys. 116, 3985 共2002兲. 42C. Bloch, Nucl. Phys. 6, 329 共1958兲. 43J. des Cloizeaux, Nucl. Phys. 20, 321 共1960兲. 44See, for instance, W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes 共Cambridge University Press, Cambridge, 1986兲. 45J. P. Malrieu, P. Durand, and J. P. Daudey, J. Phys. A 18, 809 共1985兲. 46C. J. Calzado, J. P. Malrieu, and J. F. Sanz, J. Phys. Chem. A 102, 3659 共1998兲. 47J. F. Sanz, C. J. Calzado, and A. Márquez, Int. J. Quantum Chem. 76, 458 共2000兲. 48C. J. Calzado, J. F. Sanz, and J. P. Malrieu, J. Chem. Phys. 112, 5158 共2000兲. 49C. J. Calzado and J. P. Malrieu, Chem. Phys. Lett. 317, 404 共2000兲. EFFECT OF ON-SITE COULOMB REPULSION TERM U…PHYSICAL REVIEW B 77, 045118 共2008兲 045118-9