Features of structure and phase transitions in pure uranium and U–Mo alloys: atomistic simulation

We study structural properties of cubic and tetragonal phases of U-Mo alloys using atomistic simulations: molecular dynamics and density functional theory. For pure uranium and U-Mo alloys at low temperatures we observe body-centered tetragonal (bct) structure, which is similar to the metastable γ°-phase found in the experiments. At higher temperatures bct structure transforms to a quasi body-centered cubic (q-bcc) phase that exhibits cubic symmetry just on the scale of several interatomic spacings or when averaged over time. Instantaneous pair distribution function (PDF) differs from PDF for the time-averaged atomic coordinates corresponding to the bcc lattice. The local positions of uranium atoms in q-bcc lattice correspond to the bct structure, which is energetically favourable due to formation of short U-U bonds. Transition from bct to q-bcc could be considered as ferro-to paraelastic transition of order-disorder type. The temperature of transition depends on Mo concentration. For pure uranium it is equal to about 700 K, which is well below than the upper boundary of the stability region of the α-U phase. Due to this reason, bct phase is observed only in uranium alloys containing metals with low solubility in α-U.


Introduction
Pure uranium has three different allotropes: a low-temperature orthorhombic α-U, a hightemperature body-centered cubic (bcc) γ-U, and an intermediate tetragonal β-U observed in a very small interval of pressures and temperatures.
To improve mechanical properties, irradiation and corrosion resistance of uranium at room temperature while maintaining the high density, uranium is frequently alloyed with other metals. For example, molybdenum exhibits high solubility in γ-U and allows to stabilize γ-like phase up to room temperature. Solid solution based on γ-U alloyed with Mo shows stable swelling behavior in U-Mo fuels, also it has high thermal conductivity, low thermal expansion, and high melting point [1][2][3].
Ab initio methods have been applied to investigate the structure and thermodynamics of various phases of pure U [13][14][15][16] and U alloys [2,17]. These calculations revealed the energetic hierarchy of the uranium structures and predicted the existence of the body-centered tetragonal (bct) structure that is more stable than the bcc lattice [13,16,[18][19][20]. The energy difference between bcc and bct lattices is approximately 0.1 eV/atom. Ab initio calculations [14][15][16] of the elastic constant C ′ = C 11 − C 12 have shown that C ′ is negative for bcc phase, and this structure is mechanically unstable with respect to the tetragonal distortion.
However, to our knowledge, bct phase have not been registered experimentally in pure uranium at ambient pressure. Conversely, experimental studies of the U-Mo alloys observe a tetragonally distorted bcc phase [4][5][6][7][8]10]. Atomistic simulations [21] can be used to shed some light on the nature of high-temperature bcc and bct phases of U and U-Mo alloys.
At the same time, molecular dynamics (MD) simulations with different interatomic potentials predict dissimilar structures for these phases [20,22,23]. So one should be sure that the model applied agrees with the experiments and ab initio simulations. However, static instability of the bcc lattice does not allow to apply ab initio simulations; thus, a lot of properties usually extracted from the static calculations become unavailable.
Parameters of the potential [23] were found via the force-matching method [24,25] based on the ab initio calculations. So this model accurately describes the fitting data computed in the framework of the density functional theory (DFT). Nevertheless, accuracy of description of bcc γ-U by the DFT with currently available pseudopotentials is still debated [26].
In this work we analyse structural properties of cubic and tetragonally distorted phases by means of atomistic simulations. Also we pay a special attention to comparison with the existing experiments.

Computational technique
Molecular dynamics is a powerful tool for study of physical properties of matter at the atomistic level. In this work we consider structure of U-Mo alloys using molecular dynamics simulation with the novel interatomic ADP potential [23]. This potential contains an angular-dependent term and was parameterized using the force-matching method [24,25]. For pure uranium we additionally made tests with EAM potential [20] and MEAM potential [22]. All calculations were performed with the LAMMPS code [27]. Dynamics of atoms was visualized by the Atom-Eye program [28], where it was necessary.
DFT calculations have been performed with the code VASP 5.2 [29]. Typical structures contain approximately 54-128 atoms in a simulation box with periodic boundary conditions. The Brillouin zone is sampled with the 2×2×2 Monkhorst-Pack k-point mesh for static relaxation, while calculations with just Γ-point are used for DFT-MD runs, since relatively long MD trajectories (up to 100 ps) are required. The cut-off energy of a plane-wave basis set is 330 eV. We use projector augmented wave pseudopotentials included in VASP package and the exchangecorrelation functional within generalized-gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof. Fourteen electrons 6s 2 6p 6 5f 3 6d 1 7s 2 for uranium are taken into account as valence electrons.

