Quantum Chemistry of Solids / 23-Modeling and LCAO Calculations of Point Defects in Crystals

10.1 Symmetry and Models of Defective Crystals

419

In the SCM (CCM) models two finite energies are compared when determining the defect-formation energy: they are the total energies per unit supercell of the two crystals (perfect and defective). In view of the large energies involved, it is important that the same cell be used for the perfect and defective crystals to ensure cancelation of errors.

As an example, we consider the convergence of the LCAO HF AE supercell calculations of Ca and Ba atom subsitution in bulk MgO, the details of the computational scheme used can be found in [684]. In Table 10.2 the substitution energies are given as a function of the supercell size for nonrelaxed and relaxed defective lattice. The supercells of 8, 16 and 32 atoms correspond to the host crystal fcc lattice vector transformations, giving simple, face-centered and body-centered cubic lattices of supercells with L=4, 8, 16 (8, 16 and 32 atoms), respectively(see Sect. 10.1.2). Table 10.2 shows that in an unrelaxed lattice the increasing of the supercell from 8 to 32 atoms changes the defect-formation energy from 0.11 eV (Ca substitution) and 0.02 eV (Be substitution). It is also seen that for the largest supercell the Ca substitutional energy reduces from 7.09 eV to 6.43 eV when the first star of neighbors (six O ions) is allowed to relax; a further energy gain of 0.19 eV is obtained after allowing the second neighbors (12 Mg ions) to relax, the total relaxation energy, 0.84 eV, being substantial. In the case of the Be substitution, the relaxation energy, 0.27 eV, is considerably smaller. In the Be case the displacements are smaller and in the opposite direction to those for the Ca substitution, but the pattern is otherwise qualitatively siimilar. Thus, quite small unit cells appear to be able to describe the e ect of electronic rearrangement. In contrast, Table 10.2 indicates that nuclear-positional relaxation is much more long ranged, at least for fully ionic compounds such as MgO.

Table 10.2. Ca and Be HF substitutional energies ∆Ef (in eV) in MgO as a function of the size of the supercell. M is the number of atoms in the supercell, ∆R is the variation (in ˚A) of the distance between the defect and its first (6O) and second (12Mg) [684]

			Ca				Be
				∆R				∆R
M	∆Ef	First(6O)		Second(12Mg)	∆Ef	First(6O)		Second(12Mg)
8	7.20	—		—	–3.81	—		—
16	7.12	—		—	–3.80	—		—
32	7.09	—		—	–3.83	—		—
16	6.56	0.073		—	–3.93	–0.044		—
32	6.43	0.084		—	–4.00	–0.053		—
32	6.25	0.097		0.031	–4.10	–0.070		–0.022

The dependence of the defect-calculation results on the supercell size in MgO was obtained for oxygen atom removal energies in DFT PW calculations [685]: increasing of the supercell from 8 to 64 atoms changed the defect-formation energy from 10.656 eV to 10.568 eV in unrelaxed lattice and from 10.553 eV to 10.547 eV for the relaxed lattice. In this case, the e ect of relaxation appeared to be completely negligible.

The convergence of the results with the supercell size is strongly dependent on the point-defect nature. Ca and Be substitutions of Mg atom in MgO are examples

420 10 Modeling and LCAO Calculations of Point Defects in Crystals

of isoelectronic substitution, giving small impurity-atom displacement from the substituted atom site. However, in the case of the MO : X point defect (M=Mg, Ca, Sr; X = H, Li, Na, K) this displacement becomes large due the di erence between the host and substitution atoms. The di erence in the number of valence electrons in M and X atoms (two and one, respectively) is called hole trapping. Ionizing radiation produces a variety of trapped hole centers in metal oxides, both alkaline-earth oxides and perovskite-type compounds.

