Derivation of Transferable Pair Potentials and the Calculation of

Subscriber access provided by University of Sunderland

C: Physical Processes in Nanomaterials and Nanostructures

Derivation of Transferable Pair Potentials and the Calculation of Intrinsic Defect Properties for Xenotime. Geoffrey L Cutts, Joseph A Hriljac, and Mark Stuart David Read J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b06978 • Publication Date (Web): 04 Oct 2018 Downloaded from http://pubs.acs.org on October 9, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Derivation of Transferable Pair Potentials and the Calculation of Intrinsic Defect Properties for Xenotime Geoffrey L. Cutts,†,‡ Joseph A. Hriljac,† and Mark S. D. Read∗,† †School of Chemistry, University of Birmingham, Edgbaston, Birmingham, B15 2TT, U.K. ‡Present address: Diamond Light Source Ltd., Diamond House, Harwell Science & Innovation Campus, Didcot, Oxfordshire,OX11 0DE. [email protected] E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract A new force field has been empirically derived that is transferable across the YPO4 , Y2 O3 and P2 O5 phases, utilising a reverse Monte Carlo (RMC) method. This method employs a simulated annealing technique with a logarithmic quench to fit potential parameters to observed crystallographic and mechanical properties, producing a force field suitable for simulating radiation damage events in an atomistic molecular dynamics regime. These potentials are used to investigate the defect properties of xenotime, where a wide range of intrinsic defects including Schottky, Schottky-like, Frenkel pairs and anti-site defects have been investigated, both at infinite dilution and as defect clusters. A common feature in the lowest energy defect configurations was the presence of polymerised phosphate tetrahedra, forming P2 O7 units. The trend in the formation energies for the Frenkel pair defects at infinite dilution was in good agreement with previously published simulations. However, the binding energy associated with the aggregation of point defects was found to have a profound impact on the defect formation energies, significantly lowering the formation energy of the phosphorous Frenkel pair in particular. The intrinsic defect calculations presented here have been compared with previous work in zircon, to gain insight into differences that may contribute to the disparity in the radiation resistance of the two minerals.

1

Introduction

Of the waste produced over the lifetime of a nuclear power plant (including decommissioning), only a very small proportion (∼ 0.1%) is categorised as high level waste (HLW). 1 However, this waste fraction contains ∼ 95% of the radioactivity. 2 One of the main contributors to this waste originates from the reprocessing of spent nuclear fuel (SNF). The HLW contains a complex mixture of radionuclides, including fission products (e.g. Cs, Sr, Nb, Y and Zr etc.) 3–5 and trans-uranics such as Pu, Am and Cm. 6,7 Although the waste can be stored at the processing facility on a temporary basis while the short lived radioisotopes decay, a 2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


more long term strategy is required for dealing with the long lived radioisotopes. One of the most promising solutions is the construction of a geological disposal facility (GDF), using a multi-barrier approach. Here, the waste is immobilised into a waste form that can be stored more easily and safely than the liquid waste. 8 Vitreous waste forms such as borosilicate glass have been studied extensively and the process of immobilizing the waste well established. Despite the capacity for high waste loading and the ability to accommodate a wide range of waste elements, these waste forms are inherently metastable. Indeed, their stability has been questioned over geological timescales, an important consideration since any instability may result in de-vitrification and possible release of the waste elements. Ceramic waste forms offer a possible alternative with many natural analogues existing, demonstrating good durability over geological timescales for minerals containing reasonably high concentrations of actinides and fission products. 9 Orthophosphate and orthosilicate minerals with the general formula ABO4 , such as monazite (LaPO4 ) and zircon (ZrSiO4 ) have received increased interest as ceramic waste forms. It is well documented that zircon is regularly found to be fully or partially amorphised due to the accumulation of damage from radioactive decay events (metamict) and is a common mineral used in geochronology. 10–16 Xenotime (YPO4 ) is another accessory mineral in plutonic rocks 17 , crystallising in the tetragonal space group I41 /amd 18 , isostructural with zircon, Figure 1. In contrast to zircon, natural samples are rarely found to be metamict. One metric often used to compare the resistance of a material to radiation damage is the critical temperature of annealing (TC ). The high TC value for zircon (∼1100 K) compared with that of xenotime (∼512 K) is indicative of the increased resistance of xenotime to the accumulation of radiation damage. 19 Whilst zircon has been extensively studied both experimentally 20–25 and computationally 26–30 , xenotime has received much less attention. 31–33 As both zircon and xenotime crystallise in the same space group, they offer an excellent opportunity to investigate the inherent differences in the behaviour of phosphate and silicate minerals to radiation damage. First, the intrinsic defect properties of xenotime must be understood, with relevance to applica3



Page 4 of 36

tions including nuclear engineering and the development of electronic devices. 34 This work aims to derive a robust potential set for xenotime (transferable across YPO4 , Y2 O3 and P2 O5 phases) prior to a rigorous investigation of the defect chemistry. This is important to further the understanding of the disparity in the radiation resistance of the two minerals.

(a)

(b)

Figure 1: Crystal structure of xenotime viewed down the c axis (a) and a slice in the bc plane (b).

2 2.1

Computational Methodology Interatomic potential model

The simulations reported here are based upon atomistic modelling techniques in which the total energy of the system, U , is, in concept, calculated as an expansion in terms of the contributions to the potential energy from all atoms in the system, over all length scales. In practice, this expansion is truncated at the three-body term, φijk , and contributions from pairwise interactions, φij , are truncated at a distance , rmax . The force field developed in this work utilises pair potentials of the Buckingham and Morse forms defined in Equations 1 and 2, respectively. Here qi and qj represent the ionic charges, rij is the interatomic separation and all other constants are adjustable parameters which may be fitted empirically. qi qj φij (rij ) = + Aij exp rij 4

−rij ρij


−

Cij 6 rij

(1)

Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


φij (rij ) =

q i qj + Dij (1 − exp (−αij (rij − r0 )))2 − 1 rij

(2)

A harmonic three-body potential has been applied to help constrain the geometry of the tetrahedral phosphate units, Equation 3, where (θijk − θ0 ) is the deviation from the equilibrium angle and kijk is the force constant which may also be fitted empirically.

φijk (θijk ) =

kijk (θijk − θ0 )2 2

(3)

Once a suitable force field is defined, the geometry is optimised by altering the lattice parameters and atomic positions to minimise the forces acting upon the atoms and therefore the total energy of the system. Properties of the material such as elastic constants are subsequently calculated from the minimum energy structure. All simulations performed in this work were carried out in the static regime using the General Utilities Lattice Program (GULP) 35 and optimised at constant pressure.

2.2

Defect Calculations

