**2. Materials and Methods**

Based on the formalism of EAM, several modified versions including the MEAM and ADP potentials were developed by the addition of angular dependent terms to increase the description accuracy. Similarly, for the present study, the new angular-dependent interatomic potential calculated atomic energy contribution *Ei* by the following formula:

$$E\_i = \frac{1}{2} \sum\_{j \neq i} \phi\_{ij} \left( r\_{ij} \right) + F(\overline{\rho}\_i) - \vartheta\_i \tag{1}$$

where

$$\begin{array}{c} \boldsymbol{\rho}\_{lj}(\boldsymbol{r}\_{lj}) = \sum\_{i=0}^{3} a\_{i} (\boldsymbol{r}\_{c} - \boldsymbol{r}\_{lj})^{3+l}, \\ \boldsymbol{F}(\boldsymbol{\overline{\boldsymbol{\sigma}}}\_{l}) = -\sqrt{\boldsymbol{\overline{\boldsymbol{\sigma}}}\_{l}}, \boldsymbol{\overline{\boldsymbol{\sigma}}}\_{l} = \boldsymbol{\Sigma}\_{l \neq i} \boldsymbol{\rho}\_{lj}(\boldsymbol{r}\_{lj}), \\ \boldsymbol{\rho}\_{lj}(\boldsymbol{r}\_{lj}) = \boldsymbol{a}^{2} \left(\boldsymbol{r}\_{c} - \boldsymbol{r}\_{lj}\right)^{4} \exp\left[-\boldsymbol{\beta} \left(\frac{\boldsymbol{r}\_{ij}}{r\_{0}} - \boldsymbol{1}\right)^{2}\right], \ \boldsymbol{\vartheta}\_{l} = \boldsymbol{\theta}\_{l}^{\rm u} + \boldsymbol{\theta}\_{l}^{\rm v} + \boldsymbol{\theta}\_{l}^{\rm w}. \end{array} \tag{2}$$

It could be indicated by the first two terms in Equation (1) that the main part of the atomic energy maintains the form of the widely used EAM potential. To be exact, *φij* and

*F*(*ρ<sup>i</sup>* ) describe the atomic interactions of the electrostatic repulsion between ion cores and the attraction induced by embedding atoms in the electron atmosphere provided by their neighbor atoms, respectively. In calculation of *Ei* of the *i*-th atom, summation is over its *j*-th neighbor atom within the cutoff distance *r*c, which was set to 4.7 Å in the present study. In particular, *r*<sup>0</sup> and *rij* denotes the first neighbor distance in a reference lattice at equilibrium and the distance between the *i*-th center atom and its *j*-th neighbor atom, respectively. Due to the inherent sphere symmetry in the potential form, EAM poorly describes some peculiar atomic behaviors in systems of U alloys, such as the stability of the alpha phase with low symmetry at room temperature. In the present work, the description of the atomic energy was modified with the addition of a new angular-dependent term *ϑi*, which has three components given as the following:

$$\begin{split} \theta\_{l}^{\text{au}} &= \sum\_{\mathfrak{a}} \left( \sum\_{j \neq i} \Psi\_{j}^{\text{au}} \left( r\_{ij} \right) \frac{r\_{ij}^{\text{a}}}{r\_{ij}} \right)^{2}, \\ \theta\_{l}^{\text{v}} &= \sum\_{\mathfrak{a}, \mathfrak{d}} \left( \sum\_{j \neq i} \Psi\_{j}^{\text{v}} \left( r\_{lj} \right) \frac{r\_{il}^{\text{a}}}{r\_{lj}^{\text{v}}} \frac{r\_{lj}^{\text{b}}}{r\_{lj}} \right)^{2} - \frac{1}{3} \left( \sum\_{j \neq i} \Psi\_{j}^{\text{v}} \left( r\_{lj} \right) \right)^{2}, \\ \theta\_{l}^{\text{av}} &= \sum\_{\mathfrak{a}, \mathfrak{d}, \mathfrak{e}} \left( \sum\_{j \neq i} \Psi\_{j}^{\text{v}} \left( r\_{lj} \right) \frac{r\_{il}^{\text{b}}}{r\_{lj}^{\text{v}}} \frac{r\_{lj}^{\text{b}}}{r\_{lj}} \frac{r\_{lj}^{\text{v}}}{r\_{lj}} \right)^{2} - \frac{3}{\mathfrak{e}} \sum\_{\mathfrak{a}} \left( \sum\_{j \neq i} \Psi\_{j}^{\text{v}\ell} \left( r\_{lj} \right) \frac{r\_{lj}^{\text{v}}}{r\_{lj}} \right)^{2} \end{split} \tag{3}$$