As a result of the impurity-atom displacement the point defect in MgO has axial symmetry along one of the fourth-order symmetry axes of the cubic lattice. This lowering of the crystalline symmetry from cubic to tetragonal generates a dipole moment within the cell. A long-range dipolar defect–defect interaction then originates among defects in neighboring cells, and larger supercells are needed to obtain the converged results. As an example, the defect formation energy of the MgO : Li center in UHF calculations (for the Li only relaxation taken into account) changes with the supercell increase are the following [568]:4.82 eV(8), 4.99eV(16), 5.06eV(32), 5.10eV(64). Here, the number of atoms in the supercell is given in brackets. When all the atoms in the supercell are allowed to relax, the defect-formation energy changes (in the calculations of the largest supercell of 64 atoms) from 5.43 eV (unrelaxed) to 4.15 eV (fully relaxed). The comparison with the data given in Table 10.2, demonstrates that for the hole-trapping center the convergence is slower than for isoelectronic substitution and the geometry optimization plays a crucial role for the hole-trapping center: the di erence in the defect-formation energy for fully relaxed and nonrelaxed largest supercells is 0.84 eV (Ca substitution), 0.27 eV (Be substitution) and 1.28 eV for the hole-trapping center. The atomic relaxations themselves are especially large for the Li atom (about 0.3 ˚A). It was found [568] that in the fully relaxed case, the convergence of the defect-formation energy is much faster (0.03 eV in going from a supercell of 8 atoms to one of 64 atoms) that shows that structural relaxation is an e ective mechanism to screen and minimize long-range electrostatic interactions induced by the dipolar nature of the defect center.

When the alkali-metal ion replaces an alkaline-earth cation, it relaxes from the perfect lattice position toward the oxygen ion (O2) along the axial direction, which brings a formal +2 charge; the electron hole localizes at the opposite oxygen (O1), which in turn relaxes away from the X monovalent ion. However, this relaxation is periodically repeated, and it would be desirable to compare its role for the supercell calculations at the center and symmetry points of the BZ.

The supercell calculations give important information concerning the charge redistribution in the defective crystal. In particular, the UHF electronic structure of the trapped-hole defect MgO : [Li]0 indicates the localization of the unpaired electron at the O1 atom [568]. However, the degree of localization is quantified by the spin moment of O1 and changes significantly when di erent Hamiltonians are considered: the spin moment decreases from 0.98 at the UHF level, to 0.41 for B3LYP, to about 0.1 with other DFT Hamiltonians [568]. The di erent degree of localization produced by DFT methods has important consequences for the atomic relaxation. Magnetic coupling constants determined by EPR and ENDOR techniques permit a direct comparison with experimental data to be made. In the particular case of the Li defect, the agreement is reasonable for the UHF result, where the hole is localized at O1. For

10.1 Symmetry and Models of Defective Crystals

421

the other Hamiltonians, the disagreement increases in parallel with the delocalization of the hole [568].

In PW calculations of defective crystals the use of PBC is the only possibility in the model choice. Apart from the problem of getting rid of interactions between defects by considering su ciently large cells or by introducing corrections by averaging over the sampled k-points [686], the supercell technique seems to provide a natural “universal” reference for the neutral point defects in crystals: the two structures that are compared satisfy the same boundary conditions, and the two wavefunctions (or density matrices) di er only in the vicinity of the defect in each supercell [687]. This is, however, no longer true when charged defects are considered. The artificial longrange Coulomb interactions between the periodic charged defect images introduce the divergence in the energy. The problem of the neutrality of the periodic array of charged defects is solved by superimposing a neutralizing charge uniformly spread in each supercell (a uniform, neutralizing “jellium” background charge is introduced to restore the electroneutrality of the unit cell).

It is also necessary to correct the defect-formation energy using an estimate of the dielectric response of the host system to the periodic ensemble of neutralized defects. While the self-energy of this charge becomes vanishingly small with increasing supercell size, this is not true for its interaction with the host crystal. Polarization of the crystalline medium by the local charge could be taken into account by allowing full relaxation of nuclear positions with very large supercell sizes. This is impossible in practice, and the corresponding finite contribution to the stabilization of the defect must be evaluated separately. With “reasonable” supercell sizes, it is necessary to correct for the interaction between charged defects. This correction can be performed classically, using the experimental ε static dielectric constant [688]. However, using the experimental ε value makes the model inconsistent, not only because of the semiempirical character of this assumption, but also because use of the static dielectric constant would imply full relaxation of electrons and nuclei in the field of the defect charge, while use of the optical dielectric constant would require that no nuclear relaxation takes place [687]. To calculate the charged point defect formation energy in SCM some other correction schemes are suggested, [689–692]. These schemes were applied for the supercell calculations of point defects in semiconductors: tetrahedral

