Determination of Vacancy Formation Energies in Binary UZr Alloys Using Special Quasirandom Structure Methods

The use of predictive models to examine defect production and migration in metallic systems requires a thorough understanding of the energetics of defect formation and migration. In fully miscible alloys, atomistic properties will all have a range of values that are heavily dependent on local atomic configurations. In this work we have used the atomistic simulation tool Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) to investigate the impact of first nearest neighbor configuration on vacancy formation energies at 0 K in γ-U-Zr alloys of varying Zr concentrations. The properties of randomly generated alloy microstructures were also compared with those produced as special quasi-random structures (SQS) using the “mcsqs” code within the Alloy Theoretic Automated Toolkit. Results have confirmed that local configuration can have a significant impact on measured properties and must be considered when characterizing miscible alloy systems. Results also indicated that the generation method of the random structure (i.e., via random species assignment or a method of enforced randomness) does not result in a measurable difference in average vacancy formation energies in miscible U-Zr systems.


INTRODUCTION
The use of predictive models to examine defect production and migration in metallic systems requires a thorough understanding of the energetics of defect formation and migration. In fully miscible alloys, atomistic properties are heavily dependent on local atomic configurations, and can vary significantly across the volume of even a homogeneous alloy mixture. One such atomistic property is the vacancy formation energy (E v ), which is a very important material property that is frequently required for continuum and Monte Carlo simulations of materials (Soisson, 2006;Badillo et al., 2015;Piochaud et al., 2016). It plays a significant role in governing the diffusion and formation of point defects in materials, and must be parameterized accurately for the specific material system to allow for representative simulations. In alloy systems, such as the U-Zr alloy system examined in this work, the vacancy formation energy is a function of local atomic arrangement and must be tracked separately for the two species, since continuum methods such as phase field or Monte Carlo methods must track which species is removed to form the defect. The E v for elemental uranium has been explored extensively in both experiments and computational calculations, with a mono-vacancy formation energy of 1.20 ± 0.25 eV being reported by Matter et al. (1980). E v in Zr metal has not been studied as extensively, with different authors reporting values over a large range of energies (1.68 eV by Korzhavyi et al. (1999) and 2.1 eV by Le Bacq et al. (1999). Additional reported values for E v in elemental uranium and zirconium are provided in Tables 1, 2 in section Results and Discussion.
For this work, the Modified Embedded Atom Method (MEAM) based interatomic potential for U-Zr developed by Moore et al. (2015) was used. This potential was developed to simulate c-U-Zr alloys at operational temperatures. The uranium portion of this potential was designed to recreate the properties of c-phase uranium at high temperatures, but was also fit to the elastic properties and lattice constant of c-phase uranium at 0 K. However, the c-phase of pure uranium is unstable at 0 K, and must be stabilized through the application of an external pressure. Due to this mechanical instability, calculations of properties by other common simulation methods such as density functional theory (DFT) or ab initio calculations cannot be performed, due to their reliance on the presence of an equilibrium state at 0 K (Moore et al., 2015). However, through the use of a volume or pressure constraint in molecular dynamics simulations, simulations of pure c-phase uranium can be performed. It is also important to note that the presence of a small atom percentage of an alloying element can stabilize the c-phase of uranium, meaning that all of the studied U-Zr alloy compositions in this work were stable without requiring additional constraints.
It is common practice in molecular dynamics (MD) modeling of metallic alloys to populate the lattice structure via random atom placement. Bulk material properties of crystals generated in this way can closely mimic the properties of real crystals; however, due to the finite nature of any modeled supercell, the random population of lattice sites may deviate from "perfect" randomness on a local scale. To address this issue, Wei et al. (1990) proposed a method of generating "special quasirandom structures" (SQS) that would more closely mimic the configuration average of an infinite perfectly random A 1−x B x alloy, and can be generated in sizes as small as N 8. SQS "supercells" can be generated, and then combined into a larger lattices if bulk properties needed to be tested. The process of generating SQS is provided in Wei et al. (1990), where multisite correlation functions are used to determine if a randomly generated microstructure mimics the ensemble averages that would be expected for an infinite random alloy. This is done to ensure that there are not regions within the randomly generated microstructure that statistically would be considered not random, thereby ensuring that properties determined using SQS are representative of a homogeneously mixed alloy. In their original work, Wei et al. (1990) found that a supercell size of N 8 was sufficient to mimic the first and second nearest neighbor correlation functions for a binary alloy, and N 64 was sufficient to mimic the correlation functions for the third nearest neighbors. For this work, SQS's of 2000 atoms and 250 atoms were developed using the stochastic SQS generation system that was implemented in the "mcsqs" code by van de Walle et al. (2013). The 250 atom structures were then repeated to produce a single system containing 2,000 atoms.
For any alloy, there exists a distribution of possible vacancy formation energies due to the species arrangement and thermal oscillations of the atoms in the neighboring sites (Zhao et al., 2016). For this work, to examine the impact of different configurations of nearest neighbor sites on vacancy formation energy, all 256 possible nearest neighbor configurations for a binary bcc alloy were examined for each atom percentage and vacancy type tested. However, not all configurations are equally likely to occur in all alloys with differing atom percentages (e.g., all eight nearest neighbor lattice positions of a vacancy in a U-10% Zr alloys being filled by Zr atoms, etc.) In a perfectly random binary alloy, the probability of an atom of species A inhabiting any given lattice position can be defined by the concentration of that species in the bulk lattice, C A , while the probability of species B inhabiting the same lattice position can be defined as the concentration of species B in the lattice, C B . For a number of lattice positions z, the probability that i lattice positions are inhabited by species A and z − i lattice positions are inhabited by species B can be defined as This probability refers to the probability of i species A atoms and z − i species B atoms inhabiting specific lattice positions. For each species distribution in z lattice positions, with i atoms of species A and z − i atoms of species B, the total number of configurations with that species distribution can be described by the following expression: Finally, the expression that defines the probability of having a specific species distribution in z lattice positions in a binary alloy can be expressed as the product of the probability of having the

Vizoso and Deo
Vacancies in U-Zr Alloys species distribution in a particular configuration and the total number of possible configurations for that species distribution, resulting in a binomial probability distribution (Kim, 1984), In a BCC structure, there are eight first nearest neighbor sites and 6 s nearest neighbor sites. For a binary BCC alloy, there then exist 256 unique first nearest neighbor configurations and 104 unique second nearest neighbor configurations, combining for a total of 26,624 unique local configurations around a single atomic position, which in a binary system can be one of two different species. To account for this the results of all possible configurations were compared with results for the configurations with atom percentages that were reasonable when compared to the atom percentage of the alloy as a whole. All MD simulations were performed with a supercell consisting of 10 × 10 × 10 unit cells held in an NVE ensemble.

Interatomic Potential of Uranium-Zirconium
Molecular dynamics simulations model interactions between atoms via interatomic potential functions that are used to calculate the force on any atom in the system caused by the other N−1 atoms in the system. A primary advantage of interatomic potential functions over first principle calculations is that potentials can be used for the calculation of larger scale atomic properties at temperature, as well as the ability to perform calculations on systems with many atoms. This is particularly advantageous for materials such as c-uranium, which is mechanically unstable at 0 K (Moore et al., 2015). For this work, a modified embedded atom method, or MEAM, potential was used. MEAM potentials are modified versions of embedded atom method potentials (EAM), which are semi-empirical potentials that express cohesive energy in terms of embedding energies, with the atoms in the system being embedded in an electron gas. The difference between MEAM and EAM potentials is that MEAM potentials account for angular forces, which are present in certain materials that form directional bonds (Beeler et al., 2012). Full descriptions and comparisons of EAM and MEAM potentials have been completed in other works (Beeler et al., 2012;Foiles and Baskes, 2012), and only a short description of the theory will be included here. A more complete description of MEAM potentials is also provided in the Supplementary Material, paraphrased from Beeler et al. (2012). The general equation for the interatomic potential for both EAM and MEAM is shown by Eq. 4.
In this equation, E represents the energy of the atoms i, ϕ ij represents the pair potential between atoms i and j, ρ is the effective electron density, and F i is the embedding function, which introduces many-body effects. For EAM potentials, the background electron density is taken to be the linear superposition of the spherically averaged atomic electron densities (ρ 0 ), whereas in MEAM, angular effects are considered (Beeler et al., 2012). Vacancy formation energies can be calculated for a bulk lattice using the following equation, where E (n−1) is the energy of the system with a single vacancy, E (n) is the energy of the perfect lattice, and n is the total number of atoms in the perfect lattice (Beeler et al., 2012). This equation was used due to the emphasis of this work on the impact of local atomic configuration and alloy composition on the value of the vacancy formation energy with respect to irradiation, rather than desiring a global minimum value for use in determining the concentration of thermal vacancies Zhang and Sluiter (2015), Ruban (2016).

Special Quasirandom Structures Generation Procedures
For the generation of SQS, three input files were used with the "mcsqs" code from the Alloy Theoretic Automated Toolkit (ATAT) (van de Walle and Asta, 2002;van de Walle et al., 2002;Hart and Forcade, 2008;van de Walle et al., 2013;Liu et al., 2017). The first input file (by default rndstr.in, but this can be changed with the −1 filename option in the command line) defines the structure of the alloy and the ratio of the components. This input file must be made by the user. The second input file (by default clusters.out, but this can be changed with the −cf filename option in the command line) tells the code what clusters must be considered while generating the SQS. This input file can be generated by code by using the command where −2 and −3 define the ranges of pairs and triplets, while in the same directory as rndstr.in. For the purpose of this work, the range of pairs was set to 6 Angstroms, and the range of triplets was set to 8 Angstroms. The third input file used was sqscell. out, which defines the supercell that "mcsqs" will use while generating the SQS. This input file is required only if the user requires a supercell of a specific shape. For the purpose of this work, a cubic supercell was used, to allow for ease of comparison between the SQS structure and the cubic structures that LAMMPS generates. To generate the SQS, the command mcsqs − n . . . − rc must be used, where −n defines the number of atoms that will be placed within the supercell, and −rc tells the code to use the supercell defined in sqscell.out. For the code to work, the value for −n must be an integer, and the number of atoms must match with the size of the supercell provided in sqscell.out.
To run multiple ATAT sequences in parallel, the command −ip ... can be added to the previous command, with an integer following the equals sign. This integer will be appended to the names of the output files produced by the code. The produced SQS will be found in the output file bestsqs(i). out, which will have the dimensions of the unit cell and the supercell, and then the xyz location and type of each atom in the supercell in units of the specified lattice parameter. The file bestcorr.out tracks the correlations of the best SQS that has been found. For the purpose of this work, the code was stopped after the first SQS was generated.

Random Structures
All of the random structures used for this work were produced by Plimpton (1995) LAMMPS. To produce a random structure, the dimensions of the supercell must be given in terms of unit cells, and the lattice parameter of the unit cell must be provided. Also, the structure of the unit cell must be specified (bcc, fcc, etc. . .). LAMMPS will produce a homogeneous structure with the specified dimensions and structure. To create a random alloy, the set command is used, which allows the user to specify what percentage of atoms should be replaced by a specified atom type. Each alloy was minimized to energy (unitless) and force (eV/Angstrom) stopping tolerances of 1.0e−15. After running the initial minimization on the as-generated crystal, the atom at the center of the supercell was removed, and the alloy was again run for one time-step, then minimized to the same tolerance limits. For each random supercell generated by LAMMPS, the calculation of the vacancy formation energy was performed for both zirconium and uranium vacancies by setting the atom type of the atom that was to be removed to the desired element. To prevent the mechanical instability of pure uranium at 0 K from influencing the results of simulations performed with pure uranium, the energy minimization step was only performed on the first nearest neighbors and second nearest neighbors of the vacancy.

Special Quasirandom Structures
LAMMPS has the ability to import supercell structures through the use of the read_data command. After the supercell structure is read into the LAMMPS, the same process for the determination of the vacancy formation energy was used. For the 2,000 atom SQS structures, a MATLAB script was used that determined the atom type of the atom at the center of the supercell, and found the xyz coordinates of the closest atom to the center of the supercell that was of the other type. For the 250 SQS structures, the atom at the center of the supercell was removed, and random atoms within the structure were chosen until an atom of the different type was found. Then, two runs were performed for each SQS, one where the atom at the center of the supercell was removed, and then a run where the atom of the opposite type was removed. In this way, the structure of the SQS was preserved.

RESULTS AND DISCUSSION
Using LAMMPS, vacancy formation energies were calculated for the random structures and the SQS by minimizing the energy of the lattice as it was produced, then removing an atom from the center of the structure, minimizing the energy, and then taking the difference between the energies before and after the removal of the atom, as described by Eq. 5 in the Section Methods. All of these simulations were performed at a temperature of 0 K. Figures 1, 2 compare the vacancy formation energies for the removal of zirconium and uranium atoms, respectively, for structures that were generated using the SQS method and by random atom placement. SQS can not be generated for pure metals, so only a single LAMMPS simulation was performed to determine the E v of uranium and zirconium in pure uranium and zirconium metals. Tabulated results for the E v in the pure metals can be seen in Tables 1, 2, compared with results from other works. The value of 1.34 eV obtained for the uranium E v in pure uranium at 0 K agrees well with the value of 1.38 eV produced by Beeler et al. using a MEAM potential for c-uranium in LAMMPS (Beeler et al., 2013).
For the E v of zirconium, the value seems to be fairly constant with changing atom percent of uranium, with any variation in the lines being within two standard deviations. The value of 2.074 eV obtained for the zirconium E v in pure zirconium matches very well with the results of 2.07 and 2.10 eV attained by Willaime and Massobrio (1991), Le Bacq et al. (1999), and Moore et al. (2015), respectively.
The results for the vacancy formation energies of the 2000 atom SQS structures and the random structures were in good agreement. The E v for the 250 atom SQS structures tended to be higher than those from the random and the 2,000 atom SQS structures at higher uranium percentages; however, this difference never exceeded more than one standard deviations, indicating that it was not a true difference in the properties between the structures. One important note is that random alloys produced by LAMMPS do not have exact atom percentages, with the atom percentages for each run being randomly distributed around the specified percent, with a maximum deviation of a few percentage points. Atom percentages for SQS structures are exact.
To examine the impact of local configuration on the vacancy formation energy, all possible nearest neighbor configurations were tested for the U-Zr alloy system for both uranium and zirconium vacancies. It has been shown in previous studies that the vacancy formation energy in binary alloys is highly dependent on the species distribution in both the first and second nearest neighbor sites (del Rio et al., 2011). For the simulations performed for this work only the first nearest neighbor configurations were varied, while keeping the second nearest neighbor configurations as they were generated. Therefore, the impact of second nearest neighbor configurations on the Frontiers in Materials | www.frontiersin.org July 2021 | Volume 8 | Article 692660 vacancy formation energy cannot be determined from the collected data.
It is important to note that not all local atomic configurations are equally likely in an alloy with a set concentration. For example based on the probabilities provided in Table 3, a U-10%Zr alloy will have a 96.2% chance that there will be three or fewer zirconium atoms in nearest neighbor sites, and a 99.55% chance that there will be four or fewer zirconium atoms in nearest neighbor positions. Therefore, when considering the impact that different configurations have on the average vacancy formation energy for a given material composition, it is valid to sample configurations based on their probability of occuring. Figures 3, 4 show the vacancy formation energy as a function of local atomic composition and configuration calculated for a zirconium and uranium vacancy in U-10%Zr alloy, respectively. FIGURE 1 | Plots of the vacancy formation energies for the removal of zirconium calculated for SQS and random structures of U-Zr alloy ranging from 0% to 100% U.
FIGURE 2 | Plots of the vacancy formation energies for the removal of uranium calculated for SQS and random structures of U-Zr alloy ranging from 0% to 100% U.
Frontiers in Materials | www.frontiersin.org July 2021 | Volume 8 | Article 692660 5 Each data point represents the E v for a particular nearest neighbor configuration with the number of Zr atoms specified along the x-axis. The plots of the E v results for the 256 configurations at different atom percentages can be seen in Figures 5, 6. All of the simulations for Figures 5, 6 were performed in randomly generated supercells with the nearest neighbors adjusted manually.
By examining Figures 3, 4, it can be seen that the energy required to remove either a uranium or zirconium atom from a U-10%Zr alloy generally increases as the total number of Zr atoms in nearest neighbor positions increases. In the case of removal of Zr, there is a slight drop in E v as the number of zirconium atoms in nearest neighbor sites increases from five atoms to eight atoms, while the trend in the removal of uranium consistently increases with increasing zirconium atoms in nearest neighbor positions. This trend could be interpreted as an increase in the vacancy formation energy near positions where the local atomic concentrations deviate significantly from the composition of the alloy. This exact trend is observed in Figure 5, where the U-10%Zr, U-20%Zr, U-80%Zr, and U-90%Zr alloy compositions all show the trend of increasing E v as the local composition differs significantly from that of the alloy. For the middle concentrations (U-30% Zr to U-70%Zr), very similar trends are observed where there is a slight decreasing trend in E v with an increasing number of zirconium nearest neighbors. In the case of the dependence of E v on local atomic composition for the removal of uranium atoms, the same trend is observed for all alloy compositions, where the energy required to form a vacancy by removing a uranium atom always increases as the number of zirconium nearest neighbors increases.
Determining specific trends in how the local arrangement of the species (i.e., species clumped together or distributed evenly across the nearest neighbor sites) impacts the formation energy of vacancies is difficult, particularly since this will also depend on the composition and arrangement of the FIGURE 3 | Plot of the vacancy formation energy for removing a zirconium atom from a U-10% Zr alloy. Each data point represents the vacancy formation energy for a particular configuration for the number of Zr atoms in the nearest neighbor sites as specified by the x-axis position.
FIGURE 4 | Plot of the vacancy formation energy for removing a uranium atom from a U-10% Zr alloy. Each data point represents the vacancy formation energy for a particular configuration for the number of Zr atoms in the nearest neighbor sites as specified by the x-axis position.
Frontiers in Materials | www.frontiersin.org July 2021 | Volume 8 | Article 692660 second nearest neighbor atoms, which were not manipulated in this study, as well as the alloy composition. Examination of specific configurations by hand was performed for a series of initial alloy compositions and configurations, and it was determined that E v is higher when the different species are grouped together into adjacent nearest neighbor sites, and that E v decreases as the species are spread out as much as possible.

CONCLUSIONS
The formation energies of uranium and zirconium vacancies in U-Zr binary bcc alloys were determined using MD simulations with a MEAM potential. Non-linear dependence on atom percent was observed for both uranium and zirconium vacancies. Good agreement was found for the E v of both pure uranium and pure zirconium with the works of Matter et al. (1980), Beeler et al. (2010), and Beeler et al. (2013) and for uranium and Willaime and Massobrio (1991), Le Bacq et al. (1999), and Moore et al. (2015) for zirconium. The impact of crystal generation method on E v was examined by generating the alloys both as random alloys and as SQS. It was found that, for the supercell size used in our MD simulations, there is no statistically significant difference between the vacancy formation energies in the crystals generated by the two methods. This indicates that the supercell size used in our simulations was sufficiently large to mimic a truly random structure. Additional work could be performed to determine the supercell size where a random alloy ceases to perform similarly to SQS structures, which are intended to mimic perfectly random systems.
The effect of local atomic configuration on E v was examined by varying the configuration of atoms in the nearest neighbor sites around a vacancy in U-Zr alloys. It was observed that for both uranium and zirconium, the configuration of nearest neighbor atoms did have an impact on the E v . The removal energy of uranium was found to increase with increasing Zr nearest neighbors up until an atom fraction of 90% Zr, where no strong dependence on nearest neighbor configuration was found. The removal energy of zirconium was found to increase with increasing Zr nearest neighbors in low atom percent Zr alloys, and decrease with increasing Zr nearest neighbors in alloys with Zr atom percentages higher than percent. It was observed that E v increased when the nearest neighbors were arranged in a manner that grouped atoms of the same species together, while E v was found to decrease as the species were evenly distributed across the nearest neighbor sites.
With a better understanding of the E v in c-U-Zr alloys, additional work can be performed with the examination of defect diffusion via vacancies, as well as the formation of defects due to radiation damage events. This work could also support the development or improvement of conventional interatomic potentials for alloys, or as a possible data source for machine learning interatomic potentials or machine learnining-based simulation techniques. Further insight into the relationship between the E v and local atomic configuration could be explored through the use of DFT simulations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CD conceptualized and initiated the ideas in this work, DV performed calculations and made modifications to initial FIGURE 6 | Configuration E v of uranium vacancies for different atom percentages.
Frontiers in Materials | www.frontiersin.org July 2021 | Volume 8 | Article 692660 8 ideas based on results and consultation with CD. DV wrote the initial draft and CD edited and revised the paper draft.

FUNDING
CD would like to acknowledge funding from the G.W. Woodruff Fellowship. CD and DV acknowledge support from Department of Energy Nuclear Energy UNiveristy Progranms Office.