Once the geometry of the pure bulk lattice has been optimised and relaxations inherent in the potential model accounted for, lattice defects are introduced. The Mott-Littleton 36 approach is adopted to calculate defect energies efficiently, taking into account the lattice relaxation that occurs about point defects or defect clusters. Here, the ions within a spherical region centred about the defect site, which exhibit the majority of the lattice relaxations, are allowed to relax to zero strain explicitly. This is referred to as region I. Region II comprises all remaining atoms outside of region I, extending to infinity, and is subdivided into regions IIa and IIb. In this region, the displacements of lattice ions are assumed to occur as a result of the polarisation due to the effective charge of the defect. The contribution from region IIb is calculated implicitly whilst region IIa acts as a buffer region in which the energy contributions for each atom are calculated explicitly, ensuring smooth convergence between the explicitly summed region I and the continuum region IIb. Radii for the spherical regions 5



I and IIa are determined via convergence testing and selected so as to maximise the accuracy of the calculated defect energy, whilst attempting to minimise the computation time. The defect energy is defined as the difference between the total energy of the defective system and that of the perfect lattice.

3

Results and Discussion

3.1

Potential Models

To the knowledge of the authors there are four sets of previously derived interatomic potentials for the xenotime system. 31–33 The suitability of these potentials for the simulation of radiation damage has been analysed and it was immediately clear that only one set of potentials (Gao et al. 31 ) is transferable across the YPO4 , Y2 O3 and P2 O5 phases, based on the charges of the constituent ions. This is essential to model reliably phase decomposition and segregation phenomena that may occur upon irradiation. These potentials reproduce the crystallographic and mechanical properties of YPO4 and Y2 O3 to a high degree, but are unable to reproduce accurately the crystallographic properties of P2 O5 , as summarised in Tables 2–4, respectively. As such derivation of a new self-consistent force field was required.

3.2

Derivation of Interatomic Potentials

A Reverse Monte Carlo (RMC) fitting procedure was developed to fit a range of potential parameters within multiple phases, simultaneously. The pseudo-code for the procedure is listed in Figure 2 A cost function is minimised by randomly moving around the parameter space, with moves accepted or rejected with a probability dictated by the difference between the quality of the fit for both the previous, LSq_Old, and current, LSq_New, cycles. The cost function, LSqw , is calculated as the weighted sum of the squared relative differences to n observed properties across N phases following Equation 4. The squared relative difference for each observable value, yj , for each phase is weighted by a factor βj and the 6


Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33

RMCFitting ( c y c l e s ) c y c l e := 0 s t e p S i z e := 0 . 0 0 5 p a r a m e t e r s := s t a r t i n g p o t e n t i a l p a r a m e t e r s CALL GULP( p a r a m e t e r s ) LSq_Old := g e t F i t Q u a l i t y ( p a r a m e t e r s ) p a r a m e t e r s _ B e s t := p a r a m e t e r s LSq_Best := LSq_Old WHILE c y c l e ≤ c y c l e s : o p t i m i z a t i o n C r i t e r i a := FALSE WHILE o p t i m i z a t i o n C r i t e r i a = FALSE rand := a random p o t e n t i a l p a r a m e t e r step_max := rand ∗ s t e p S i z e s t e p := a random number s u c h t h a t −step_max ≤ s t e p ≤ step_max and s t e p 6= 0 rand_new := rand + s t e p IF rand_min ≤ rand_new ≤ rand_max o p t i m i z a t i o n C r i t e r i a := TRUE IF o p t i m i z a t i o n C r i t e r i a := TRUE parameters_new := updated p o t e n t i a l p a r a m e t e r s CALL GULP( parameters_new ) IF o p t i m i z a t i o n c r i t e r i a a r e not met o p t i m i z a t i o n C r i t e r i a := FALSE LSq_New := g e t F i t Q u a l i t y ( parameters_new ) A c c e p t := c h e c k A c c e p t a n c e ( LSq_Old , LSq_New ) IF A c c e p t = TRUE p a r a m e t e r s := parameters_new LSq_Old := LSq_New IF LSq_Old < LSq_Best LSq_Best := LSq_Old p a r a m e t e r s _ B e s t := p a r a m e t e r s INCREMENT c y c l e RETURN p a r a m e t e r s _ B e s t END

Figure 2: Fitting routine pseudo-code total function weighted for each phase by a factor of αi . This allows for particular properties to be weighted more heavily based on the user’s priorities and confidence in the accuracy of the experimental values, ultimately biasing the fit towards the phases of primary interest. Properties such as unit cell parameters, atomic positions, elastic constants and the bulk modulus are calculated using routines embodied within the GULP code, subsequent to successful geometry optimisation.

LSqw =

N X i

αi

n X j

βj

|y − yj | yj

2

(4)

A proportion of ‘bad’ moves about the parameter space are accepted in an attempt to avoid trapping within local minima, whilst searching for the global minimum, Equation 5. Here, the value of LSqw(P c) is calculated using the same weighting factors as discussed previously and is equivalent to the cost function for a fit with a uniform percentage difference of P c % across all observable properties used within the fitting. Simulated annealing regimes were included to provide a balance between exploration and exploitation of the parameter 7



Page 8 of 36

space. This systematically lowers the value of P c with each successive cycle, Figure 3. Both linear and logarithmic quenching routines were made available. The value of P clin for cycle Ni of N cycles for the linear quench, P clin , is described simply following Equations 6 and 7 by defining a maximum and a minimum percentage difference, P cmax and P cmin , respectively. The value of P clog for the i’th cycle following the logarithmic quench is calculated following Equations 8 and 9. Here, an additional variable, P cmid , is defined as the value of P clog to be achieved when Ni = N/2. This allows the first half of the search to be biased towards exploration with a higher percentage of ‘bad’ moves being accepted, whilst allowing the latter half of the search to be focussed on exploitation. This effect is illustrated in Figure 3 for the logarithmic regime.

−∆LSqw LSqw(P c)

P (x0 ) = min 1, exp δP c =

!!

(P cmax − P cmin ) (N − 1)

P clin = P cmax − (Ni − 1)δP c

(5) (6) (7)

2 (ln(P cmid ) − ln(P cmax )) δP c(−N + 2)

(8)

P clog = P cmax exp (klog (P clin − P cmax ))

(9)

klog =

Optimisation criteria are defined to avoid searching the parameter space in areas which produce models with unphysical properties. Along with the adjustable parameters outlined previously for the Buckingham and Morse potentials, the force constant for three-body terms and a scaling factor for the ionic charges were also included in the fitting procedure. This factor scales the formal charges uniformly across all species to maintain charge neutrality for all phases. The parameter space was constrained for each parameter by setting upper and lower limits, as defined by the user. If the proposed candidate for a cycle was within the allowed search area then the model was checked against further optimization conditions: 8


Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