– Si [690], diamond [693], InP [692], SiC [691] and hexagonal –AlN [694], GaN [695]. LCAO calculations of the charged point defects in metal oxides are made mainly in the molecular-cluster model, considered in the next section. As we already noted

PW molecular-cluster calculations are impossible as use of the PW basis requires the periodicity of the structure under consideration.

10.1.4 Molecular-cluster Models of Defective Solids

The molecular-cluster model (MCM) of the defective crystal is the simplest and most direct approach. As was already noted, MCM is obtained by cutting out in the crystal some fraction of atoms consisting of the point defect and several spheres of nearest neighbors, followed by an embedding of this cluster into the field of the surrounding crystal and (or) by saturating with the pseudoatoms the dangling bonds on the cluster surface. The embedding pseudopotential can also be used (see Sect. 8.2.3).

The main attractive features of a cluster model are the following: the simplicity of the mathematical formulation; the possibility of direct transfer of computational

422 10 Modeling and LCAO Calculations of Point Defects in Crystals

schemes, worked out in quantum chemistry of molecules; and applicability to almost all types of solids and point defects in them. A reasonable choice of cluster is possible when well-localized point defects are considered [696].

There are well-known [294] di culties of the cluster model connected with changes of host-crystal symmetry, pseudoatom choice at the cluster boundaries and the necessity to consider nonstoichiometric (charged) clusters. The validity of the results obtained within the cluster approach is often questionable due to some serious drawbacks: the strong unpredictable influence of the cluster shape and size on the results of calculations; the di culties arising when relating the cluster one-electron energy levels with the band energies of the perfect crystal; and the appearance of “pseudosurface” states in the one-electron energy spectrum of the cluster. The appearance of the “pseudosurface” states is accompanied by the unrealistic distortion of electronic density on the cluster boundaries and can lead to artefactual resonances between the “pseudosurface” and the point defect one-electron states.

In order to make the cluster model more realistic, a large number of embedding schemes have been developed. Most of them are oriented on special types of crystal (ionic, simple covalent) or are based on the use of special approximations in the Hamiltonian operator (local exchange, tight-binding approximations). The other group of cluster models is based on adding the e ective-potential operator (the embedding operator) to the unperturbed-cluster Hartree–Fock operator, giving the proper oneelectron solution for the perfect crystal. When the cluster shape and size are chosen, this embedding operator can be determined and used in the following perturbedcluster calculations. This so-called embedded-cluster approach is su ciently general and can be applied to di erent types of solids.

The only characteristic that to a certain extent depends on the nature of the solid and local center in it is the minimal appropriate cluster size. This size is limited by the extent of localization of one-electron states in the perturbed and nonperturbed solid. At the same time the extent of the defect-potential localization is not so critical if the cluster is chosen in such a way that the long-range polarization outside the cluster can be taken into account.

The use of localized orbitals for the cluster calculation is an e cient approach for defective crystals. To connect the perfect crystal localized orbitals and molecular cluster one-electron states the molecular cluster having the shape of a supercell was considered [699]. Such a cluster di ers from the cyclic cluster by the absence of PBC introduction for the one-electron states. Evidently, the molecular cluster chosen is neutral and stoichiometric but its point symmetry can be lower than that of the cyclic cluster. Let the localized orthogonal crystalline orbitals (Wannier functions) be defined for the infinite crystal composed of supercells. The corresponding BZ is L-times reduced (the supercell is supposed to consist of L primitive unit cells). The Wannier functions Wn(r − Ai) are now introduced for the supercells with the trans-

lation vectors and satisfy the following equation:
ˆ				εn (Al − Ai)	(10.4)
HWn(r − Ai) =
				A
ˆ				l
Here, H is the one-electron Hamiltonian of the crystal,
1
εn(Al) =	√		k	En(k ) exp(ikAl)	(10.5)
		N

10.1 Symmetry and Models of Defective Crystals

423

are the Fourier coe cients of the energy En(k ) of the nth band. The range of k summation, numbering of energy bands and the coe cients εn(Al) in (10.5) correspond to the supercell chosen, N is the number of supercells in the macrocrystal. Let

