{'fraction_non_alphanumeric': 0.025010734220695576, 'fraction_numerical': 0.011753971661657365, 'mean_word_length': 4.78484942564421, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ab Initio Construction of Interatomic Potentials for Uranium Dioxide Across all Interatomic Distances

We provide a methodology for generating interatomic potentials for use in classical molecular dynamics simulations of atomistic phenomena occurring at energy scales ranging from lattice vibrations to crystal defects to high energy collisions. A rigorous method to objectively determine the shape of an interatomic potential over all length scales is introduced by building upon a charged-ion generalization of the well-known Ziegler-Biersack-Littmark universal potential that provides the short- and long-range limiting behavior of the potential. At intermediate ranges the potential is smoothly adjusted by fitting to ab initio data. Our formalism provides a complete description of the interatomic potentials that can be used at any energy scale, and thus, eliminates the inherent ambiguity of splining different potentials generated to study different kinds of atomic materials behavior. We exemplify the method by developing rigid-ion potentials for uranium dioxide interactions under conditions ranging from thermodynamic equilibrium to very high atomic energy collisions relevant for fission events.

Molecular Dynamics (MD) simulations provide a convenient tool for studying the dynamics of large atomic ensembles, provided that the dynamics of interest is not of a duration that makes the simulation time impractical. There are many examples of beautiful applications of how MD has been able to provide detailed insight into the dynamics and statistics of materials behavior. One contemporary interest is the field of high energy radiation damage of crystalline materials, such as nuclear fuel. A core material of interest in this context is uranium dioxide, and one of the important aspects of the interest in this material is to understand the evolution and statistics of atomic displacement cascades due to high energy radiation 1 . Classical Molecular Dynamics is ideally suited for this kind of study since it strikes a fine balance between being coarse enough to simulate the spatial scale necessary to represent the extent of a damage cascade due to, e.g., a 100 MeV atomic collision with being detailed enough to retain the atomic structure of the material. The high energy range of the potential (short range in atomic separation) is consistent with the well-accepted Ziegler-Biersack-Littmark (ZBL) universal pair potentials 2 , which treat close range atomic interactions as screened Coulomb forces between the nuclei. However, the complexity of the true interatomic interactions cannot be fully represented in an efficient manner by a simple classical functional form. Thus, one needs to develop a set of essential interaction features that are necessary for a given application. This is a particularly challenging exercise for radiation damage simulations due to the disparate scales of energies involved. Therefore, interatomic potentials, suitable for this purpose, are typically constructed by smoothly joining different types of interactions. At medium to long range distances, a traditional potential (e.g., Buckingham, electrostatic etc.) fitted to a variety of thermodynamic and structural data is used. At short-ranges, accurate potentials are developed by fitting to ab initio data. The ZBL universal potential is one such very popular pair-potential developed by Ziegler et al. in the 1980's, as a generic function of the atomic numbers of the species involved. 2 Although each type of interaction is directly determined through fitting, the determination of a suitable spline that smoothly joins these two pieces is a highly non-unique process. Splining leads to an inherent ambiguity in the behavior of the complete potential, since the exact cutoff distances and the spline's algebraic form will have consequences for how large the "cores" of the atomic interactions are and how the potential behaves in the region of transition. This, in turn, will have a significant impact on the ion trajectory and damage production one is ultimately interested in. Such ambiguity is illustrated in Figure 1. Of course, a complete description of the dynamics of the system will also require the inclusion of a suitable electronic stopping model 2,3,4,5,6,7 .
In the present communication we introduce a rigorous method to unambiguously determine the form of the potential over all distances using ab initio data. Although our approach is applicable to any material system, we illustrate it using the important example of the nuclear fuel UO 2 , which has been the focus of several detailed computational investigations, both through ab initio 8 and MD 9,10,11 methods, due to its critical importance in the nuclear industry.
Our approach is to first generalize the universal ZBL potential to include charged ions that behaves correctly in both short range and long range limits. The advantage of building upon the ZBL formalism is that our potential automatically inherits the well-tested ability of the ZBL potential to describe high-energy scattering phenomena associated with the short-range behavior of the potential. This generalized potential smoothly interpolates between these two regimes over a physically-motivated length scale that is based on atomic orbital sizes instead of necessitating a user-specified transition radius for the electrostatic interactions at short ranges to prevent double-counting. The only component that remains to be determined by fitting is then the medium-range energy contribution associated with chemical bonding. As this contribution is only significant over a relatively small range of distances it is possible to introduce lower and upper cutoffs, where this contribution must smoothly vanish. Importantly, these physically motivated cutoffs can be fitted to ab initio data and are therefore no longer arbitrary; unlike in the current practice.
The ZBL potential 2 , which properly accounts for the screening of nuclear charge by the electronic clouds as a function of interatomic distances, is built on considering two interacting spherically symmetric rigid electron clouds. In this spirit, we also consider two interacting spherically symmetric rigid electron clouds with electron densities determined ab initio and fitted to a sum of Slater functions. 12 This approximation is valid beyond distances where electron clouds overlap and chemical bonds form. For short distances (i.e., distances less than chemical bond lengths) we obtain the energy as a function of distance between any two atoms through first order perturbation theory. It was demonstrated by Ziegler et al. that more sophisticated self-consistent field calculations incorporating the distortion of electronic clouds did not lead to any significant differences in the resulting interatomic potential at short distances. 2 Thus, since we employ the same electronic density and the same screening function (ratio of the actual atomic potential at some radius to the potential caused by an unscreened nucleus) as used by Ziegler et al., we recover the ZBL potential at these short distances, as we will show later.
We consider two spherically symmetric charge densities per unit volume ρ 1 (r) and ρ 2 (r), with central point charges of Z 1 e and Z 2 e respectively, ρ 1 (r) and ρ 2 (r) being normalized to equal (Z 1 + q 1 )e and (Z 2 + q 2 )e respectively. We further think of the point charge Z 1 e as being made of two point charges, (Z 1 + q 1 )e and (−q 1 )e; Z 2 e is similarly decomposed into (Z 2 + q 2 )e and (−q 2 )e. q 1 e and q 2 e here denote the net ionic charges on atoms 1 and 2 respectively. We make this decomposition so that the Coulombic interaction term naturally arises in the expression for the net interaction potential V (r), where ZBL Z1+q1,Z2+q2 (r) denotes the ZBL form of interaction between two neutral atoms having atomic numbers Z 1 + q 1 and Z 2 + q 2 , but using the screening length for Z 1 and Z 2 . This is because the screening length is governed by the electronic charge distribution close to the nucleus and not far away from it. Since the extra charges q 1 and q 2 have been added to the valence shells of neutral atoms with atomic numbers Z 1 and Z 2 , we do not change the screening length. This was found critical in recovering the standard short-range ZBL potential. t 1 denotes the interaction between the point charge −q 1 e and the (Z 2 + q 2 )e point charge plus electron cloud system of the atom 2, and is given by while t 2 is the converse remaining point-atom interaction, expressed similarly. Eq.(1) is correct for atomic separations smaller than the 'chemical bond length' (where it recovers the original neutral atom ZBL) and for very large atomic separations as well. In the latter case, only the 2nd term in Eq.(1) survives as we show below.
The task now is the determination of ρ O (r) and ρ U (r). It is here important to notice that the potential for very small interatomic separations is only as good as the ZBL form (see Eq. (1)). Thus, we use the charge densities employed for ZBL, which are primarily Hartree-Fock-Slater atomic distributions for most of the atomic pairs. We fit the numerical data for charge density used by Ziegler et al. to a sum of Slater functions. While the density ρ(r) is known to be a monotonic decreasing function of radial distance for all atoms, the graph of 4πr 2 ρ(r) exhibits a number of peaks (see Figure 2) corresponding to atomic shells. To ensure the best possible accuracy, we fit to 4πr 2 ρ(r) because this is the quantity entering Eq.(1). We find that 2 and 4 Slater functions are sufficient to capture the behavior of 4πr 2 ρ(r) for neutral Oxygen and Uranium atoms (see Figure 2), i.e., Fitting the ZBL charge density with Slater functions does not exactly ensure that the areas under the 4πr 2 ρ(r) curves for Oxygen and Uranium are 8 and 92 respectively. The Slater function fits therefore needed a slight modification in their pre-factors. In addition, the pre-factors in the above two equations as reported in Table I have been multiplied by 10/8 for Oxygen and 88/92 for Uranium because we are interested in the electronic cloud of the ionic species and not of the neutral atoms. We also experimented with more sophisticated corrections (e.g., using noble gas densities instead) but this did not change the results by more than the intrinsic accuracy of the ZBL potential. By using Eqs.
(3) and (4) in Eq.(1) and performing the integrations, we obtain the following pair potentials for Oxygen-Oxygen, Uranium-Uranium and Oxygen-Uranium respectively: where we have We illustrate in Figure 3 how closely the potentials given in Eqs. (5) - (7) match, for small r, the neutral atom ZBL, and for large r, the relevant Coulombic interaction.
As for any empirical potential, there is an intermediate distance range for which the interactions follow neither a ZBL nor a purely Coulombic form. A correction term is thus needed for this regime. We find this correction term by fitting to an extensive database of GGA+U ab initio calculations on UO 2 . GGA+U is known to provide electronic and magnetic behaviors of UO 2 that are consistent with experiments 13 , and a correct treatment of the localized and strongly correlated 5f electrons of Uranium. 14 Our ab initio calculations also take into account the experimentally observed noncollinear antiferromagnetic magnetic moment ordering and the Oxygen cage distortion in UO 2 . 15 Therefore, in addition to capturing correct elastic and defect properties, our potential also covers a much more vast energy landscape in the material due to the richness of the ab initio data used. We fit to a database obtained by GGA+U calculations with the projector-augmented-wave method implemented in the VASP 16 package. In the GGA+U approximation, the spin-polarized GGA potential is supplemented by a Hubbard-like term to account for the strongly correlated 5f orbitals. 17 We use the rotationally invariant approach to GGA+U due to Dudarev et al. 18 , wherein the parameter U-J is set to 3.99 eV. This is the generally accepted value for this parameter to reproduce the correct band structure for UO 2 . 19 The magnetic moments were allowed to be fully non-collinear. The ab initio database so obtained comprises: (i) isochoric relaxed runs on a 12 atom unit cell which was isometrically contracted and expanded by various amounts (i.e., equation of state calculations wherein each data point was calculated under constraint of constant cell volume), and for which an energy cutoff of 500 eV and a 8×8×8 k-point grid were taken; k-point convergence was ascertained before choosing this value for the k-point grid. The cell was allowed to relax in shape but not in size. Ionic relaxations were carried out until residual forces less than 0.01 eV/Åwere achieved.
(ii) static (i.e., no ionic relaxation) runs on 96 atom 2×2×2 supercell in which one atom at a time (i.e., Oxygen or Uranium) was perturbed from its equilibrium position by varying distances in different directions. Energy cutoff was 500 eV. After performing convergence studies on the k-point grid, gamma point only version of VASP was found to be satisfactorily accurate for this. Note that any interactions between atoms and their periodic images do not systematically bias the fit of the potential because the same supercell geometry is used in both the ab initio and the empirical potential energy calculations.
(iii) 96 atom 2×2×2 supercell for the formation energies of three kinds of stoichiometric defects, namely Oxygen Frenkel pair, Uranium Frenkel pair and Schottky trio. The vacancies and the interstitials were taken as far from each other as the supercell would allow. The details of the calculations are the same as that for case (ii) above. Correct prediction of defect energies has been given great importance in generating interatomic potentials for cascade simulations.
(iv) first-order transition states in a 2×2×2 supercell for the migration energy of Oxygen and Uranium vacancy and interstitial. Nudged Elastic Band(NEB) method 20 in conjunction with the climbing image method 21 for determination of saddle point energy, as implemented in VASP, was used for this.
With the ab initio database so generated, we now fit the final potential forms as follows.
The long-range Coulomb terms in the interaction here were calculated through the standard Ewald summation technique. The upper cut-offs for all terms except the Coulombic in Eq. (10) may be chosen as per the availability of computational resources. As seen from Eq. (10) we now have absolutely no splines for the U +4 -U +4 interaction, reflecting the fact that no chemical bonding takes place. There are splines in the other two interactions but these are now unambiguously determined since the respective cut-offs are not imposed but instead determined through fitting. The splines maintain continuity through the second derivatives of the potential and the specific form of the splines can be uniquely recovered from these conditions. Since they have been fitted to accurate ab initio data, these splines do not introduce any spurious wriggles in the potential. For the interaction between two Oxygen ions, the potential has one (and only one) minimum at r min = 2.28Å, as may be seen from Figure 4. We have thus minimized the unphysical features of all of the interatomic potentials in UO 2 , which were demonstrated for the particular case of Oxygen-Oxygen interaction in Figure 1.
The downhill simplex method of Nelder-Mead was used to carry out the fitting. 24 The fitting involved minimizing the sum of the squares of the differences between the ab initio energies and the energies predicted by the potential for all the classes of data points as detailed above. The package GULP 25 was used for energy calculations and for atomic positions optimization. The quality of fit for the equation of state data and the perturbed atom data can be seen in Figure 5. The ab initio/experimental and predicted defect formation/migration energies are compared in Table II,  while Table III lists the predicted ground state lattice parameter and other elastic properties as compared with the corresponding experimental values 10,26 (extrapolated accordingly) and with values obtained with the Morelon et al. potential 23 . The agreement is very satisfactory.
As a final validation of the developed potential we considered various dynamic properties by performing MD simulations in a constant number, pressure and temperature (NPT) ensemble comprising 6×6×6 unit cells. The system was equilibrated for 10.0 ps, while production runs were carried out for 100.0 ps, with a time step of 0.001 ps and sampling every 0.05 ps. The fluorite structure remained stable during all the runs we performed, up to temperatures of 2500 K. The properties we considered are the variation in the lattice parameter and the enthalpy as functions of the temperature. These are compared in Figure 6 with the corresponding experimental data. 26 The quality is similar to what is given by the previous potentials 23,27,28 , as tabulated in the work by Morelon et al. 23 To summarize, we have shown a methodology for developing an interatomic pair potential such that it is appropriate for all relevant interatomic separations, without the need for any ambiguous splines. Splining between regions of different characteristics is not just an inconvenience in the implementation of a potential in MD simulations, but also introduces an uncertainty regarding which distances, and by which functions, one realizes the splines. The potential we have focused on in this presentation, the nuclear fuel UO 2 , is just one example of a model material in which very relevant materials physics depends on accurate and reliable interactions over many orders of magnitude, and we have obtained the first complete description that allows for direct simulations of damage cascades due to high energy radiation effects. The potential has been generated based on a slight revision of the ZBL universal potential to account for ionic materials with the intermediate interatomic distances fitted to a broad database of ab initio structural energies. In view of these qualities, we expect it to be a very reliable potential for studying displacement cascades in UO 2 .      26 The scatter in the experimental values is also shown (solid lines). (b) Enthalpy variation using potential from current work (open circles) compared with corresponding experimental values (solid line). 26 The scatter in the experimental values was less than 1 percent and is thus not shown.