On the accurate prediction of the optical absorption energy of F-centers in MgO from explicitly correlated ab initio cluster model calculations

centers that are near to each other as experimentally observed. However, the absolute value of the excitation energy is in error by ;1 eV or ;20%. Increasing the basis set reduces the calculated excitation energy for the allowed transition, reaching a value of 5.44 eV for the F center, only 9% in error with respect to experiment. Improving the basis set does not result in a better value of the excitation energy of the charged F 1 center. Attempts to improve the calculated result by geometry optimization of the region near the oxygen vacancy, enlarging the cluster model, improving the primitive Gaussian set, or enlarging the auxiliary basis set centered on the vacancy failed to further reduce the error. It is concluded that much larger basis sets are required to predict excitation energies of electrons trapped at oxygen vacancies in ionic oxides with accuracy of or better than 0.4 eV. © 2001 American Institute of Physics. @DOI: 10.1063/1.1381011 #


I. INTRODUCTION
The theoretical study of electronic excited states is a necessary step for understanding and interpreting electronic spectra. Therefore, theoretical methods have been developed to obtain accurate information about spectroscopic properties of molecular systems ͑see, for instance, Ref. 1 and references therein͒. Recent advances in electronic structure theory [2][3][4][5][6][7][8] have permitted us to reproduce the electronic spectra of a large number of small to medium sized molecules with accuracy that is comparable to experiment. 9,10 This impressive achievement permits us to interpret the experimental spectra and to assign the corresponding peaks in a rather univocal manner. Unfortunately, there are no equivalent approaches in solid state physics, and the study of excited states in extended systems is still in its infancy. 11 Among the different formalisms that have been proposed to study excited states in extended systems, we mention the efficient approach suggested by Rohlfing and Louie. 12 This approach is based on the use of the two-particle Green's function and provides an accurate description of excitonic effects in semiconductors and simple ionic insulators such as LiF.
In spite of the difficulties encountered when attempting to properly describe excited states in extended systems, there is quite a large number of situations in which accurate results can be obtained. The common feature of these cases arises from the local character of the electronic excitations, and therefore these systems can be studied by means of a cluster model approach and accurate configuration interaction wave functions. It is interesting to point out that such localized excited states are encountered in a broad range of situations with excitation energies ranging from ϳ0.001 eV for magnetic problems, to ϳ200 eV for core-level excitations and ionization energies. 13 Several examples of accurate predic-tions of the spin excitations governing the magnetic coupling constants of wide gap insulators have been reported recently. 14 -28 Likewise, the study of d→d excitations in oxides and related ionic systems has been largely improved from the earlier qualitative predictions 29-31 now reaching almost experimental accuracy, [32][33][34] and when necessary, fully relativistic calculations of the optical spectra are now possible. 35,36 Finally, the accurate prediction of core level spectra-including the fine details of satellites, and within the Dirac-Fock relativistic framework-completes this overview on the high energy range. 37,38 The cases described above are not exclusive to pure solids; impurities and defects can also be treated with the same level of accuracy, and this cluster model approach has permitted us to solve important problems. This is the case with the recent study on the optical spectra of the neutral Ag͑0͒ impurity in a KCl host crystal, 39,40 the assignment of several spectroscopic features of SiO 2 as due to well-defined point defects [41][42][43][44] and the description of core excitons in several oxides. 45 The configuration interaction cluster model approach has also been applied to characterize the spectroscopic features of oxygen vacancies-usually referred to as F-centers-in bulk and surface MgO. 46,47 These studies corroborate the experimental assignment for bulk F-centers and attempt to assign the main features of surface F-centers. Unfortunately, the accuracy reached in this particular system is less satisfactory than that reached in the cases discussed above. In fact, it is well established that optical spectra of O deficient MgO presents a characteristic band on the optical spectra at ϳ5 eV that is associated to the presence of bulk F-centers. 48 Moreover, Chen, Williams, and Sibley 48 have been able to decompose the band into two components due to F ϩ centers-one single electron trapped in the O vacancy-at 4.96 eV, and to neutral F centers-two electrons trapped in the O vacancy-at 5.03 eV. However, the theoretical calculations reported in previous works lead to an estimate of 6.0 and 5.8 eV, for bulk F and F ϩ centers with too large an error with respect to the experimental values. 46,47 Attempts to improve the result by employing a more extended cluster model and by taking into account possible geometry relaxation effects were unsuccessful. It was concluded that the remaining difference with respect to experiment was due to basis sets limitations. This is an important point that must be solved before attempting to reach a quantitative understanding of other important features of these materials, such as luminescence, or to extend the procedure to more complex materials. In this work a systematic study is presented aimed to establish the influence of the basis sets on the accurate computation of F-center excitation energies by means of explicitly correlated wave functions. It will be shown that accurate results can also be obtained in this class of point defects, but at quite a high computational cost.

