Potential Model for Actinide Oxides and their Solid Solutions
We have developed a new potential model describing a range of actinide oxides which includes manybody effects to improve the description of their thermophysical properties. The model has been developed to describe the following oxides and as of v1.1 of the model, their solid solutions: CeO_{2}, ThO_{2}, UO_{2}, NpO_{2}, PuO_{2}, AmO_{2} and CmO_{2}.
 Enquiries:
Version History
Versions  Description  Relevant Papers 

v1.11  10/11/2014  Explanation of error function for short range cutoff on EAM  [1],[2] 
v1.1  08/10/2014  Model extended to allow simulation of mixed oxide (MOX) systems  [1],[2] 
v1.0  09/01/2014  Initial potential model as published in J. Phys. Condens. Matter  [1] 
 The model was originally published in:
 The extension of the model to mixed oxide solid solutions was published in:
We would be grateful if you could cite these papers if you make use, or derive, any work from the materials contained here.
Currently, this page hosts descriptions of the model suitable for the GULP and LAMMPS codes, we will be providing tabulation files for DL_POLY in the near future. Furthermore, we will be releasing the EAM tabulation python library used during the development of the model in the next few weeks. If you haven any questions please get in touch.
Overview
 Downloads
 Model Description and Parameters
 Pairwise Component
 Embedded Atom Component
 Mixed CationCation Interactions
 Model Validation
 Properties of Single Actinide Oxide Systems
 Mixed Oxide System (U,Th)O_{2}
 Using the Model in Common Simulation Codes
 Using the Model in LAMMPS
 Using the Model in GULP
Downloads
File  Description 

GULP_MO2.lib  Gulp .lib containing potential parameters for all actinide oxides. 
CeThUNpPuAmCmO.eam.alloy  LAMMPS tabulation for all actinide oxides for use with pair_style eam/alloy

GULP_Example_UO2.zip  Example GULP minimisation for UO_{2} 
GULP_UThO2_example.zip  Example Mixed Oxide GULP minimisation 
LAMMPS_example.zip  Example LAMMPS equilibration 
LAMMPS_UThO2_example.zip  Example LAMMPS equilibration U, Th mixed oxide 
Model Description and Parameters
Our approach used the embedded atom method to introduce manybody effects to the more conventional Buckingham + Morse pairwise description. This approach is advantageous as it allows the inclusion of manybody interactions in MD simulations without the tricky massless entities required by the shell model.
The actinide oxide potential set has two distinct components: pairwise and manybody. The energy, \(E_i\), of atom \(i\) with respect to all other atoms can be written as follows:
\[ E_{\textit{i}} = \color{orange} \frac{1}{2}\sum_\textit{j}\phi_{\alpha\beta}(r_{ij}) \color{magenta} G_\alpha \sqrt{\sum_\textit{j}\sigma_{\beta}(r_{\textit{ij}})} \]
Pairwise Component
The first, pairwise term describes short range cationoxygen pair interactions. This uses the conventional BuckinghamMorse potential with long range Coulombic interactions also included:
\[ \definecolor{buck}{RGB}{255,0,0} \definecolor{coulomb}{RGB}{0,0,255} \definecolor{morse}{RGB}{0,128,0} \phi_{\alpha\beta}(r_{\textit{ij}}) = \color{coulomb} \frac{q_{\alpha}q_{\beta}}{4\pi \epsilon_0 r_{\textit{ij}}} \color{black} + \color{buck} A_{\alpha\beta}\exp \left(\frac{r_{\textit{ij}}}{\rho_{\alpha\beta}}\right)\frac{C_{\alpha\beta}}{r_{\textit{ij}}^6} \color{black} + \color{morse} D_{\alpha\beta}[\exp(2\gamma_{\alpha\beta}(r_{\textit{ij}}r_0))2\exp(\gamma_{\alpha\beta}(r_{\textit{ij}}r_0))] \]
For Coulombic interactions the charges of the ions are nonformal such that \(q_\alpha=Z^{eff}_\alpha e \), however, they are proportional to their formal charges ensuring the system is charge neutral (for tetravalent cations \(Z^{eff}_\alpha\) = 2.2208 and for oxygen \(Z^{eff}_\alpha\) = –1.1104 ). \(A_{\alpha\beta}\), \(\rho_{\alpha\beta}\), \(C_{\alpha\beta}\), \(D_{\alpha\beta}\), \(\gamma_{\alpha\beta}\) and \(r_0\) are empirical parameters that describe the Buckingham and Morse potentials. For cationcation and oxygenoxygen interactions the covalent Morse term is not required (i.e. D=0). The pairwise parameters are now summarised:
Interaction  \(\phi_B (r_{ij})\)  \(\phi_M (r_{ij})\)  

