Handling radioactive elements is one of the research areas where modeling often does a better job than actual experiments. This is due to the fact that radioactivity is a pretty well known and well described piece of science that, at least in principle, is straightforward to model. On the other hand, experiments on radioactive materials carry an obvious risk and therefore have to be limited to minimimum, according to ALARA principle: (risk) As Low As Reasonably Achievable. To model radioactive behavior only one could employ Monte Carlo — but what if we are interested in purely mechanical or chemical properties of an element or its compound? In the first case Molecular Dynamics (MD) comes at hand and, as for the latter one, a first principles modelling will do a better job.
In this project, we are going to study the mechanical properties of uranium oxide, also known as urania, using classical Molecular Dynamics simulations. Urania crystal has a simple cubic structure with uranium occupying the corners and oxygen placed in more tetragonal positions of the primitive cell. The ratio of one to the other is 1:2, so the chemical formula of the oxide is UO2. In our model, the crystal dimensions are 15x10x5, in terms of the unit vector, so there are 3000 uranium atoms and 6000 oxygen atoms, totalling to 9000. Though this number might seem small, it is enough to provide a resonable micromechanical moduli.
The very core of each MD simulation is not solving the underlying Newton's equations of motion but rather the force field. In this term we understand all parametrization (specific to the studied system) that includes: partial charges, energies and bond lengths of covalent bonds (plus angles and dihedrals), as well as the non-bonded interactions or Van der Waals forces. There is also an extra layer of settings that are, among others: cut-off distances, correction for dispersion forces, functions used for bond dispersion (doesn't have to be harmonic), long-range Coulomb evaluation, or artificial fixed. In this project, the force field parameters are going to be taken from the paper of Arima et al. (2005). We set the uranium and oxygen to have partial charges of +4 and -2, respectively. The two types of atoms interact with each other solely through electrostatic and Van der Waals forces implemented here by Buckingham potential. Thus, the energy V of an interaction between two atoms A and B at distance r is:
V = - k qA qB r -1 + Q e - r/r0 - C r -6
where qA and qB are partial charges, r0 is equilibrium distance and k, Q, and C are other parameters.
Now, we set up the simulation for 100'000 time steps, one equivalent to 1 fs (10-15 s). To make the simulation more realistic, we apply a Nose-Hoover thermostat and a Berendsen barostat. The first one keeps the temperature oscillating around 298 K and the other one is supposed to hold the pressure of 1 bar independently for each direction. The action of these two and a natural movement of the lattice result in thermal and volume fluctuations. The fluctuations, in turn, can provide us with some thermodynamic properties, such as heat capacity, thermal expansion coefficient, and compressibility (inverse of bulk modulus). This amazing capability has been first described in Parrinello and Rahman paper in 1982. By getting averages (< . >) and variances (< .2 >-< . >2) of internal energy U and volume V, the following relations can be obtained:
The figure below shows the fluctuations of temperature and volume — these will be different for different material. Although difficult so state by only looking at seemingly chaotic time series, it can be caught by the formula above. In our case, we obtain the following values: 162 GPa of bulk modulus, 4.55e-05 1/K of thermal expansion coefficient, and 410 J/(kg K) of heat capacity — all well within experimental range.
In the subsequent step the point defects are created by removing randomly atoms in batches to gradually get to 2% defecting level. The drop is quasi-linear down to approximately 202 GPa. Voids make the average interactions between atoms weaker and overall structure less coherent.
Obviously, the point defects present in a crystal deteriorate its mechanical properties. In a nanoscale it can be conveniently studied with Molecular Dynamics. At higher scales, crystal grain boundaries have to be taken into account — this kind of dynamics is often researched with Monte Carlo simulations.
If you would like to know more details about this topic, feel free to contact us or check the provided links and references.
author: Dr. Karol Kulasinski
M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, Oxford Science Publications (1987).
T. Arima, S. Yamasaki, Y. Inagaki, K. Idemitsu, Evaluation of thermal properties of UO2 and PuO2 by equilibrium molecular dynamics simulations from 300 to 2000K, Journal of Alloys and Compounds 400, (2005), pp. 43-50.
M. Parrinello and A. Rahman, Strain fluctuations and elastic constants, J. Chem. Phys. 76 (1982).
K. Kulasinski, Physical and Mechanical Aspects of Moisture Adsorption in Wood Biopolymers Investigated with Atomistic Simulations, ETH Zurich, 2015.
Gromacs: B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theor. Comput. 4, 435 (2008)
Wikipedia: Uranium dioxide