Simulation of uranium β-phase
According to the experimental studies β-U is stable in a temperature range between 945 K and 1045 K, at pressure lower than 3 GPa [30,31]. This uranium phase has a complicated structure compared to the other uranium allotropes: it was determined to be tetragonal with 30 atoms in the unit cell [32]. At the same time, attempts of identification of the exact space group corresponding to such complex atomic structure led to the controversy. It was suggested that the atomic arrangement can be possibly described by three different space groups: P4 2 Figure 1. Illustration of the atomic arrangement in the model of β-U described by space group No. 102, together with the corresponding G(r) computed with the ADP potential at 900 K and zero pressure.
During the last decades experimental and computational efforts have been made to determine the most possible space group. For example, in 1988, analyses of the time-of-flight neutron diffraction data reported by Lawson et al [33] led them to conclusion that the symmetry of β-U should be described by space group No. 136.
Ab initio simulations of β-phase structures performed by Li et al [34] at zero temperature also resume that application of space group No. 136 gives the most credible results. It is also important to note that one of the conclusions made in [34] concerns the fact that the structures described by space groups No. 136 and No. 118 can be identical for β-U. Structure and elastic constants of the P4 2 /mnm (No. 136) have been than additionally tested in [16]. Nevertheless, the first-principles computations [35] of energy per atom and the phonon properties revealed that the β-U model defined by space group No. 102 turned out to be more stable (energetically and dynamically) than the model corresponding to space group No. 136. Summarizing information reported about the β-U it can be said that to date there is no unequivocal opinion about its structure.
In this work we perform MD simulations of β-U at finite temperatures. For this purpose we use two interatomic potentials designed previously for investigation of pure uranium (EAM potential [20]) and its alloys (ADP potential [23]). Originally, any information about β-U have not been used during the fitting of these potentials, therefore, it can be a challenge for them to provide a model of this complex phase. Two different space groups suggested for β-U are considered: one is P4 2 /mnm (No. 136), and the other is P4 2 nm (No. 102). Figure 1 illustrates atomic arrangement (projected on a plane) and the pair distribution function G(r) corresponding to the actual MD simulation of β-U described by space group No. 102.
A series of simulations of β-U was carried out for different temperatures (in a range between 230 K and 950 K) and zero pressure. MD simulation box contained 7500 atoms with periodic boundary conditions along x, y, and z directions. The number of atoms, volume of the simulation box and the total energy were kept constant during the computation run. Time-length of each simulation was set equal to 350 ps.
We see that the ADP potential predicts existence of the stable structures of P4 2 /mnm and P4 2 nm β-U almost in all range of the studied temperatures. The results also allow to conclude that the area of stability seems to be larger for the case of P4 2 nm structure-mechanical instabilities are observed only at T > 940 K, while P4 2 /mnm starts to become unstable at T > 900 K. Thereby, the temperature of the upper stability limit in the MD model of β-U is about 100 K lower than the experimental temperature of β-U ←→ γ-U transition. The G(r) computed for the more stable β-U is plotted in figure 1. We should note that it is quite similar to the G(r) obtained for the q-bcc structure (see next section). It is important to compare the equilibrium lattice parameters of the simulated β-U with the existing results of the experimental studies. Lattice constants of pure β-U measured at high temperatures (955 K and 1030 K) have been reported previously in [33]. Also one can take into account the data obtained for β-U stabilized at the room temperature by a small addition of Cr [36]. From the experiments [36] it was found that a = 10.590Å and c = 5.634Å for T = 300 K, while the MD simulations with the ADP potential gives a = 10.594Å and c = 5.583Å at the same T . As for the high temperatures, experimental constants are equal to a = 10.7589Å and c = 5.6531Å for T = 950 K; MD simulations at the upper temperature limit provides a = 10.706Å and c = 5.642Å for T = 940 K. According to the comparison given above it is possible to conclude that the lattice constants of the simulated β-U are close to that obtained from the different experimental investigations. Here we should also note that the ADP potential predicts the same temperature dependence of the lattice constants for both structures of β-U studied in this work.
At the same time, the simulation attempts made using the EAM potential do not lead to successful results. The EAM potential does not provide stable structure of β-U in the studied temperature range. As opposed to the results reported for the ADP potential, both structures (No. 102 and No. 136) are found to be unstable, they deteriorate completely over a time ∼ 1 ps.