• All valid candidates must achieve geometry optimisation under the criteria defined within the GULP package. • The computation time for each candidate must not exceed a set limit to prevent excessive run times for unphysical models which fail to converge during the geometry optimisation stage. • The elastic constants must meet the symmetry–specific necessary and sufficient conditions outlined by Mouhat et al. 37 This procedure was iterated over a defined number of cycles and the set of parameters producing the model with the smallest value of LSqw returned. This process was subsequently repeated for multiple runs to attempt to locate the global minimum.

9



(a)

Page 10 of 36

(b)

Figure 3: Comparison of the linear and logarithmic annealing regimes (a) and the effect of the logarithmic regime on the acceptance probability (b) . In practice, the equilibrium bond length for the P-O Morse potential and ideal bond angles for three-body terms were fixed. The three-body term relating to the P-O-O angle within the phosphate tetrahedra originates from Gao et al., with all other parameters allowed to vary within constraints. The resulting force field parameters are summarised in Table 1 and a comparison of the fit to the observed properties for this force field with that of Gao et al. for the YPO4 , Y2 O3 and P2 O5 phases is summarised in Tables 2–4, respectively. The force field derived here reproduces both the crystallographic and mechanical properties of all three phases to an acceptable degree of accuracy, with only a small sacrifice in accuracy to achieve transferability across these key phases. This transferability enhances the robustness of this force field, with the potentials adequately describing the interactions over a wider range of interatomic separations.

3.3

Modifications to potentials

The force field described in the previous section considers the forces acting upon the atoms around the equilibrium bond lengths successfully. However, at shorter interatomic separations the atoms experience unphysical attractive forces aptly known as the ‘Buckingham

10


Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Force field parameters Buckingham Y1.11 · · · O−0.74 O−0.74 · · · O−0.74 Morse P1.85 · · · O−0.74

Aij eV 27025.72941 1458.23941 Dij eV 3.51305

ρij Å 0.22225 0.33460 αij Å−2 3.76122

Three-body harmonic

kijk

θ0 o

P1.85 · · · O−0.74 · · · O−0.74 O−0.74 · · · P1.85 · · · P1.85

eV rad−2 3.5000 20.9326

109.47 135.58

Cij Å6 188.27260 99.41669 r0 Å6 1.54000 rmax (ij) Å 1.8 1.8

rmax (ik) Å 1.8 1.8

rmax Å 10.0 10.0 rmax Å 10.0 rmax (jk) Å 3.2 3.2

catastrophe’. This becomes especially troublesome when modelling radiation damage within a molecular dynamics regime in which atoms are likely to come into much closer proximity to each other than would be the case for the perfect lattice. Similarly, the grid search method for isolating defect configurations discussed in section 3.4 may encounter a similar issue, with interstitial or relaxed atoms coming into close proximity with other lattice atoms. To overcome this, the offending potentials are splined to the strongly repulsive ZBL potential 43 , Figure 4. The potentials were splined using a fifth order exponential polynomial, ensuring continuity of the first and second derivatives at the cut-offs, summarised in Table 5. The total potential was then also subjected to a linear shift to ensure that the potential was zero at the potential cut-off in both value and first derivative. The version of GULP used in this study did not support the ZBL functional form and as such the splined potential was approximated using a least squares fit to the linear sum of a third order exponential polynomial and a Morse potential, Equation 10. Here all terms were treated as adjustable parameters and the resulting potential is illustrated in Figure 5. The potential parameters are listed in Table 6. This modification to the force field made only a marginal impact on the overall fit of the equilibrium properties to the observed experimental

11



Page 12 of 36

Table 2: Comparison of predicted and observed properties for YPO4 Property

Observed

Lattice Parameters (Å) (a) a 6.8947 c 6.0276 Elastic Constants (GPa) (b) C11 220.0 C12 55.0 C13 86.0 C33 332.0 C44 64.6 C66 17.3 Interatomic Separation (Å) (a) Y· · · O 2.31 Y· · · O 2.38 P· · · O 1.54 Bulk Modulus (G Pa) (c) B0 144 Lattice Enthalpy (eV mol−1 )

∆%

This study Calculated

−0.29 0.39

6.8747 6.0511 282.1 44.0 84.7 278.6 66.4 14.9

28.23 −19.99 −1.51 −16.09 2.74 −13.94

2.32 2.42 1.52

0.26 1.80 −1.27

144

0.00

Gao et al. Calculated 6.7658 6.2349 319.5 56.2 97.7 284.6 76.6 24.4 2.29 2.51 1.49 164

∆%

−1.87 3.44 45.23 2.19 13.65 −14.29 18.56 40.87 −0.86 5.63 −3.13 13.97

−57.99 (a)

Ref. 18 , (b)

Ref. 38 , (c)

Ref. 39

Figure 4: Splined potential for the Y-O interaction with the Buckingham, fifth order exponential polynomial and ZBL regions highlighted in red, green and black respectively.

12


Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Comparison of predicted and observed properties for Y2 O3 Property

Observed

Lattice Parameters (Å) (a) a 10.604 Interatomic Separation (Å) (a) Y· · · O 2.29 Y· · · O 2.27 Y· · · O 2.34 Y· · · O 2.36 Bulk Modulus (G Pa) (b) B0 146 Lattice Enthalpy (eV mol−1 )

∆%


Gao et al. Calculated

∆%

10.513

−0.86

10.596

−0.07

2.27 2.26 2.29 2.25

−1.05 −0.27 −2.10 −4.80

2.29 2.28 2.25 2.31

−0.13 0.44 −4.00 −1.96

−19.64

1470.86

117 −27.02

(a)

Ref. 40 , (b)

Ref. 41

Table 4: Comparison of predicted and observed properties for P2 O5 (ortho) Property

Observed

Lattice Parameters (Å) (a) a 16.314 b 8.115 c 5.265 Interatomic Separation (Å) (a) P· · · O 1.44 P· · · O 1.56 P· · · O 1.58 P· · · O 1.58 Lattice Enthalpy (eV mol−1 )

∆%


Gao et al. Calculated

∆%

16.435 8.164 5.260

0.74 0.60 −0.10

16.563 8.364 4.727

1.52 3.07 −10.22

1.51 1.53 1.53 1.53

4.69 −2.16 −3.29 −3.45

1.46 1.50 1.52 1.52

1.72 −3.87 −3.68 −3.99

−85.58 (a)

Ref. 42

data for each phase. These modified pair potentials were therefore deemed suitable for performing defect calculations.

13



Page 14 of 36

Table 5: Cut-offs used for ZBL spline Interaction · · · O−0.74 O−0.74 · · · O−0.74

rmin (Å) 0.24 0.40

Y1.11

rmax (Å) 1.25 1.37

3 2 + Dij (1 − exp (−αij (rij − r0 )))2 − 1 (10) + B3 rij φij = Aij exp B0 + B1 rij + B2 rij