\(\alpha  \beta\)  \(A_{\alpha\beta}\) / eV  \(\rho_{\alpha\beta}\) / Å  \(D_{\alpha\beta}\) / eV  \(\gamma_{\alpha\beta}\) / Å^{1}  \(r_0\) / Å  
CeCe  18600  0.2664  0.0       
ThTh  18600  0.2884  0.0       
UU  18600  0.2747  0.0       
NpNp  18600  0.2692  0.0       
PuPu  18600  0.2637  0.0       
AmAm  18600  0.2609  0.0       
CmCm  18600  0.2609  0.0       
CeO  351.341  0.3805  0.0  0.7193  1.869  2.356 
ThO  315.544  0.3959  0.0  0.6261  1.860  2.498 
UO  448.779  0.3878  0.0  0.6608  2.058  2.381 
NpO  360.436  0.3830  0.0  0.7638  1.840  2.368 
PuO  377.395  0.3793  0.0  0.7019  1.980  2.346 
AmO  364.546  0.3774  0.0  0.7467  1.955  2.325 
CmO  356.083  0.3775  0.0  0.7613  2.081  2.311 
OO  830.283  0.3529  3.8843       
Embedded Atom Component
The second term uses the Embedded Atom Method to introduce a subtle manybody perturbation to the conventional pairwise description. This manybody dependency is achieved by summing a set of pairwise functions, \(\mathrm{\sigma_\beta(r_{ij})}\), between atom \(i\) and its surrounding atoms. \(\mathrm{\sigma_\beta(r_{ij})}\) is proportional to the 8th power of the interionic separation with \(n_\beta\) being the constant of proportionality and is defined by the species of the surrounding atom \(j\):
\[\sigma_{\beta}(r_{\textit{ij}}) =\left( \frac{n_{\beta}}{r_{\textit{ij}}^8} \right)\frac{1}{2}\left(1+\mathrm{erf}\left(20(r1.5)\right) \right)\]
\(\mathrm{\Sigma\sigma_\beta(r_{ij})}\) is then passed through a nonlinear embedding function. The embedding function used means that the manybody energy perturbation is proportional to the square root of \(\mathrm{\Sigma\sigma_\beta(r_{ij})}\) with \(G_\alpha\) being the constant of proportionality and is defined by the species of atom \(i\). Due to the unrealistic energies that arise from the EAM at short distances an error function is applied at 1.5 Å that reduces the EAM contribution gradually. This ensures that there is no discontinuity near the bottom of the potential well that would arise from an abrupt cut off.
Using the EAM and pairwise parameters, the figure below shows the interaction energy as a function of interionic separation for PuO, ThO and UO interactions where the actinide is coordinated by one, two or three oxygen anions (OO interactions are omitted). This plot reaffirms both the manybody nature of the potential and that the error function works well as a means of tapering the EAM at short distances.
Species  \(G_\alpha\) / eVÅ^{1.5}  \(n_\beta\) / Å^{5} 

Ce  0.308  1556.803 
Th  1.185  1742.622 
U  1.806  3450.995 
Np  0.343  1796.945 
Pu  1.231  1456.773 
Am  0.333  1631.091 
Cm  0.494  1503.704 
O  0.690  106.856 
Mixed CationCation Interactions
Although the OO interactions have been kept consistent across the actinide oxide systems, it is still necessary to develop pair potentials that describe the interaction between two dissimilar cations (e.g. a UTh potential). Similarly to the cationcation interactions defined previously, the Buckingham A term is fixed at 18600 eV for all pairs of cations, whilst C parameters are set to zero. However, the mixed cation \(\rho_{\alpha\beta}\) parameters are defined in terms of the self cationcation \(\rho_{\alpha\alpha}\) and \(\rho_{\beta\beta}\) terms using the following relationship:
\[\rho_{\alpha\beta}=(\rho_{\alpha\alpha}\rho_{\beta\beta})^{\frac{1}{2}}\]
\(\rho_{\alpha\beta}\)  Ce  Th  U  Np  Pu  Am  Cm 