Structure of pure uranium at high temperature
In order to study the structure of γ-U we used pair distribution functions G(r), which can be easily calculated from the molecular dynamics. Contrary to the simulations, extracting PDF from the experiment is not convenient. However, there are methods for it based on the x-ray powder diffraction and on the extended x-ray absorption fine structure (EXAFS) spectroscopy that probes the element specific local structure [37]. Unfortunately, the experiments of such kind have not been found for uranium alloys.
Structure of crystals are succesfully revealed from the more widespread x-ray diffraction (XRD) techniques based on the analysis of the dependence of scattered x-rays intensity on scattering angle for single crystals or powders. In order to "guess" a structure from the XRD pattern one have to solve the inversed problem: it is necessary to calculate theoretical XRD pattern of the probable structure A and compare it to an experimental pattern.
Theoretical XRD pattern can be calculated in two steps. At the first step one have to calculate structure factors for a primitive cell F = N j=1 f j (θ) exp(−2iπk · r j ), intensities I = AF F * and positions θ = arcsin(λ|k|/2) of peaks at zero temperature. Here k is the change of x-rays wave vector, r j is the position of each atom in a primitive cell, f j are atomic scattering factors, θ is the scattering angle of diffraction and λ is the wavelength of x-rays. At the second step, calculated peaks must be multiplied by Debye-Waller factors DWF = exp(−1/3 u 2 |k| 2 ), where u 2 is the average square of atom displacement (the latter is invalid for correlated motion of atoms).
We have calculated XRD patterns from our MD simulations with the help of the additional package [38] available in LAMMPS. Calculating XRD patterns from the MD snapshots provides different technique to take into account thermal motion of atoms, including the effects of correlated displacements. Unfortunately, this type of calculation may be performed only with the MD simulation, as it requires large number of atoms in a simulation box.
Instantaneous pair distribution functions G(r) for γ-phase at 1000 K and zero presssure are presented in figure 2. One can see results obtained from DFT-MD calculations and from the MD with two potentials: ADP [23] and MEAM [22]. EAM potential gives the results that are similar to the ones computed with ADP potential. We also calculated pair distribution functions from the time-averaged positions for each atom R (it may be denoted like G R (r)).
For calculations with the MEAM potential, both G(r) and time-averaged G R (r) represent pair distribution of atoms in a simple bcc lattice, the only effect of time averaging is a reduction of peaks' width. One could see that in DFT-MD and in MD with ADP potential, G(r) and  Figure 2. Typical instantaneous G(r) and G R (r) from the time-averaged positions of atoms in the γ-U (at T = 1000 K) described by various models: a) ADP potential [23]; b) MEAM potential [22]; c) DFT-MD. Vertical black lines mark the maximum of G(r) corresponding to the perfect bcc lattice. Inline image (d) shows XRD patterns calculated for γ-U by ADP [23] and MEAM [22] potentials.
G R (r) are strongly different from each other, so the most probable instantaneous positions of atoms are not identical to a time-averaged arrangement. Moreover, comparison of the peaks of the first and the second coordination spheres tells us that G R (r) is much closer to bcc structure than G(r): the intensity of the first peak should be larger in bcc due to a larger number of neighbours in the first shell (8), than in the second (6). Thus, the high-temperature phase has cubic symmetry only in average over time or on the scale of several interatomic spacings. Detailed analysis of the snapshots shows that the local positions of uranium atoms in the γ-phase correspond to the structure which is similar to body-centered tetragonal lattice with displacements of the central atoms (as it will be discussed in more details in the next section). The latter could be interpreted as the formation of 4 short U-U contacts, quite similar to a structure of α-U. The length of these "bonds" is approximately 2.7Å, which is quite close to the estimates (2.87-2.91Å) obtained from the experimental XRD data [6,7].
One could note that the first peak in the instantaneous G(r) is more pronounced for DFT-MD calculations, than for the classical MD. Moreover, averaging to an ideal bcc structure proceeds slower in DFT-MD: G R (r) curve presented in figure 2a (ADP potential) is obtained for averaging interval of 0.5 ps, while in figure 2c (DFT) it is averaged over 10 ps. So, we can conclude that correlations in atomic motion (connected with the displacements of atoms from the base centers of the cubes) in a high-temperature γ-U are more pronounced in DFT-MD.
For DFT-MD and ADP potential it is more correct to denote the high-temperature state of pure uranium as q-bcc (or quasi-bcc). Such designation allows to distinguish the properties of such phase from classical bcc, which is reproduced by MEAM-potential. There is no difference between q-bcc and bcc structures after time-averaging of atom positions during several picoseconds. This fact means that it is hard to distinguish these structures from the x-ray diffraction. Indeed, calculated XRD patterns (figure 2d) of structure for both potentials are practically the same. The calculated XRD patterns are identical to classical bcc structure (for instance, see experimental works [10,39]). This fact explains why q-bcc structure have not been revealed in the experiments for pure uranium.