Figure 5: Comparison of the ZBL splined and exponential polynomial + Morse potential functions. Table 6: Modified potential parameters Exponential Polynomial 1.11 Y · · · O−0.74 O−0.74 · · · O−0.74 Morse Y1.11 · · · O−0.74 O−0.74 · · · O−0.74

Aij eV −3.937 × 10−6 3.232

B0 16.316 9.358

B1 Å−1 7.995 −25.539 αij Å−2 2.143 1.404

Dij eV 0.444 0.017

14


B2 Å−2 −8.617 55.926

B3 Å−3 0.497 −52.763 r0 Å 2.405 3.782

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.4

Intrinsic Defects

When simulating defect configurations it is simple to define vacancy sites based on the positions of the lattice atoms within the unit cell. However, the position of prospective interstitial atoms is much more ambiguous. In this study, the location of interstitial atoms as both point defects and as Frenkel pair clusters were thoroughly investigated by utilising a 3D grid set on the unit cell which subdivided the unit cell into multiple small regions of equal volume. The unit cell was subdivided into 15625 regions (25 × 25 × 25) for the point defect calculations and 1000 regions (10 × 10 × 10) for the defect cluster calculations. An interstitial atom was placed at the centre of each region sequentially and allowed to relax freely. The vacancy position was fixed throughout the defect cluster calculations. An initial search was performed using region I and IIa sizes of 7 and 14 A, respectively, centred on the geometric centre of the interstitial and vacancy sites. The unique defect configurations were filtered using geometric transformations, with a new set of vacancy and interstitial sites identified for each configuration based on a radius method. Here, the atomic positions of the lattice atoms before optimisation are used as nodes. After geometric relaxation, the unique defect configurations were filtered using geometric transformations (i.e. reflection and rotation). Subsequently a new set of vacancy and interstitial positions were defined for each configuration using a radius based method. When multiple atoms were associated with a single node, only the closest atom claimed the node and all other atoms were classed as interstitials. Once these new vacancy and interstitial sites had been identified for each configuration, further defect calculations were performed where convergence was achieved using region I and IIa sizes of 17 and 34 A, respectively. The new defect centre for the defect calculation was taken as the geometric centre of the new defect sites. 3.4.1

Schottky and Schottky-like defects

Defect energies were initially calculated for isolated point defects and vacancy defect clusters, summarised in Table 7. These can be combined following Equations 11– 13 to calculate the 15



Page 16 of 36

relative formation energies for Schottky and Schottky-like defects, summarised in Table 8. Throughout this work defects contained within braces denote a defect cluster. Those not grouped within braces are point defects at infinite dilution.

000 00000 •• × × * YY + P× P + 4OO ) VY + VP + 4VO + YPO4

(11)

000 •• × * 2YY + 3O× O ) 2VY + 3VO + Y 2 O3

(12)

00000 •• × * 2P× P + 5OO ) 2VP + 5VO + P2 O5

(13)

The formation energy calculated from the point defects at infinite dilution do not take into account the binding energy associated with the aggregation of point defects. As such the binding energy for defect clusters can be calculated as the difference between the formation energy of the Schottky defect cluster and the formation energy at infinite dilution. The lowest energy Schottky and Schottky-like clusters are illustrated in Figure 6. Whilst the YPO4 Schottky and P2 O5 Schottky-like defect clusters experience minimal relaxation in the region surrounding the defects, the Y2 O3 Schottky-like cluster causes considerable deformation to the lattice. Here, localised phosphate units become polymerised forming both P3 O10 and P2 O7 units. These deformations to the lattice are manifest in the relatively small binding energy for this cluster. At infinite dilution the Y2 O3 Schottky-like defect is found to be the lowest energy of the three Schottky and Schottky-like defects investigated, with a formation energy of 20.34 eV mol−1 . However, the clustering of the vacancy point defects for the YPO4 Schottky defect results in a significant binding energy, reducing the formation energy to 6.05 eV mol−1 . Once the binding energy is taken into consideration, it can be seen that the YPO4 Schottky cluster is the lowest energy of the three Schottky and Schottky-like defect mechanisms.

16


Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 7: Vacancy defect formation energies Defect V000 Y 00000 VP V•• 00000 O •• V + 4V 000 P 00000 O •• VY + VP + 4VO 000 2VY + 3V•• O 00000 2VP + 5V•• O N.B. Braces define defects contained within a defect cluster.

Defect Formation Energy eV 14.71 48.25 5.98 53.25 64.04 38.52 95.74

Table 8: Formation and binding energies calculated for the Schottky and Schottky-like defects Schottky Defect

YPO4 Y2 O3

P2 O5

Defect 00000 •• V000 Y +VP + 4VO 000 00000 •• V + VP + 4VO Y000 •• VY + V00000 P + 4VO 2V000 + 3V•• O Y •• 2V000 + 3V Y O •• 2V00000 + 5V P O •• 00000 V00000 + V + V + 4V•• P O P O 00000 2VP + 5V•• O

17


Formation Energy eV mol−1 28.89 9.96 6.05 20.34 11.48 40.82 21.89 10.14

Binding Energy eV mol−1 − −18.94 −22.84 − −8.86 − −18.94 −30.68


(a)

(b)

Page 18 of 36

(c)

Figure 6: Lowest energy Schottky and Schottky-like defect clusters for (a) YPO4 , (b) Y2 O3 and (c) P2 O5 . Vacancy sites are highlighted in orange.

3.4.2

Frenkel pair defect configurations

The formation energies for the three possible Frenkel pair defects can be calculated from the defect energies of the constituent vacancy and interstitial defects at infinite dilution or using defect cluster calculations following Equations 14– 16.

••• × * YY ) V000 Y + Yi

(14)

••••• * 00000 P× P ) VP + Pi

(15)

00 * •• O× O ) VO + Oi

(16)