where

$$\psi\_j^t(r) = c\_t \exp\left(-\frac{\left(r - d\_t\right)^2}{2\lambda\_t^2}\right), t = \mathbf{u}\_\prime \mathbf{v}\_\prime \mathbf{w}.\tag{4}$$

In Equation (3), superscripts *α*, *β*, *γ* = 1, 2, 3 refer to the Cartesian components of position vectors. *ϑ<sup>u</sup> <sup>i</sup>* , *<sup>ϑ</sup><sup>v</sup> <sup>i</sup>* , and *<sup>ϑ</sup><sup>w</sup> <sup>i</sup>* denote the angular-dependent components of the first, second, and third order, with *ψ<sup>t</sup> <sup>j</sup>* as basis function in the form of the normal distribution function. In the present study, it was assumed that all the three angular-dependent components are required in the fitting process for the description of U and the first two ones for Mo.

A total of 15 parameters were adjusted in the fitting process, which include *ai*, *α*, *β*, *ct*, *dt*, *λ<sup>t</sup>* listed in Table 1. In particular, the angular-dependent parameters of the third order for Mo were kept to zero. During the fitting, an error function was optimized. Macroscopic properties including the cohesive energy, lattice parameters, and elastic constants of certain lattices were calculated with a set of fitting parameters and compared with reference data from the literature or calculations by the First Principle (FP) method. The framework of the fitting procedure was implemented in Python, where symbolic computation and optimization are supported by Sympy and Scipy packages. In particular, in order to improve the description of point defects, several lattices with certain point defects were included in the reference data and allowed to relax during the calculations of the formation energies.

In order to capture atomic behaviors within equilibrium distances, *φij* and *ρij rij* within the inner cutoff were further modified. Interpolations are implemented for *ρij rij* to approach a constant and for *φij* to connect smoothly with the well-known Ziegler-Biersack-Littmark (ZBL) potential [17]. In particular, special attention was paid to the intermediate repulsive range in development of the potential for further employment in cascade simulations. Traditionally speaking, a smooth interpolation is used between the pairwise part of many-body potentials in the near-equilibrium range and the ZBL potential in the short range. Moreover, in 2016, Stoller et al. [18] pointed out that no accepted standard method had been developed in the interpolation process and that how the force fields are linked can sensibly influence the results of cascade simulation. Similar to the modifications implemented by Stoller, the present study also corrected the potential with calculation results from density functional theory (DFT) as a benchmark. The DFT calculations were performed using the Vienna ab initio simulation package (VASP) [19] with projector augmented wave pseudo-potentials. A kinetic energy cutoff of 400 eV and a 4 × 4 × 4 grid with the Monkhorst−Pack method were employed. The ADP calculations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [20]. A supercell containing 3 × 3 × 3 conventional bcc cells under equilibrium

(*a*<sup>U</sup> = 3.48 Å) was chosen for calculations of energetic benchmark and then an atom in the supercell was translated in several equidistant steps toward its neighbors along the direction of <110> and <100>.


**Table 1.** Fitted parameters of ADP for U–Mo system.

For simplicity in description of energetic variance, the present study used Δ*E*<sup>t</sup> as the following:

$$
\Delta E\_{\text{f}} = E\_{\text{transulated}} - E\_{\text{perfect}} \tag{5}
$$