Ce  0.2664  0.2772  0.2705  0.2678  0.2651  0.2637  0.2637 
Th    0.2884  0.2815  0.2786  0.2758  0.2743  0.2742 
U      0.2747  0.2719  0.2691  0.2677  0.2677 
Np        0.2692  0.2664  0.2650  0.2650 
Pu          0.2637  0.2623  0.2623 
Am            0.2609  0.2609 
Cm              0.2609 
Model Validation
The following section compares the model’s predictions with experimental data.
Properties of Single Actinide Oxide Systems
The potential was successfully fitted for all actinide oxides to experimental data of thermal expansion via MD and elastic constants via static calculations. The agreement obtained between model predictions and experimental data is now given.
Thermal Expansion:
Comparison of experimental data (solid lines) with the data points from MD simulations obtained using the potential model for CeO_{2}, ThO_{2}, UO_{2}, NpO_{2}, PuO_{2}, AmO_{2} and CmO_{2}.
Bulk Modulus:
Comparison of experimental data for UO_{2} bulk modulus (open circles) with the data points from MD simulations (red halfcircles) using the new potential.
Specific Heat Capacity:
Change in enthalpy relative to 300K as a function of temperature calculated using the new potential. Polynomials have been fitted to the data for PuO_{2} ThO_{2} and UO_{2}, shown in green, purple and red respectively. The polynomial form is based on that used for the experimental data, also shown in black. The derivative of the enthalpy gives the specific heat and is shown alongside the enthalpy data.
Mixed Oxide System (U,Th)O_{2}
The results of simulations on the (U,Th)O_{2} system are given here to demonstrate the predictive ability of this model.
Thermal Expansion:
The lattice parameter as function of composition observes Vegard’s Law below 2000 K. Above this temperature, however, there is a ‘bump’ in lattice parameter associated with oxygen disorder.
Specific Heat Capacity:
The increase in specific heat capacity as function of temperature is approximately linear below 2000 K. Above this temperature, however, there is a ‘bump’ in specific heat capacity associated with the enthalpy required to generate high temperature oxygen disorder.
Oxygen Diffusivity:
The change in gradient of the diffusivity as a function of temperature (left) indicates a change in the oxygen diffusion mechanism as the system undergoes the superionic transition. At low temperatures it is apparent that the solid solutions exhibit enhanced oxygen diffusivity relative to the end member systems.
Oxygen Defect Energies:
Unlike the end members, the solid solution oxygen sites that have a wide range of possible coordination environments due to the random distribution of U and Th on the cation sublattice. The peaks in vacancy formation enthalpy are associated with the 5 possible first nearest neighbour coordination environments of the oxygen sites  i.e. from lowest to highest enthalpy they are surrounded by 4, 3, 2, 1, and 0 uranium cations (or 0, 1, 2, 3 and 4 thorium cations).
All of the solid solutions contain a significant proportion of defect enthalpies that are lower in energy that for pure UO_{2}. This contributes to the enhanced oxygen diffusivity.
Using the Model in Common Simulation Codes
Using the Model in LAMMPS
The following examples show how to use the relevant files within the Downloads table with the LAMMPS code.
Single Oxide Example
The LAMMPS_example.zip contains an example equilibration for a UO_{2} supercell with the fluorite structure. See below for a thorough description of how the model is specified to LAMMPS.
To run the example, unzip the archive, enter its directory and type:
lammps in EAM_equil.lmpin log EAM_equil.lmpout
Mixed Oxide Example
In order to show how the MOX potential model can be used in LAMMPS a working example has been provided. LAMMPS_UThO2_example.zip provides an example where a 5×5×5 (U_{0.5},Th_{0.5})O_{2} supercell is energy minimised then equilibrated under NPT conditions at 300K for 50ps. This should give a good idea of how other mixed oxide systems can be simulated by modifying the .lmpstruct
file.
Running the Example
Unzip the LAMMPS_UThO2_example.zip archive.
In your terminal enter the example directory.

Run LAMMPS:
lammps in EAM_NPT.lmpin log EAM_NPT.lmpout
Details
Both the EAM and pairwise interactions of the actinide oxide potential described below are tabulated for use in LAMMPS using the eam/alloy
pair_style, as described in the LAMMPS manual. The pair_coeff
command is used to assign the elements within a tabulated EAM file to LAMMPS species IDs. The following example shows how the CeThUNpPuAmCmO.eam.alloy
table file is used in LAMMPS.
pair_style hybrid/overlay coul/long ${SR_CUTOFF} eam/alloy
pair_coeff * * coul/long
pair_coeff * * eam/alloy CeThUNpPuAmCmO.eam.alloy O Ce Th U Np Pu Am Cm
Notes:
 It is important to make sure the order in which the species appear in the
pair_coeff
command is consistent with the IDs assigned in the.lmpstruct
file. 
pair_style hybrid/overlay
is used to combinecoul/long
with the EAM model (eam/alloy
), this tells LAMMPS to calculate the electrostatic interactions between ions (NB: akspace_style
must also be defined).  Ion charges of 2.2208 for uranium and –1.1104 for oxygen should be
specified. This can be performed either in a LAMMPS file read in
using
read_data
or using theset charge
directive. The latter approach is used in our example.  The variable
${SR_CUTOFF}
is used to define the cutoff parameter forcoul/long
.
Using the Model in GULP
The following examples demonstrate how to use the files provided in the Downloads table to perform simulations in the GULP code.
A single GULP library files is provided for all the actinide oxides and their solid solution. This can be downloaded using the link in the downloads table.
Single Oxide Example
An example showing how to energy minimise a single fluorite UO_{2} unitcell can be found in GULP_Example_UO2.zip.
Running the Example
Unzip GULP_Example_UO2.zip.
From the terminal enter the example directory

Run the provided GULP file:
gulp < GULP_UO2_fluorite_example.gin
Mixed Oxide Example:
An example showing how to energy minimise a single fluorite UO_{2}
unitcell with a Th cation substituted at one of the uranium sites using the GULP_MO2.lib
file is provided in the
GULP_UThO2_example.zip.
Running the Example
Unzip GULP_UThO2_example.zip.
From the terminal enter the example directory

Run the provided GULP file:
gulp < GULP_Th_in_UO2_example.gin
Note: The other actinide oxide systems can be studied by changing the structure inGULP_Th_in_UO2_example.gin
without changing the library
directive.