Yttrium interstitial point defects In total, 11 unique yttrium interstitial point defects were identified ranging from −5.90 to −3.49 eV. The lowest energy of these consisted of a six coordinate yttrium site located at the centre of the channel parallel to the c axis, edge sharing with one and corner sharing with four phosphate tetrahedra, Figure 7. This defect is comparable with that identified by Gao et al. Phosphorus interstitial point defects A total of 18 unique configurations were isolated for the phosphorus interstitial point defect, with defect energies ranging from −27.41 to 3.94 eV. The lowest energy isolated phosphorus interstitial was found to occupy a four 18


Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Lowest energy configuration for an yttrium interstitial point defect. coordinate site down the channel parallel to the c axis, Figure 8(a), with a defect energy of −26.81 eV. Here, the phosphorous interstitial is corner sharing with four neighbouring phosphate tetrahedra to form a P5 O16 unit. This is in agreement with Gao et al. and is comparable with the lowest energy silicon interstitial site in zircon, determined by Crocombette. 26 The presence of a secondary phosphorous Frenkel pair defect can be seen to lower the formation energy of the interstitial, with two such configurations illustrated in Figure 8(b) and (c) , both with a defect energy of −27.41 eV. In both cases the interstitial ion causes polymerisation of the local phosphate tetrahedra. The interstitial atom in the configuration illustrated in Figure 8(b) is corner sharing with three phosphate tetrahedra to form a P4 O13 unit. A second smaller P2 O7 chain is also formed, resulting from the polymerisation of two neighbouring phosphate units. The residual oxygen atom from this polymerisation is highlighted in blue and corresponds to the non-bridging oxygen atom coordinated to the phosphorus interstitial ion. In contrast, the interstitial ion in the configuration illustrated in Figure 8(c) causes the formation of only one polymerised P5 O16 chain. The defect configuration illustrated in Figure 8(d) shows a split interstitial centred on an yttrium lattice site. Here, the yttrium and phosphorous atoms are split diagonally into neighbouring channels at distances of 1.95 and 2.18 A from the yttrium vacancy site, respectively. The yttrium atom adopts a six-coordinate site similar to that discussed in the previous section whilst the phosphorus interstitial is corner sharing with four phosphate tetrahedra, forming a P5 O16 unit. The defect energy calculated for this defect is −27.11 eV, which is only 0.3 eV greater than the energy of the most favourable configurations. 19



(a)

(b)

Page 20 of 36

(c)

(d)

Figure 8: Phosphorus interstitial point defect configurations. Residual oxygen atoms and vacancy sites are highlighted in blue and orange, respectively. Oxygen interstitial point defects A total of 17 unique oxygen interstitial point defects were isolated with defect energies ranging from 1.69 to 12.21 eV. The lowest energy defect is found to be located in the channel running parallel to the c axis, bridging between two yttrium lattice atoms, Figure 9(a). A split interstitial was identified with the oxygen atoms at ca 1.03 A from the vacancy site, Figure 9(b). The interstitial is split parallel to the a axis with both oxygen atoms bridging between a phosphorus and an yttrium atom, increasing their coordination to five and nine, respectively. The defect energy calculated for this configuration is 2.95 eV. This higher energy split interstitial is comparable with the lowest energy oxygen interstitial configuration determined by Crocombette in zircon 26 . However, a lower energy split interstitial was isolated, Figure 9(c), where the oxygen interstitials are split along the (101) vector at a distance of ca 1.15 A from the vacancy site. Here, the split interstitial causes a rotation of the phosphate tetrahedron, maintaining the coordination of the phosphorous and yttrium atoms. This stabilises the split interstitial, with a defect energy of 1.82 eV. This defect is similar to that previously identified by Gao et al., although it was not found to be 20


Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the lowest energy defect configuration for an oxygen interstitial in this study.

(a)

(b)

(c)

Figure 9: Oxygen interstitial point defect configurations. Vacancy sites are highlighted in orange.

Yttrium Frenkel pair clusters A total of 22 unique configurations were identified containing an yttrium Frenkel pair cluster, with defect energies ranging from 4.50 to 10.16 eV. The lowest energy configuration containing an isolated Frenkel pair consisted of a six coordinate yttrium interstitial within the channel parallel to the c axis, at a distance of 5.46 A from the vacancy site, Figure 10(a). The defect energy for this configuration is 6.53 eV. A similar, lower energy defect configuration is illustrated in Figure 10(b), with a defect energy of 4.93 eV. Here the configuration is stabilised by the polymerisation of neighbouring phosphate tetrahedra, with the residual oxygen atom remaining coordinated to the yttrium interstitial. The lowest energy configuration consisted of two Frenkel pairs with a defect energy of 4.50 eV, Figure 10(c). Here, yttrium interstitials are edge sharing, adopting sixcoordinate sites in adjacent channels. They are symmetrical about the centre of the cluster, at distances of 1.93 and 3.10 A from the vacancy sites. Phosphorus Frenkel pair clusters A total of 25 unique defect configurations were identified with defect energies ranging from 3.02 to 29.30 eV. The lowest energy configuration for a phosphorus Frenkel pair is illustrated in Figure 11(a). Here, the phosphorus atom is shifted within the ac plane to form a polymerised P2 O7 unit. This configuration is consider21



(a)

(b)

Page 22 of 36

(c)

Figure 10: Yttrium Frenkel pair cluster defect configurations. Residual oxygen atoms and vacancy sites are highlighted in blue and orange, respectively. ably more energetically favourable than the defect configuration illustrated in Figure 11(b), in which the phosphorous interstitial is located in the channel. Here, the interstitial is corner sharing with two phosphate tetrahedra, and is coordinated to two non-bridging oxygen atoms, forming a P3 O10 polymerised unit. The defect energy calculated for this configuration is 6.02 eV.

(a)

(b)

Figure 11: Phosphorus Frenkel pair cluster defect configurations. Residual oxygen atoms and vacancy sites are highlighted in blue and orange, respectively.

Oxygen Frenkel pair clusters In total 12 unique configurations were identified containing an oxygen Frenkel pair cluster, with defect energies ranging from 3.85 to 9.15 eV. In 22


Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


general, the lower energy configurations are produced by either rotations or distortions of phosphate tetrahedra, where the oxygen interstitial remained coordinated to the original phosphorous atom. The lowest energy configuration is illustrated in Figure 12(a) where the oxygen interstitial is removed from a vacancy site in the ac plane and relocated into the bc plane. A similar configuration is illustrated in Figure 12(b) where the oxygen interstitial is located in the channel. The defect energies for these configurations are calculated to be 3.85 and 3.96 eV, respectively. In the defect configurations where the oxygen atom was removed from the coordination sphere of a phosphate unit, the resulting three-coordinate phosphorous atom would polymerise with a neighbouring phosphate unit forming a P2 O7 unit with a vacancy on the residual oxygen site. The lowest energy of these configurations is illustrated in Figure 12(c), with a defect energy of 6.33 eV. Here, the oxygen interstitial is bridging between two nine coordinate yttrium atoms, at a distance of 3.63 A from the vacancy site. Finally, a low energy defect configuration is identified consisting of three oxygen Frenkel pairs resulting from the rotation of a phosphate tetrahedron of ca 60° about one of the P · · · O bonds, Figure 12(d). The defect energy for this configuration is calculated to be 6.01 eV. Table 9: Interstitial defect formation energies at infinite dilution Defect

Defect Formation Energy eV −5.90 −26.81 −27.41 −27.11 1.69 1.82

Y••• i P••••• 00000 i VP + 2P••••• i 000 VY + Y••• + P••••• i i O00 •• i 00 VO + 2Oi

3.4.3

Frenkel pair formation energies

Bound and unbound Frenkel defect formation energies have been calculated with the associated binding energies, listed in Table 10. Simple Frenkel defects have been considered along 23



(a)

(b)

Page 24 of 36

(c)

(d)