Enmin, Enmax be the minimum and maximum values of the energy in the nth band.

Introducing ∆En = 21 Enmax − Enmin and En = Enmin + ∆En we get

k

1

N

(En(k ) − En) exp(ik Al) ≤

1

N

max

En(k )

En

exp(ik Al)

≤ N ·

·

−

||

{|

|} ≤

The inequality obtained ≤ max (|En(k ) − Enmin − ∆E|3

≤ ∆En

k

1

N

(En(k ) − En) exp(ik Al) ≤ ∆En

(10.6)

is fulfilled for any choice of the supercell, defining the translation vectors Al, the

range of k summation and the band numbering. For any choice of the supercell

exp(ik Al) = N · δAl 0. Using (10.6) we have the relation

k

|ε(Al) − EnδAl 0| ≤ ∆En	(10.7)
that can be rewritten as
|ε(0) − En| ≤ ∆En, |ε(Al)| ≤ ∆En	(10.8)
Hence
ε(Al) → 0, ε(0) → En
∆En → 0	(10.9)

On increasing the size of the supercell one can approach the limits in (10.9) as closely as desired. When the size of the supercell increases, the number of energy bands in some fixed energy interval ∆E increases proportionally, and the width ∆En of the nth band decreases. From (10.8) it follows that the narrower the energy bands of the perfect crystal (in the usual band numbering for the primitive cell), the smaller will be the unit cell for which one can neglect on the right side of (10.4) all the terms with Al =Ai and consider the approximate equation

ˆ		A	W (r		A	)	(10.10)
HW (r	−			−
!			i) = εn(0)!		i

This equation is written for the Wannier function of the ith supercell. Due to (10.9) εn(0) can be considered as the average energy of the nth band (in the supercell classification of crystalline states). In fact, (10.10) corresponds to the molecular cluster having the shape of the supercell. The Bloch functions of the macrocrystal can be constructed from the cluster one-electron states:

424 10 Modeling and LCAO Calculations of Point Defects in Crystals

	1	l	!
Ψn(k, r) =	N R		exp(ikRl)W (r − Al)	(10.11)

This consideration is the foundation of the approach originally called “the cyclically embedded cluster” model [697]: the crystal is divided into the sets of regions chosen in the shape of a Wigner–Seitz supercell of the perfect crystal. The central region A contains the defect with its surrounding, the other regions together form the rest of the crystal. We refer the reader to [697] for the mathematical details of this approach. We note only that in the limit when the region A is the whole crystal this model gives the traditional HF LCAO equations for a defective crystal considered as a large molecule. In the other limit, when the defect potential is equal to zero, the cyclic cluster of the perfect crystal is obtained. The cyclically embedded cluster model is relatively simple in implementation and has been used for point defects in ionic [698] and covalent [697] crystals. The drawback of the model is that it does not take into account the polarization of the crystal by the perturbative potential of the defect.

The consideration above using a Wannier function for the supercell shows that the molecular-cluster model is more appropriate for the solids with the relatively narrow (1–3 eV) upper valence bands and su ciently large energy gap. For ionic insulating crystals these conditions are fulfilled and the convergence of the results with increasing cluster size is rapid. The results of recent molecular-cluster calculations [432] confirmed these conclusions made about thirty years before [699]. As an example, we mention the localized orbitals approach used to introduce the embeddedcluster model of α-quartz [700]. The localized orbitals are constructed by applying various localization methods to canonical HF orbitals calculated for a succession of finite molecular clusters of increased size with appropriate boundary conditions. The calculated orbitals span the same occupied Fock space as the canonical HF solutions, but have the advantage of reflecting the true chemical nature of the bonding in the system.

In strongly ionic systems, where coupling between the cluster and the background is mainly by classical Coulomb interactions, the long-range forces from the background can be accounted for by a Madelung term. A higher level of approximation applies a shell model for the environment, i.e. places a set of ions with rigid positive cores and polarizable, negatively charged shells around the cluster. The parameters of this shell model can be determined, e.g. , from self-consistent calculations on model systems. The charges in the embedding point charge model or in the shell model must be determined self-consistently with charges in the cluster [294]. Except for ionic solids, the interaction between the cluster and the background is quantum mechanical in nature, i.e. the electronic equation is to be solved.

A rigorous mathematical basis for the separation of cluster and background is provided by Green-function (GF) techniques. There are essentially two possible ways: the perturbed-crystal and the perturbed-cluster approaches. In the former approach it is assumed that the defect introduces a perturbation relative to the perfect crystal, localized to the cluster. If the GF of the perfect crystal is known in any localized basis representation, and the matrix elements of the perturbation can be constructed, the GF of the perturbed crystal can be calculated [294]. The perturbed-crystal approach looks conceptually like the cleanest and most sophisticated solution as far as the description of the crystalline background is concerned. The problems arise more with the cluster. The perturbation should be constructed self-consistently and at the same

10.1 Symmetry and Models of Defective Crystals

425

level as the crystal potential. The perfect-crystal problem can easily be solved on a plane-wave basis but then the solutions must be transformed into Wannier functions in order to construct the perturbation. Alternatively, a localized basis set can be used for the perfect crystal in the first place. The perturbed-crystal model is most often used in metals where the screening e ect of electrons is strong and the extent of atomic displacements is smaller.

The perturbed-cluster approach [701, 702] is formulated in the LCAO approximation and is based on the following sequence of steps, [703]: a) subdivide the entire defect system into a molecular cluster (C), containing the defect, and an external region (D), the indented crystal; b) calculate the wavefunction for the molecular cluster in the field of the indented crystal; c) correct the cluster solution in order to allow for the propagation of the wavefunction into the indented crystal while generating the density matrix of the defect system. Steps “b” and “c” are repeated to self-consistency. The corrective terms in step “c” are evaluated by assuming that the density of states projected onto the indented crystal is the same as in the perfect host crystal (fundamental approximation).

