Home
Search
Collections
Journals
About
Contact us
My IOPscience
An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 J. Phys.: Condens. Matter 26 105501 (http://iopscience.iop.org/0953-8984/26/10/105501) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 155.97.178.73 This content was downloaded on 01/07/2014 at 15:50
Please note that terms and conditions apply.
:
Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 229501 (1pp)
doi:10.1088/0953-8984/26/22/229501
Corrigendum: An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach (2014 J. Phys.: Condens. Matter 26 105501) A E Thompson, B Meredig and C Wolverton Department of Materials Science and Engineering, Northwestern University, Evanston, IL 60208, USA Received 28 April 2014 Accepted for publication 28 April 2014 Published 14 May 2014
We would like to add the following acknowledgments to our paper: The authors would like to acknowledge funding from the Department of Energy under NERI-C grant no. DE-FG07-07ID1489, and under grant no. DE-F602-07ER46433. BM also acknowledges the National Defense Science and Engineering Graduate Fellowship.
0953-8984/14/229501+1$33.00
1
© 2014 IOP Publishing Ltd Printed in the UK
Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 105501 (9pp)
doi:10.1088/0953-8984/26/10/105501
An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach Alexander E Thompson, Bryce Meredig and C Wolverton Department of Materials Science and Engineering, Northwestern University, Evanston, IL 60208, USA E-mail: [email protected] Received 30 August 2013, revised 4 November 2013 Accepted for publication 3 January 2014 Published 19 February 2014
Abstract
We have created an improved xenon interatomic potential for use with existing UO2 potentials. This potential was fit to density functional theory calculations with the Hubbard U correction (DFT + U ) using a genetic algorithm approach called iterative potential refinement (IPR). We examine the defect energetics of the IPR-fitted xenon interatomic potential as well as other, previously published xenon potentials. We compare these potentials to DFT + U derived energetics for a series of xenon defects in a variety of incorporation sites (large, intermediate, and small vacant sites). We find the existing xenon potentials overestimate the energy needed to add a xenon atom to a wide set of defect sites representing a range of incorporation sites, including failing to correctly rank the energetics of the small incorporation site defects (xenon in an interstitial and xenon in a uranium site neighboring uranium in an interstitial). These failures are due to problematic descriptions of Xe–O and/or Xe–U interactions of the previous xenon potentials. These failures are corrected by our newly created xenon potential: our IPR-generated potential gives good agreement with DFT + U calculations to which it was not fitted, such as xenon in an interstitial (small incorporation site) and xenon in a double Schottky defect cluster (large incorporation site). Finally, we note that IPR is very flexible and can be applied to a wide variety of potential forms and materials systems, including metals and EAM potentials. Keywords: density functional theory, potential fitting, molecular dynamics, UO2 (Some figures may appear in colour only in the online journal)
1. Introduction
dynamics simulations [2, 3, 10–16]. There have been many interatomic potentials parameterized for UO2 [12, 13], but only a few for xenon in UO2 , despite the crucial role of xenon in the failure of nuclear fuels. In this paper we consider two recently-developed empirical Xe–UO2 models as well as our own potential, generated using a new genetic algorithm-based approach which we call iterative potential refinement (IPR). The first xenon potential we consider was proposed by Geng et al [17] who developed Xe–O and Xe–U interactions to use with the UO2 potential of Basak et al [18] in order to examine the properties of xenon bubbles in UO2 . For simplicity, we will refer to the potentials
There is much interest in the behavior of xenon in uranium dioxide (UO2 ) due to the fact that xenon atoms, which are created as a consequence of uranium fission, cluster into bubbles that can swell the host UO2 material, degrading its properties and causing it to fail. There have been many studies of xenon in UO2 using density functional theory (DFT) [1–9], which can achieve good defect energetics, but DFT is too computationally demanding for large systems such as xenon bubbles in UO2 . Interatomic potentials, on the other hand, enable calculations of larger systems and long molecular 0953-8984/14/105501+09$33.00
1
c 2014 IOP Publishing Ltd
Printed in the UK
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Table 1. Potential parameters of uranium dioxide and xenon potentials. The xenon is neutrally charged in the potentials.
Uranium dioxide–xenon potentials Basak A (eV) U–U O–O U–O U O
294.64 1633.0051 693.6487 Charge (e) +2.4 −1.2
Buckingham ρ (Å) c (eV Å6 ) 0.32702 0.32702 0.32702
Geng
Born–Mayer A (eV) ρ (Å)
Xe–U Xe–O
4887.7
Morelon
Buckingham A (eV) ρ (Å)
U O
11272.6 566.498 Charge (e) +3.227252 −1.613626
0 3.94879 0
De (eV)
Morse a (Å−1 )
r0 (Å)
0.57719
1.65
2.369
Lennard-Jones E (eV) σ (Å)
0.415
0.1363 0.42056
Chartier
Born–Mayer A (eV) ρ (Å)
Xe–U Xe–O
478875.11 94.489
0.0099
2.5
c (eV Å6 )
Buckingham four-range rmin1 (Å) cut1 (Å) rmin2 (Å)
cut2 (Å)
0
2.6
134 0
1.2
2.1
0.179 0.58
Oxides are often described well with Buckingham potentials: ri j c UiBuck = A exp − (1) − 6 j ρ ri j
by the name of the first author. The Xe–O potential was fit to XeO3 in a hypothetical (i.e., non-ground-state) CuAu3 structure, and the Xe–U was fit to the same Born–Meyer form as Jackson et al [19]. The parameters of Jackson et al were used as a starting point and were refined with ab initio calculations. However, instead of fitting to xenon in a complete UO2 lattice, the interaction was fit to a ‘virtual structure’. This virtual structure is constructed from a 12 atom unit cell of UO2 in the fluorite structure, but with three of the UO2 formula units removed. The remaining UO2 formula unit and xenon atom are aligned in the [111] direction. The Basak potential was parameterized against thermal expansion of UO2 and achieves good high temperature properties such as thermal expansion and isothermal compressibility [18]. The second xenon potential was proposed by Chartier et al [20] to use in conjunction with the UO2 potential of Morelon et al [21]. The Chartier potential was parameterized by fitting to DFT calculations of four xenon-containing defects (xenon in the following sites: interstitial, uranium vacancy, oxygen vacancy, and a Schottky defect cluster). The Morelon potential was fit to experimental defect and migration energies of UO2 [21]. The Chartier potential treats high-energy defects such as a xenon atom in an interstitial (Xei ) accurately, but overestimates the defect energy of xenon in small vacancy clusters, like xenon in a Schottky defect (XeSD) [21]. We have previously used these potentials to study different diffusion pathways of xenon in tetravacancies in UO2 [22]. Table 1 shows the parameterizations of these various UO2 –Xe potentials. These potentials take on various forms.
where r is the separation between atom I and atom j. A variation on the Buckingham is sometimes used to avoid an unphysical attraction at short separations between atoms by introducing different functional forms at discrete separations in the Buckingham four-range potential: ri j A exp − if rmin1 < ri j ≤ cut1 ρ fifth order polyomial if cut < r ≤ r 1 ij min2 UiBuck4 = (2) j third order polyomial if r < r ≤ cut min2 i j 2 C if cut2 < ri j ≤ rmax . − r 6 ij The Morse potential incorporates some covalent behavior in the interactions between atoms: UiMorse = De [1 − exp(−a(ri j − r0 ))]. j
(3)
The xenon potentials are fairly simple and include the Lennard-Jones potential: " 6 # 12 σ σ UiLJ − (4) j = 4ε ri j ri j as well as the Born–Mayer potential: Born–Mayer
Ui j 2
r . = A exp − ρ
(5)
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
We motivate the need for new, improved Xe–UO2 potentials by first describing some shortcomings of the existing Geng–Basak and Chartier–Morelon Xe–UO2 potentials, both in terms of their fitting procedures and their quantitative defect energetics. The fitting procedures for both the Geng and Chartier potentials have certain peculiarities. For the Geng potential, the Xe–O interaction was fit to a XeO3 crystal in a hypothetical CuAu3 structure, which is not likely to be representative of xenon in UO2 since oxygen will have significantly different chemical environments and xenon exists in different charge states in the two cases. Experimentally, XeO3 ¯ spacegroup), does not exist in the CuAu3 structure (Pm 3m but instead in the P21 21 21 space group [23]. The P21 21 21 structure of XeO3 is metastable and has a large, positive formation enthalpy [24]. The Chartier potential, although it was fit to DFT calculations, does not accurately reproduce the DFT energetics for one of the configurations to which it was fit (XeSD). For calculating xenon bubbles, the accuracy of defect energetics for defects such as the XeSD may be important. To correct the deficiencies of existing potentials, we turn to an alternative approach: We have developed a method, called iterative potential refinement, which uses a genetic algorithm to fit interatomic potentials of arbitrary form to energies, forces, and stresses from DFT calculations. We have used IPR, combined with ab initio molecular dynamics calculations of xenon in a Schottky defect in UO2 , to create a new interatomic potential for xenon in UO2 . We fit a new xenon potential for use with existing UO2 potentials from several geometric snapshots of the aforementioned molecular dynamics calculation. In order to determine whether this new potential is more accurate than the Geng or Chartier potentials, we perform a detailed comparison of the potentials across a range of defects. IPR is a very general approach to creating new potentials and we have performed a separate study fitting a UO2 potential for accurate phonons and defect energetics [25]. Because there is no systematic study of the accuracy of xenon interatomic potentials for UO2 , we compare defect energetics from three xenon interatomic potentials (Geng, Chartier, and our present IPR-generated) paired with two different UO2 potentials (Basak and Morelon) to DFT + U values for a variety of defects, representing different incorporation sites (large, intermediate, and small vacant sites). When adding xenon to these different sites, the Geng potential overestimates the defect energy for every size of incorporation site tested. The Chartier potential performs better than the Geng potential for small incorporation sites, but overestimates intermediate and large incorporation sites. The IPR-generated potential outperforms the Geng and Chartier potentials across all tested defects and has good agreement with DFT + U across incorporation sites of small to large sizes, even though the potential was only fit to a single defect configuration.
potentials as well as density functional theory (DFT). For DFT, we used spin-polarized calculations as implemented in the Vienna Ab Initio Simulation Package (VASP) [26–29]. We use the generalized gradient approximation of Perdew and Wang [30] with the Hubbard +U correction (DFT + U ) for strongly correlated systems as described by Dudarev [31]. Our energy cutoff is 500 eV with a 2 × 2 × 2 gamma centered kpoint grid. For the perfect UO2 lattice, we use the ‘U-ramping’ method described by Meredig et al [32] to ensure low energy orbital occupations were achieved in the DFT + U calculation. The defect calculations are performed without U -ramping at a value of U = 3.99 eV, starting at the bulk, distorted UO2 geometry. We have previously shown that this methodology produces accurate defect energetics, as described in [33]. We used 20 snapshots of molecular dynamics (MD) calculations as input to our genetic algorithm. Ab initio MD was performed with a 1 × 1 × 1 k-point grid at 1000 K, a Nos´e thermostat [33] in the NVT ensemble, and a 2.0 fs time step for 0.5 ps after equilibration. It should be noted that the DFT + U calculations contain a distortion in the oxygen sublattice that is not present in any of the interatomic potentials that we test. However, this distortion is relatively small. None of the potentials contain this distortion, so the small geometrical difference due to the distortion will be present for all of the potentials and will not significantly affect the results. For these calculations, we use a 96 atom cell (2 × 2 × 2 unit cells) supercell. 2.2. Interatomic potentials
The GULP code [34] is used for all defect and NEB calculations with interatomic potentials. We utilize two sets of UO2 /Xe potentials: Basak [18]/Geng [17] and Morelon [21]/ Chartier [20]. In some cases we perform calculations with different combinations of these potentials (e.g., Basak UO2 potential with Chartier xenon potential or Morelon UO2 potential with Geng xenon potential). In order to achieve a direct comparison to our DFT + U calculations, we perform our interatomic calculations using the same supercell size (96 atoms) as DFT. 3. Iterative potential refinement, a novel genetic algorithm/density functional theory approach for generating interatomic potentials
We describe here our approach for automated fitting of interatomic potentials to density functional theory calculations. Our iterative potential refinement (IPR) method is a genetic algorithm (GA)-based approach for fitting an interatomic potential to the results of first-principles calculations. The concept of fitting an interatomic potential to first-principles data has been well-explored, and dates at least to the force-matching method of Ercolessi et al [35]. In our case, we use the first-principles energies, forces, and stresses as targets, and attempt to create an interatomic potential that best reproduces those input data. In other words, we employ a GA (i.e., a global optimization technique) to minimize a cost function consisting of error relative to first-principles calculations, similar to the method proposed by Brommer and G¨ahler [36]. On startup,
2. Methodology 2.1. DFT
We calculate defects in uranium dioxide, including xenon in interstitials and xenon in vacancy clusters, using interatomic 3
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
IPR creates a population of random-guess potential forms, and applies GA operations of selection, crossover, mutation, and replacement to iteratively create improved interatomic potentials. It terminates either when the error has reached an acceptable threshold, or when a maximum number of iterations has been performed. IPR is, at its core, a relatively standard real-numberencoded GA [37]. In such a scheme, instead of the traditional GA representation of solutions as binary strings, each solution is simply a real-number array of the potential parameters being optimized. As a simple example, if we wished to fit a Lennard-Jones potential of the form of equation (4), where ε and σ are adjustable parameters, the IPR population would consist initially of n random combinations Uk = [ε, σ ], where k is the index of the potential, within some specified search space. Each candidate potential Vi is then used to calculate a series of input structures whose first-principles energies, forces, and stresses are known. The resulting interatomic potential outputs are compared to the reference first-principles values to calculate an error, and this error is the fitness function that IPR minimizes. As noted above, the GA relies on evolutionary operators such as selection, crossover, mutation, and replacement to produce better interatomic potentials over time. Selection describes the process by which we choose ‘parent’ potentials from among our population of n that will be used to generate new ‘child’ potentials during crossover. IPR performs tournament selection, in which two potentials are drawn at random from the population, and the potential with the better fitness value of the two becomes one crossover parent. IPR’s crossover operator mixes two parent potentials Uk and Ul to create a child potential Uchild that has characteristics of each parent: Uchild = w · Uk + (1 − w) · Ul ,
4.1. Analysis of previous potentials
Geng potential: fitting to-X eO 3 . First, we examine the Geng potential to find how accurately it represents xenon in UO2 . In order to better understand these compounds we used DFT to find the relative stabilities of these two structures. We performed four calculations: an O2 molecule (with the O2 molecular energy correction of [38]), gaseous xenon, XeO3 in the CuAu3 structure, and XeO3 in the P21 21 21 space group. The calculated formation energy of the experimentally observed P21 21 21 structure is +3.6 eV/formula unit. The formation energy of XeO3 in the hypothetical CuAu3 structure is even more significantly positive, at +6.6 eV/formula unit. The calculated energy of the P21 21 21 structure is close to the experimentally measured formation energy of +4.2 eV/formula unit [24]. These positive formation energies for both phases means that these XeO3 crystals are strongly thermodynamically unstable and the most stable states for xenon and oxygen are their respective elemental phases. Intuitively, neither structure seems appropriate for fitting a Xe–O interaction for xenon defects in UO2 . In addition to these structures being thermodynamically unstable, any compound of XeO3 will not be appropriate for fitting interactions of xenon in UO2 . While XeO3 -type interactions might be appropriate for modeling xenon and oxygen in a compound of only these two elements, Geng et al use them to model the interaction of xenon with oxygen ions in UO2 . The net result of the fitting procedure for the Geng potential is that the potential describes a very short equilibrium separation between xenon and oxygen (2.8 Å). This short interaction can be traced back to the distance between xenon and oxygen in the CuAu3 structure (2.6 Å). We can use this separation distance along with the xenon–xenon separation (2.4 Å) from the Geng potential to approximate the atomic radii of xenon and oxygen: the Geng potential describes an effective O2− radius of just 0.4 Å, whereas the empirical value for an oxygen ion is 1.24 Å [39]. Thus, we might expect that any situation in which xenon and oxygen come into close proximity with one another will not be well described by the Geng potential. Geng potential: small incorporation sites. We assess the accuracy of the Geng potential by calculating the energetics of inserting a single xenon atom into the UO2 lattice. One possible way to insert xenon into the lattice of UO2 without removing any uranium or oxygen atoms is to insert xenon into the octahedral interstitial position (Xei ) where xenon would strongly repel the neighboring oxygen atoms, causing a strain penalty. Another choice would be to substitute xenon in for a uranium atom and place that uranium atom in a neighboring interstitial site (XeU + Ui ). In addition to xenon repelling the neighboring oxygen atoms and causing a strain penalty, like the Xei case, XeU + Ui should also incur an electrostatic penalty with uranium moved from its ideal fluorite position, (the XeU will have an effective charge of −4e while the Ui would have an effective charge of 4e). Therefore, we would expect that the XeU + Ui configuration should have a larger defect energy than the Xei configuration. We have calculated the incorporation energies of these defects from
where w is a random number between 0 and 1. Mutation simply replaces child parameter values with new, entirely random values with some low probability pmutation (in this work, 0.1) during crossover. Selection is the process by which we replace poorer potentials with hopefully-improved children; here, one GA generation (i.e., iteration) consists of generating and evaluating child potentials to replace the bottom 75% of the population. Performing many such generations gives improved interatomic potentials as the GA ‘learns’ the solution space and locates optima. Here, we use IPR, along with DFT + U molecular dynamics simulations, to parameterize a new xenon potential using the Basak as the base UO2 potential. 4. Results
We begin by benchmarking the Basak/Geng and Morelon/ Chartier potentials against DFT calculations of xenoncontaining defects in UO2 . We consider all combinations of the xenon and UO2 potentials and we refer to each possible potential combination by, first, the base UO2 potential, and second, the add-on xenon potential: Basak/Geng, Morelon/Chartier, Basak/Chartier, and Morelon/Geng. Firstly, we examine the fitting and form of the Geng and Chartier xenon potentials and their ability to reproduce defect energetics from DFT + U calculations. 4
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
distance (0.006 eV), despite xenon being closer to oxygen than to uranium for both defect configurations. This unphysically small Xe–O interaction in UO2 is a direct consequence of the choice to fit to XeO3 , which contains short Xe–O bonds, as discussed above. These small Xe–O interactions mean that even though xenon neighbors eight oxygen atoms in both defects we are considering, they will not contribute significantly to the defect energy for either case. Next, we consider the xenon–uranium interactions. The xenon–uranium distances differ for the two defect sites. The Xei –U distance is 2.74 Å, while the XeU –U distance is 3.87 Å. For the Geng potential, the Xei –U and XeU –U the interaction energies at those distances are 6.63 eV and 0.44 eV, respectively, showing that the Xe–U interaction is fairly long ranged. Although the oxygen atoms are closer to xenon than the uranium atoms, the repulsive interaction between uranium and xenon is much stronger than the interaction between xenon and oxygen, showing that the Xe–U interaction is too strong compared to the Xe–O interaction. The large Xe–U interaction for Xei translates into a larger defect energy than for XeU + Ui , despite the latter having an unfavorable electrostatic penalty. Because the Xe–O interactions are too small and the Xe–U interactions are too large, the Geng potential gives an energy of XeU + Ui that is lower than the energy of Xei . The fitting approach for the Geng potential has created a qualitative failure for the ordering of these defect energies. As discussed above, the Chartier potential accurately describes small incorporation sites such as xenon in an interstitial, but it overestimates the energy of the large incorporation site of the XeSDs [21]. Because the Chartier potential does not accurately reproduce the XeSD incorporation energies, it should be tested for other xenon in defects in large incorporation sites.
Figure 1. The interatomic potential functions of energy versus distance from the Geng, Chartier, and IPR-generated potentials. The interatomic potentials are generally characterized by repulsive interactions only. The Geng potential has the smallest separation for Xe–O and the largest separation for Xe–U. The Chartier Xe–O potential has a low curvature compared to the other xenon interactions, which are treated more like hard spheres. The vertical lines represent the bond distances for atoms that are nearest neighbors to xenon in Xei and XeU + Ui .
5. New Xe–O and Xe–U interatomic potential fitted by IPR
As we have discussed above, the existing Geng and Chartier potentials have shortcomings in describing defects, so we have created a new set of parameters for the Xe–O and Xe–U interactions (leaving the U–U, O–O, and U–O interactions equal to the values given by Basak et al) using iterative potential refinement. IPR works by fitting to DFT + U energies, forces, and stresses with a genetic algorithm. For the IPR-generated Xe–O and Xe–U potentials, we consider only the repulsive interactions of a Born–Mayer potential. For the IPR fitting procedure. we used snapshots from a single DFT + U molecular dynamics calculation of xenon in a Schottky defect cluster with the oxygen vacancies at nearest neighbor positions with respect to one another (XeSD1). The IPR-generated potential was not fit to any other defects. We also examined the effect of using the Morelon potential as the base UO2 potential, and the results were very similar to those from the Basak potential. Table 2 shows the result of the IPR fitting procedure. We can examine the form of the IPR potential compared to the Geng and Chartier potentials. Figure 1 shows the form of the Basak, Chartier, and our IPR-generated potentials as a function of atomic separation. The Geng Xe–O potential
DFT + U and the Basak/Geng potential. From the DFT + U calculations, we find that the energetics follow our intuition with the incorporation energy of XeU + Ui having a larger defect energy than Xei (11.3 eV versus 9.7 eV). However, the Basak/Geng potential finds there is a much larger penalty to create these defects compared to DFT + U and that XeU + Ui is lower in energy than Xei (17.0 eV versus 22.7 eV). For comparison, the Morelon/Chartier potential finds that the two defect configurations are nearly degenerate in energy (11.8 eV for both defects). To find why the interatomic potential energies are not in agreement with DFT + U for these defects, we can examine the energetic components of the potential (Xe–O and Xe–U pairs). For this test, uranium and oxygen ions are placed at their equilibrium, fluorite positions. For both defect cases, the distance between xenon and its eight oxygen neighbors is the same (2.36 Å, for a0 = 5.47 Å). Counter to our intuition, the Geng Xe–O potential has almost no interaction at this 5
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Finally, we calculate the defect energetics of a single large incorporation site, xenon in a double Schottky defect cluster (Xe2 · SD). It is important to note that our IPR potential was fit only to one of the test configurations (XeSD1). For these defects, we calculate the incorporation energy (the energy required to add xenon to an existing defect site). These defect energetics are defined as follows: Xenon incorporation energy (1E Inc. Xe ). This quantity is the energy of incorporating a xenon atom into the octahedral interstitial site of UO2 .
Table 2. Potential parameters for IPR-generated Xe–O and Xe–U
potentials. The form of the potential is of equation (5). The U–O, O–O, and U–U parameters are those of Basak et al [18].
Xe–U Xe–O
A (eV)
ρ (Å)
6607 300 95 913
0.178 0.210
differs significantly from the Chartier and IPR potentials in terms of the effective size for xenon, a result of fitting to XeO3 . The Chartier Xe–O potential differs from the other xenon interactions in the ‘hardness’ (the curvature of the potential near its minimum) of the Xe–O interactions. The difference in hardness may be a source of error in the Chartier potential. In contrast, using the IPR algorithm, we find that the Xe–O and Xe–U IPR potentials have similar, hard curvatures. The Geng and Chartier potentials describe very different Xe–U separations, and the Xe–U separation described by our IPR potential falls between them. Next, we create a quantitative comparison of all possible combinations of UO2 potentials and xenon potentials for a set of defects in a variety of incorporation sites to determine how well each combination of the potentials are able to reproduce DFT + U energetics.
N UO2 1E Inc. − E Xe , Xe = E Xei − N E
(6)
N is the energy of xenon in the interstitial position where E Xe i of a supercell with N UO2 formula units, E N is the energy per formula unit of bulk UO2 , and E Xe is the energy of a single, isolated xenon atom (for the interatomic potentials considered, this quantity is zero). Similarly, we can calculate the incorporation energy of XeU + Ui . Inc. ): XeU + Ui incorporation energy (1E Xe N UO2 1E Inc. − E Xe . XeU +Ui = E XeU ,Ui − N E
(7)
Inc. Schottky defect-xenon atom incorporation energy (E Xe−SD ). This quantity represents the energy to add a xenon atom into an existing Schottky defect cluster.
5.1. Benchmarking interatomic potentials with DFT + U
Now we benchmark the Geng, Chartier, and our IPR-generated interatomic potentials against DFT + U for various xenoncontaining defects. We can determine how accurate the interatomic potentials are by comparison to the energetics of relatively small defects that are practical to calculate using DFT + U . We have created a series of physically relevant benchmarks to compare all three xenon potentials on equal footing. We place xenon in different incorporation sites in UO2 (small, intermediate, and large vacant sites). For small incorporation site, we used the following defects: xenon in an interstitial (Xei ), xenon in an interstitial without relaxing the ideal fluorite UO2 geometry, and xenon on a uranium site with uranium in a neighboring interstitial (XeU + Ui ). The intermediate-sized incorporation sites are xenon in the three configurations of Schottky defect cluster (XeSD1, XeSD2, and XeSD3, in order of increasing oxygen vacancy separation).
N −1 N −1 1E Inc. − E Xe , Xe−SD = E Xe − E
(8)
N −1 where E Xe is the is the energy of a supercell with N sites for UO2 formula units that is missing one UO2 formula unit (the Schottky defect cluster) plus a xenon atom, and E N −1 represents the energy of a supercell with an empty SD. We can also define the incorporation energy for xenon in a double Schottky defect N −2 − E N −2 − E Xe . 1E Xe−2:SDInc. = E Xe
(9)
We then calculated these incorporation energies with DFT + U as well as combinations of the UO2 potentials (Basak and Morelon) with xenon potentials (Geng, Chartier, and IPR). The energies of these xenon-containing defects calculated are shown in table 3. The IPR potentials give the best
Table 3. Comparison of xenon point defects of a 96 atom cell of UO2 calculated with DFT + U and various interatomic potentials. The
Basak and Morelon potentials are used for the U–O, U–U, and O–O interactions, while the Geng, Chartier, and IPR potentials are used for Xe–O and Xe–U interactions. The IPR potential is fitted to a DFT + U molecular dynamics calculation of xenon in SD1. The IPR potentials give the best agreement with DFT calculations, even in configurations which it was not directly fit. Bold energies correspond to the closest to the DFT + U energies. Incorporation energies (eV)
DFT + U
Basak/Geng
Basak/Chartier
Morelon/Geng
Morelon/Chartier
Basak/IPR
Morelon/IPR
Xei (relaxed) Xei (static) XeU + Ui Xe in SD1 Xe in SD2 Xe in SD3 Xe in two SDs
9.73 15.64 11.33 1.06 1.83 1.94 0.27
22.74 41.45 17.03 4.49 4.57 4.60 5.08
12.85 14.61 16.95 4.80 5.26 5.72 5.03
21.72 41.46 17.45 4.57 4.63 4.67 2.68
11.76 14.61 11.76 4.10 4.43 4.61 2.65
9.67 18.70 11.90 0.60 0.87 1.17 0.32
7.96 18.70 9.87 0.39 0.56 0.62 0.14
6
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Figure 2. The percentage difference between the formation energy of DFT + U and interatomic potentials for a series of defects. The xenon is placed in small incorporation sites which have large defect formation energies (Xei , Xei without geometry optimization, and XeU + Ui ), intermediate-sized incorporation sites (xenon in three Schottky defect configurations: XeSD1, XeSD2, and XeSD3, in order of increasing VO –VO separation), and large incorporation sites (xenon in a double Schottky defect). The IPR potential performs better than either Chartier or Morelon, independent of the UO2 potential used. As size of the incorporation site increases, Chartier and Morelon perform much worse.
Figure 3. The relaxation energy error for xenon defects in different
incorporation sites as described by equation (6). The defects were calculated at their DFT + U geometries without geometry optimization. The difference in defect energy between interatomic potential calculations at the DFT + U geometry and DFT + U calculations is a measure of how well the interatomic potential calculations match the DFT + U geometry. Both the Basak/IPR and Morelon/IPR combinations have much smaller differences from DFT + U than the Geng or Chartier containing calculations, even for defects in which IPR was not fit.
agreement with DFT + U calculations, even in configurations to which it was not directly fit. As we have previously shown, there is qualitative disagreement between DFT + U calculations and the Basak/Geng and Morelon/Chartier interatomic potentials for the energetic ordering of XeU + Ui and Xei . The IPR potentials (as well as the Basak/Chartier potential combination), agree with the DFT + U ordering for those defects. For comparison, table 4 shows the defect energetics
calculated with the interatomic potentials obtained using a larger supercell (2592 atoms) for better convergence. Figure 2 shows the percentage error of the interatomic potential defect energies compared to DFT + U calculations. For calculations with the IPR potential, the choice of UO2 potential does not have a large impact on the defect energetics; the IPR potential gives better agreement with DFT + U when 7
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Table 4. Comparison of xenon point defects of a 2592 atom cell of UO2 calculated with various interatomic potentials. Because of the larger supercell size, the image effects are negligible and these calculations are well converged.
Incorporation energies (eV)
Basak/Geng
Basak/Chartier
Morelon/Geng
Morelon/Chartier
Basak/IPR
Morelon/IPR
Xei XeU + Ui Xe in SD1 Xe in SD2 Xe in SD3 Xe in two SDs
22.61 16.93 4.75 4.84 4.94 4.08
12.83 17.13 5.14 6.05 6.55 3.16
22.20 17.62 4.72 4.76 4.81 4.07
11.84 11.84 4.19 4.60 4.78 2.69
9.62 11.86 0.90 1.51 1.79 0.18
8.15 10.29 0.44 0.66 0.70 0.13
paired with the Basak UO2 potential (which was used for the UO2 interactions during the IPR fitting), but also performs well with the Morelon UO2 potential. The IPR potential (whether paired with Basak or Morelon) is comparable with or outperforms the Geng and Chartier potentials across all incorporation sites. The Basak/Geng, Basak/Chartier, and Morelon/Geng potentials have high error for the Xei and static Xei and low error for the XeU + Ui ; however, with the error near 50%, the XeU + Ui energy difference represents 5.7 eV deviation from DFT + U . The Morelon/Chartier potential has good agreement with DFT + U for the three small incorporation sites. All of the Geng and Chartier derived calculations have a large error for the intermediate-sized incorporation sites of the XeSDs (>100%) and even larger errors for double Schottky defects, which are essentially very small bubbles (∼900% for Chartier and ∼1700% for Geng). The IPR potential has the most accurate energetics for the double Schottky defect by a large margin, even though it was not fit to this geometry, which suggests that the IPR potential is much more transferable to a wide range of geometries than the other potentials. When considering all of the defects calculated, IPR potentials find the best agreement with DFT + U . The deviation of defect energies from DFT values is just one metric we can use to validate the interatomic potentials. We can also compare the geometries of the DFT + U and interatomic potential calculations, albeit indirectly. We calculate the defect energetics from the potentials at the DFT + U geometries and compare to the DFT + U defect energies. We normalize the difference in defect energies to the DFT + U Relax defect energy to find a relaxation energy error, E DFT Geometry , which is defined as Relax E DFT Geometry =
Pot DFT+U ) (E DFT Geo − E , E DFT+U
and intermediate-sized incorporation sites calculated with the IPR potential have markedly smaller relaxation energy error compared to Geng and Chartier. For the double Schottky defect, the differences among the three xenon potentials are smaller, but the IPR-generated potential still performs the best. Overall, the IPR potential (paired with Basak or Morelon) performs best not only on energetics, but also on geometries. 6. Conclusions
Using iterative potential refinement, or IPR, we have created a new xenon interatomic potential for use with the Morelon and Basak UO2 potentials. This IPR-generated potential was fit by a genetic algorithm to snapshots of DFT + U molecular dynamics of xenon in a Schottky defect cluster. We have compared defect energetics of xenon in a variety of incorporation sites (small, intermediate, and large vacant sites) using combinations of xenon potentials (Geng, Chartier, and IPR-generated) and UO2 potentials (Basak and Morelon) to DFT + U calculated energetics. The IPR-generated potential, despite only being fit to a single defect configuration in a intermediate-sized incorporation site, gives accurate energetics for small and large incorporation sites as well, pointing to its predictive power. The existing xenon potentials overestimate the energetics of xenon incorporation across all sizes of incorporation sites tested. The existing potentials do not yield accurate defect energetics due, in part, to their fitting methodologies. The Geng potential describes the interaction between oxygen and xenon as too short, and the Chartier potential finds that the xenon–oxygen interaction is not ‘hard’ enough. The IPR-generated potential performs better than both the Geng and Chartier potentials across a range of defects, including geometries to which it was not fit such as xenon in an interstitial and xenon in a double Schottky defect, while adding no additional fitting parameters compared to the Chartier or Geng potentials. The IPR-generated potential will facilitate calculations of xenon bubbles in a variety of morphologies and sizes.
(10)
Pot where E DFT Geo is the defect energy as calculated with an interatomic potential as described in equations (2)–(5), but calculated at the DFT geometry and E DFT+U is defect energy from DFT + U . Figure 3 shows relaxation energy error described by equation (10) for xenon defects in different incorporation sites. As with the defect energy comparison above, the IPR potentials perform the best, even for defect configurations for which IPR was not fit. However, in contrast to the defect energy comparison above, the relaxation energy error does not depend strongly on incorporation sites. The small
References [1] Zhou F and Ozolins V 2009 Phys. Rev. B 80 125127 [2] Tiwary P, van de Walle A and Grønbech-Jensen N 2009 Phys. Rev. B 80 174302 [3] Tiwary P, van de Walle A, Jeon B and Grønbech-Jensen N 2011 Phys. Rev. B 83 094104 8
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
[4] Hanken B E, Stanek C R, Grønbech-Jensen N and Asta M 2011 Phys. Rev. B 84 85131 [5] Freyss M, Vergnet N and Petit T 2006 J. Nucl. Mater. 352 144 [6] Petit T 2002 Study of the atomic rare gas behavior by Ab Initio calculations Fission Gas Behaviour in Water Reactor Fuels (Cadarache: OECD Publishing) [7] Yun Y, Kim H, Kim H and Park K 2008 J. Nucl. Mater. 378 40 [8] Yun Y, Eriksson O and Oppeneer P M 2009 J. Nucl. Mater. 385 510 [9] Yun Y, Eriksson O, Oppeneer P M, Kim H and Park K 2009 J. Nucl. Mater. 385 364 [10] Jeon B, Asta M, Valone S M and Grønbech-Jensen N 2010 Nucl. Instrum. Methods Phys. Res. B 268 2688 [11] Grimes R W and Catlow C R A 1991 Phil. Trans. R. Soc. A 335 609 [12] Govers K, Lemehov S, Hou M and Verwerft M 2007 J. Nucl. Mater. 366 161 [13] Govers K, Lemehov S, Hou M and Verwerft M 2008 J. Nucl. Mater. 376 66 [14] Govers K, Lemehov S E and Verwerft M 2010 J. Nucl. Mater. 405 252 [15] Arima T, Yamasaki S, Inagaki Y and Idemitsu K 2005 J. Alloys Compounds 400 43 [16] Ball R G J and Grimes R W 1992 J. Nucl. Mater. 188 216 [17] Geng H Y, Chen Y, Kaneta Y and Kinoshita M 2008 J. Alloys Compounds 457 465 [18] Basak C B, Sengupta A K and Kamath H S 2003 J. Alloys Compounds 360 210 [19] Jackson R A and Catlow C R A 1985 J. Nucl. Mater. 127 161 [20] Chartier A, Van Brutzel L and Freyss M 2010 Phys. Rev. B 81 174111
[21] Morelon N-D, Ghaleb D, Delaye J-M and Van Brutzel L 2003 Phil. Mag. 83 1533 [22] Thompson A and Wolverton C 2013 Phys. Rev. B 87 104105 [23] Templeton D H, Zalkin A, Forrester J D and Williamson S M 1963 J. Am. Chem. Soc. 85 817 [24] Bard A J, Parsons R and Jordan J 1985 Chemistry IU of P and A. Standard Potentials in Aqueous Solution (Boca Raton, FL: CRC Press) [25] Thompson A E, Meredig B, Stan M and Wolverton C 2014 J. Nucl. Mater. 446 155–62 [26] Kresse G and Hafner J 1993 Phys. Rev. B 47 558 [27] Kresse G 1993 Thesis, Technische Universit¨at Wien, Vienna, Austria [28] Kresse G and Furthm¨uller J 1996 Comput. Mater. Sci. 6 15 [29] Kresse G and Furthm¨uller J 1996 Phys. Rev. B 54 11169 [30] Perdew J P and Ziesche P (ed) 1991 Electron. Struct. Solids vol. 11 (Berlin: Akademie Verlag) [31] Dudarev S L, Botton G A, Savrasov S Y, Humphreys C J and Sutton A P 1998 Phys. Rev. B 57 1505 [32] Meredig B, Thompson A, Hansen H A, Wolverton C and van de Walle A 2010 Phys. Rev. B 82 195128 [33] Thompson A and Wolverton C 2011 Phys. Rev. B 84 134111 [34] Gale J D and Rohl A L 2003 Mol. Simul. 29 291 [35] Ercolessi F and Adams J B 1994 Europhys. Lett. (EPL) 26 583 [36] Brommer P and G¨ahler F 2007 Modelling Simul. Mater. Sci. Eng. 15 295 [37] Michalewicz Z, Janikow C Z and Krawczyk J B 1992 Comput. & Math. Appl. 23 83 [38] Wang L, Maxisch T and Ceder G 2006 Phys. Rev. B 73 195107 [39] Shannon R D 1976 Acta. Cryst. A 32 751
9
Search
Collections
Journals
About
Contact us
My IOPscience
An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 J. Phys.: Condens. Matter 26 105501 (http://iopscience.iop.org/0953-8984/26/10/105501) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 155.97.178.73 This content was downloaded on 01/07/2014 at 15:50
Please note that terms and conditions apply.
:
Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 229501 (1pp)
doi:10.1088/0953-8984/26/22/229501
Corrigendum: An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach (2014 J. Phys.: Condens. Matter 26 105501) A E Thompson, B Meredig and C Wolverton Department of Materials Science and Engineering, Northwestern University, Evanston, IL 60208, USA Received 28 April 2014 Accepted for publication 28 April 2014 Published 14 May 2014
We would like to add the following acknowledgments to our paper: The authors would like to acknowledge funding from the Department of Energy under NERI-C grant no. DE-FG07-07ID1489, and under grant no. DE-F602-07ER46433. BM also acknowledges the National Defense Science and Engineering Graduate Fellowship.
0953-8984/14/229501+1$33.00
1
© 2014 IOP Publishing Ltd Printed in the UK
Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 105501 (9pp)
doi:10.1088/0953-8984/26/10/105501
An improved interatomic potential for xenon in UO2: a combined density functional theory/genetic algorithm approach Alexander E Thompson, Bryce Meredig and C Wolverton Department of Materials Science and Engineering, Northwestern University, Evanston, IL 60208, USA E-mail: [email protected] Received 30 August 2013, revised 4 November 2013 Accepted for publication 3 January 2014 Published 19 February 2014
Abstract
We have created an improved xenon interatomic potential for use with existing UO2 potentials. This potential was fit to density functional theory calculations with the Hubbard U correction (DFT + U ) using a genetic algorithm approach called iterative potential refinement (IPR). We examine the defect energetics of the IPR-fitted xenon interatomic potential as well as other, previously published xenon potentials. We compare these potentials to DFT + U derived energetics for a series of xenon defects in a variety of incorporation sites (large, intermediate, and small vacant sites). We find the existing xenon potentials overestimate the energy needed to add a xenon atom to a wide set of defect sites representing a range of incorporation sites, including failing to correctly rank the energetics of the small incorporation site defects (xenon in an interstitial and xenon in a uranium site neighboring uranium in an interstitial). These failures are due to problematic descriptions of Xe–O and/or Xe–U interactions of the previous xenon potentials. These failures are corrected by our newly created xenon potential: our IPR-generated potential gives good agreement with DFT + U calculations to which it was not fitted, such as xenon in an interstitial (small incorporation site) and xenon in a double Schottky defect cluster (large incorporation site). Finally, we note that IPR is very flexible and can be applied to a wide variety of potential forms and materials systems, including metals and EAM potentials. Keywords: density functional theory, potential fitting, molecular dynamics, UO2 (Some figures may appear in colour only in the online journal)
1. Introduction
dynamics simulations [2, 3, 10–16]. There have been many interatomic potentials parameterized for UO2 [12, 13], but only a few for xenon in UO2 , despite the crucial role of xenon in the failure of nuclear fuels. In this paper we consider two recently-developed empirical Xe–UO2 models as well as our own potential, generated using a new genetic algorithm-based approach which we call iterative potential refinement (IPR). The first xenon potential we consider was proposed by Geng et al [17] who developed Xe–O and Xe–U interactions to use with the UO2 potential of Basak et al [18] in order to examine the properties of xenon bubbles in UO2 . For simplicity, we will refer to the potentials
There is much interest in the behavior of xenon in uranium dioxide (UO2 ) due to the fact that xenon atoms, which are created as a consequence of uranium fission, cluster into bubbles that can swell the host UO2 material, degrading its properties and causing it to fail. There have been many studies of xenon in UO2 using density functional theory (DFT) [1–9], which can achieve good defect energetics, but DFT is too computationally demanding for large systems such as xenon bubbles in UO2 . Interatomic potentials, on the other hand, enable calculations of larger systems and long molecular 0953-8984/14/105501+09$33.00
1
c 2014 IOP Publishing Ltd
Printed in the UK
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Table 1. Potential parameters of uranium dioxide and xenon potentials. The xenon is neutrally charged in the potentials.
Uranium dioxide–xenon potentials Basak A (eV) U–U O–O U–O U O
294.64 1633.0051 693.6487 Charge (e) +2.4 −1.2
Buckingham ρ (Å) c (eV Å6 ) 0.32702 0.32702 0.32702
Geng
Born–Mayer A (eV) ρ (Å)
Xe–U Xe–O
4887.7
Morelon
Buckingham A (eV) ρ (Å)
U O
11272.6 566.498 Charge (e) +3.227252 −1.613626
0 3.94879 0
De (eV)
Morse a (Å−1 )
r0 (Å)
0.57719
1.65
2.369
Lennard-Jones E (eV) σ (Å)
0.415
0.1363 0.42056
Chartier
Born–Mayer A (eV) ρ (Å)
Xe–U Xe–O
478875.11 94.489
0.0099
2.5
c (eV Å6 )
Buckingham four-range rmin1 (Å) cut1 (Å) rmin2 (Å)
cut2 (Å)
0
2.6
134 0
1.2
2.1
0.179 0.58
Oxides are often described well with Buckingham potentials: ri j c UiBuck = A exp − (1) − 6 j ρ ri j
by the name of the first author. The Xe–O potential was fit to XeO3 in a hypothetical (i.e., non-ground-state) CuAu3 structure, and the Xe–U was fit to the same Born–Meyer form as Jackson et al [19]. The parameters of Jackson et al were used as a starting point and were refined with ab initio calculations. However, instead of fitting to xenon in a complete UO2 lattice, the interaction was fit to a ‘virtual structure’. This virtual structure is constructed from a 12 atom unit cell of UO2 in the fluorite structure, but with three of the UO2 formula units removed. The remaining UO2 formula unit and xenon atom are aligned in the [111] direction. The Basak potential was parameterized against thermal expansion of UO2 and achieves good high temperature properties such as thermal expansion and isothermal compressibility [18]. The second xenon potential was proposed by Chartier et al [20] to use in conjunction with the UO2 potential of Morelon et al [21]. The Chartier potential was parameterized by fitting to DFT calculations of four xenon-containing defects (xenon in the following sites: interstitial, uranium vacancy, oxygen vacancy, and a Schottky defect cluster). The Morelon potential was fit to experimental defect and migration energies of UO2 [21]. The Chartier potential treats high-energy defects such as a xenon atom in an interstitial (Xei ) accurately, but overestimates the defect energy of xenon in small vacancy clusters, like xenon in a Schottky defect (XeSD) [21]. We have previously used these potentials to study different diffusion pathways of xenon in tetravacancies in UO2 [22]. Table 1 shows the parameterizations of these various UO2 –Xe potentials. These potentials take on various forms.
where r is the separation between atom I and atom j. A variation on the Buckingham is sometimes used to avoid an unphysical attraction at short separations between atoms by introducing different functional forms at discrete separations in the Buckingham four-range potential: ri j A exp − if rmin1 < ri j ≤ cut1 ρ fifth order polyomial if cut < r ≤ r 1 ij min2 UiBuck4 = (2) j third order polyomial if r < r ≤ cut min2 i j 2 C if cut2 < ri j ≤ rmax . − r 6 ij The Morse potential incorporates some covalent behavior in the interactions between atoms: UiMorse = De [1 − exp(−a(ri j − r0 ))]. j
(3)
The xenon potentials are fairly simple and include the Lennard-Jones potential: " 6 # 12 σ σ UiLJ − (4) j = 4ε ri j ri j as well as the Born–Mayer potential: Born–Mayer
Ui j 2
r . = A exp − ρ
(5)
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
We motivate the need for new, improved Xe–UO2 potentials by first describing some shortcomings of the existing Geng–Basak and Chartier–Morelon Xe–UO2 potentials, both in terms of their fitting procedures and their quantitative defect energetics. The fitting procedures for both the Geng and Chartier potentials have certain peculiarities. For the Geng potential, the Xe–O interaction was fit to a XeO3 crystal in a hypothetical CuAu3 structure, which is not likely to be representative of xenon in UO2 since oxygen will have significantly different chemical environments and xenon exists in different charge states in the two cases. Experimentally, XeO3 ¯ spacegroup), does not exist in the CuAu3 structure (Pm 3m but instead in the P21 21 21 space group [23]. The P21 21 21 structure of XeO3 is metastable and has a large, positive formation enthalpy [24]. The Chartier potential, although it was fit to DFT calculations, does not accurately reproduce the DFT energetics for one of the configurations to which it was fit (XeSD). For calculating xenon bubbles, the accuracy of defect energetics for defects such as the XeSD may be important. To correct the deficiencies of existing potentials, we turn to an alternative approach: We have developed a method, called iterative potential refinement, which uses a genetic algorithm to fit interatomic potentials of arbitrary form to energies, forces, and stresses from DFT calculations. We have used IPR, combined with ab initio molecular dynamics calculations of xenon in a Schottky defect in UO2 , to create a new interatomic potential for xenon in UO2 . We fit a new xenon potential for use with existing UO2 potentials from several geometric snapshots of the aforementioned molecular dynamics calculation. In order to determine whether this new potential is more accurate than the Geng or Chartier potentials, we perform a detailed comparison of the potentials across a range of defects. IPR is a very general approach to creating new potentials and we have performed a separate study fitting a UO2 potential for accurate phonons and defect energetics [25]. Because there is no systematic study of the accuracy of xenon interatomic potentials for UO2 , we compare defect energetics from three xenon interatomic potentials (Geng, Chartier, and our present IPR-generated) paired with two different UO2 potentials (Basak and Morelon) to DFT + U values for a variety of defects, representing different incorporation sites (large, intermediate, and small vacant sites). When adding xenon to these different sites, the Geng potential overestimates the defect energy for every size of incorporation site tested. The Chartier potential performs better than the Geng potential for small incorporation sites, but overestimates intermediate and large incorporation sites. The IPR-generated potential outperforms the Geng and Chartier potentials across all tested defects and has good agreement with DFT + U across incorporation sites of small to large sizes, even though the potential was only fit to a single defect configuration.
potentials as well as density functional theory (DFT). For DFT, we used spin-polarized calculations as implemented in the Vienna Ab Initio Simulation Package (VASP) [26–29]. We use the generalized gradient approximation of Perdew and Wang [30] with the Hubbard +U correction (DFT + U ) for strongly correlated systems as described by Dudarev [31]. Our energy cutoff is 500 eV with a 2 × 2 × 2 gamma centered kpoint grid. For the perfect UO2 lattice, we use the ‘U-ramping’ method described by Meredig et al [32] to ensure low energy orbital occupations were achieved in the DFT + U calculation. The defect calculations are performed without U -ramping at a value of U = 3.99 eV, starting at the bulk, distorted UO2 geometry. We have previously shown that this methodology produces accurate defect energetics, as described in [33]. We used 20 snapshots of molecular dynamics (MD) calculations as input to our genetic algorithm. Ab initio MD was performed with a 1 × 1 × 1 k-point grid at 1000 K, a Nos´e thermostat [33] in the NVT ensemble, and a 2.0 fs time step for 0.5 ps after equilibration. It should be noted that the DFT + U calculations contain a distortion in the oxygen sublattice that is not present in any of the interatomic potentials that we test. However, this distortion is relatively small. None of the potentials contain this distortion, so the small geometrical difference due to the distortion will be present for all of the potentials and will not significantly affect the results. For these calculations, we use a 96 atom cell (2 × 2 × 2 unit cells) supercell. 2.2. Interatomic potentials
The GULP code [34] is used for all defect and NEB calculations with interatomic potentials. We utilize two sets of UO2 /Xe potentials: Basak [18]/Geng [17] and Morelon [21]/ Chartier [20]. In some cases we perform calculations with different combinations of these potentials (e.g., Basak UO2 potential with Chartier xenon potential or Morelon UO2 potential with Geng xenon potential). In order to achieve a direct comparison to our DFT + U calculations, we perform our interatomic calculations using the same supercell size (96 atoms) as DFT. 3. Iterative potential refinement, a novel genetic algorithm/density functional theory approach for generating interatomic potentials
We describe here our approach for automated fitting of interatomic potentials to density functional theory calculations. Our iterative potential refinement (IPR) method is a genetic algorithm (GA)-based approach for fitting an interatomic potential to the results of first-principles calculations. The concept of fitting an interatomic potential to first-principles data has been well-explored, and dates at least to the force-matching method of Ercolessi et al [35]. In our case, we use the first-principles energies, forces, and stresses as targets, and attempt to create an interatomic potential that best reproduces those input data. In other words, we employ a GA (i.e., a global optimization technique) to minimize a cost function consisting of error relative to first-principles calculations, similar to the method proposed by Brommer and G¨ahler [36]. On startup,
2. Methodology 2.1. DFT
We calculate defects in uranium dioxide, including xenon in interstitials and xenon in vacancy clusters, using interatomic 3
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
IPR creates a population of random-guess potential forms, and applies GA operations of selection, crossover, mutation, and replacement to iteratively create improved interatomic potentials. It terminates either when the error has reached an acceptable threshold, or when a maximum number of iterations has been performed. IPR is, at its core, a relatively standard real-numberencoded GA [37]. In such a scheme, instead of the traditional GA representation of solutions as binary strings, each solution is simply a real-number array of the potential parameters being optimized. As a simple example, if we wished to fit a Lennard-Jones potential of the form of equation (4), where ε and σ are adjustable parameters, the IPR population would consist initially of n random combinations Uk = [ε, σ ], where k is the index of the potential, within some specified search space. Each candidate potential Vi is then used to calculate a series of input structures whose first-principles energies, forces, and stresses are known. The resulting interatomic potential outputs are compared to the reference first-principles values to calculate an error, and this error is the fitness function that IPR minimizes. As noted above, the GA relies on evolutionary operators such as selection, crossover, mutation, and replacement to produce better interatomic potentials over time. Selection describes the process by which we choose ‘parent’ potentials from among our population of n that will be used to generate new ‘child’ potentials during crossover. IPR performs tournament selection, in which two potentials are drawn at random from the population, and the potential with the better fitness value of the two becomes one crossover parent. IPR’s crossover operator mixes two parent potentials Uk and Ul to create a child potential Uchild that has characteristics of each parent: Uchild = w · Uk + (1 − w) · Ul ,
4.1. Analysis of previous potentials
Geng potential: fitting to-X eO 3 . First, we examine the Geng potential to find how accurately it represents xenon in UO2 . In order to better understand these compounds we used DFT to find the relative stabilities of these two structures. We performed four calculations: an O2 molecule (with the O2 molecular energy correction of [38]), gaseous xenon, XeO3 in the CuAu3 structure, and XeO3 in the P21 21 21 space group. The calculated formation energy of the experimentally observed P21 21 21 structure is +3.6 eV/formula unit. The formation energy of XeO3 in the hypothetical CuAu3 structure is even more significantly positive, at +6.6 eV/formula unit. The calculated energy of the P21 21 21 structure is close to the experimentally measured formation energy of +4.2 eV/formula unit [24]. These positive formation energies for both phases means that these XeO3 crystals are strongly thermodynamically unstable and the most stable states for xenon and oxygen are their respective elemental phases. Intuitively, neither structure seems appropriate for fitting a Xe–O interaction for xenon defects in UO2 . In addition to these structures being thermodynamically unstable, any compound of XeO3 will not be appropriate for fitting interactions of xenon in UO2 . While XeO3 -type interactions might be appropriate for modeling xenon and oxygen in a compound of only these two elements, Geng et al use them to model the interaction of xenon with oxygen ions in UO2 . The net result of the fitting procedure for the Geng potential is that the potential describes a very short equilibrium separation between xenon and oxygen (2.8 Å). This short interaction can be traced back to the distance between xenon and oxygen in the CuAu3 structure (2.6 Å). We can use this separation distance along with the xenon–xenon separation (2.4 Å) from the Geng potential to approximate the atomic radii of xenon and oxygen: the Geng potential describes an effective O2− radius of just 0.4 Å, whereas the empirical value for an oxygen ion is 1.24 Å [39]. Thus, we might expect that any situation in which xenon and oxygen come into close proximity with one another will not be well described by the Geng potential. Geng potential: small incorporation sites. We assess the accuracy of the Geng potential by calculating the energetics of inserting a single xenon atom into the UO2 lattice. One possible way to insert xenon into the lattice of UO2 without removing any uranium or oxygen atoms is to insert xenon into the octahedral interstitial position (Xei ) where xenon would strongly repel the neighboring oxygen atoms, causing a strain penalty. Another choice would be to substitute xenon in for a uranium atom and place that uranium atom in a neighboring interstitial site (XeU + Ui ). In addition to xenon repelling the neighboring oxygen atoms and causing a strain penalty, like the Xei case, XeU + Ui should also incur an electrostatic penalty with uranium moved from its ideal fluorite position, (the XeU will have an effective charge of −4e while the Ui would have an effective charge of 4e). Therefore, we would expect that the XeU + Ui configuration should have a larger defect energy than the Xei configuration. We have calculated the incorporation energies of these defects from
where w is a random number between 0 and 1. Mutation simply replaces child parameter values with new, entirely random values with some low probability pmutation (in this work, 0.1) during crossover. Selection is the process by which we replace poorer potentials with hopefully-improved children; here, one GA generation (i.e., iteration) consists of generating and evaluating child potentials to replace the bottom 75% of the population. Performing many such generations gives improved interatomic potentials as the GA ‘learns’ the solution space and locates optima. Here, we use IPR, along with DFT + U molecular dynamics simulations, to parameterize a new xenon potential using the Basak as the base UO2 potential. 4. Results
We begin by benchmarking the Basak/Geng and Morelon/ Chartier potentials against DFT calculations of xenoncontaining defects in UO2 . We consider all combinations of the xenon and UO2 potentials and we refer to each possible potential combination by, first, the base UO2 potential, and second, the add-on xenon potential: Basak/Geng, Morelon/Chartier, Basak/Chartier, and Morelon/Geng. Firstly, we examine the fitting and form of the Geng and Chartier xenon potentials and their ability to reproduce defect energetics from DFT + U calculations. 4
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
distance (0.006 eV), despite xenon being closer to oxygen than to uranium for both defect configurations. This unphysically small Xe–O interaction in UO2 is a direct consequence of the choice to fit to XeO3 , which contains short Xe–O bonds, as discussed above. These small Xe–O interactions mean that even though xenon neighbors eight oxygen atoms in both defects we are considering, they will not contribute significantly to the defect energy for either case. Next, we consider the xenon–uranium interactions. The xenon–uranium distances differ for the two defect sites. The Xei –U distance is 2.74 Å, while the XeU –U distance is 3.87 Å. For the Geng potential, the Xei –U and XeU –U the interaction energies at those distances are 6.63 eV and 0.44 eV, respectively, showing that the Xe–U interaction is fairly long ranged. Although the oxygen atoms are closer to xenon than the uranium atoms, the repulsive interaction between uranium and xenon is much stronger than the interaction between xenon and oxygen, showing that the Xe–U interaction is too strong compared to the Xe–O interaction. The large Xe–U interaction for Xei translates into a larger defect energy than for XeU + Ui , despite the latter having an unfavorable electrostatic penalty. Because the Xe–O interactions are too small and the Xe–U interactions are too large, the Geng potential gives an energy of XeU + Ui that is lower than the energy of Xei . The fitting approach for the Geng potential has created a qualitative failure for the ordering of these defect energies. As discussed above, the Chartier potential accurately describes small incorporation sites such as xenon in an interstitial, but it overestimates the energy of the large incorporation site of the XeSDs [21]. Because the Chartier potential does not accurately reproduce the XeSD incorporation energies, it should be tested for other xenon in defects in large incorporation sites.
Figure 1. The interatomic potential functions of energy versus distance from the Geng, Chartier, and IPR-generated potentials. The interatomic potentials are generally characterized by repulsive interactions only. The Geng potential has the smallest separation for Xe–O and the largest separation for Xe–U. The Chartier Xe–O potential has a low curvature compared to the other xenon interactions, which are treated more like hard spheres. The vertical lines represent the bond distances for atoms that are nearest neighbors to xenon in Xei and XeU + Ui .
5. New Xe–O and Xe–U interatomic potential fitted by IPR
As we have discussed above, the existing Geng and Chartier potentials have shortcomings in describing defects, so we have created a new set of parameters for the Xe–O and Xe–U interactions (leaving the U–U, O–O, and U–O interactions equal to the values given by Basak et al) using iterative potential refinement. IPR works by fitting to DFT + U energies, forces, and stresses with a genetic algorithm. For the IPR-generated Xe–O and Xe–U potentials, we consider only the repulsive interactions of a Born–Mayer potential. For the IPR fitting procedure. we used snapshots from a single DFT + U molecular dynamics calculation of xenon in a Schottky defect cluster with the oxygen vacancies at nearest neighbor positions with respect to one another (XeSD1). The IPR-generated potential was not fit to any other defects. We also examined the effect of using the Morelon potential as the base UO2 potential, and the results were very similar to those from the Basak potential. Table 2 shows the result of the IPR fitting procedure. We can examine the form of the IPR potential compared to the Geng and Chartier potentials. Figure 1 shows the form of the Basak, Chartier, and our IPR-generated potentials as a function of atomic separation. The Geng Xe–O potential
DFT + U and the Basak/Geng potential. From the DFT + U calculations, we find that the energetics follow our intuition with the incorporation energy of XeU + Ui having a larger defect energy than Xei (11.3 eV versus 9.7 eV). However, the Basak/Geng potential finds there is a much larger penalty to create these defects compared to DFT + U and that XeU + Ui is lower in energy than Xei (17.0 eV versus 22.7 eV). For comparison, the Morelon/Chartier potential finds that the two defect configurations are nearly degenerate in energy (11.8 eV for both defects). To find why the interatomic potential energies are not in agreement with DFT + U for these defects, we can examine the energetic components of the potential (Xe–O and Xe–U pairs). For this test, uranium and oxygen ions are placed at their equilibrium, fluorite positions. For both defect cases, the distance between xenon and its eight oxygen neighbors is the same (2.36 Å, for a0 = 5.47 Å). Counter to our intuition, the Geng Xe–O potential has almost no interaction at this 5
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Finally, we calculate the defect energetics of a single large incorporation site, xenon in a double Schottky defect cluster (Xe2 · SD). It is important to note that our IPR potential was fit only to one of the test configurations (XeSD1). For these defects, we calculate the incorporation energy (the energy required to add xenon to an existing defect site). These defect energetics are defined as follows: Xenon incorporation energy (1E Inc. Xe ). This quantity is the energy of incorporating a xenon atom into the octahedral interstitial site of UO2 .
Table 2. Potential parameters for IPR-generated Xe–O and Xe–U
potentials. The form of the potential is of equation (5). The U–O, O–O, and U–U parameters are those of Basak et al [18].
Xe–U Xe–O
A (eV)
ρ (Å)
6607 300 95 913
0.178 0.210
differs significantly from the Chartier and IPR potentials in terms of the effective size for xenon, a result of fitting to XeO3 . The Chartier Xe–O potential differs from the other xenon interactions in the ‘hardness’ (the curvature of the potential near its minimum) of the Xe–O interactions. The difference in hardness may be a source of error in the Chartier potential. In contrast, using the IPR algorithm, we find that the Xe–O and Xe–U IPR potentials have similar, hard curvatures. The Geng and Chartier potentials describe very different Xe–U separations, and the Xe–U separation described by our IPR potential falls between them. Next, we create a quantitative comparison of all possible combinations of UO2 potentials and xenon potentials for a set of defects in a variety of incorporation sites to determine how well each combination of the potentials are able to reproduce DFT + U energetics.
N UO2 1E Inc. − E Xe , Xe = E Xei − N E
(6)
N is the energy of xenon in the interstitial position where E Xe i of a supercell with N UO2 formula units, E N is the energy per formula unit of bulk UO2 , and E Xe is the energy of a single, isolated xenon atom (for the interatomic potentials considered, this quantity is zero). Similarly, we can calculate the incorporation energy of XeU + Ui . Inc. ): XeU + Ui incorporation energy (1E Xe N UO2 1E Inc. − E Xe . XeU +Ui = E XeU ,Ui − N E
(7)
Inc. Schottky defect-xenon atom incorporation energy (E Xe−SD ). This quantity represents the energy to add a xenon atom into an existing Schottky defect cluster.
5.1. Benchmarking interatomic potentials with DFT + U
Now we benchmark the Geng, Chartier, and our IPR-generated interatomic potentials against DFT + U for various xenoncontaining defects. We can determine how accurate the interatomic potentials are by comparison to the energetics of relatively small defects that are practical to calculate using DFT + U . We have created a series of physically relevant benchmarks to compare all three xenon potentials on equal footing. We place xenon in different incorporation sites in UO2 (small, intermediate, and large vacant sites). For small incorporation site, we used the following defects: xenon in an interstitial (Xei ), xenon in an interstitial without relaxing the ideal fluorite UO2 geometry, and xenon on a uranium site with uranium in a neighboring interstitial (XeU + Ui ). The intermediate-sized incorporation sites are xenon in the three configurations of Schottky defect cluster (XeSD1, XeSD2, and XeSD3, in order of increasing oxygen vacancy separation).
N −1 N −1 1E Inc. − E Xe , Xe−SD = E Xe − E
(8)
N −1 where E Xe is the is the energy of a supercell with N sites for UO2 formula units that is missing one UO2 formula unit (the Schottky defect cluster) plus a xenon atom, and E N −1 represents the energy of a supercell with an empty SD. We can also define the incorporation energy for xenon in a double Schottky defect N −2 − E N −2 − E Xe . 1E Xe−2:SDInc. = E Xe
(9)
We then calculated these incorporation energies with DFT + U as well as combinations of the UO2 potentials (Basak and Morelon) with xenon potentials (Geng, Chartier, and IPR). The energies of these xenon-containing defects calculated are shown in table 3. The IPR potentials give the best
Table 3. Comparison of xenon point defects of a 96 atom cell of UO2 calculated with DFT + U and various interatomic potentials. The
Basak and Morelon potentials are used for the U–O, U–U, and O–O interactions, while the Geng, Chartier, and IPR potentials are used for Xe–O and Xe–U interactions. The IPR potential is fitted to a DFT + U molecular dynamics calculation of xenon in SD1. The IPR potentials give the best agreement with DFT calculations, even in configurations which it was not directly fit. Bold energies correspond to the closest to the DFT + U energies. Incorporation energies (eV)
DFT + U
Basak/Geng
Basak/Chartier
Morelon/Geng
Morelon/Chartier
Basak/IPR
Morelon/IPR
Xei (relaxed) Xei (static) XeU + Ui Xe in SD1 Xe in SD2 Xe in SD3 Xe in two SDs
9.73 15.64 11.33 1.06 1.83 1.94 0.27
22.74 41.45 17.03 4.49 4.57 4.60 5.08
12.85 14.61 16.95 4.80 5.26 5.72 5.03
21.72 41.46 17.45 4.57 4.63 4.67 2.68
11.76 14.61 11.76 4.10 4.43 4.61 2.65
9.67 18.70 11.90 0.60 0.87 1.17 0.32
7.96 18.70 9.87 0.39 0.56 0.62 0.14
6
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Figure 2. The percentage difference between the formation energy of DFT + U and interatomic potentials for a series of defects. The xenon is placed in small incorporation sites which have large defect formation energies (Xei , Xei without geometry optimization, and XeU + Ui ), intermediate-sized incorporation sites (xenon in three Schottky defect configurations: XeSD1, XeSD2, and XeSD3, in order of increasing VO –VO separation), and large incorporation sites (xenon in a double Schottky defect). The IPR potential performs better than either Chartier or Morelon, independent of the UO2 potential used. As size of the incorporation site increases, Chartier and Morelon perform much worse.
Figure 3. The relaxation energy error for xenon defects in different
incorporation sites as described by equation (6). The defects were calculated at their DFT + U geometries without geometry optimization. The difference in defect energy between interatomic potential calculations at the DFT + U geometry and DFT + U calculations is a measure of how well the interatomic potential calculations match the DFT + U geometry. Both the Basak/IPR and Morelon/IPR combinations have much smaller differences from DFT + U than the Geng or Chartier containing calculations, even for defects in which IPR was not fit.
agreement with DFT + U calculations, even in configurations to which it was not directly fit. As we have previously shown, there is qualitative disagreement between DFT + U calculations and the Basak/Geng and Morelon/Chartier interatomic potentials for the energetic ordering of XeU + Ui and Xei . The IPR potentials (as well as the Basak/Chartier potential combination), agree with the DFT + U ordering for those defects. For comparison, table 4 shows the defect energetics
calculated with the interatomic potentials obtained using a larger supercell (2592 atoms) for better convergence. Figure 2 shows the percentage error of the interatomic potential defect energies compared to DFT + U calculations. For calculations with the IPR potential, the choice of UO2 potential does not have a large impact on the defect energetics; the IPR potential gives better agreement with DFT + U when 7
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
Table 4. Comparison of xenon point defects of a 2592 atom cell of UO2 calculated with various interatomic potentials. Because of the larger supercell size, the image effects are negligible and these calculations are well converged.
Incorporation energies (eV)
Basak/Geng
Basak/Chartier
Morelon/Geng
Morelon/Chartier
Basak/IPR
Morelon/IPR
Xei XeU + Ui Xe in SD1 Xe in SD2 Xe in SD3 Xe in two SDs
22.61 16.93 4.75 4.84 4.94 4.08
12.83 17.13 5.14 6.05 6.55 3.16
22.20 17.62 4.72 4.76 4.81 4.07
11.84 11.84 4.19 4.60 4.78 2.69
9.62 11.86 0.90 1.51 1.79 0.18
8.15 10.29 0.44 0.66 0.70 0.13
paired with the Basak UO2 potential (which was used for the UO2 interactions during the IPR fitting), but also performs well with the Morelon UO2 potential. The IPR potential (whether paired with Basak or Morelon) is comparable with or outperforms the Geng and Chartier potentials across all incorporation sites. The Basak/Geng, Basak/Chartier, and Morelon/Geng potentials have high error for the Xei and static Xei and low error for the XeU + Ui ; however, with the error near 50%, the XeU + Ui energy difference represents 5.7 eV deviation from DFT + U . The Morelon/Chartier potential has good agreement with DFT + U for the three small incorporation sites. All of the Geng and Chartier derived calculations have a large error for the intermediate-sized incorporation sites of the XeSDs (>100%) and even larger errors for double Schottky defects, which are essentially very small bubbles (∼900% for Chartier and ∼1700% for Geng). The IPR potential has the most accurate energetics for the double Schottky defect by a large margin, even though it was not fit to this geometry, which suggests that the IPR potential is much more transferable to a wide range of geometries than the other potentials. When considering all of the defects calculated, IPR potentials find the best agreement with DFT + U . The deviation of defect energies from DFT values is just one metric we can use to validate the interatomic potentials. We can also compare the geometries of the DFT + U and interatomic potential calculations, albeit indirectly. We calculate the defect energetics from the potentials at the DFT + U geometries and compare to the DFT + U defect energies. We normalize the difference in defect energies to the DFT + U Relax defect energy to find a relaxation energy error, E DFT Geometry , which is defined as Relax E DFT Geometry =
Pot DFT+U ) (E DFT Geo − E , E DFT+U
and intermediate-sized incorporation sites calculated with the IPR potential have markedly smaller relaxation energy error compared to Geng and Chartier. For the double Schottky defect, the differences among the three xenon potentials are smaller, but the IPR-generated potential still performs the best. Overall, the IPR potential (paired with Basak or Morelon) performs best not only on energetics, but also on geometries. 6. Conclusions
Using iterative potential refinement, or IPR, we have created a new xenon interatomic potential for use with the Morelon and Basak UO2 potentials. This IPR-generated potential was fit by a genetic algorithm to snapshots of DFT + U molecular dynamics of xenon in a Schottky defect cluster. We have compared defect energetics of xenon in a variety of incorporation sites (small, intermediate, and large vacant sites) using combinations of xenon potentials (Geng, Chartier, and IPR-generated) and UO2 potentials (Basak and Morelon) to DFT + U calculated energetics. The IPR-generated potential, despite only being fit to a single defect configuration in a intermediate-sized incorporation site, gives accurate energetics for small and large incorporation sites as well, pointing to its predictive power. The existing xenon potentials overestimate the energetics of xenon incorporation across all sizes of incorporation sites tested. The existing potentials do not yield accurate defect energetics due, in part, to their fitting methodologies. The Geng potential describes the interaction between oxygen and xenon as too short, and the Chartier potential finds that the xenon–oxygen interaction is not ‘hard’ enough. The IPR-generated potential performs better than both the Geng and Chartier potentials across a range of defects, including geometries to which it was not fit such as xenon in an interstitial and xenon in a double Schottky defect, while adding no additional fitting parameters compared to the Chartier or Geng potentials. The IPR-generated potential will facilitate calculations of xenon bubbles in a variety of morphologies and sizes.
(10)
Pot where E DFT Geo is the defect energy as calculated with an interatomic potential as described in equations (2)–(5), but calculated at the DFT geometry and E DFT+U is defect energy from DFT + U . Figure 3 shows relaxation energy error described by equation (10) for xenon defects in different incorporation sites. As with the defect energy comparison above, the IPR potentials perform the best, even for defect configurations for which IPR was not fit. However, in contrast to the defect energy comparison above, the relaxation energy error does not depend strongly on incorporation sites. The small
References [1] Zhou F and Ozolins V 2009 Phys. Rev. B 80 125127 [2] Tiwary P, van de Walle A and Grønbech-Jensen N 2009 Phys. Rev. B 80 174302 [3] Tiwary P, van de Walle A, Jeon B and Grønbech-Jensen N 2011 Phys. Rev. B 83 094104 8
J. Phys.: Condens. Matter 26 (2014) 105501
A E Thompson et al
[4] Hanken B E, Stanek C R, Grønbech-Jensen N and Asta M 2011 Phys. Rev. B 84 85131 [5] Freyss M, Vergnet N and Petit T 2006 J. Nucl. Mater. 352 144 [6] Petit T 2002 Study of the atomic rare gas behavior by Ab Initio calculations Fission Gas Behaviour in Water Reactor Fuels (Cadarache: OECD Publishing) [7] Yun Y, Kim H, Kim H and Park K 2008 J. Nucl. Mater. 378 40 [8] Yun Y, Eriksson O and Oppeneer P M 2009 J. Nucl. Mater. 385 510 [9] Yun Y, Eriksson O, Oppeneer P M, Kim H and Park K 2009 J. Nucl. Mater. 385 364 [10] Jeon B, Asta M, Valone S M and Grønbech-Jensen N 2010 Nucl. Instrum. Methods Phys. Res. B 268 2688 [11] Grimes R W and Catlow C R A 1991 Phil. Trans. R. Soc. A 335 609 [12] Govers K, Lemehov S, Hou M and Verwerft M 2007 J. Nucl. Mater. 366 161 [13] Govers K, Lemehov S, Hou M and Verwerft M 2008 J. Nucl. Mater. 376 66 [14] Govers K, Lemehov S E and Verwerft M 2010 J. Nucl. Mater. 405 252 [15] Arima T, Yamasaki S, Inagaki Y and Idemitsu K 2005 J. Alloys Compounds 400 43 [16] Ball R G J and Grimes R W 1992 J. Nucl. Mater. 188 216 [17] Geng H Y, Chen Y, Kaneta Y and Kinoshita M 2008 J. Alloys Compounds 457 465 [18] Basak C B, Sengupta A K and Kamath H S 2003 J. Alloys Compounds 360 210 [19] Jackson R A and Catlow C R A 1985 J. Nucl. Mater. 127 161 [20] Chartier A, Van Brutzel L and Freyss M 2010 Phys. Rev. B 81 174111
[21] Morelon N-D, Ghaleb D, Delaye J-M and Van Brutzel L 2003 Phil. Mag. 83 1533 [22] Thompson A and Wolverton C 2013 Phys. Rev. B 87 104105 [23] Templeton D H, Zalkin A, Forrester J D and Williamson S M 1963 J. Am. Chem. Soc. 85 817 [24] Bard A J, Parsons R and Jordan J 1985 Chemistry IU of P and A. Standard Potentials in Aqueous Solution (Boca Raton, FL: CRC Press) [25] Thompson A E, Meredig B, Stan M and Wolverton C 2014 J. Nucl. Mater. 446 155–62 [26] Kresse G and Hafner J 1993 Phys. Rev. B 47 558 [27] Kresse G 1993 Thesis, Technische Universit¨at Wien, Vienna, Austria [28] Kresse G and Furthm¨uller J 1996 Comput. Mater. Sci. 6 15 [29] Kresse G and Furthm¨uller J 1996 Phys. Rev. B 54 11169 [30] Perdew J P and Ziesche P (ed) 1991 Electron. Struct. Solids vol. 11 (Berlin: Akademie Verlag) [31] Dudarev S L, Botton G A, Savrasov S Y, Humphreys C J and Sutton A P 1998 Phys. Rev. B 57 1505 [32] Meredig B, Thompson A, Hansen H A, Wolverton C and van de Walle A 2010 Phys. Rev. B 82 195128 [33] Thompson A and Wolverton C 2011 Phys. Rev. B 84 134111 [34] Gale J D and Rohl A L 2003 Mol. Simul. 29 291 [35] Ercolessi F and Adams J B 1994 Europhys. Lett. (EPL) 26 583 [36] Brommer P and G¨ahler F 2007 Modelling Simul. Mater. Sci. Eng. 15 295 [37] Michalewicz Z, Janikow C Z and Krawczyk J B 1992 Comput. & Math. Appl. 23 83 [38] Wang L, Maxisch T and Ceder G 2006 Phys. Rev. B 73 195107 [39] Shannon R D 1976 Acta. Cryst. A 32 751
9