Figure 12: Phosphorus Frenkel pair cluster defect configurations. Vacancy sites are highlighted in orange. with some more complex Frenkel pair clusters where multiple Frenkel pairs have aggregated together. It is important to note that it is energetically favourable for the vacancy and interstitial defects to cluster together for all three of the isolated Frenkel pairs, in particular for the phosphorus Frenkel pair. This trend extends to the more complex aggregated Frenkel pairs with the phosphorous-containing defect clusters having larger binding energies. In general, the formation energies calculated for isolated Frenkel pairs at infinite dilution for yttrium, phosphorus and oxygen are in good agreement with the work of Gao et al., with the phosphorus and oxygen Frenkel pairs being the highest and lowest energy, respectively, Table 11. The binding energy associated with the aggregation of the point defects has been shown to play a significant role in lowering the formation energy of defect configurations. When this is taken into account, the phosphorus Frenkel pair cluster becomes the lowest 24


Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 10: Formation and binding energies calculated for the Frenkel defects Defect

Formation Energy eV mol−1 000 ••• V + Yi 8.81 000Y a VY + Y••• 6.53 i 00000 ••••• VP + Pi 21.44 00000 ••••• b V P + Pi 3.02 •• 00 VO + Oi 7.67 •• 00 c VO + Oi 3.85 000 ••• 2VY + 2Yi 17.62 000 ••• 2 V + Yi 13.06 000Y ••• d 2VY + 2Yi 4.50 00000 ••••• 2VP + 2Pi 42.88 00000 ••••• 2 V P + Pi 6.04 00000 ••••• 00000 20.84 VP + VP + 2Pi •• 00 2VO + 2Oi 15.34 00 + O 7.70 2 V•• O •• i 00 •• VO + VO + 2Oi 7.80 00 •• 23.01 3VO + 3Oi •• 00 3 V O + Oi 11.55 •• •• •• 00 00 VO + VO + 2Oi + VO + Oi 11.65 •• 00 e 6.01 3VO + 3Oi 00000 000 ••• ••••• V + VY + Yi + Pi 30.25 000 00000P ••• ••••• 9.55 + VY + Yi V P + Pi 000 00000 ••• ••••• V + VY + Yi + Pi 21.14 P00000 ••• ••••• f 000 4.93 VP + VY + Yi + Pi 00000 •• ••••• 00 V + VO + P i +O 29.11 00000P •• i 00 ••••• V + Pi + VO + Oi 6.87 P00000 ••••• 00 g VP + V•• + P + O 6.33 i i O a Figure 10(a),b Figure 11(a),c Figure 12(a),d Figure 10(c) e Figure 12(d),f Figure 10(b),g Figure 12(c),

Binding Energy eV mol−1 − −2.28 − −18.42 − −3.82 − −4.56 −13.12 − −36.84 −22.04 − −7.64 −7.54 − −11.46 −11.36 −17.00 − −20.70 −9.11 −25.32 − −22.24 −22.78

energy of the three species, with a formation energy of 3.02 eV. This defect configuration results from a small displacement of the phosphorus atom in the ac or bc plane, where it becomes polymerised with a neighbouring phosphate tetrahedron to form a P2 O7 unit. Similarly, the binding energy for the defect cluster consisting of two localised yttrium Frenkel 25



Page 26 of 36

••• pairs, 2V000 Y + 2Yi , is sufficiently large to reduce the formation energy of the cluster to less

than that of a single yttrium Frenkel pair cluster. In general, all phosphorus Frenkel pairs resulted in the polymerisation of multiple phosphate units, frequently forming P2 O7 units preferentially. Oxygen Frenkel pairs also often introduced polymerisation of local phosphate units when the oxygen interstitial is located outside of the nearest neighbour distance of the phosphorus atom to which it was originally coordinated. This residual three-coordinate phosphorus atom subsequently polymerises with a neighbouring phosphate tetrahedron to form a P2 O7 unit. The lowest energy oxygen Frenkel pair clusters were, however, found to be those resulting from simple rotations or distortions of phosphate tetrahedra. With respect to oxygen split interstitials, very little difference is found between the formation energies of the 00 configuration where two oxygen Frenkel pair clusters are at infinite dilution, 2 {V•• O + Oi }, •• 00 and for the split interstitial, V•• O + {VO + 2Oi }, with formation energies of 7.70 and 7.80 eV,

respectively. It is noteworthy that when comparing the trend in formation energies for the three possible Frenkel pairs in xenotime and zircon, Table 11, the Frenkel pair for the B cation at infinite dilution is unanimously found to be the least energetically favourable in xenotime. In contrast, the Frenkel pair for the A cation has the highest formation energy in zircon. 3.4.4

Anti-site Defects

Anti-site defects are simple to define in xenotime with the yttrium and phosphorous ions exchanging crystallographic sites following the Kröger Vink notation in Equation 17

•• × * 00 YY + P× P ) YP + P Y

(17)

The lowest energy anti-site cluster was identified during the search for phosphorus Frenkel pairs via the grid search method, Figure 13. Here, the yttrium and phosphorus defects sit slightly offset from the corresponding lattice sites at a distance of ca 1.10 A. This results in a six-coordinate yttrium anti-site and a four coordinate phosphorus anti-site, polymerising 26


Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 11: Comparison of Frenkel pair formation energies

Formation Energy / eV Defect

This This Gao Urusov Pruneda Crocoma 31 a 33 b 30 Work Work et al. et al. et al. betteb 26 (∞) (Cluster) AFP 8.81 6.53 7.37 9.53 13.20 24.00 BFP 21.44 3.02 16.86 40.25 10.70 22.90 OFP 7.67 3.85 4.43 10.79 7.60 7.30 General formula for orthophosphates and silicates is ABO4 . a YPO and b ZrSiO 4 4 N.B. reported defect energies are not calculated from models using a consistent forcefield. Therefore absolute values should not be compared directly.

with two neighbouring phosphate tetrahedra to form a P3 O10 unit. The formation energies for the bound and unbound anti-site defects are listed in Table 12, along with the associated binding energy. The formation energy for the anti-site defect at infinite dilution is much higher than that of the anti-site defect in zircon, 8.5 eV, calculated by Pruneda et al. 30 The large binding energy associated with the cluster reduces the formation energy to a more comparable value.

Figure 13: Anti-site defect cluster found using the grid search methodology. Vacancy sites are highlighted in orange.

27



Page 28 of 36

Table 12: Formation and binding energies calculated for the anti-site defects

Anti-site Defect

Formation Energy / eV

Binding Energy / eV

Y00 + P•• Y 00P YP + P•• Y

21.93 7.60

− −14.33

4

Conclusion