where *E*translated is the energy of the configuration with atomic translation and *E*perfect is the energy of the original configuration with perfect symmetry. The energetics of the supercell as a function of the ratio between the distance of the nearing dimers and the lattice constant are shown in Figure 1. It could be seen that a reasonable reproduction of energetics variance was obtained in the intermediate range, especially for U-U and U-Mo dimers. Moreover, it should be noted that in cases of dimers with different elemental components, the choice of which atom is moved toward the other produces similar but actually different energetics variances. The difference depends mainly on the distinct interactions between the translated atom and its atomic environment. Good reproduction was achieved in the present study, indicating reasonable accuracy of the refitted potential in the intermediate range. Meanwhile, sensible deviations could be seen at the outer end of the intermediate range, especially for Mo-Mo dimers, which are induced by the difficulty of resolving the energetic uplift at the outer end of the intermediate range while leaving the description by the ADP of near-equilibrium atomic interactions intact. Nonetheless, even under the deviations, the general curve shapes of energetics variance were reserved to a large extent, indicating a good description of dynamic properties during atomic interactions in the intermediate range.

The newly developed potential was then applied in MD simulations of displacement cascades in U-Mo alloys in the present study. Simulation boxes containing 60 × 60 × 60 bcc unit cells were set up under periodic boundary conditions (PBC) with up to 648,000 atoms included. Alloying elements were randomly distributed in initial lattices before simulations. A total of four proportions of alloying elements, that is, 0 wt.%, 2.5 wt.%, 5 wt.%, 7.5 wt.%, and 10 wt.%, were chosen to evaluate the effect of the content of the alloying element on displacement cascades. Given the high symmetry of the bcc lattice, a total of 15 evenly distributed directions were sampled in the fundamental orientation zone, which is shown in Figure 2. The initial kinetic energies of the primary knock-on atoms (PKA) was set to 5 keV for all the simulations. For each set of initial conditions, simulations were carried out up to eight times to reduce statistical error in the following analysis.

**Figure 1.** The energetics of the supercell as a function of the ratio between the distance of the nearing dimers and the lattice constant of bcc U. In particular, nearing atom pairs include (**a**) U–U in <110>, (**b**) U–Mo in <110>, (**c**) Mo–U in <110>, (**d**) Mo–Mo in <110>, (**e**) U–U in <100>, (**f**) U–Mo in <100>, (**g**) Mo–U in <100> and (**h**) Mo–Mo in <100>. The black and red dash lines represent calculation results from DFT and MD using ADP. In the configuration diagrams, blue and orange circles represent U and Mo atoms, respectively, and the original positions of displaced atoms are indicated by dotted circles.

**Figure 2.** Angular distribution of launch directions sampled in the fundamental orientation zone of bcc lattice.

Prior to each set of cascade simulations, the simulation block was equilibrated under 600 K and 0 GPa in an isothermal-isobaric (NPT) ensemble until the temperature and pressure reached a stable level. The initial structures for cascade simulations were chosen randomly from those under equilibrium in a series of time steps. During the cascade simulations, the simulation block was divided into three layers. In the outmost layer, the atoms within 2 Å from the boundary were maintained as static all along the simulation to avoid excessive overall displacement of the simulation block. Atoms located at a distance of 2 to 8 Å from the boundary were partitioned into the second layer. In this layer, a thermostat was applied by rescaling the velocities of the atoms to maintain a stable level of temperature and to absorb the excess kinetic energy from the cascade. The rest of atoms in the core region were simulated with the constant NVE ensemble and a variable time step in a range of 0.01 to 1 fs was used in all simulations, which depends on the maximum of atom displacements between consecutive steps.

To determine and record defect sites during simulations, some native commands in LAMMPS were invoked, which construct a voronoi tessellation dividing each atom into an exclusive cell at the start of each simulation. In particular, at equal intervals of timesteps, voronoi cells with occupancy not equal to 1 were recorded with their position coordinates, while those with occupancy equal to 0 and larger than 1 were marked as vacancy and interstitial sites, respectively. In cluster analysis, adjacent defect sites within a cutoff radius of 4.1 Å are classified into the same cluster group and the size of a cluster is defined as the total count of defect sites within the cluster. By using other commands in LAMMPS, atom displacements and element types were also recorded. The program OVITO [21] was utilized for data extraction and statistical analysis.
