The different roles of Pu-oxide overlayers in the hydrogenation of Pu-metal: An ab initio molecular dynamics study based on vdW-DFT+U
Based on the van der Waals density functional theory (vdW-DFT)+U scheme, we carry out the ab initio molecular dynamics (AIMD) study of the interaction dynamics for H impingement against the stoichiometric PuO(111), the reduced PuO(111), and the stoichiometric -PuO(111) surfaces. The hydrogen molecular physisorption states, which can not be captured by pure DFT+U method, are obtained by employing the vdW-DFT+U scheme. We show that except for the weak physisorption, PuO(111) surfaces are so difficult of access that almost all of the H molecules will bounce back to the vacuum when their initial kinetic energies are not sufficient. Although the dissociative adsorption of H on PuO(111) surfaces is found to be very exothermic, the collision-induced dissociation barriers of H are calculated to be as high as eV and eV for stoichiometric and reduced PuO surfaces, respectively. Unlike PuO, our AIMD study directly reveals that the hydrogen molecules can penetrate into -PuO(111) surface and diffuse easily due to the % native O vacancies located along the 111 diagonals of -PuO matrix. By examining the temperature effect and the internal vibrational excitations of H, we provide a detailed insight into the interaction dynamics of H in -PuO. The optimum pathways for hydrogen penetration and diffusion, the corresponding energy barriers ( eV and eV, respectively) and rate constants are systematically calculated. Overall, our study fairly reveals the different interaction mechanisms between H and Pu-oxide surfaces, which have strong implications to the interpretation of experimental observations.
pacs:68.35.Ja, 71.15.Pd, 71.27.+a, 82.65.Pa
Plutonium (Pu) is one kind of strategically vital fissile material owning attractive nuclear properties for both energy production and nuclear explosives. From basic point of view, Pu is the most complex element in the periodic table showing such intricate properties as the complicated Pu- states Savrasov2000 (); Petit2003 (); Moore2009 (); Zhu2013 (), which are very sensitive to the slight variations of the surrounding physical and chemical environments. Thus, since it was first isolated and the samples became available, Pu has always been challenging the boundaries of knowledge of the fundamental science and attracting extensive researching attentions. However, during the preparation/handling of Pu-sample and the long-term storage of Pu-device, the reactive Pu metal will easily suffer from the environment-dependent chemical corrosions, among which the surface oxidation is to some extent unavoidable because of the ubiquitous oxygen in the realistic environment. Also, the surface hydrogenation is catastrophic for Pu and need to be avoided since the Pu-hydride product can violently catalyze the oxidation (by a factor of ) and even induce the pyrophoricity of Pu. Actually, the surface corrosion and oxidation of Pu are often treated as equivalent topics since the readily generated surface layers of Pu oxides prominently influence other corrosion courses of Pu metal such as the first step of Pu-hydrogenation named the induction period. Moreover, given the high radioactivity and toxicity of Pu, the easily dispersed Pu-hydride and Pu-oxide undoubtedly bring more hazards to the personnel and the environment during handling and storage of Pu metal. Therefore, understanding the effects of Pu-oxide surface layers in the chemical corrosion (especially the hydrogenation) of Pu metal is crucial for the plutonium science, the environment, and the nuclear energy sectors. However, it is also full of challenges that beset experimentalists.
As a matter of fact, the oxidized surface layers are directly involved in the corrosion process of Pu metal. In many technological applications of Pu oxides such as in nuclear fuel cycle, it is always desirable to study the surface properties and surface reactivity of Pu oxides. Many experiments Colmenares1975 (); Has2000 (); Butterfield2006 () have been devoted to studying the oxidation mechanisms of Pu metal under various atmospheres and temperatures, and eventually come to a main conclusion that the multivalent nature of Pu presents an active response to the varying environment, in other words, it is difficult to stabilize Pu at one oxidation state when the surrounding conditions (such as relative humidity, oxygen content, and the temperature) change. As a result, the complex Pu-O phase diagram Wriedt1990 () has not been completely described yet, let along the detailed surface conditions of such many kinds of Pu oxides. When exposed to oxygen-rich air, metallic Pu surface oxidizes to a protective dioxide (PuO) layer Has2000 (); Butterfield2006 (); Butterfield2004 () and a thin sesquioxide (PuO) layer existing in between. Under special aqueous condition, the hydration reaction on PuO surfaces (PuOHOPuOH) has been reported to generate higher compound PuO (x) HaschkeSci2000 (). Whereas, it has been theoretically proved to be strongly endothermic Korzhavyi2004 () and the products (PuO) have been experimentally found to be chemically unstable at elevated temperatures Gouder2007 (). Thus, PuO is the stable oxide as the compound of choice for the long-term storage of Pu and as the nuclear fuels. However, under oxygen-lean conditions (in the vacuum or inert gas), the PuO-layer undergoes the thermodynamically driven auto-reduction (AR) to PuO, including the low-temperature phase of -PuO in cubic structure () and the high-temperature -PuO in a hexagonal structure (). More recently, sub-stoichiometric sesquioxide PuO ( to ) has been obtained from further AR of PuO layer on Pu surface GARC2011 (). All of those studies have shown us the complex surface conditions of oxide-coated Pu, and made us believe that they are one of the controlling factors for the subsequent corrosion of Pu. As pointed out hereinbefore, for the safe handling and storage of Pu-based materials, it is a very important topic to control or avoid the catastrophic Pu-hydrogenation. Whereas, early experiments have found that the length of induction period in Pu-hydrogenation is mainly determined by the component and configuration of the preexisting Pu-oxide overlayer, namely, compared with the PuO overlayers under ambient conditions, the PuO under ultra high vacuum (UHV) conditions can promote the Pu-hydrogenation reaction by markedly shortening the induction period Has2000 (). Since then, experimentalists have invested persistent efforts to better understand the roles of Pu-oxide overlayers and to quantify the controlling factors/parameters in detail MORRALL2007 (), aiming to establish a predictive model of Pu-hydrogenation. More recently, Dinh et al. Dinh2011 () have investigated the role of cubic -PuO in the hydrogenation of PuO-coated Pu by chemically and mechanically altering the PuO surface on Pu, and demonstrated that induction periods are consistently eliminated at conditions that produce catalytic -PuO. They also pointed out that more experimental efforts are surely needed to clarify such unclear factors for the unpredictable induction period as the solubility and diffusion of hydrogen species in Pu oxides. McGillivray et al. McGillivray2011 () have determined hydrogen pressure () dependence of two early steps (namely, induction and nucleation) of Pu-hydrogenation reaction by synthesizing a well reproducible PuO overlayer. They found that the induction time and the nucleation rate vary inversely and linearly with , respectively. From their standpoint, the diffusion barrier properties of hydrogen in PuO are not sensitive to the temperature (cooling from to ) and the oxygen content (varying from vacuum to O-rich) during the PuO synthesis. Although all experimental studies so far have definitely demonstrated the distinct roles of PuO and -PuO overlayers, such details as the interaction mechanisms of molecular H with Pu-oxide surfaces and the diffusion behaviors of hydrogen species in Pu oxides have not been clarified, and demand further study from a very fundamental viewpoint. However, given the complexity and hazard listed above, it is extraordinarily difficult to experimentally clarify what happens in the early stages when the H molecules are either physisorbed or chemisorbed, and subsequently the H either dissociate to create reactive species or diffuse into the Pu-oxide overlayers. Thus, the atomic simulations are highly required for the in-depth understanding of the interaction dynamics of H molecule with Pu-oxide overlayers.
In order to address the above-mentioned issues, theoretical schemes have to reasonably describe the ground state properties of Pu oxides at first, and then guarantee that the ab initio molecular dynamics (AIMD) simulations of the interactions between Pu-oxide surfaces and H molecules are methodologically reliable and computationally viable. However, the conventional density-functional theory (DFT) underestimates the strong on-site Coulomb repulsion of the electrons and incorrectly describes PuO as a metal Boettger2002 () instead of a Mott insulator reported by experiments McNeilly1964 (). Fortunately, several beyond-DFT approaches, including the DFT+ Dudarev1998 (), the self-interaction corrected LDA Petit2003 (); Petit2010 (), the hybrid density functional of (Heyd, Scuseria, and Enzerhof) HSE Prodan2005 (), and LDA+dynamical mean-field theory (DMFT) Yin2008 (), have been developed to correct the failures of conventional DFT in calculations of Pu oxides. The effective modification of the pure DFT by the DFT+ formalism has been confirmed widely in studies of physical properties of Pu-oxides (PuO, -PuO and -PuO) SunJCP2008 (); SunCPB2008 (); ZhangP2010 (); Shi2010 (); Shi2012 (); Andersson2009 (); Jomard2008 (); Jomard2011 (); SunJNM2012 (); SunPLA2012 (), such as the insulating ground states, the lattice parameters, the bulk modulus, the phonon spectra and density of states (DOS), the surface stability and chemical activity, and the redox energetics. Note that within the DFT+ framework, both PuO and PuO are described as the antiferromagnetic (AFM) insulator. The calculated AFM ground state of PuO is actually consistent with experimental measurements McCart1981 (); Wulff1988 (). However, unlike the affirmative AFM PuO and UO, PuO is peculiar due to its complex ground-state multiplet, so that strictly speaking its paramagnetic ground state has not been well understood and an AFM exchange Santini1999 (); Colarieti2002 () between Pu ions in PuO is still used to explain the discrepancy between neutron-scattering experiments Kern1999 () and magnetic susceptibility measurements Raphael1968 (), which is beyond the scope of our current work. In spite of all the published papers and increasing attentions focusing on Pu oxides, our theoretical understanding on the basic mechanisms for the hydrogenation of Pu-oxide coated Pu is, to put it mildly, very poor. Especially, little is known regarding the interaction behaviors of H on Pu-oxide surfaces, which is in sharp contrast to the depth and comprehensiveness of researches conducted upon the surface reaction mechanisms of ceria VICARIO2006 (); CHEN2007 (); WATKINS2007 (); ALAM2011 (); Yang2010 (); Paier2013 (). As far as we are aware, Wu et al. Andersson2009 (); Jomard2008 () were the only ones who have studied the interaction between gaseous molecule (HO) and PuO(100)/PuO(110) surfaces, whereas in their literature both the strong correlation effect of Pu- electrons and the van der Waals force between polar HO molecule and PuO surfaces were not taken into account. Recently, Jomard and Bottin Jomard2011 () have via DFT (PBE)+ calculations discussed the influence of electronic correlations on thermodynamic stability of PuO surfaces and pointed out that the O-terminated PuO(111) is the most stable surface and DFT+ method can well describe the electronic structure of the insulating PuO surfaces, which is a prerequisite to study surface reactivity of PuO. In order to thoroughly understand the surface properties of PuO, we have discussed the effects of thickness and O-vacancy on surface stability and chemical activity of PuO based on the DFT+ formalism with the spin-orbit coupling (SOC) effect included, and found that under O-rich conditions, the stoichiometric O-terminated PuO(111) surface is thermodynamically stable and chemically inert, while the O-reducing condition facilitates the formation of on-surface O-vacancy, which is expected to enhance the chemical activities of PuO(111) surface SunJNM2012 ().
As a continuation of our theoretical studies on Pu-surface corrosion SunJNM2012 (); SunPLA2012 (), in this paper, we study the interaction dynamics of H molecule with Pu-oxide surfaces to clarify the different roles of PuO and -PuO overlayers in the hydrogenation of Pu-metal. Specifically, based on the systematic simulations within the vdW-DFT+ framework, we address the following issues: (i) the interaction behaviors such as adsorption, collision, and dissociation of a cold or hot H on stoichiometric and reduced PuO(111) surfaces; (ii) the interaction between H and stoichiometric -PuO(111) surface and the kinetic energy-dependent penetration of H into -PuO(111) overlayer; (iii) the adsorption and diffusion behaviors of H in the -PuO matrix, and the interaction between the H of vibrational excitation and the -PuO during the diffusion. The rest of this paper is organized as follows. The details of our calculations are described in Sec. II. In Sec. III we present and discuss the results. In Sec. IV, we summarize our main conclusions.
Ii Details of Calculation
The spin-polarized vdW-DFT+ calculations are performed using the Vienna ab initio simulation package (VASP) Kresse1996 () with the projector-augmented wave (PAW) method Blo1994 () and a plane-wave cutoff energy of eV. The Pu-67665, O-22, and H-1 electrons are treated as valence electrons, while the remaining electrons were kept frozen as core states in the PAW method. The electronic exchange and correlation are treated within the generalized-gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) functional PERDEW1996 (), and the strong on-site Coulomb repulsion among the Pu-5 electrons is described with the DFT+ scheme formulated by Dudarev Dudarev1998 (). Thus, in this study a Hubbard -like term (=, i.e. the difference between the Coulomb and exchange parameters, hereinafter referred to as ) is added to the PBE functional. In the previous DFT+ studies of PuO, the choices of and parameters are based on an overall agreement between available experimental data and theoretical results as regarding the basic physical properties of the ground state, and the combination of = and = eV has been proved to be reasonable by our previous calculations SunJCP2008 (); SunCPB2008 (); ZhangP2010 (); Shi2010 (); SunJNM2012 (); SunPLA2012 () of PuO and PuO. Furthermore, the DFT+ method has been shown to give substantially similar results to those obtained by hybrid functionals Prodan2005 () with partial exact exchange. Note that the computational demand of DFT+ is close to that of conventional LDA or GGA functionals and much less than that of hybrid functionals, which are computationally too time-consuming to do the AIMD simulations of such surface systems as studied in this paper.
The Pu-oxide surfaces are modeled by finite-sized periodic slab supercells, consisting of a number of oxide layers infinite in and directions and separated in the direction by a vacuum of Å, which is found to sufficiently minimize the interaction between adjacent slabs. Following our previous DFT+ studies of PuO surfaces SunJNM2012 (), in this work we first calculate the surface energy of the stoichiometric -PuO slabs. This quantity is written as
where is the total energy of the fully relaxed -PuO surface slab, is the energy of the reference bulk -PuO with the same number of atoms, and the denominator is the total area of both surfaces of the slab with a finite thickness. Based on extensive test calculations, the of a slab with varying thickness is found to be converged within meV/Å. According to Tasker’s general criteria Tasker1979 () on surface stability of an ionic crystal, the O-terminated -PuO(111) slab, consisting of successive and electrically neutral “O-Pu-O” blocks along the direction, has the “Tasker-type-II” surfaces and its surface energy is calculated to be eV/Å. The -PuO(110) is the typical “type-I” surface stacked with identical neutral planes fulfilling the -PuO stoichiometry with its of eV/Å. The polar “type-III” (001) surface, consisting of oppositely charged planes, is in this work modeled as two oxygen-terminated surfaces with 50% oxygen vacancies to fulfill the stoichiometric formula and to quench the dipole moment of the repeated unit in direction, and its is found to be the largest one ( eV/Å). Thus, for both PuO and -PuO, the O-terminated (111) surface is most stable and is briefly named as the (111) surface if not mentioned differently in this study.
Considering the auto-reduction (AR) of PuO to -PuO driven by the oxygen reducing/poor conditions, in this study, we therefore mainly study three model surfaces of Pu-oxides, corresponding to three typical cases of varying oxygen conditions, namely, stoichiometric/ideal PuO(111) surface (oxygen rich condition), reduced/defective PuO(111) surface with % O vacancies (oxygen reducing), and stoichiometric -PuO(111) surface (oxygen poor). The top views of stoichiometric and reduced PuO(111) surfaces and stoichiometric -PuO(111) surface used in the calculations are shown in Fig. 1. The Brillouin zone (BZ) integration is performed using the Monkhorst-Pack (MP) -point mesh Monkhorst1976 (). Specifically, for both stoichiometric and reduced PuO(111) surfaces, the two-dimensional (2D) cell with a -point MP grid is employed in the static calculations, and a -point MP grid is used in AIMD simulations. Since the ideal stoichiometric -PuO is structurally similar to the 222 fluorite PuO () supercell containing % O vacancy located in the 16 () sites, the surface area of unit -PuO(111) cell is equal to that of the PuO(111) cell. Thus, a -point MP grid is used in static calculations and it is restricted to the -point in AIMD simulations. The slab of -PuO(111) consists of “O–Pu–O” blocks along the direction, including totally plutonium and oxygen atoms, and the -PuO(111) surface model contains “O–Pu–O” blocks with Pu and O atoms in the supercell. In the spin-polarized static calculations of Pu-oxide slabs with the AFM orders of the Pu-sublattice set to be in a simple “ ” alternative manner along the direction, the two “O–Pu–O” oxide bottom blocks are fixed to their calculated bulk positions ( Å and Å) and all other atoms are fully relaxed until the Hellmann-Feynman forces become less than eV/Å. The dipole-dipole interactions along the direction perpendicular to the slab are also corrected, as implemented in the VASP code Makov1995 (); Neugebauer1992 ().
Due to the “Tasker-type-II” character of PuO(111) and -PuO(111) surfaces, the non-coplanar O-anion and Pu-cation may generate surface dipole moment along the direction, especially when the surface O-vacancies exist. In order to detailedly discuss the multiple interaction mechanisms between H and Pu-oxide surfaces, van der Waals (vdW) dispersion forces upon H had better be considered. In this study, we employ the non-local vdW-DFT+ approach Dion2004 (); Klime2010 () to calculate the static potential-energy pathway (PEP) along which a cold H approaches to Pu-oxide surfaces, and to simulate the dynamic interaction between a hot H (with certain initial kinetic energy) and the Pu-oxide surfaces. Moreover, we expect that the static PEP calculations and the AIMD simulations can complement each other to achieve a thorough understanding of the energetics factor and dynamics process determining the distinct roles of PuO and -PuO overlayers in Pu hydrogenation. From the PEP results, we can analyze the binding strength between Pu-oxide and H molecule via the adsorption energy , which is given by
Here, , and are the total energies of the adsorption system, the pure Pu-oxide system, and the H molecule without zero-point vibration, respectively.
The AIMD simulations under the Born-Oppenheimer approximation are performed using the Verlet algorithm with a time step of fs within the micro canonical ensemble (NVE). With the NVE-AIMD, the very fast impingement of H molecule against the Pu-oxide surfaces is simulated to uncover the possible molecular dissociation and estimate the high dissociation barriers. Based on the harmonic approximation and the vibrational energy written as
the ground-state zero-point vibration () and internal vibrational excitation ( or ) of H molecule with a calculated harmonic frequency ( THz) are also considered. In addition, the canonical ensemble (NVT) at several prefixed temperatures through a Nosé thermostat Nose1984 () has also been employed as an alternative approach to discuss the high temperature effect on the state or behavior of H molecule in -PuO. However, the short time scale that can be simulated by the AIMD based on vdW-DFT+ is a major limitation for such rare events as the penetration and diffusion of H molecule in this study, which usually take place on a much longer timescale (of milliseconds or more) than the atomic vibrational timescale of femtoseconds. Here, we attempt to depict such rare events by using the harmonic transition state theory (hTST) with the rate constant given by
where and are the normal mode vibrational frequencies of the initial and transition states, respectively. Here, the energy barrier is defined as the energy difference between the initial state (IS) and transition state (TS), where the TS structures for relevant processes are located by employing the climbing image nudged elastic band (CI-NEB) method HENKELMAN2000 () using at least three images along each pathway. The vibrational frequencies are obtained by diagonalization of the force-constant matrix in Cartesian coordinates (Hessian), where the force constants are obtained from finite differences of the forces with atomic displacements of Å. Then, we can further use the frequencies to characterize whether an optimized stationary point is an initial/minimum state without imaginary frequency or a TS with only one imaginary frequency.
Iii Results and discussion
Since the surface interaction with molecular monomer is very site-specific, as shown in Fig. 1, in this study we consider four high-symmetry target sites (T1, H1, F1, and B1) for a H molecule approaching the stoichiometric PuO(111) surface, and also four sites (T1, T2, H1, and H2) on the reduced PuO(111) surface with an artificial surface oxygen vacancy marked by a solid green circle in Fig. 1. For the stoichiometric -PuO(111) surface shown in Fig. 1, where the visible O vacancies on the first/outmost, the second/subsurface, and the third/deep oxygen layers are marked by big solid, small solid, and small dashed green circles, respectively, we consider three adsorption sites (T1, T3, and F2). Due to the % O vacancies existing inherently in -PuO, one can see that T3 is actually a special site, under which two superposed O vacancies seem to bore a hole in the -PuO(111) slab and may act as an entrance for the smallest diatomic molecule H with its bond-length of only Å.
The H molecule is initially placed over different Pu-oxide surface sites with a height of or Å and the initial orientation at each site is set to be along the , , and axes, respectively. Altogether, we construct high-symmetry initial geometries of H on three Pu-oxide surfaces and name them by adding a corresponding postfix—T1- for example, among which ‘-’ structures are presented in Fig. 1. In addition, we also generate some low-symmetry initial geometries by rotating H with small angles, and find that low-symmetry geometries can relax into the high-symmetry ones after geometry optimization. Thus, in the static PEP calculations, the H-H bond is fixed to be either parellel or perpendicular (i.e., the ‘-’ cases) to the (111) surface when the H molecule drops step by step from its initial location. Since a H molecule on a solid surface has degrees of freedom (-DOF) [see the inset to Fig. 3(e), i.e., bond-length , height , azimuth angles ( and ) of H, and projected positions (, ) of mass center], so that in the static PEP calculations we consider in practice -DOF of H except for . However, during the AIMD simulations of the H molecule with a given kinetic energy impingement against the initially resting/cold Pu-oxide surface, we let both the molecule and the surface layers relax freely.
iii.1 The interaction between H and PuO(111) surfaces
In this subsection, we focus on the interaction mechanisms between H and PuO(111) surfaces, namely, stoichiometric/ideal PuO(111) and reduced/defective PuO(111) surfaces (see Fig. 1).
First of all, the effects of vdW forces upon H molecule from the ideal-PuO(111) surface are probed by an overall comparison between vdW-DFT+ and DFT+ calculated PEPs. We find that the pure DFT+ calculations usually can not capture the physisorption state of H (with the adsorption energy of H underestimated to less than meV); not but what the vdW-DFT+ scheme is able to do that. The improved description of H physisorption by vdW-DFT+ is particularly so when diatomic molecule H is perpendicular to the ideal-PuO(111) surface, for example, the F1- PEPs presented in Fig. 2(a). When H drops from F1- of = Å, the DFT+ PEP indicates nearly an onefold force of repulsion existing between H molecule and ideal-PuO(111) surface, since the variation of relative energy as a function of can be well fitted by . However, the vdW-DFT+ PEP shows the twofold interaction of both repulsion and attraction between H and ideal-PuO(111) surface, and thus the physisorption state of H ( Å and eV) is well discerned mainly due to the vdW attractive force upon the H from the ‘Tasker-type-II’ PuO(111) surface, which is not entirely considered in the pure DFT+ method. We note that the DFT+ calculated of H with the T1- geometry on the ideal-CeO(111) surface is only meV CHEN2007 (), which is also notably underestimated. Since the surface O-vacancy can enhance the chemical activity of the inert ideal-PuO(111) surface by efficiently lowering the surface work function SunJNM2012 (), some DFT+ PEPs seem to be able to give the weak physisorption states of H on reduced-PuO(111) with of only meV more or less, but even so they can not really describe the interaction between H and reduced-PuO(111) because of the missing vdW force in the treatment. Figure 2(b) shows as an example the T2- PEP at reduced-PuO(111). As a matter of fact, the surface O-vacancy will induce the local polarization of ideal-PuO(111) surface and produce a stronger vdW attraction for H. We can see from Fig. 2(b) that the vdW-DFT+ calculated physisorption energy of H near the O-vacancy (namely, Å over T2- site) on reduced PuO(111) surface is about eV instead of meV from pure DFT+, and is somewhat higher of meV than the cases of H on ideal-PuO(111) such as F1- in Fig. 2(a). Therefore, in order to reasonably describe the multiple interaction between H and Pu-oxide surfaces, we choose the vdW-DFT+ scheme instead of pure DFT+ method in the following static PEP calculations and AIMD simulations.
Figure 3(a)-(d) present all the vdW-DFT+ PEPs of H approaching the ideal-PuO(111) surface considered in this work. One can see that H molecules first undergo a weak physisorption state with eV when Å, and then it seems to be hard for H to get further close ( Å) to the inert surface, as the dominant force of repulsion between H and ideal-PuO(111) increases rapidly and notably raises the total energy of these systems, especially when Å. We also find that many H molecules can diverge away from their initial pathways driven by the lateral components of the repulsive force. These surface migrations of H are labeled by their approximate directions in Figs. 3(a), 3(b), and 3(d). Much to our surprise, in Fig. 3(a) H molecules with the T1- and T1- initial geometries can first overcome an energy barrier of eV at Å to migrate to H1 site, and then move towards F1 site. Among all these PEPs with ranging from to Å, none but the H1- PEP presents the dissociative adsorption of H molecule when Å [see Fig. 3(b)], with one H atom binding to the subsurface Pu-cation and the other leaving over the surface, and the corresponding energy barrier is much higher than eV, which indicates that the dissociative adsorption of H on ideal-PuO(111) surface should be an extremely endothermic and rare event. However, after artificially breaking the H-H bond, the adsorption energy of H atoms binding to the surface O anions is calculated to be eV, which is really a very exothermic reaction to hydroxylate the PuO(111) surface with H atoms. So that the dissociative adsorption should be the most stable and favorable state of H molecule on ideal-PuO(111) surface, which is also found to be similar to the strongly exothermic dissociation of H on the ideal-CeO(111) surface VICARIO2006 (); CHEN2007 (); WATKINS2007 () with the of H atoms calculated by the DFT+ to be about eV. And the two OH groups formed on the oxide surface can help to reduce the nearby Pu/Ce cations to state.
Given the existing high barrier, the most controversial issue about the dissociative adsorption of H actually rests with the reliable method to effectively simulate this extremely rare event and reasonably estimate the corresponding dissociation barriers. Recently, Alam et al. ALAM2011 () have simulated the interaction process of hot H molecule with the ideal-CeO(111) and (110) surfaces by using the ultra accelerated quantum chemical molecular dynamic (UA-QCMD) method, which is based on the tight-binding quantum chemistry and the classical MD programs. In their UA-QCMD simulations, the mean velocity of several km/s was given to the H molecule. Here, we perform comprehensive NVE-AIMD simulations to fulfil this task. From Fig. 3(e) we can see that if its initial kinetic energy is set at eV, the initially nonrotating and nonvibrating H from T1- of = Å will bounce back to the vacuum soon and even can not touch the ideal-PuO(111) surface. When the initial of H is less than eV, the diffraction behaviors of molecular H seem to be similar to the case presented in Fig. 3(e) without reference to its initial geometries. By continually elevating the initial of H up to eV, we find that most of H molecules parallelly hitting the ideal-PuO(111) surface can be broken into two hot H atoms owning high velocities, which are expected to be captured by the surface O anions and eventually form two surface hydroxyl groups, i.e., the above-mentioned stable dissociative adsorption state of H. Whereas, among all AIMD trajectories, the -dependent dissociation behavior of H from T1- is an interesting exception and is thus presented in Fig. 3(f). It shows in Fig. 3(f) that (i) when the initial eV, H can not touch the surface O-anion and bounce straightly without arousing the molecular zero-point vibration, according to the inset of Fig. 7(b); (ii) when its reaches eV, the internal vibrational excitation of H molecule emerges [ in Eq. (3)] after the head-on collision with the O-anion; (iii) until reaches eV, molecular H is dissociated by its violent collision with surface O-anion. However, the products are not two hot H atoms but a hot HO molecule released to the vacuum and an O-vacancy left on the surface. This surface reduction process has also been reported by the above-mentioned UA-QCMD simulations of high-energy colliding H molecule on CeO(111) and (110) surfaces ALAM2011 (). The energy barrier of collision induced dissociation process in this study is estimated to be eV based on the electronic free energy of this system in Fig. 3(f), which is close to the bond energy ( eV) of a free-standing H molecule and indicates that the inert ideal-PuO(111) surface can not act as the catalyst for H dissociation. Since moisture can strongly enhance the further oxidation of Pu metal Colmenares1975 (); Has2000 (), one can see that both the surface hydroxylation of PuO by H and the interaction between HO and PuO surface are important topics, which will be investigated in our next work. In order to search out the possible lowest dissociation barrier, different incident conditions (namely, the initial and the angles of incidence in some trajectories) are considered. Figure 3(g) presents two limiting cases. When is eV, the AIMD trajectory of H started from the horizontally shifted T1- geometry (with one H-atom right above the surface O-anion) turns out to be the most possible dissociative adsorption path with the lowest barrier of eV. However, the trajectory started from F1- indicates that H will never dissociate even if we would further elevate the in the NVE-AIMD simulation. Thus, the collision-dissociation behavior of H is very site specific on ideal-PuO(111) surface and depends heavily on the incident condition of H. Overall, all the AIMD results well support the static PEPs and leads to a conclusion that the ideal-PuO(111) surface is really difficult of access for H molecule.
As the format of Fig. 3, Fig. 4 presents the interaction behaviors of H molecule with the reduced-PuO(111) surface, based on all static PEP results and some representative AIMD trajectories plotted in the left- and right-hand columns, respectively. From the vdW-DFT+ PEPs, we can see that the molecular physisorption states are ubiquitous for H on reduced-PuO(111) surface when Å, and the physisorption states of H over the surface O-vacancy (i.e., T2 site) seem to be more stable than others, so that when further depressing the nearby adsorbed H molecules, most of them tend to stride over certain barriers (less than eV) to migrate to T2 site. Note that H seems to get much closer to the reduced-PuO(111) surface than to the ideal-PuO(111) surface mainly because of the 25% O-vacancy on the former. However, when Å, it becomes more and more difficult for H to get further closer to the surface and all these PEPs with ranging from to Å do not present the dissociative adsorption of H, which is to some extent similar to the behavior of H on ideal-PuO(111) surface.
Figure 4(e) shows the total kinetic energy of the system and the bond-length of H with its initial of eV along AIMD trajectory started from T2- of = Å. We can see that during the initial fs, the vdW force of attraction causes the steady augmentation of , and then the increasing force of repulsion exceeds the vdW force and pushes H back to the vacuum. Meanwhile, the zero-point vibration of H is aroused. Comparing Fig. 3(e) with Fig. 4(e), we can see that the interaction of H with ideal- and reduced-PuO(111) surfaces are actually unalike, and the surface O-vacancy modulates the interaction behavior between H and PuO(111) surfaces by changing the surface atomic and electronic structure properties. Based on enough NVE-AIMD trajectories of H impinging against the reduced-PuO(111) surface, we also try to reveal the possible dissociation behavior and search out the lowest dissociation barrier of H on reduced-PuO(111) surface. We find that several contact collisions between H and reduced-PuO(111) can not occur until the initial of H reaches eV, which is consistent with most of the PEPs in the left-hand column of Fig. 4. When the of H reaches eV, the collision induced dissociation of H will take place to yield two hot H atoms in most AIMD trajectories. Whereas, the trajectory started from T1- still make an interesting exception, which is revealed in Fig. 4(f). From Fig. 4(f), we can see that when eV, after the quick migration [from T1 to T2 site driven by the repulsive force, see Fig. 4(a)] and the subsequent collision (at T2 site), the rebounded H molecule seems to be at the second vibrational excitation state [ in Eq. (3)]. When reaches eV, without the lateral migration, the H molecule straight hits the O-anion target to break its bond, and soon two hot H atoms pull the O-anion out of the surface to yield a hot HO molecule released to the vacuum. Upon that, when starting from T1-, the collision-dissociation dynamics of H on reduced-PuO(111) surface is similar to that occurred on the ideal-PuO(111) surface. However, the corresponding dissociation barriers are quite different, namely, and eV for H on reduced- and ideal-PuO(111) surfaces, respectively. For the dissociation possibility of H (with eV) on reduced-PuO(111) surface, Fig. 4(g) presents two limiting cases as follows: (i) The AIMD trajectory started from the horizontally shifted T1- geometry [as in Fig. 3(g)] gives the lowest energy barrier of eV, much lower than the lowest barrier eV in Fig. 3(g); (ii) Whereas, the trajectory started from T2- (or T2-, not shown here) indicates that H will never dissociate even with a higher . Given the incidental migrations of H from other sites to T2 (see the PEP results in Fig. 4), we can conclude that the reduced-PuO(111) surface is also difficult of access for H molecule, let along the collision induced dissociation.
If H is barely able to react with PuO(111) surfaces, as the smallest diatomic molecule, is there any possibility for it to penetrate vertically into the PuO(111) overlayers? As the continuations of the vdW-DFT+ PEPs in Fig. 2(a) and 2(b), the PEPs for H penetration into ideal- and reduced-PuO(111) surfaces from the corresponding physisorption states are presented in Fig. 5(a) and 5(b), respectively. We can see from Fig. 5(a) that the relative energy of the system increases fleetly during the penetration of H molecule into ideal-PuO(111) surface, wherein two insets show atomic structures of two reference points (namely, and Å). Since H touches the surface ( Å), its bond ( Å) is stretched notably, and once the whole H enters the surface ( Å), the H-H bond is snapped ( Å) with one H atom binding to the subsurface Pu-cation and the other segregating out of the surface. Upon that, the penetration of H molecule into ideal-PuO(111) surface is too endothermic to occur under ambient conditions because the stoichiometric PuO is so close-grained that the H needs to overcome a high energy barrier of more than eV to be broken before squashing in the PuO overlayer. For the penetration of H molecule into reduced-PuO(111) surface, we can see from Fig. 5(b) that the PEP curve of T2- is stepwise, namely, the relative energy first increases rapidly during molecular entry into the surface ( Å) with the penetration barrier of eV, then increases slowly for the subsurface migration of H ( Å), and ultimately increases sharply when H approaches and binds to the Pu-cation. Once the H reacts with the Pu-cation below, it will dissociate (see the lower inset), and soon the upper H-atom will be segregated out and stay at the site of surface O-vacancy binding with three Pu cations nearby (see the upper inset). Before molecular dissociation, interestingly, the H-H bond is first stretched slightly while H is hauled step by step into the surface, and then relaxes back to the equilibrium state after H entry into the surface, which differs from the variation of in Fig. 5(a) and indicates that just one surface O-vacancy can let the smallest H molecule enter the PuO surface layer. But the further penetration of H in PuO overlayer ( Å) is found to be very difficult since H meets with the close-grained PuO(111) subsurface layer subsequently.
Thus, according to our current comprehensive study of the interaction properties between H molecule and PuO(111) surfaces, it is found that the close-grained PuO overlayer can efficiently prevent the penetration and diffusion of hydrogen (into the Pu-oxide film), and protect the underlying Pu-metal from hydrogenation.
iii.2 The interaction between H and -PuO(111) surface
In this subsection, we turn to study the interaction mechanisms between H and the stoichiometric -PuO(111) surface. Note that in addition to the formation of % oxygen vacancies, there is an interesting volume expansion (% in theory and % in experiment) during the PuO¡úŠÁ-PuO isostructure reduction mainly due to the strong on-site Coulomb repulsion among the Pu- electrons SunPLA2012 (). Thus, compared to the close-grained and smooth PuO(111) surface, the -PuO(111) seems to be a porous slab with % native O vacancies located on every O-layer, and its surface is not smooth any more.
Four representative PEPs of H molecule approaching -PuO(111) surface are shown in Fig. 6(a). One can see that when H molecule is dropped from T1- and F2- geometries, the corresponding PEPs are similar to those PEPs of H approaching ideal-PuO(111) surface [see Figs. 3(a)-(d)], and thus the T1 and F2 sites on -PuO(111) surface are also very difficult of access. However, when starting from T3- or T3-, after the weak physisorption states, H molecule can further be pulled into the -PuO(111) overlayer via the surface O-vacancy, which is similar to the case of H entry into the reduced-PuO(111) surface [see Fig. 5(b)]. The energy barrier of H molecule penetration into -PuO reported by the PEPs is eV less than the barrier of eV for H entry the reduced-PuO(111). In order to accurately obtain the energy barrier and the rate constant of H penetration, we first use the CI-NEB method HENKELMAN2000 () with seven images to search out the transition-state (TS) from the minimum energy pathway (MEP), and then calculate the vibrational frequencies of H molecule at the initial physisorption-state (IS) and the TS. The MEP result and the orientations of H molecule during the penetration process are plotted in Fig. 6(b). One can see that there are two TSs for H molecule first entry (TS-1) and then diffusion (TS-2) in the -PuO overlayer, respectively. And H at TS-1 almost horizontally locates Å below the surface-layer (with the penetration barrier of eV) but H at TS-2 tends to be perpendicular to the surface (with the diffusion barrier of eV), indicating that H seems to be able to rotate freely in the tunnel made of the vertically distributed O vacancies.
Based on the vibrational frequencies of H listed in Table I, we can see that H at the IS vibrates, nearly like a free molecule, with mainly a symmetric stretching vibration ( THz). Whereas, H at both TSs vibrates somewhat complicatedly under the influence of -PuO matrix and has one imaginary frequency. At the room temperature of K, the of H is estimated to be /s, which is really small mainly due to the relatively high penetration barrier ( eV) of H entry the stoichiometric -PuO(111) overlayer. However, in practice the Pu-oxide overlayer on Pu-metal is really far from being an ideal (defect-free) PuO or PuO, and other defects such as microcracks or some forms of impurity inclusions in the overlayer can prominently promote the probability of H penetration. Therefore, it is certainly one challenging task to predict and control the length of induction period in Pu-hydrogenation. Here, at the atomic and molecular level, we just focus on understanding the different roles of ideal PuO and -PuO overlayers in Pu-hydrogenation, and at this point our theoretical results have well revealed the major difference between PuO and -PuO, namely, PuO is airtight but -PuO is porous for the smallest H molecule.
Furthermore, the AIMD simulations are also performed to depict the interaction dynamics for H impingement against the -PuO(111) surface. Figure 6(c) shows the representative NVE-AIMD trajectories started from the T1- and T3- with Å, where the evolutions of height and bond-length of H (upper panel), the relative electronic free-energy , and the kinetic energy of the system (lower panel) are shown. We can see that even when the initial reaches eV, H seems to be unable to touch the T1 site on -PuO(111) surface ( Å) and soon bounces back with the zero-point vibration, which is very unlike the dissociation of H (with of eV) on reduced-PuO(111) surface [see Fig. 4(f)] but is like the diffraction behavior on ideal-PuO(111) surface [see Fig. 3(f)]. The reasons are that (i) the -PuO(111) surface is not smooth with many protrusions such as the T1 site, and actually H can contact the O-anion at T1 site; and (ii) on the other hand it will take much more energy to dig an O-anion out of the surface since further reduction of -PuO is very difficult SunPLA2012 (). However, the H molecule with an initial of eV can stride over an energy barrier of eV to penetrate via T3 site into the -PuO(111) overlayer, and then diffuses easily along the devious alleyway made of the continuous distributed O vacancies. During the H penetration and diffusion, the internal vibration of H molecule is fitful and is especially so after the first collision with Pu-cation ( fs). Although the inhalation of a H molecule in the -PuO(111) slab is endothermic (see of the system), the hot H molecule can release the redundant kinetic energy and change the direction of movement via successive collisions with -PuO matrix, and as a result H seems to be neither escape out of nor dissociate in -PuO(111) during the simulation process of ps.
In order to gain detailed insights into the adsorption states and the diffusion behaviors of H molecule in the -PuO(111) matrix, we first carry out systematic total-energy calculations to search for favorable sites for H to stay at. Then, based on the AIMD simulations within canonical ensemble (NVT-AIMD), we discuss the temperature effect on the state or behavior of H in -PuO. Figure 8(a) shows four relatively favorable molecular mass-center sites A (=-) for H adsorption, and the adsorption energies of H at these four sites are listed in Table II, where H molecule can rotates slightly around the four sites and the A site seems to be a little more favorable with a positive (endothermic) adsorption energy of eV. Subsequent NVT-AIMD simulations at three temperatures (, , and K) show that the molecules just moves around their initial adsorption sites with the zero-point vibration and do not diffuse or dissociate within the isothermal period of ps. The atomic adsorption energy is calculated to be endothermic ( eV/atom), which indicates that the dissociative adsorption of H in -PuO should be more endothermic than the molecular adsorption. Thus, in contrary to the slightly exothermic adsorption of H atom to form OH group in bulk PuO, H molecule and even H atom can not easily react with the -PuO matrix since the reduction of -PuO is much more difficult than that of PuO SunPLA2012 (). The exothermic adsorption of H atom in bulk PuO is also found to be similar to that of H atom in CeO SOHLBERG2001 (), which can reduce the cation in both fluorite-structured oxides.
According to Table II, the adsorption properties of H molecules at four sites seem to be close to each other, and their temperature-dependent behaviors in -PuO are also very similar. For the convenience of depiction, we just plot in Fig. 7(a) the atomic- and electronic-structure results of the A-site adsorption system. We can see that the H molecule actually stays in suspension at A site, which is near the location of an O-anion in the isostructure PuO supercell, i.e., the native O-vacancy site in -PuO. The charge density distribution in Fig. 7(a) indicates that there is no evidence of chemical binding of H molecule to the nearby Pu-cation or O-anion, and there seems to be enough space for the movement or diffusion of H molecule in the -PuO matrix. Furthermore, in the NVE-AIMD, we first consider the internal vibrational excitation of H molecule at A site [with in Eq. (3)] and then let the system relax freely, as is shown in Fig. 7(b). We find that during the initial fs, the H vibrates drastically at A site. Then, this hot H molecule touches the nearby Pu-cation at fs, and meanwhile transforms most of the vibration energy into the translation and rotation energy. Soon after that, H migrates quickly to the A site and transfers its kinetic energy towards the vibration of the -PuO matrix via the successive collisions with the matrix within fs. Eventually ( fs), the H molecule at A site comes back to the zero-point vibration state after releasing most of its kinetic energy. The displacement of H as a function time ( ps) plotted in Fig. 7(b) indicates that the local movement (AA) of H is easy but the long-distance diffusion through the four sites seems to be very time-consuming, which can not be directly studied by the AIMD simulations. Overall, this interesting example reveals that the H molecule tends to keep its ground state with the zero-point vibration.
Figure 8(a) shows that the four sites in a crystal cell of -PuO are actually connectible by the O-vacancy distributed continuously in the -PuO matrix, so that if given enough time, H can penetrate through the -PuO overlayer and the diffusion of H in -PuO seems to be continuous in a stepwise way. By using the CI-NEB method, we search out four TSs [with their sites labeled by T (=-) in Fig. 8(a)] from the calculated minimum energy pathways (MEPs) between adjacent sites, which are plotted in Fig. 8(b). We can see from Fig. 8 that (i) the whole diffusion pathway is a close loop and consists of four short routes in the periodical -PuO matrix; and (ii) when diffusing along the linear routes (from A to A and from A to A, and vice versa), the H can keep parallel and vertical to the routes and the corresponding lowest diffusion barriers are eV (AA) and eV (AA); whereas (iii) when migrating along the meandering routes (A A and A A). The corresponding rightward diffusion barriers are doubled to be and eV because H needs to rotate during the migration courses and swerve at the transition states (T and T). Furthermore, we calculate the vibrational frequencies of H at four adsorption sites and four transition states, with the major stretching vibration frequencies listed in Table II. Then, following Eq. (4), we figure out the rate constant () of H diffusion in bulk -PuO and list the calculated results of at room temperature ( K) in Table II. We can see that the of H diffusion from a metastable state to a relative stable state is always larger than the of the reverse diffusion, and the probabilities for H diffusion along linear routes (A A and A A) are about five or six orders of magnitude larger than the case of H migration along meandering routes, indicating that the diffusion kinetics of H in -PuO is mainly determined by the thermodynamical parameter diffusion barrier.
Assuming sequential first-order kinetics, the overall diffusion rate for H in the -PuO crystal cell can be simply estimated as , where is the rate constant of the th step diffusion. At K, the rightward (AA) of H in -PuO is /s, and the leftward one is /s, indicating that the H migration along meandering routes is the rate-determining step for the diffusion kinetics of H in -PuO. Considering the small of /s for H penetration into -PuO(111) overlayer, the diffusion of H in bulk -PuO seems to be much easier since the is at least ten orders of magnitude larger than the . Thus, under oxygen-poor conditions, the course of H entry the -PuO overlayer is the final rate-determining step for H penetration through the -PuO to reach the Pu-metal. In order to quantificationally determine the rate constants of H penetration, the practical surface configuration of -PuO overlayer is a critical factor. But the complex surface conditions of oxide-coated Pu are undoubtedly not as simple as the stoichiometric -PuO(111) surface considered in this work, so that we just take the first step towards understanding the hydrogenation of Pu-oxide coated Pu metal.
In summary, based on the vdW-DFT+ scheme, we have explored the different roles of the Pu-oxide overlayers in the hydrogenation of Pu underlayer. Three model surfaces of Pu-oxides are considered by varying oxygen conditions. We have found that the physisorption state of H on PuO(111) surfaces can not be captured by pure DFT+ calculations until the van der Waals correction is taken into account. We have shown that the surface O-vacancy induces a local polarization and produce a stronger van der Waals attraction for H, thus the physisorption energy of H near the O-vacancy on the reduced PuO(111) surface becomes higher by meV than that on the ideal PuO(111) surface. However, the physisorbed H molecules can not further get close to PuO(111) surfaces. In agreement with the static PEP results, our AIMD simulations have shown that the H molecule will bounce back to the vacuum when its initial kinetic energy is not high enough. Until the initial reaches eV ( eV), the H molecule parallel to stoichiometric (reduced) PuO(111) surface can be broken. Interestingly, these collision-induced dissociation behaviors of H have been found to be very site-specific. For instance, two H atoms can chemically bind with two O-anions to form hydroxyl groups and can also react with one O-anion to produce a HO molecule. Although the dissociative adsorption of H is very exothermic and can reduce PuO(111) surfaces, it should be an extremely rare event due to the high energy barriers ( and eV for H on the ideal- and reduced-PuO(111) surfaces, respectively). Therefore, the close-grained PuO overlayer acts as a diffusion barrier to control the permeation of hydrogen and lengthen the initiation time of Pu-hydrogenation.
As a product of isostructure reduction of PuO, -PuO has % oxygen vacancies located along 111 diagonals and thus seems to be a porous slab. Our AIMD study have directly revealed that the H molecule can overcome an barrier of eV and directly penetrate into the -PuO(111) film via the O vacancies. By examining the temperature effect and the internal vibrational excitations of H, we have also provided a detailed insight into the interaction dynamics between H and -PuO. In -PuO, the hot H molecule prefers to release its energy via successive collisions and come back to its ground state with zero-point vibration. The diffusion behavior of H in -PuO is also systematically investigated by searching out the minimum diffusion paths and calculating the diffusion rate constants of H. Our results are consistent with the general experimental observations, and come to the conclusion that the PuO overlayer can hinder the hydrogen penetration, provided the isostructure reduction of PuO to -PuO is efficiently suppressed.
Acknowledgements.This work was supported by the Foundations for Development of Science and Technology of China Academy of Engineering Physics under Grants Nos. 2010B0301048, 2011A0301016, and 909-07, and partially supported by the NSFC under Grants No. 11004012, No. 51071032, No. 11205019. and No. 11275032.
- (1) S. Y. Savrasov and G. Kotliar, Phys. Rev. Lett. 84, 3670 (2000).
- (2) L. Petit, A. Svane, Z. Szotek, W. M. Temmerman, Science 301, 498 (2003).
- (3) K. T. Moore and G. van der Laan, Rev. Mod. Phys. 81, 235 (2009).
- (4) J. X. Zhu, R. C. Albers, K. Haule, G. Kotliar and J. M. Wills, Nat. Commun. 4, 2644 (2013).
- (5) C. A. Colmenares, Prog. Solid St. Chem. 9, 139 (1975).
- (6) J. M. Haschke, Los Alamos Science 26, 253 (2000).
- (7) M. T. Butterfield, T. Durakiewicz, I. D. Prodan, G. E. Scuseria, E. Guziewicz, J. A. Sordo, K. N. Kudin, R. L. Martin, J. J. Joyce, A. J. Arko, K. S. Graham, D. P. Moore, and L. A. Morales, Surf. Sci. 600, 1637 (2006).
- (8) H. A. Wriedt, Bull. All. Ph. Dia. 11(2), 184-202 (1990).
- (9) M. T. Butterfield, T. Durakiewicz, E. Guziewicz, J. Joyce, A. Arko, K. Graham, D. Moore, and L. Morales, Surf. Sci. 571, 74 (2004).
- (10) J. M. Haschke, T. H. Allen, and L. A. Morales, Science 287, 285 (2000).
- (11) P. A. Korzhavyi, L. Vitos, D. A. Andersson and B. Johansson, Nature Mater. 3, 225 (2004).
- (12) T. Gouder, A. Seibert, L. Havela, and J. Rebizant, Surf. Sci. 601, L77 (2007).
- (13) H. G. García Flores, P. Roussel, D. P. Moore, D. L. Pugmire, Surf. Sci. 605, 314 (2011).
- (14) P. Morrall, S. Tull, J. Glascott, P. Roussel, J. Alloys Comp. 444-445, 352-355 (2007).
- (15) L. N. Dinh, J. M. Haschke, C. K. Saw, P. G. Allen, W. McLean II, J. Nucl. Mater. 408, 171 (2011).
- (16) G. W. McGillivray, J. P. Knowles, I. M. Findlay, M. J. Dawes, J. Nucl. Mater. 412, 35 (2011).
- (17) C. E. Boettger and A. K. Ray, Int. J. Quantum Chem. 90, 1470 (2002).
- (18) C. McNeilly, J. Nucl. Mater. 11, 53 (1964).
- (19) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
- (20) L. Petit, A. Svane, Z. Szotek, W. M. Temmerman, and G. M. Stocks, Phys. Rev. B 81, 045108 (2010).
- (21) I. D. Prodan, G. E. Scuseria, J. A. Sordo, K. N. Kudin, and R. L. Martin, J. Chem. Phys. 123, 014703 (2005).
- (22) Q. Yin and S. Y. Savrasov, Phys. Rev. Lett. 100, 225504 (2008).
- (23) B. Sun, P. Zhang, and X.-G. Zhao, J. Chem. Phys. 128, 084705 (2008).
- (24) B. Sun, and P. Zhang, Chin. Phys. B 18, 1364 (2008).
- (25) P. Zhang, B.-T. Wang, and X.-G. Zhao, Phys. Rev. B 82, 144110 (2010).
- (26) G. Jomard, B. Amadon, F. Bottin, and M. Torrent, Phys. Rev. B 78, 075125 (2008).
- (27) D. A. Andersson, J. Lezama, B. P. Uberuaga, C. Deo, and S. D. Conradson, Phys. Rev. B 79, 024110 (2009).
- (28) H. Shi, M. Chu, and P. Zhang, J. Nucl. Mater. 400, 151 (2010).
- (29) H. Shi, and P. Zhang, J. Nucl. Mater. 420, 159 (2012).
- (30) G. Jomard, and F. Bottin, Phys. Rev. B 84, 195469 (2011).
- (31) B. Sun, H. Liu, H. Song, G. Zhang, H. Zheng, X.-G. Zhao, and P. Zhang, J. Nucl. Mater. 426, 139 (2012).
- (32) B. Sun, H. Liu, H. Song, G. Zhang, H. Zheng, X.-G. Zhao, and P. Zhang, Phys. Lett. A 376, 2672 (2012).
- (33) B. McCart, G. H. Lander, and A. T. Aldred, J. Chem. Phys. 74, 5263 (1981).
- (34) M. Wulff and G. H. Lander, J. Chem. Phys. 89, 3295 (1988).
- (35) P. Santini, R. Lémanski, and P. Erdös, Adv. Phys. 48, 537 (1999).
- (36) M. Colarieti-Tosti, O. Eriksson, L. Nordström, J. Wills, and M. S. S. Brooks, Phys. Rev. B 65, 195102 (2002).
- (37) S. Kern, R. A. Robinson, H. Nakotte, G. H. Lander, B. Cort, P. Watson, and F. A. Vigil, Phys. Rev. B 59, 104 (1999).
- (38) G. Raphael and R. Lallement, Solid State Commun. 6, 383 (1968).
- (39) G. Vicario, G. Balducci, S. Fabris, S. Gironcoli, and S. Baroni, J. Phys. Chem. B 110, 19380 (2006).
- (40) H.-T. Chen, Y. M. Choi, M. Liu, and M. C. Lin, ChemPhysChem 8, 849 (2007).
- (41) M. B. Watkins, A. S. Foster, and A. L. Shluger, J. Phys. Chem. C 111, 15337 (2007).
- (42) M. K. Alam, F. Ahmed, R. Miura, A. Suzuki, H. Tsuboi, N. Hatakeyama, A. Endou, H. Takaba, M. Kubo, and A. Miyamoto, Catal. Today 164, 9 (2011).
- (43) Z. Yang, Q. Wang, S. Wei, D. Ma, and Q. Sun, J. Phys. Chem. C 114, 14891 (2010).
- (44) J. Paier, C. Penschke, and J. Sauer, Chem. Rev. 113, 3949 (2013).
- (45) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
- (46) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
- (47) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (48) P. W. Tasker, Solid State Phys. 12, 4977 (1979).
- (49) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
- (50) G. Makov, M. C. Payne, Phys. Rev. B 51, 4014 (1995).
- (51) J. Neugebauer, M. Scheffler, Phys. Rev. B 46, 16067 (1992).
- (52) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
- (53) J. Klimeš, D. R. Bowler, and A. Michaelides, J. Phys.: Condens. Matter 22, 022201 (2010); J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).
- (54) S. Nosé, J. Chem. Phys. 81, 511 (1984); S. Nosé, Prog. Theor. Phys. Suppl. 103, 1 (1991).
- (55) G. Henkelman, B. Uberuaga, H. Jonsson, J. Chem. Phys. 113, 9901 (2000).
- (56) K. Sohlberg, S. T. Pantelides, and S. J. Pennycool, J. Am. Chem. Soc. 123, 6609 (2001).