A new self-consistent force field has been developed across the YPO4 , Y2 O3 and P2 O5 phases using a fitting methodology based upon a reverse Monte Carlo procedure. These potentials successfully reproduce the crystallographic and mechanical properties of these phases, in relatively good agreement with experimental data. This new force field has been used to calculate the formation energies for a range of intrinsic defects in xenotime including Schottky, Schottky-like, Frenkel pair and anti-site defects. The binding energy associated with the aggregation of point defects into clusters is shown to play a significant role in lowering the formation energy of defect configurations. Though the Frenkel pair defects at infinite dilution are generally in good agreement with the previous work of Gao et al., this study has identified that, once the binding energy is taken into account, the phosphorus Frenkel pair cluster is the lowest energy Frenkel pair, forming a P2 O7 polymerised unit in the ac or bc plane. The polymerisation of phosphate units was found to be a key feature of the defect configurations for both phosphorus and oxygen Frenkel pairs, with the lower energy configurations generally forming P2 O7 units and maintaining a coordination number of four for the phosphorus atoms. In contrast, the lowest energy yttrium Frenkel pair configurations generally resulted in an yttrium interstitial with a reduced coordination number of six. Intrinsic differences have been identified between the lowest energy defect configuration in this study for xenotime and the equivalent defects in zircon outlined by Crocombette and Pruneda et al. In zircon, an oxygen split interstitial, resulting in a five coordinate silicon 28


Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atom, was identified as the lowest energy oxygen interstitial. The equivalent split interstitial calculated in this study was much higher in energy than an oxygen interstitial located in the c channel, bridging between two yttrium atoms. However, a similar lower energy split interstitial was identified for xenotime where the phosphorus atom remained in a four coordinate geometry. This may be indicative of the relative stability of the five coordinate silicon and phosphorus species. In general the phosphorus Frenkel pair at infinite dilution was found to be the least energetically favourable cation Frenkel pair in this study and those of Gao et al. 31 and Urusov et al. 33 , whereas in contrast, the zirconium Frenkel pair was found to have the highest formation energy in zircon by both Crocombette 26 and Punedaet al. 30 This again indicates another inherent difference in the intrinsic defect properties of the two isostructural phases. Finally, the formation energy for anti-site defects at infinite dilution are calculated to be much larger than those calculated for zircon by Pruneda et al., likely owing to the difference in the relative charges of the Zr, Si, Y and P atoms, respectively. The transferability of these potentials across the YPO4 , Y2 O3 and P2 O5 phases along with the good agreement between the observed and calculated crystallographic and mechanical properties makes them ideal for simulating radiation damage events in the xenotime system within a molecular dynamics regime. N.B. employing the fitting method described above, further parameterisation of charges and short-range pair potentials for systems such as UO2 and PuO2 would enable a consistent force field to fully model the effect of incorporating such aliovalent dopant species as extrinsic defects into xenotime. An additional point to consider when simulating radiation damage with molecular dynamic techniques is the possibility that energy would not be conserved during the simulation if an O-P-O bond breaks, rendering the three body term inapplicable. Although not of concern to the static defect simulations reported in this paper, a solution would be to implement a variant of the three-body potential that decays smoothly to zero as the P-O distance increases, to conserve energy. We have shown in this work that the novel potential set for xenotime, derived empirically using in-house reverse Monte Carlo techniques, is robust for the investigation of the defect 29



Page 30 of 36

chemistry. Our results further the understanding of radiation damage in this material, which is essential for long term assessment of the safety and surety of this material for use in the immobilisation of High Level Radioactive Waste.

5

Acknowledgements

This work was funded by the School of Chemistry at the University of Birmingham. The calculations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service, which provides High Performance Computing to the University’s research community 1 . GLC thanks the School of Chemistry for the provision of a studentship.

References (1) Abbott, D. Limits to growth: Can nuclear power supply the world’s needs? Bull. At. Sci. 2012, 68, 23–32. (2) Radioactive

wastes

in

the

UK:

A

summary

of

the

2016

inventory.

https://ukinventory.nda.gov.uk/wp-content/uploads/sites/18/2017/03/ High-Level-Summary-UK-Radwaste-Inventory-2016.pdf. (3) Kleykamp, H.; Paschoal, J. O.; Pejsa, R.; Thummler, F. Composition and Structure of Fission-Product Precipitates in Irradiated Oxide Fuels - Correlation With Phase Studies In The Mo-Ru-Rh-Pd and BaO-UO2-ZrO2-MoO2 Systems. J. Nucl. Mater. 1985, 130, 426–433. (4) Sari, C.; Schumacher, G.; Walker, C. T. Solubility and Migration of Fission-Product Barium in Oxide Fuel. J. Nucl. Mater. 1979, 79, 255–259. 1

http://www.birmingham.ac.uk/bear

30


Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(5) Romberger, K. A.; Baes, C. F.; Stone, H. H. Phase Equilibrium Studies in UO2-ZrO2 System. J. Inorg. Nucl. Chem. 1967, 29, 1619–1630. (6) Markin, T. L.; Street, R. S.; Crouch, E. C. Uranium-Cerium-Oxygen Ternary Phase Diagram. J. Inorg. Nucl. Chem. 1970, 32, 59–&. (7) Heyn, F. A.; Aten, A. H. W.; Bakker, C. J. Transmutation of uranium and thorium by neutrons. Nature 1939, 143, 516–517. (8) Donald, I. W.; Metcalfe, B. L.; Taylor, R. N. J. The immobilization of high level radioactive wastes using ceramics and glasses. J. Mater. Sci. 1997, 32, 5851–5887. (9) Weber, W. J.; Ewing, R. C.; Angell, C. A.; Arnold, G. W.; Cormack, A. N.; Delaye, J. M.; Griscom, D. L.; Hobbs, L. W.; Navrotsky, A.; Price, D. L. et al. Radiation effects in glasses used for immobilization of high-level waste and plutonium disposition. J. Mater. Res. 1997, 12, 1946–1978. (10) Pupin, J. P. Zircon And Granite Petrology. Contrib. Mineral Petr. 1980, 73, 207–220. (11) Nasdala, L.; Wenzel, M.; Vavra, G.; Irmer, G.; Wenzel, T.; Kober, B. Metamictisation of natural zircon: accumulation versus thermal annealing of radioactivity-induced damage. Contrib. Mineral Petr. 2001, 141, 125–144. (12) Murakami, T.; Chakoumakos, B. C.; Ewing, R. C.; Lumpkin, G. R.; Weber, W. J. Alpha-Decay Event Damage In Zircon. Am. Mineral. 1991, 76, 1510–1532. (13) Rubatto, D.; Muentener, O.; Barnhoorn, A.; Gregory, C. Dissolution-reprecipitation of zircon at low-temperature, high-pressure conditions (Lanzo Massif, Italy). Am. Mineral. 2008, 93, 1519–1529. (14) Mezger, K.; Krogstad, E. J. Interpretation of discordant U-Pb zircon ages: An evaluation. J. Metamorph. Geol. 1997, 15, 127–140.