II. CLUSTER MODEL AND COMPUTATIONAL DETAILS
Following the strategy used in previous works, 46,47 a cluster model was used to represent the oxygen vacancyhereafter referred to as Vac-in bulk MgO. The cluster model is similar in size to that used in our previous work, but more accurately described, i.e., all electrons in the cluster atoms are explicitly considered, a more efficient embedding technique is used, and rather large basis sets are employed. Hence, the cluster model contains the oxygen vacancy, its six-nearest neighbor Mg 2ϩ cations, the twelve secondnearest neighbor anions, and the eight cations in the third shell. To avoid an artificial polarization of the anions' electronic density they have been fully coordinated by the addition of an extra shell of Mg 2ϩ ions represented by ab initio model potentials ͑AIMP͒. [49][50][51][52] Finally, an array of 292 optimized point charges 53 was added to the cluster model to account for the crystal Madelung potential. The total cluster contains 343 centers and may be described as ͑Vac͒ (Mg) 6 (O) 12 (Mg) 8 (AIMP Mg 2ϩ ) 24 (Point Charges) 292 . Local relaxation effects have been shown to be of importance in surface F-centers, but rather negligible-less than 0.1 eV in the excitation energies-for the bulk F-centers. 46,47 Therefore, the ideal bulk geometry was always preserved. Finally, a series of calculations was also carried out in a larger model, explicitly including the fourth shell of anions, the final cluster model being ͑Vac͒ (Mg) 6 (O) 12 (Mg) 8 (O) 6 (AIMP Mg 2ϩ ) 30 (Point Charges) 280 . These two clusters will be referred to as clusters A and B, respectively.
Since the main goal of this work is to investigate whether the optical spectroscopic features of bulk MgO F-centers can be accurately reproduced by means of a configuration interaction cluster model approach, it is necessary to avoid possible sources of errors. Therefore, all electrons on the ͑Vac͒ (Mg) 6 (O) 12 (Mg) 8 region of cluster A, and the ͑Vac͒ (Mg) 6 (O) 12 (Mg) 8 (O) 6 region of cluster B have been explicitly considered. The molecular orbitals used to construct the many electron wave functions for the ground and excited states have been expressed as a linear combination of atomic natural orbital, ANO, Gaussian-type functions, 54,55 which are specially designed to accurately describe electron correlation effects while having a rather compact form. The ANO basis sets are derived from a large primitive set and contracted using a general scheme in which all primitive functions contribute to each contracted basis function. The Mg and O primitive basis sets are (13s8p3d) and (10s6 p3d), respectively. 56 The adequacy of these basis sets has been shown in a recent study of core level excitons in MgO. 45 Starting from the primitive sets described above, up to five different contracted basis sets of increasing size have been employed; the number of ANO basis functions included in each shell of ions is reported in Table I. These five different basis sets systematically improve the description of the quantum mechanical region. For cluster A the five sets were used, whereas only Basis 1 and Basis 4 ͑Table I͒ were employed within the larger cluster B. This is enough to check the stability of the results with respect to the number of ions included in the quantum mechanical region. It is worth pointing out that since these basis sets are based on the ANO scheme they provide an effective representation of the atomic orbitals. An auxiliary basis consisting of 3s, 2p, and 1d uncontracted Gaussian functions centered at the vacancy have been added mainly to help to characterize the different excited states. Exponents for this auxiliary basis have been optimized for the vacancy ground state within a cluster model and have been reported in previous works. 46,47 Second-order multiconfigurational ab initio calculations have been carried out for the lowest-lying electronic states of F and F ϩ centers of bulk MgO. For the F-centers the electronic ground state is a 1 A 1g with two electrons in the outermost orbital of a 1g symmetry, which indeed is largely dominated by the vacancy basis functions. The low-lying excited states are obtained by single excitation from the a 1g to the first virtual orbital which is of t 1u symmetry, which has a clearly marked vacancy character. Both, singlet and triplet coupling of the unpaired electrons in the excited configuration have been considered; the resulting electronic states are 1 T 1u and 3 T 1u . The local character of ground and excited states will be established in the forthcoming discussion ͑vide infra͒. For the F ϩ center the same orbitals are involved, but there is only one electron in the a 1g orbital in the ground state that is excited to the t 1u orbital, the corresponding electronic states being 2 A 1g and 2 T 1u , respectively. Notice that in this case Complete Active Space Self Consistent Field ͑CASSCF͒ is equivalent to restricted open shell Hartree-Fock. The N-electron wave functions are obtained by first defining a complete active space, CAS, with two active electrons in the a 1g and t 1u orbitals for the F center, and one electron in the same active orbitals for the F ϩ center. ͑CASSCF͒ wave functions for the corresponding electronic states were first obtained and used as zero order wave functions for a subsequent second-order perturbation calculation within the CASPT2 method developed by Andersson and co-workers, 2,3 and widely used in spectroscopic problems in molecules 9,10 and cluster models of extended systems. 13,40,45 This method can be applied to single and multideterminantal wave functions, and reduces to the well-known second-order Møller-Plesset perturbation scheme 57 to electron correlation energy for closed shell zero-order wave functions. 58 It is worth pointing out that the CAS defined above is rather limited, in fact, it is the minimum possible to define the electronic states of interest. Consequently, the CASSCF excitation energies are meaningless and will not be reported here.
Let us just point out that CASSCF excitation energies are always about 0.8 -1.0 eV larger than the CASPT2 values where all Mg and O 2s and 2p electrons are explicitly correlated.
A further technical point deserves some additional comments. In the present calculations only D 2h symmetry is imposed. For the excited states this can lead to a broken symmetry solution where the excited electron is localized in one component of the three-fold degenerate t 1u instead of being equally distributed among the three degenerate components. This is a consequence of using symmetry adapted molecular orbitals. Since each component of the t 1u manifold reduces to a different irreducible representation of the D 2h point group, it is not possible to average orbitals. Running calculations imposing either the real O h symmetry or no symmetry at all permits to check the importance of breaking the symmetry on the calculated excitation energies. Technically, this has to be carried out by explicitly removing all symmetry constraints and averaging all three orbitals of the t 1u manifold in a CASSCF calculation where the active space contains four vacancy orbitals, the a 1g and the three components of the t 1u manifold, and involves two active electrons. In this way the three a 1g 1 t 1u 1 solutions are degenerate. Unfortunately, within the present basis sets and the number of atoms being quantum mechanically described, CASPT2 calculations with no symmetry become too expensive and have not been carried out. Nevertheless, at the CASSCF level, it is found that the broken symmetry solution is as much as 0.1 eV lower than the solution which keeps the real O h symmetry.
All calculations have been carried out using the MOLCAS 4 suite of programs 59 implemented in HP J282 and J2240 workstations.