With respect to standard molecular-cluster techniques, this approach has some attractive features: explicit reference is made to the HF LCAO periodic solution for the unperturbed (or perfect) host crystal. In particular, the self-embedding-consistent condition is satisfied, that is, in the absence of defects, the electronic structure in the cluster region coincides with that of the perfect host crystal; there is no need to saturate dangling bonds; the geometric constraints and the Madelung field of the environment are automatically included. With respect to the supercell technique, this approach does not present the problem of interaction between defects in di erent supercells, allows a more flexible definition of the cluster subspace, and permits the study of charged defects. The perturbed-cluster approach is implemented in the computer code EMBED01 [703] and applied in the calculations of the point defects both in the bulk crystal, [704] and on the surface [705]. The di culties of this approach are connected with the lattice-relaxation calculations.

Concluding this section we compare the results obtained in the supercell and perturbed-cluster HF LCAO treatment for the cation vacancy in MgO [706]. The supercells with 16, 32 and 64 atoms have been considered. The perturbed-cluster was chosen so that the most important e ects of electronic rearrangement and ionic relaxation occur within the cluster, i.e. the cluster included su ciently large number of oxygen ions, which are the most polarizable species (perturbed-clusters of 7, 13 and 25 atoms were considered). For the largest supercell (64 atoms) and the largest perturbed-cluster (25 atoms) the cation-vacancy-formation energy was obtained (in eV) as 1.58 (1.48) and 1.29 (1.19), respectively. In brackets are given the formation energies after relaxation of first-neighbor oxygens; this relaxation is close in both models – 0.36 a.u and 0.37 a.u for the supercell and perturbed-cluster, respectively. This agreement can be called close.

The majority of the point-defect calculations are now made with the use of the supercell model as it allows use of the computer codes developed for the perfectcrystals calculations. The examples of supercell LCAO calculations of some point defects in the metal oxides are given in the next sections.

426 10 Modeling and LCAO Calculations of Point Defects in Crystals

10.2 Point Defects in Binary Oxides

10.2.1 Oxygen Interstitials in Magnesium Oxide:

Supercell LCAO Calculations

Calculations on interstitial oxygen atoms Oi in MgO crystal are of interest because of the possibility of their formation as a primary product during radiation damage.