31



(15) Koksal, S.; Goncuoglu, M. C.; Toksoy-Koksal, F.; Moller, A.; Kemnitz, H. Zircon typologies and internal structures as petrogenetic indicators in contrasting granitoid types from central Anatolia, Turkey. Mineral. Petrol. 2008, 93, 185–211. (16) Geisler, T.; Pidgeon, R. T.; van Bronswijk, W.; Kurtz, R. Transport of uranium, thorium, and lead in metamict zircon under low-temperature hydrothermal conditions. Chem. Geol. 2002, 191, 141–154, 11th Meeting of the European Union of Geosciences (EUG), STRASBOURG, FRANCE, APR 12, 2001. (17) Meldrum, A.; Boatner, L. A.; Ewing, R. C. Displacive radiation effects in the monaziteand zircon-structure orthophosphates. Phys. Rev. B 1997, 56, 13805–13814. (18) Ni, Y. X.; Hughes, J. M.; Mariano, A. N. Crystal-Chemistry Of The Monazite And Xenotime Structures. Am. Mineral. 1995, 80, 21–26. (19) Ewing, R. C.; Meldrum, A.; Wang, L. M.; Wang, S. X. Transformation Processes in Minerals; Reviews in Mineralogy & Geochemistry; 2000; Vol. 39; pp 319–361. (20) Meldrum, A.; Boatner, L. A.; Weber, W. J.; Ewing, R. C. Radiation damage in zircon and Monazite. (vol 62, pg 2509, 1998). Geochim. Cosmochim. Acta. 1999, 63, 2693. (21) Meldrum, A.; Zinkle, S. J.; Boatner, L. A.; Ewing, R. C. Heavy-ion irradiation effects in the ABO(4) orthosilicates: Decomposition, amorphization, and recrystallization. Phys. Rev. B 1999, 59, 3981–3992. (22) Zhang, M.; Salje, E. K. H.; Farnan, I.; Graeme-Barber, A.; Daniel, P.; Ewing, R. C.; Clark, A. M.; Leroux, H. Metamictization of zircon: Raman spectroscopic study. J. Phys.: Condens. Matter 2000, 12, 1915–1925. (23) Farnan, I.; Salje, E. K. H. The degree and nature of radiation damage in zircon observed by Si-29 nuclear magnetic resonance. J. Appl. Phys. 2001, 89, 2084–2090.

32


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Ashbrook, S. E.; Farnan, I. Solid-state O-17 nuclear magnetic resonance spectroscopy without isotopic enrichment: direct detection of bridging oxygen in radiation damaged zircon. Solid State Nucl. Mag. 2004, 26, 105–112. (25) Farnan, I.; Cho, H.; Weber, W. J. Quantification of actinide alpha-radiation damage in minerals and ceramics. Nature 2007, 445, 190–193. (26) Crocombette, J. P. Theoretical study of point defects in crystalline zircon. Phys. Chem. Miner. 1999, 27, 138–143. (27) Williford, R. E.; Begg, B. D.; Weber, W. J.; Hess, N. J. Computer simulation of Pu3+ and Pu4+ substitutions in zircon. J. Nucl. Mater. 2000, 278, 207–211. (28) Moreira, P. A. F. P.; Devanathan, R.; Yu, J.; Weber, W. J. Molecular-dynamics simulation of threshold displacement energies in zircon. Nucl. Instrum. Methods Phys. Res., Sect. B 2009, 267, 3431–3436. (29) Grechanovsky, A. E. Computer Modelling of Autoradiation Damages to ZrSiO4 Zircon, YPO4 Xenotime, La2Zr2O7 Pyrochlore and CaZrTi2O7 Zirconolite. Metallofiz. Nov. Tekh. 2011, 33, 287–296. (30) Pruneda, J. M.; Artacho, E. Energetics of intrinsic point defects in ZrSiO4. Phys. Rev. B 2005, 71, 094113–1 – 094113–7. (31) Gao, F.; Xiao, H. Y.; Zhou, Y. G.; Devanathan, R.; Hu, S. Y.; Li, Y. L.; Sun, X.; Khaleel, M. A. Ab initio study of defect properties in YPO4. Comput. Mater. Sci. 2012, 54, 170–175. (32) Rabone, J. A. L.; De Leeuw, N. H. Interatomic potential models for natural apatite crystals: Incorporating strontium and the lanthanides. J. Comput. Chem. 2006, 27, 253–266.

33



(33) Urusov, V. S.; Grechanovsky, A. E.; Eremin, N. N. Radiation Resistance of the Xenotime YPO4 from the Computer Simulation Data. Glass Phys. Chem. 2012, 38, 55–62. (34) Kaminskii, A. A.; Bettinelli, M.; Speghini, A.; Rhee, H.; Eichler, H. J.; Mariotto, G. Tetragonal YPO4 - a novel SRS-active crystal. Laser Phys. Lett. 2008, 5, 367–374. (35) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simulat. 2003, 29, 291–341. (36) Mott, N. F.; Littleton, M. J. Conduction in Polar Crystals .1. Electrolytic Conduction in Solid Salts (Reprinted From, Vol 34, PG 485, 1938). J. Chem. Soc., Faraday Trans. II 1989, 85, 565–579. (37) Mouhat, F.; Coudert, F.-X. Necessary and sufficient elastic stability conditions in various crystal systems. Phys. Rev. B 2014, 90, 224104–1 – 224104–4. (38) Mogilevsky, P.; Zaretsky, E. B.; Parthasarathy, T. A.; Meisenkothen, F. Composition, lattice parameters, and room temperature elastic constants of natural single crystal xenotime from Novo Horizonte. Phys. Chem. Miner. 2006, 33, 691–698. (39) Lacomba-Perales, R.; Errandonea, D.; Meng, Y.; Bettinelli, M. High-pressure stability and compressibility of APO(4) (A = La, Nd, Eu, Gd, Er, and Y) orthophosphates: An x-ray diffraction study using synchrotron radiation. Phys. Rev. B 2010, 81, 064113–1 – 064113–9. (40) Paton, M. G.; Maslen, E. N. A Refinement of Crystal Structure of Yttria. Appl. Catal. 1965, 19, 307–310. (41) McClure, J. P. PhD Thesis. Ph.D. Thesis, University of Nevada 2009, (42) Arbib, E. H.; Elouadi, B.; Chaminade, J. P.; Darriet, J. New refinement of the crystal structure of o-P2O5. J. Solid State Chem. 1996, 127, 350–353.

34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(43) Ziegler, J. F.; Ziegler, M. D.; Biersack, J. P. . Nucl. Instrum. Methods Phys. Res., Sect. B 2010, 268, 1818–1823.

35



Graphical TOC Entry

36


Page 36 of 36

Derivation of Transferable Pair Potentials and the Calculation of

Recommend Documents