III. RESULTS AND DISCUSSION
The main purpose of the present work is to establish the computational requirements, especially the quality of the basis set, that are needed for an accurate description of the F and F ϩ center excitation energies in MgO. To this purpose, a cluster model representation of MgO is used. However, representing an extended system by means of an embedded cluster model is not free of limitations and one may wonder whether such an approach is physically reasonable. There-fore, a necessary first step involves defining a suitable procedure to check the adequacy of the models defined in the previous section. In previous works it has been suggested that a finite model is adequate because the electron density of the electrons trapped at the oxygen vacancy are localized precisely around this vacancy. For the ground state, this has been illustrated by means of density plots 60 whereas for the excited states, the extent of , the orbital to which the electron is promoted, can been measured as ͱ ͗r 2 ()͘, the square root of the expectation value of the square of the r operator for this orbital. 61 This technique has permitted us to validate the use of a cluster model approach in rather different physical situations. For instance, it permits us to differentiate the delocalized character of the MgO conduction band as opposite to the localized character of the one electron state resulting from a transition from the Mg(2p) core level to the conduction band. 62 It has also been applied to show the localized character of the low-lying excited states arising from a single excitation from the trapped electrons in MgO, 46 and the local character of the core level excitons in Al 2 O 3 and SiO 2 . 45 In the present work we use the same procedure to show that extending the basis set on the cluster atoms does not change the physical description of the resulting electronic state. Results on Table II show that changing for a more extended basis set, Basis 4, the extent of the a 1g orbital occupied by the trapped electron͑s͒ is well localized inside the cluster. This a 1g orbital appears to be rather more localized in the excited state as expected from the reduction from double to single occupancy. It is interesting to remark that the size of this orbital is almost the same when considering the F ϩ center where, again, the a 1g orbital is singly occupied. Likewise, the extent of the t 1u orbital is well localized within the atoms that are quantum mechanically treated. Nevertheless, the t 1u orbital size differs from singlet to triplet even if it is singly occupied in both cases. The smaller extent of this orbital in the triplet state is a consequence of the inclusion of the Fermi hole already in the zero-order, CASSCF wave function. Finally, we note that enlarging the basis set does not increase the size of this orbital. On the contrary, with Basis 4 ͱ ͗r 2 ()͘Ϸ5.7 a.u., whereas a value of 7.0 a.u. was obtained using a considerably smaller basis set and pseudopotentials to describe the cluster atoms. The more contracted description found here is a direct consequence of the use of the atomic-like ANO basis set.
Having established the adequacy of the present approach, we now turn our attention to the excitation energies and their dependence with respect to the basis set ͑Table III͒. Overall, the excitation energies decrease with increasing basis sets and seem to converge to the experimental result. The quality of the smallest basis set ͑Basis 1͒ is similar to that used in previous work in which pseudopotentials were use to describe the effect of the core electrons of the cluster atoms 46 and, as expected, the present calculated excitation energies are remarkably close to those reported in Ref. 46. For the largest basis set the excitation energy of the F center is 5.44 eV, only 0.4 eV ͑ϳ9%͒ in error with respect to the experimental value of 5.03 eV. Unfortunately, the error corresponding to the charged F ϩ center is still too high and indicates that the basis set requirements for this center are more stringent. This is not completely unexpected, because the presence of this charged point defect will polarize the electronic structure of its surrounding to a larger extent than that of the neutral F center. For the F center, the calculated allowed excitation energy is roughly on the error method of CASPT2, but on the high error bar, CASPT2 excitation energies are often within 0.1-0.2 eV of the experimental result.
In order to investigate the possible sources of error, several benchmark calculations have been carried out. First, we comment on the effect of geometry optimization of the atomic positions in the region near the point defect. In a previous work, 46 it was found that geometry optimization of the MgO F center does not introduce significant changes on the geometry, and consequently the excitations change only to a little extent. For the F center, and using Basis 4, the same result is again obtained: the d(Vac-Mg) slightly decreases and therefore, slightly higher energies are expected. In fact, geometry optimization increases the excitation energy by less than 0.1 eV and its effect is opposite to that of increasing the basis set size. However, for the F ϩ center the optimization enlarges the d(Vac-Mg) distance around 2-3%, and this excitation appears at slightly lower energies. Nevertheless, the changes induced by geometry optimization are too small to account for the remaining error with respect to the experimental value. Another potential source of error arises from differential correlation effects at a higher order of perturbation. To investigate this effect we computed the excitation energy using cluster A and Basis 4, but enlarging the active space, i.e., improving the zero-order wave function for the second-order treatment. Including six active electrons and six active orbitals in the CASSCF calculation leads to a CASPT2 excitation energy of 5.64 eV to be compared with 5.59 eV which corresponds to the minimum CAS used throughout this work. Enlarging the CAS further to six active electrons and eight active orbitals predicts an excitation en-ergy of 5.65 eV. This is a very small effect and strongly suggests that correlation effects are not at the origin of the discrepancy with respect to the experiment. The only additional effect that one may wish to consider before enlarging the basis set even further is the cluster size. Table IV presents a set of results obtained within Basis 4 and clusters A and B. The inclusion of a fourth shell of quantum mechanical ions induces changes in the excitation energy of 0.05 eV only, and therefore cluster size does not appear to be the source of error.
Finally, we considered a further extent of the atomic basis. On the one hand, we substituted the vacancy basis for a (3s2 p1d) ANO basis set for oxygen. Using cluster A and Basis 5 for the remaining cluster atoms, the F center allowed excitation energy becomes 5.48 eV, only 0.04 eV higher than the value reported in Table III. On the other hand, we had considered using a larger primitive set for the ANO basis for the quantum mechanical region of cluster A. This larger ANO set is specially designed for cases in which differential electronic correlation effects are important. The discussion above suggests that this is not the case, and calculated results do indeed support this statement. For Basis 5, the effect of improving the primitive set only reduces the allowed transition of the F ϩ center from 5.95 to 5.92 eV.