Base centered tetragonal phase in pure uranium and U-Mo
Attempt to create γ-phase (in bcc or q-bcc structure) of U-Mo alloy with random distribution of Mo leads to formation of bct, which is similar to a metastable γ 0 -phase [21]. The martensitic phase transformation γ-phase → bct-phase takes place at low T and at low concentration of Mo. Besides, structure of high-temperature phase reflects features of bct structure. This circumstance is responsible for existence of the q-bcc instead of bcc structure in pure uranium at high temperature. Pair distribution functions G(r) are used to investigate structure details. Figure 3 shows G(r) for four cases: two U-U distributions for U-7.5 at.% Mo at T = 200 K (bct-phase) and at T = 700 K (γ-phase); U-U and Mo-Mo distributions for γ-phase of U-40at.%Mo at T = 200 K. Analysis of G(r) together with the visualization of atomic positions allow us to explore the alloy structures. The lattice parameters of bct-phase fulfill the following condition: a = b > c. Central atom in the basic cell is displaced, the displacement is parallel to c-parameter and equal to about 0.55Å (from center of the basic cell). This direction (c-axis) may be named as "anisotropy direction". In this case, the local surrounding of U atom (in region of two coordination spheres) is similar to α-phase of uranium. Such bct structure can be obtained from bcc state with minimal atomic reorganization and has lower energy, compared with bcc structure. It should be noticed that the existence of molybdenum in the system prevents such displacement. Moreover, pairs of Mo-Mo atoms are mostly "stabilization centers" of bcc lattice. Such pairs of atoms, especially located along the [111] direction, prevent the displacement of the atomic planes.
The simulated bct phase is not completely identical to γ 0 -phase. It is generally considered that the "anisotropy directions" in γ 0 -phase are oriented in opposite directions in neighboring basic cells [6,8] (i.e. the displacements have ordering which is similar to the ordering of magnetic moments in an antiferromagnetic material). In this case γ 0 -phase may be named by antiferrodisplacement phase. Contrariwise, the simulated bct phase is ferro-displacement structure with uniform orientation of the displacements of basic cell atoms. The ADP potential may describe both ordering types of "anisotropy directions" in bct phase. However, our model describes that antiferro-displacement type of bct structure is less stable than ferro-displacement structure (the difference between energies of these phases is about 0.005 eV/atom). Despite this circumstance the simulated bct phase has features, which are similar to the properties of γ 0 -phase (tetragonal symmetry, stability area on T -x diagram, mechanical hysteresis). The question about similarity of these phases requires further investigation and will be considered in future works.
The dependence of averaged bct lattice parameters on Mo concentration is shown in figure 4. The figure contains the simulation results together with the experimental data [10,40]. The data are given for room temperature. The simulation results qualitatively agree with the experiments. The calculated lattice parameters correspond to all simulation cell without specification of local atomic arrangement. For high Mo concentration the "a = b = c" condition is correct on average, but local anisotropy for uranium subsystem also takes place.
Primarily, the anisotropy is produced by displacement of central atom in the basic cell. Magnitude of this displacement is mostly independent from such conditions as temperature or composition. Anisotropy disappears at high concentrations of molybdenum, and alloy   corresponds to a bcc lattice on the average. However, the local arrangement of the uranium atoms continues to form the tetragonal phase with the central atom displacement in the entire concentration range (calculations were carried out up to 70 atomic percent of molybdenum). One can say that at high concentrations of molybdenum long-range correlations between the  Figure 5. Calculated lines of phase stability on T -x diagram for various structures: 1-α-phase turns into γ-phase at heating; 2-ferro-displacement bct phase turns into γ-phase at heating; 3-γ 0 -phase turns into ferro-displacement bct phase at heating.
atoms of uranium is missing, but the nearest neighbors in several interatomic spacings try to reproduce bct lattice. The same process takes place when the temperature increases (instead of molybdenum concentration). An analogy exists between the bct → q-bcc transition studied in this work and the ferroelectric-paraelectric and ferroelastic-paraelastic transitions [41]. At a low level of disorder (corresponding to low values of temperature or Mo concentration), the system exhibits longrange correlations in the "anisotropy direction", and formation of bct lattice takes place. At a high level of disorder, long-range correlation in the "anisotropy direction" is disturbed, and formation of γ-phase with q-bcc structure occurs. In this terminology, γ-phase is similar to the paraelastic state. Moreover, analysis of G(r) and G R (r) reveals that bct-like atomic structure preserves locally in a high-temperature γ-phase. The transformation bct → q-bcc may be described as order-disorder transition. Apparently, the transition γ 0 → γ must be described by the same way. It is quite similar to a rearrangement of Li ions in the MD simulations [42] of the ferroelectric-paraelectric transition in LiNbO 3 , coupled to order-disorder dynamics as well.
All studied structures have the areas of stability on the temperature-concentration diagram (T -x diagram). Thus, the change of temperature results in spontaneous phase transition at heating (or at cooling) in the atomic simulation for the N P T -ensemble. Such phase transformation is accompanied with the volume change. During the simulation, the following condition was set as a necessary requirement: P xx = P yy = P zz = 0. Figure 5 illustrates calculated areas of phases stability at heating. The heating rate in our simulation was about 10 11 K/s. In addition to α-phase and ferro-displacement bct phase, we calculated the area of stability for antiferro-displacement bct phase (it is common concept of γ 0 -phase). Both ferro-displacement bct phase and α-phase turn into γ-phase at heating; however, antiferrodisplacement bct phase turns into bct phase with ferro-displacement ordering.