Irradiation of alkali halides (MX) leads to the formation of primary pairs of Frenkel defects, namely F centers (electrons trapped at anion vacancies) and interstitial atoms X0. The latter are chemically active and immediately form diatomic molecular-type defects X−2 (H centers), each of which is centered on one anion lattice site.

However, much less is known about interstitial defects in MgO crystals. The e - cient creation of F centers under irradiation is known and means that complementary interstitial Oi atoms should also be formed. However, the volume change due to irradiation of MgO is surprisingly small, in contrast to that accompanying the formation of H centers in alkali halides. A di usion-controlled process related to Oi defects was experimentally observed and characterized by an activation energy of 1.6 eV. This was ascribed to the destruction of some complexes. Therefore, the theoretical study of the stability and configuration of oxygen-atom interstitials in MgO is of importance in trying to understand radiation damage of this

important ceramic material.

We discuss here the results of all-electron HF and semiempirical valence-electron CNDO LCAO calculations of oxygen interstitials in MgO crystal [682]. The supercells of 8, 16 and 32 atoms (S8, S16, S32, see Sect. 10.1.2) were taken in HF calculations, in CNDO calculations the larger supercells were taken to study the convergence of the results obtained. The di erent configurations for the interstitial oxygen atom were considered, corresponding to an oxygen atom at the volume (v), face(f ) and edge (e) centers, respectively, see Fig. 10.2. The dumbbell-shaped configurations vd, f d, ed are formed when the interstitial oxygen atom Oi is moved toward the nearest-neighbor O2− lattice anion along <111>, <110> and <001> directions, respectively. The lattice oxygen anion is simultaneously displaced from its site along the same direction, the center of the dumbbell formed is at the lattice site. To settle the question of an adequate size for the supercell, the calculations for v, f , and e configurations are made for supercells containing 9, 17, and 33 atoms. The defect energy levels in the bandgap are observed to change between the 9-atom and 33-atom supercells, although the influence of supercell size on the charge distribution is relatively small. In HF LCAO calculations the size of the supercell is limited by the computational facilities available. Therefore, the convergence of the properties of these defects was investigated also by means of the semiempirical complete neglect of di erential overlap (CNDO) approximation, using a code, written by Tupitsyn [707] for the CNDO band calculations. In the CNDO approximations it was possible to use supercells containing 65, 129, and 257 atoms, obtained by enlarging all three smaller supercell translation vectors by a factor of two. It was found that both the charge distribution and the defect levels in the bandgap for the larger supercells were close to those obtained for the 33-atom supercell.

HF and CNDO calculations were done to investigate convergence with respect to supercell size. All basis atoms, except those involved in a defect, occupy lattice sites. The CNDO results show that the defect-formation energy changes by only 0.3%

10.2 Point Defects in Binary Oxides

427

on increasing the supercell size from 33 to 257 atoms. The HF energies generally change by less than 0.03% on increasing the supercell from 17 to 33 atoms. This suggests that the variations observed in the CNDO energies may be due more to inherent approximations in the method than to defect interactions and that a 33atom supercell is definitely su cient. If we were only interested in the energy of the unrelaxed crystal then an S8 supercell, containing eight primitive unit cells, would in fact be adequate. However, the number of shells of atoms around the defect that may be relaxed is dependent upon the size of the supercell, and the 33-atom supercell is the smallest that can be used to obtain quantitatively accurate results.

In Table 10.3 is given the optimized distance between the two dumbbell atoms calculated by both the HF and CNDO methods. The labels S16, S64, and S128 are used for supercells containing 16, 64, and 128 primitive unit cells, respectively. q is the charge on each of the two oxygen atoms in the dumbbell configuration and d is the optimized distance between the two O atoms in the dumbbell. ∆Ed is the energy change on forming the dumbbell from an interstitial in the rigid lattice.

Table 10.3. Results of HF and CNDO calculations for O atom interstitials in MgO (only the dumbbell oxygen atoms are relaxed), [682]

Configuration	Supercell	˚	q(e)	−∆Ed
		d(A)
ed	S16 (HF)	1.310	–1.160	—
	S16 (CNDO)	1.280	–0.944	32.6
	S64 (CNDO)	1.280	–0.944	32.6
	S128 (CNDO)	1.282	–0.944	32.6