IV. CONCLUSIONS
In this work, a systematic study of the different computational requirements that affect the accuracy of the ab initio prediction of excitation energies of F and F ϩ centers on cluster models of MgO has been presented. Using a rather limited basis set results in excitation energies of the F and F ϩ centers that are near to each other, as experimentally observed. However, the absolute value of the excitation energy is in error by ϳ1 eV or ϳ20%. Increasing the basis set reduces the calculated excitation energy for the allowed transition reaching a value of 5.44 eV for the F center, only 9% in error with respect to experiment. However, improving the quality of the basis set does not result in such a noticeable improvement of the excitation energy of the charged F ϩ center. The origin of this difference is attributed to long range polarization effects in the case of the F ϩ center. Further attempts to improve the calculated result by geometry optimization of the region near the oxygen vacancy, enlargement of the cluster model, improvement of the primitive Gaussian set, or enlargement of the auxiliary basis set centered on the vacancy failed to further reduce the error. Likewise, con-  straining the wave functions for the excited state to be symmetry adapted has only a minor effect of ϳ0.1 eV. Therefore, it is concluded that the prediction of excitation energies of electrons trapped at oxygen vacancies in ionic oxides with accuracy better than 0.4 eV requires very large basis sets, and that it is better to improve the basis set on the atoms near the oxygen vacancy than to extend further the quantum mechanical cluster.

ACKNOWLEDGMENTS
Financial support has been provided by Spanish ''Ministerio de Educación y Ciencia'' CICyT Project No. PB98-1216-CO2-01 and, in part, by ''Generalitat de Catalunya'' Grant No. 1999SGR-00040. Part of the computer time was provided by the CESCA/CEPBA supercomputer center, thanks to a research grant from the University of Barcelona.