Conclusion
The features of structure and phase transitions in pure uranium and U-Mo alloys were studied. The ADP potential reproduces properties of α-U, β-U and γ-U with good precision. For β-U the potential shows stable structure with P4 2 nm symmetry. According to the atomistic simulations, γ-U has quasi body-centered cubic structure, where local positions of atoms correspond to bodycentered tetragonal lattice with displaced central atoms. The γ-phase exhibits cubic symmetry just on the scale of several interatomic spacings or when averaged over time. This fact is the possible origin of the difficulties encountered in the description of the γ-phase of pure uranium by ab initio methods. Attempt to create γ-phase of U-Mo alloy at low temperature leads to formation of bct, which is similar to a metastable γ 0 -phase. Unfortunately, simulated bct-phase is not completely identical to γ 0 -phase. The difference between these structures is cased by various ordering of "anisotropy directions" (i.e. it is like difference between antiferromagnetic and ferromagnetic materials). γ 0 -phase may be also described by ADP potential, but this phase is less stable in comparison to simulated bct structure. The difference between energies of these phases is about 0.005 eV/atom. Despite this circumstance the simulated bct phase with ferro-displacement ordering has features which are similar to the properties of γ 0 -phase (tetragonal symmetry, displacement of central atom in basic cell, stability area on T -x diagram). The transformation bct-phase → γ-phase takes place at high T and at high concentration of Mo. This type of phase transformation is quite similar to an atoms rearrangement at ferro-to paraelastic transition of order-disorder type. Therefore, the structure of high-temperature γ-phase reflects features of bct structure (or γ 0 -phase).