Plane-wave density functional theory was used to study the properties of oxygen vacancies and interstitials, with different charge states, in MgO. The calculated properties were the relaxed configurations, the Frenkel defect formation energies and the energies of the migration barriers, and all properties were found to be strongly dependent on the defect charge state. The lowest energy configuration of the O2− interstitial was found to be the cube centre; however, the O− and O0 interstitials formed dumb-bell configurations. The Frenkel defect energies were also strongly dependent on the defect charge, with the neutral pair energy calculated to be 3 eV lower than the doubly charged Frenkel pair defect energy. The migration barriers of the oxygen vacancies were found to increase as the net charge of the oxygen vacancies decreased, which suggests that vacancies with trapped electrons are much less mobile than the F2+ vacancies modelled by classical potentials. The migration of the oxygen interstitials showed particularly interesting behaviour. The O0 interstitial was found to have a higher migration barrier than the O2− interstitial but a very low barrier (0.06 eV) was found for the O− interstitial. The results have significant implications for the reliability of classical radiation damage simulations.
Over the past three decades, there has been extensive interest in calculating the properties of defects in magnesium oxide, both because of its many important applications and because it is often considered to be a model ceramic material. Early calculations were carried out using classical interatomic potentials (Catlow 1977; Mackrodt & Stewart 1979; Catlow & Stoneham 1983; Bush et al. 1994) and the resulting defect formation energies gave reasonable agreement with experimental measurements (Oishi & Kingery 1960; Barr & Lidiard 1970; Hashimoto et al. 1972; Shirasaki & Hama 1973; Shirasaki et al. 1973; Wuensch et al. 1973; Cohen & Gordon 1976). In recent years, the classical calculations have been supplemented by ab initio calculations and these methods have indicated that the classical methods are generally reliable (Gilbert et al. 2007).
Defects in ceramic materials are of particular interest in the field of radioactive waste disposal. With the increasing realization that nuclear power represents one of the few viable options for a low-carbon economy, research into safe methods for the disposal of high-level nuclear waste is essential if the technology is to be widely deployed. The favoured method is encapsulation of the waste into radiation-resistant ceramic materials before burial in deep depositories (Baldwin et al. 2008). The encapsulating material must remain intact for several thousand years to prevent migration of the radioactive material into the environment. A detailed fundamental understanding of the evolution of radiation damage is essential for the extrapolation of experimental data to geological time scales.
Radiation events, such as the decay of encapsulated radio-nuclides, produce high-energy particles, such as alpha particles, and recoil nuclei with energies of the order of several tens of kiloelectron volts. The high-velocity nuclei and particles exchange energy with the host atoms and create regions of local dynamic disorder in the lattice. The disorder recrystallizes over a period of a few tens of picoseconds but residual defects remain, and it is the long time-scale evolution of these residual defects that will eventually modify the microstructure of the material (Trachenko et al. 2006; Zhang et al. 2008). The initial radiation event can be successfully modelled by molecular dynamics cascade simulations, but long time-scale methods, such as kinetic Monte Carlo (Kittiratanawasin et al. 2009) and rate theory (Ortiz & Caturla 2007), are required to investigate microstructure evolution.
Radiation damage simulations have traditionally employed classical methodo-logies because the length and time scales required are orders of magnitude higher than those feasible for ab initio methods. There has long been a realization, however, that important effects are neglected in classical models, in particular the role the electronic excitations. There has been a recent revival of interest in attempting to understand and include these effects, with notable progress in metallic materials (Duffy & Rutherford 2007; Duffy et al. 2008; Race et al. 2009). Electronic excitations in insulating materials are more complex than those in metals as the band gap results in a range of additional processes including the trapping and self-trapping of electrons, holes and excitons (Itoh & Stoneham 2001). The trapping of electrons and holes at defects (Chen & Abraham 1990; vacancies and interstitials) in ionic crystals is particularly relevant to radiation damage, as the trapped carriers will effectively change the net defect charge and this will have a strong effect on defect migration energies, and hence the microstructure evolution. Such effects cannot be reliably studied using classical methods. The effects of charge localization have been studied previously in a range of oxides, such as TiO2 (Deskins & Dupuis 2007), SiO2 and HfO2 (Shluger et al. 2009).
Vacancy and interstitial migration energies in MgO have been calculated for a range of empirical potentials; however, these are, necessarily, limited to defects with the full formal charge. There have been a limited number of studies using ab initio techniques but these have mainly focused on the doubly charged O (F2+) and Mg (V2−) vacancies (Gilbert et al. 2007) in order to validate the classical models. The neutral O vacancy (F0) was calculated using plane-wave density functional theory (DFT) by Carrasco et al. (2004). The migration energies of the singly charged O vacancy (F+) and the O− interstitials have not previously been studied using ab initio methods. In this paper, DFT is used to calculate geometry and migration energies of oxygen vacancies and interstitials in MgO, with a range of charge states, and results are compared with classical potentials, semi-empirical and other ab initio calculations where these are available.
Trapped electrons and holes will also affect the defect formation energies in ionic crystals. Electrons are known to localize on O vacancies in MgO and holes can localize on O interstitials. Holes may localize on Mg vacancies but this localization is too weak to be described by the plane-wave DFT methodology employed in this paper and, in any case, the weak localization is unlikely to have a significant effect on radiation damage. In this paper, we focus on oxygen defects, because of the stronger charge localization. We investigate the effect of charge localization on oxygen vacancy and interstitial formation energies in MgO. We calculate defect formation energies for O vacancies with zero (F2+), one (F+) and two (F0) trapped electrons and O interstitials with zero, one and two trapped holes. These energies are used to calculate Frenkel defect energies for neutral, singly charged and doubly charged defect pairs along with the corresponding migration barrier energies for the vacancies and interstitials. Frenkel pairs are created by radiation events in crystals and it is the creation and migration of these defects that dominate the long time-scale evolution of the microstructure.
The defect formation energies and the migration energies were calculated using DFT with the Vienna ab initio simulation programme (VASP; Blöchl 1994; Kresse & Joubert 1999) plane-wave code. The generalized gradient approximation (GGA) with Perdew and Wang (PW91) exchange correlation potentials (Perdew et al. 1992), projected augmented waves (PAW) and a plane-wave basis set cut-off of 700 eV were used with Γ point sampling. A simulation cell with 3×3×3 unit cells of MgO (216 ions) was constructed and the cell was relaxed until the energy difference converged to 0.01 eV, which resulted in an MgO lattice parameter of 4.14 Å. This compares favourably with the (room temperature) experimental value of 4.24 Å(Schmahl & Eikerling 1968). A vacancy was created in the cell by removing an Mg or an O atom. The number of electrons was adjusted to give the required charge state of the defect. A geometry optimization, with the cell parameters fixed, was performed until the forces converged to 0.01 eV Å−1. Interstitial energies were calculated by adding an O or an Mg atom to the cell in two distinct configurations and adjusting the number of electrons to obtain the required charge state. In one configuration, the interstitial configuration was initiated with the extra atom in the centre of a cube of atoms, and in the other a dumb-bell interstitial, oriented in the 〈110〉 direction, was created by adding the interstitial atom close to a lattice atom.
The migration energies were calculated using the standard nudged elastic band (NEB) functionality of VASP (Mills et al. 1995). The initial and final configurations of the defect migration path were optimized and used as the end points for the NEB calculations. A spring constant of 5 eV Å−2 was used with five images along the path. Each image was optimized until the forces converged to 0.01 eV Å−1.
Mg vacancies in different charge states are not considered here because standard DFT methods describe hole localization on Mg vacancies incorrectly. In the singly charged (V−) vacancy, the hole is known to localize on one neighbouring O ion and the two holes localize on O ions on opposite sides of the neutral (V0) vacancy (Stoneham et al. 2007), whereas DFT calculations show the holes delocalized over all six O neighbours. Estimates of the migration barriers, calculated using constrained minimization, showed the Mg vacancy barriers to be independent of oxidation state, which is an anomalous result caused by the hole delocalization. The Mg2+ interstitial does not localize electrons because the conduction band of MgO is constructed from the s orbital of the magnesium ions (de Boer & de Groot 1998); therefore, the Mg interstitial is not considered here.
The neutral (Schottky and Frenkel pair) defect combination energies are calculated from the isolated point defect energies using the following formulae. The Schottky defect energy (Es) is defined as 2.1Here EvMg is the energy of the cell containing a magnesium vacancy, EvO is the energy of the cell containing an oxygen vacancy, qΔV is the potential alignment correction (Landy & Zunger 2005), NMgO is the number of MgO units in the perfect cell (108 in this case) and EP is the total energy of the perfect cell.
The Frenkel defect energy (EF) is defined as 2.2
3. Results and discussion
(a) Defect formation energies
The defect formation energies were calculated for O vacancies and O interstitials in three different charge states and the doubly charged Mg vacancy (V2−). The formation energies for the isolated defects were combined to obtain the Schottky and Frenkel defect energies for comparison with previous calculations. As only the doubly charged Mg vacancy energy could be reliably calculated, the Schottky defect energy (equation (2.1)) was calculated for the doubly charged defects only and the result is shown in table 1. This is lower than the corresponding values for the Schottky defect calculated using a range of classical potentials (Catlow et al. 1976; Mackrodt & Stewart 1979; Vocadlo et al. 1995; Busker et al. 2000; Gilbert et al. 2007). Doubly charged Schottky defect energies have previously been calculated using ab initio techniques using the Hartree–Fock embedded cluster methodology (Grimes & Catlow 1990) and periodic DFT (De Vita et al. 1992; Alfé & Gillan 2005; Gilbert et al. 2007) using 32, 128 and 180 ion cells, respectively, with local density approximation (LDA) pseudopotentials. De Vita et al. (1992) and Alfé & Gillan (2005) employed Γ-point sampling, while Gilbert et al. (2007) used a 6×6×6 k-point mesh with a localized orbital methodology. The lower value for the Schottky defect energy calculated here can be explained by the GGA methodology employed, which tends to predict lower cohesive energies (Perdew & Schmidt 2001), thus lower defect energies, than the LDA methodologies employed in previous publications.
The energies of O interstitials, with three different charge states, were calculated by relaxing from two different starting configurations for each charge state, as described in the previous section. The configurations with the lower energies were identified and these are shown in figure 1. For the O2− interstitials, the preferred interstitial site was found to be in the centre of a cube of atoms, as was seen in previous classical simulations, but other configurations were favoured by the O interstitials with fewer electrons. The O− interstitial favoured a 〈110〉 dumb-bell conformation (figure 1b) with the split interstitials separated by 1.87 Å. The neutral interstitial (O0) relaxed into a 〈111〉 dumb-bell configuration with a separation of 1.42 Å between the split interstitial ions (figure 1c). The semi-empirical calculations of Brudevoll et al. (1996) also found this to be the favoured dumb-bell configuration for the neutral interstitial. It was previously assumed by experiments (Halliburton & Kapper 1978) that the dumb-bell was orientated in the 〈110〉 direction but the calculations presented here suggest that the lowest energy configuration is the 〈111〉 dumb-bell.
A Bader analysis (Henkelman et al. 2006) was carried out for the relaxed interstitial configurations in order to investigate the charge distribution around the defect sites. The Bader analysis on the neutral split interstitial indicated that −0.8 e was transferred from the lattice ion to the neutral interstitial, leading to the strong dumb-bell formation. In the case of the singly charged interstitial, −0.2 e was transferred from the lattice O ion to the interstitial and the remaining −0.2 e from the lattice ion was distributed around the surrounding O ions, giving an equal charge on the dumb-bell atoms. The smaller electron transfer for the O− interstitial accounts for the larger dumb-bell separation, as the lower shared electron density forms a weaker covalent bond. The O2− interstitial, which is located in the centre of a cube of ions, has 0.1 e more charge than the lattice oxygen ions.
The calculated O interstitial formation energies were combined with the vacancy formation energies to calculate Frenkel energies (equation (2.2)) for three oxidation states and the results are summarized in table 2. The available classical and other ab initio results are included in the table. There is good agreement between the DFT results and the classical results for the O2− Frenkel defect, indicating that empirical potentials give a good description of the doubly charged defects in MgO. There is a wide scatter in the results obtained using empirical potentials, with values ranging from 13.57 to 15.2 eV (Mackrodt & Stewart 1979; Busker et al. 2000; Uberuaga et al. 2005). Previous DFT calculations give Frenkel defect energies of 12.17 and 10.35 eV for the oxygen and magnesium ions, respectively (Gilbert et al. 2007), using a 180-ion cell with localized orbitals and LDA pseudopotentials.
There are no previous calculations for O0 and O− Frenkel defects in MgO; however, our calculations indicate that the neutral O0 Frenkel defect has a formation energy which is 3 eV lower than the O2− Frenkel defect energy. This has significant implications for radiation damage because the number of oxygen Frenkel pairs created by a radiation event could be much higher than predicted from classical simulations.
(b) Defect migration energies: interstitials
The migration barrier energies of the three oxygen interstitials were calculated by NEB (five images), with the points fitted with a cubic spline in figure 2, and the results are summarized in table 3. The barrier for the O2− interstitial is lower than those predicted by classical simulations (Uberuaga et al. 2005) and previous LDA–DFT calculations (Brudevoll et al. 1996; Gilbert et al. 2007). The barrier is lower than the classical simulations because of the lower charge on the interstitial oxygen in this simulation (−1.5 e) than the formal charge in the classical simulation (−2 e), which results in a lower coulombic repulsion between the two dumb-bell oxygen ions at the saddle point in the DFT simulations. The saddle point for the O2− interstitial migration is the 〈111〉 dumb-bell configuration, with a split interstitial distance of 2.27 Å, as shown in figure 3a. The electron density of the interstitial does not change significantly during migration. The O0 interstitial migration energy was lower than that found previously although the same migration pathway was found in both studies. The lower barrier found in this study is due to the larger lattice constant and larger unit cell than the LDA simulations of Brudevoll et al. (1996). The migration path for the O0 interstitial was calculated to be along the diagonal of the face of the cube of ions, with the saddle point shown in figure 3c. The migration involves the rotation of the 〈111〉 dumb-bell to a 〈110〉 dumb-bell, the motion along the 〈110〉 direction, the formation of another 〈110〉 dumb-bell and, finally, rotation to the 〈111〉 dumb-bell configuration. Bader analysis shows that there is a redistribution of the electron density as the oxygen interstitial passes through the saddle point. The interstitial ion loses −0.09 e compared with the initial structure, while the lattice oxygen from the dumb-bell gains 0.44 e as the dumb-bell is broken. The extra electron density comes from the lattice ion that will become part of the dumb-bell, with the loss of −0.35 e compared with the initial configuration. This redistribution of the electron density results in both lattice oxygen ions involved in the migration process having the same number of electrons at the saddle point.
The singly charged interstitial was found to have a significantly lower migration barrier (0.06 eV) than the other charged states. In order to check the accuracy of the calculation, the number of NEB images was increased to 11, but this was found to have no effect on the barrier height. The saddle point for the O− interstitial migration path (figure 3b) has the interstitial slightly out of the 〈010〉 plane. There is a redistribution of the electronic density at the saddle point with the lattice ion from the initial dumb-bell regaining −0.16 e of the −0.2 e it initially donated to the dumb-bell. This electron density is donated from the lattice oxygen the interstitial is moving towards, with this oxygen losing −0.14 e compared with the initial configuration. There is no significant change in the electron density associated with the interstitial during migration.
(c) Defect migration energies: vacancies
The energies along the migration paths, calculated from the NEB images, are plotted for the O vacancies with three different charge states in figure 4. The activation energies for the three oxygen vacancies were calculated by fitting the points using a cubic spline and the migration barriers are summarized in table 4, along with the values calculated using other methods. We note that for each extra electron that is trapped at the vacancy the activation energy increases by approximately 1 eV. The F0 and F+ migration energies are significantly higher than the results obtained previously using semi-empirical methods, but the F0 migration energy is consistent with the previous DFT result. The F2+ migration barrier is in good agreement with experimental results obtained from the exchange rate of oxygen between the gas phase and the solid oxide using oxygen isotopes, which range between 2.42 and 2.71 eV (Oishi & Kingery 1960; Hashimoto et al. 1972; Shirasaki & Hama 1973).
The effects of charge localization on oxygen point defect properties in MgO were investigated using periodic DFT. The charge state was found to have a strong effect on the minimum energy configuration of O interstitials, with O0 and O− interstitials relaxing to dumb-bell configurations whereas the minimum energy configuration for the O2− interstitial was found to be in the centre of a cube of ions. The dumb-bell formed by the O0 interstitial resulted in more charge transfer from the lattice ion (0.8 e) and a shorter split interstitial separation (1.4 Å) than the O− dumb-bell (0.2 e and 1.9 Å). The relaxed configuration has a significant effect on the migration pathway of O interstitials and, therefore, the migration barriers for the different oxidation states. The O interstitial showed enhanced radiation diffusion, with the trapping of a hole on the O2− interstitial resulting in a surprisingly low migration barrier. The O interstitial formation energies were combined with O vacancy formation energies to calculate Frenkel defect energies for three oxidation states. It was found that the neutral Frenkel pair had a formation energy 3 eV lower than the doubly charged Frenkel pair and that the singly charged O Frenkel pair had an energy 2 eV lower than the doubly charged Frenkel pair.
The charge state was also found to have a large effect on the migration barriers for the oxygen vacancies. The F0 vacancy was found to have a migration barrier 2 eV higher than the F2+ vacancy, and the F+ vacancy was found to have a barrier height intermediate between the F0 and F2+ vacancies.
The results presented in this paper have significant implications for radiation damage simulations of MgO and other oxide materials, such as UO2. Radiation damage simulations have been traditionally carried out using empirical potentials with formal charges, thus only charged defects are created in these simulations. However, if the neutral defects are more stable, as suggested by the DFT calculations, then many more oxygen defects may be created than classical methods predict. In addition, the neutral oxygen defects will have significantly higher migration barriers than doubly charged defects; therefore, the diffusion rates to sinks and clusters will be lower. The results suggest that O defects could survive longer as isolated point defects than Mg defects, which have much weaker charge localization, and this will have a dramatic effect on the microstructure evolution of oxide materials subjected to radiation damage.
We are grateful to Prof. Alex Shluger, Dr Keith McKenna and Prof. Marshall Stoneham for useful discussion. Via their membership of the UK’s HPC Materials Chemistry Consortium, which is funded by the EPSRC (EP/F067496), this work made use of the facilities HECToR, the UK’s national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc. and NAG Ltd, and funded by the Office of Science and Technology through the EPSRC’s High End Computing Programme. We acknowledge support from the EPSRC via the DIAMOND consortium (EP/F055412/1).
One contribution to a Special Feature ‘High-performance computing in the chemistry and physics of materials’.
- Received October 7, 2010.
- Accepted February 4, 2011.
- This journal is © 2011 The Royal Society