vd	S16 (HF)	1.355	–1.070	—
	S16 (CNDO)	1.394	–0.918	7.9
	S64 (CNDO)	1.392	–0.917	7.9
	S128 (CNDO)	1.393	–0.917	7.8
f d	S16 (HF)	1.372	–1.033	—
	S16 (CNDO)	1.369	–0.929	3.2
	S64 (CNDO)	1.366	–0.928	3.2
	S128 (CNDO)	1.369	–0.928	3.3

The charge q on each of the two ions in the dumbbell is close to –1, negative charge on the dumbbell in excess of –2 is being made up by a slight deficit on the next-nearest ions. The optimized HF values for d are close to those found for a 17-atom supercell using the full-potential linear mu n-tin orbital (FPLMTO) method [708]. The optimized d values in the CNDO approximation are also close to those obtained using the HF and LMTO methods. Since the CNDO method is semiempirical these calculations were performed with all atoms except the added interstitial, or the two dumbbell atoms, fixed on their normal lattice sites. In the HF calculation, d was first optimized for a crystal with all except the defect atoms on perfect lattice sites. The charges on the two dumbbell atoms are again close to –1 and the values of d and q (the charge on each of the dumbbell atoms) in Table 10.3 are from these calculations. The optimum relaxed positions of the near-neighbor ions were then determined iteratively, and finally the value of ∆Ed was reoptimized.

428 10 Modeling and LCAO Calculations of Point Defects in Crystals

The values of d and q with the relaxed positions are in Table 10.4 together with the energy changes that accompany the formation of a dumbbell from the corresponding interstitial (ed from e, and so forth) and the further energy changes associated with relaxation of near neighbors. ∆Ed is the di erence between the total energy of the unrelaxed crystal in the dumbbell configuration and that for the unrelaxed crystal with the added O atom in the interstitial site, ∆Erel is the further energy change (in eV) accompanying the relaxation of near neighbors.

Table 10.4. Results of HF calculations of interstitial oxygen in MgO using S16 supercell, [682].

Type of dumbell	˚	q(e)	−∆Ed (eV)	−∆Erel (eV)			˚
	d (A)				Relaxations δ (A)
ed	1.319	21.001	28.88	0.12	δMg1	=	0.241
					δMg2	=	0.051
					δMg2	=	0.051
					δO	=	0.0245
					δO	=	0.1030
vd	1.383	21.049	15.94	0.27	δMg	=	0.077
					δO	=	0.1057
f d	1.324	21.010	7.15	2.59	δMg1	=	0.152
					δMg2	=	–0.010
					δO	=	0.0997

Subscripts 1,2 mean nearestand next-nearest atoms (of that species) to the dumbbell, see Fig. 10.2c. For ed, δO is the displacement parallel to the axis of the dumbbell and δO the displacement perpendicular to the axis of the dumbbell. In each case all the atoms of a symmetry-related set that lies wholly within the supercell were allowed to relax.

The importance of supercell size is shown, for example, by the calculations for vd where for S8, only the six nearest-neighbor Mg atoms may be relaxed because the nearest O atoms lie on the supercell boundary. This may be why comparable FPLMTO calculations [708] for S8 found the vd defect to be slightly more stable than f d, whereas in HF calculations the reverse is true (Table 10.4). The decrease in energy ∆Ed on forming a dumbbell in a perfect lattice from a v, f , or e-centered interstitial O atom (calculated as the di erence between the HF total energies for the two configurations) is largest (28.9 eV) for the ed center and smallest for the f d center (7.1 eV). The former result is to be expected because of the greater overlap in the e configuration. The CNDO results for ∆Ed agree with the HF results only qualitatively, but they serve to confirm convergence with increasing size of supercell, which was the reason for doing the calculations. The further energy gain on the relaxation of near neighbors, −∆Erel, was found to be small compared to that on forming the dumbbells, which evidently provides a very e ective mechanism for reducing charge overlap.

The optimized displacements are given in Table 10.4 and these are seen not to exceed 0.1–0.2 ˚A. With respect to the most stable f d relaxed dumbbell configuration the vd and ed configurations are less favorable by 1.31 and 1.86 eV, respectively. So O atoms formed in MgO during radiation damage will become trapped in face-