Discovery
This document was uploaded by user and they confirmed that they have the permission to share it. If you are author or own the copyright of this book, please report to us by using this DMCA report form. Report DMCA
Overview
Download & View Discovery as PDF for free.
More details
- Words: 3,707
- Pages: 10
Molecular Modelling at AWE
T
he reaction of hydrogen with uranium and plutonium is important to both the nuclear power and nuclear weapons communities1. Given this fact, it is surprising how poorly understood this reaction actually remains. Although experimental work can provide us with valuable insight into the reaction mechanism, it can never provide the full picture. To answer questions such as ‘which species of hydrogen diffuses to the metal surface’ requires the use of molecular modelling. Here, we will illustrate some of these molecular modelling techniques and their application to uranium hydriding.
the material towards further oxidation. Thus, metal oxides are not only important as the external surfaces of most metals, but also as materials in their own right. In order to quantify the hydrogen-uranium reactions, it is necessary to understand the mechanisms involved in gas transport through the surface oxide layers. Specifically, one requires knowledge of the species of the gas involved in diffusion through the oxide, the diffusion rate and its temperature dependence, the solubility of the gas species in the oxide and the mechanism by which it diffuses, e.g. via interstitial sites, vacancies or grain-boundaries within the oxide. Our modelling studies here at AWE, at Even ‘noble’ metals such as platinum Reading University and at Imperial develop some form of oxide coating College London are aimed at answering upon exposure to an oxygen-containing these questions using some of the atmosphere. The phase and composition methods explained below. of the surface oxide layer can greatly affect the reaction of the underlying Ab initio Calculations metal with gases, e.g. the oxide film on the surface of aluminium acts as a barrier The most desirable approach to to oxygen transport through to the modelling a system involves using a set metal surface and so rapidly passivates of physical constants, fundamental
18
equations and basic assumptions and predicting physical and chemical properties of the system in an ab initio manner without regard to experimental data. In this scenario, the rôle of experimental data is to provide a means by which modelling results may be validated, and not to predicate the modelling. This can be achieved if one uses methods based on quantum mechanics, which enable the calculation of structural, electronic and thermodynamic properties of materials and quantum molecular dynamics, which enable transport properties and other time and temperature-dependent properties to be calculated. Each of these techniques has proven to be a very powerful predictive tool when applied to small systems containing only a relatively low number of electrons. However the associated calculations rapidly become impossible when the number of electrons increases beyond approximately 500-1000. In fact, quantum molecular dynamics calculations can even become intractable problems when this number is much smaller.
Discovery • The Science & Technology Journal of AWE
In addition to the problems encountered with systems containing a relatively large number of electrons, actinides contain f electrons that are problematic in their own right for a number of reasons:
stoichiometric complexity. Thus, a calculation of the diffusion coefficient of a hydrogen atom in UO2 would
require a simulation cell that contained of the order of 300 uranium ions (27,000 electrons), which far exceeds f electrons are quite diffuse in what is technically possible on today’s nature, i.e. their associated super-computers. In fact, it is unlikely electron density is spread relatively that such calculations will be feasible evenly over a large radius from the within the next ten years, even if nucleus, making overlap and computer power increases at the rapid consequent interaction with other pace it has in the past 20 years. atoms and ions weak; partially filled f shells lead to low-lying unoccupied states and these cause convergence problems;
Molecular Modelling Calculations
when summed with the coulombic interactions between charged species the energy of a particular configuration of the ions is calculable. In systems which contain highly charged species, the van der Waals interactions represented by the pair-potentials detailed above, contribute approximately 10% of the interaction energy of the system. Thus, small errors in the pair-potentials in these systems can easily be accommodated and will not be significantly detrimental to the quality of the results2.
Atomistic simulations include interactions between the electrons and nuclei Fortunately, one can simplify the implicitly and not explicitly and are problem of calculating diffusion consequently thousands of times the f electrons can be localised on coefficients and solubility curves for quicker to perform on a system with an an individual ion or spread over atoms and molecules in materials by equivalent number of atoms. Key to several ions (itinerancy) which can setting up models of the system in the success of atomistic simulations change with charge state; which the atoms and ions are is the quality of the models for the represented by point charges which relativistic effects are very important interact with each other according to a individual ions (shell model or better, for example) and the pair-potentials set of functions or tabulated values. to actinide electronic properties. between the ions, though in highly This affects the energetic ordering This approach, sometimes called charged systems the pair-potentials molecular modelling or atomistic of the orbitals and means that simulation, requires parameters for the are of lesser importance. In systems spin-orbit coupling should be charge on each ion, or pair of charges containing uncharged atoms, such as included in the calculations since this will have a large effect on the for each ion in the shell model (Figure 1). hydrogen atoms, there are no permanent coulombic interactions, only instantaneous spin properties of the actinide atoms The charge in the rigid-ion model is ones from temporary polarisation of the or ions and magnetic properties of normally set to the formal charge on the ion, but this is not always the case, atoms by charged ions in the vicinity of the materials they constitute; the atom. In these cases the success of and using partial charges on ions can the simulations will rely heavily upon lead to unpredictable results and is many computer programs cannot the quality of the pair-potentials used another variable that needs to be handle f orbitals. defined in the model. The interactions to describe the van der Waals interactions between the atoms and ions and on the A purely ab initio or quantum mechanical between the ions and atoms in the treatment of these systems is impractical. system are defined by a set of functions shell-model parameters for the atoms themselves representing the potential energy for The accurate calculation of diffusion coefficients and solubility curves of gas each combination of ions and atoms Atomistic simulations can be at different separations; known as atoms and molecules in ceramics parameterised by a number of methods: pair-potentials since they act upon requires models for the ceramic that pairs of ions or atoms (Figure 2). These by fully quantum mechanical calculations contain of the order of 1000 atoms; of isolated systems such as pairs of pair-potentials represent the van der more if the oxide does not have a Waals interactions between species and ions or atoms; by semi-empirical simple crystal structure or has a
19
Figure 1 Core shell spring Shell model qcore q shell Rigid-Ion Model q
±
a) The rigid-ion model for ions and atoms, and b) the shell model for atoms and ions. In the rigid-ion model the charge on the point charge representing the ion is normally the formal charge on the ion, but some parameter sets do vary. b) In the shell model, the formal charge on the ion is usually the sum of the charge on the core and shell, qion = qcore + qshell
Equation 1
Etotal = E coulombic + E bond + E rotation + E torsion + E vdW + ( E H −bond )
Equation 2
qi q j e 2 r 4πε o ij j =1 i > j
N −1 N
Ecoulombic = where:
∑ ∑
calculations, such as electron-gas calculations, on isolated systems containing pairs of ions, where the electron density distribution of the ions is fixed at some unperturbed value and by fitting to bulk properties of the material to be simulated, such as bulk modulus elastic constants and lattice energy. In systems such as metal oxides the latter method is most desirable since it guarantees that the parameter set will reproduce the bulk properties of the material and thus its response to external stimuli such as isotropic and non-isotropic pressure changes and heating. In systems where no crystal structure data can exist, such as hydrogen interactions with an actinide oxide, the former two methods must be applied to the problem. Molecular modelling has been used to compute the energy of systems with the ions being described by one or other of the two models outlined above. The Force-Field method relies on the separation of energy terms into component parts, each of which can be equated to a physically meaningful quantity (equation 1). In the predominantly ionic systems under investigation in this work, the bond stretch term, Ebond, the bond
N=total number of ions qi, qj = formal charge on ion x
rotation term, Erotation, and the bond
rij = separation between ions i and j
torsion term, Etorsion, can be neglected and only the coulombic term, Ecoulombic, and van der Waals term, EvdW, need be
Equation 3
E vdW =
− r ij B A − exp ∑∑ ρ r ij6 j =1 i > j
N −1 N
20
considered. The coulombic term takes the form of the function that describes the potential energy between two charged particles, but summed over all possible pair interactions (equation 2). This coulombic term converges only very slowly as the separation between ions increases, so a mathematical
Discovery • The Science & Technology Journal of AWE
method known as the Ewald summation is used to model the interactions. This requires three terms: a real space cut-off distance; a reciprocal space cut-off and an Ewald sum constant. The mathematical formalism can be found in standard textbooks3.
Figure 2 5.0 Atomic Separation ( Å) R 4.0
3.0
The van der Waals term can take one of several forms, but the commonest of these is the Buckingham potential shown in equation 3.
D
2.0
1.0
These functions produce a potential energy curve somewhat like that plotted in Figure 2, having a potential minimum of depth, Dφ, at some separation of ions, Rφ. The repulsive part of the curve comes from the exponential term in the Buckingham function, and emanates from the repulsion of electron gas clouds around ions as they are brought together. The attractive part of the curve derives from the r-6 term and is a response to the dispersion forces sometimes called van der Waals attraction. Many properties that we wish to calculate depend up time and temperature, but molecular modelling as described above only deals with static systems. In order to calculate dynamic properties, molecular dynamics simulations must be carried out. These simulations involve solving the Newtonian equations of motion at very small time intervals (time steps) to calculate the positions of the atoms in a system as a series of ‘snapshots’ in time. Infinite systems such as crystals can be examined by using periodic boundary conditions, in which the simulation cell is repeated on all sides so that atoms leaving the cell on one side are never lost to the system because they re-enter the cell from the opposite side (Figure 3). The average
φ
0.0 φ -1.0 0.0
2.0
4.0
6.0
8.0
10.0
Potential Energy (kJmol -1 ) A plot of the Buckingham potential
Figure 3
The periodic boundary condition as applied to a simulation of UO2-x. The simulation cell is surrounded by replicas of itself such that an oxide ion leaving the cell on the left re-enters the central cell from the right hand cell
21
Figure 4
kinetic energy that the atoms in the system possess is determined by the temperature at which the simulation is carried out. Molecular dynamics simulations allow us to calculate properties of materials such as diffusion coefficients, melting temperatures and vibrational spectra.
X
Example of the Problems Addressed
The unit cell of UO2, with uranium ions represented by blue spheres, oxide ions represented by red spheres and the interstitial site indicated by a green cross
Figure 5
The reaction of hydrogen with uranium is a complex process that depends on many variables, one of which is the oxide over-layer1. The oxide that is formed on uranium adopts the fluorite structure which has an empty space at its centre, known as an interstitial site (marked by a green X in Figure 4). Though it is known that the thickness of the surface oxide layer influences the rate of reaction of hydrogen with uranium, it is not known what hydrogen species is transported through the oxide layer or the mechanism by which it diffuses. It is important to understand these factors in order to quantify the reaction of hydrogen with uranium and we are addressing these questions with our molecular modelling studies. We are investigating whether the hydrogen diffuses via interstitial sites, vacancies in the crystal or grain-boundaries of the oxide. We have carried out molecular dynamics simulations at a number of temperatures of perfect crystals of UO2 with different species of hydrogen in order to calculate the propensity of hydrogen to diffuse via the interstitial sites within the crystal (Figure 5).
Four unit cells of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. Atom diffusion via interstitial sites as indicated by curved arrow
22
We have built models of imperfect crystals of UO2 containing oxide vacancies
Discovery • The Science & Technology Journal of AWE
and hydrogen species and calculated the diffusion coefficients using molecular dynamics simulations at a range of temperatures (Figure 6). The lowest energy structures of 20 low Millar-index surfaces have been determined by molecular dynamics simulations at 300 K and this has revealed that some surfaces are very stable while some are very unstable as created. Those that are the most stable generally require little movement of the atoms in the top five surface layers and thus do not change much in appearance. One such example of a naturally stable surface is the [1 1 1] surface shown in Figure 7. Higher energy surfaces, such as the [4 2 0] shown in Figure 8, tend to reconstruct to minimise the surface energy. The reconstruction of the high energy surfaces can lead to the formation of microfacets, which are ‘mini-surfaces’ on top of the main surface. In Figure 8 the [4 2 0] surface has reconstructed to form [1 1 0] and [1 0 0] ‘mini-surfaces’ on top leading to the saw-tooth effect observed in this case.
Figure 6
Four unit cells of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. An arrow indicates the path of a diffusing atom through an oxide vacancy site (indicated by a green dot)
Figure 7
A number of models of surfaces have been combined to create models for grain-boundaries in UO2 and used in molecular dynamics simulations at 300K to determine which are energetically the most stable and might provide pathways for hydrogen diffusion (Figure 9). Recently we have carried out molecular dynamics simulations of helium atoms in our models for UO2 grain-boundaries and determined activation barriers from Arrhenius plots of the diffusion coefficients at different temperatures (Figure 10). These simulations allow us to place an upper limit of the rates of diffusion of hydrogen at different temperatures.
The calculated lowest energy structure of the [1 1 1] surface of UO2 at 300K with the surface running horizontally perpendicular to the page
23
Figure 8
The calculated lowest energy structure of the [4 2 0] surface of UO2 at 300K exhibiting the rumpling of microfaceting while the surface runs horizontally perpendicular to the page
Figure 9
Parameterisation of Materials The physical properties of actinide materials define their response to external stresses and strains imposed by the environment and intrinsic properties of the materials such as self-diffusion, defect energies and specific heat capacity. The oxides, carbides and hydrides of plutonium and the hydrides of uranium fall into this category. Little experimental data are available for these materials since they are difficult to prepare pure, either in a given stoichiometry or as a particular allotrope, are difficult to handle, and, of course, they are radioactive. Molecular modelling studies of these materials require physical properties to parameterise the models, but these data are unavailable for the reasons stated. Thus, they are ideal candidates for pioneering ab initio computational studies of their electronic, structural and mechanical properties. Workers at Reading University4 are using plane-wave density functional theory to calculate the ground-state structures of plutonium oxides such as PuO2, α-Pu2O3 and β-Pu2O3, and from
X
X
X
A model of the [1 1 0]-[2 2 1] grain-boundary of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. The grain-boundary runs horizontally across the page through the centre of the picture. Three green crosses indicate channels running into the page, down which diffusion can occur
24
these their mechanical properties such as bulk modulus and elastic constants, C11, C12 and C44. Very high cut-off energies are required in simulations of these systems, which are equivalent to extremely high quality basis sets in traditional ab initio calculations. This requirement for high cut-off energies means that these calculations are extraordinarily computer resource intensive and so Reading University makes use of the super-computer power at AWE. (See Box 1).
Discovery • The Science & Technology Journal of AWE
In Conclusion
they test current methodology to its limits and we require state-of-the-art Molecular modelling allows us to probe computer software to carry out our materials that are difficult or dangerous simulations. to handle and assess conditions that Our initial results suggest that both are impossible to achieve in the molecular and atomic hydrogen are laboratory. We routinely carry out immobilised in perfect crystals of UO2, simulations of radioactive materials
quantum mechanical methods and will use this parameter set in our future simulations of hydrogen solubility and diffusion in imperfect crystals and grain-boundaries of uranium and plutonium oxides. These simulations will contribute to both a better understanding of the mechanism of being trapped in the interstitial (vacant) hydrogen transport through actinide sites of the lattice. Our current studies, oxides and data for numerical models of uranium and plutonium hydriding. however, with imperfect crystals, With every simulation, we are gaining surfaces and grain-boundaries of UO2 more information about materials that indicate that helium atoms are much are technically very challenging to more mobile in these non-ideal experimentalists and theoreticians alike, environments, with activation barriers -1 to diffusion in the range 20-80 kJmol . actinide oxides.
such as plutonium and uranium oxides in order to determine diffusion coefficients of atoms, ions and molecules at temperatures up to 3000 K, which would be incredibly difficult to achieve in the laboratory. Our simulations are, however, very demanding of computational resources and require the huge memory capacity and processor We are currently defining a parameter speed of the AWE super-computers to complete successfully. Not only are our set for hydrogen interactions with actinide oxides using high level calculations computer-intensive, but
Figure 10 -21.0
He diffusion coefficients
ln D(D in m 2s -1 )1/T (10 -3 K)
-21.5 -22.0 -22.5 -23.0 -23.5 -24.0 -24.5 -25.0
0.6
0.8
1.0
1.2
Activation Energy = 25.9
1.4
kJmol -1+/-
2.5
1.6
1.8
kJmol -1
An Arrhenius plot of our helium diffusion data for the [3 2 1]-[2 2 1] grain-boundary in UO2
25
Box 1
Computer Resources at AWE
Figure 11
At AWE we have two super-computers which are constantly upgraded. A massively parallel IBM SP3 super computer5 shown in Figure 11 is equipped with 1920 Power3 II processors, 2 Terabytes (2 million Megabytes) of memory, and connected by Nighthawk II network switches provides 3 TeraFLOPS (3 x 1012 floating point operations per second) of computing power. This machine is used for calculations that can be run in parallel and require a lot of memory. Plane-wave DFT simulations of actinide oxides typically run on 40 to 250 processors over a period of 2 to 14 days and require 20 to 120 Gigabytes of memory! The 3 TeraFLOPS IBM super-computer at AWE
AWE’s beowolf super-computer shown in Figure 12 has recently been upgraded to 2.8 GHz Pentium IV Xeon CPUs delivering 1.2 TeraFLOPS with Gigabit and Myrinet networking, 552 Gigabytes of memory and 8.4 Terabytes of disk space. This machine is used for serial and parallel molecular dynamics calculations up to 32 processors requiring days or weeks of computational time.
Figure 12
The beowolf super-computer at AWE
26
Discovery • The Science & Technology Journal of AWE
Acknowledgements
References
Author Profile
The author wishes to thank Dr. Joseph 1 Glascott, AWE, Prof. Michael G.B. Drew and Dr. Mark T. Storr, University of Reading and the AWE super-computing staff of 2 Ash Vadgamma, John Dolphin, Mark Roberts, Paul Tomlinson and Pete Brown.
J. Glascott, “Hydrogen and Uranium”, Discovery, Issue 6, January 2003
3
M. P. Allen and D. J. Tildesley, “Computer Simulations of Liquids”, Clarendon Press, Oxford
4
M. G. B. Drew, M. T. Storr and D. W. Price, unpublished work.
5
R. W. Grimes, Imperial College, London University, personal communication
A. Baxter, “New Supercomputer”, Discovery, Issue 6, January 2003; J. Taylor, “History of Computing”, Discovery, Issue 3, July 2001
David can be contacted on e-mail: [email protected]
David Price David was educated at Kingshurst Comprehensive then Solihull Sixth Form College before going up to Keble College Oxford to read for a BA(Hons) in Chemistry. He stayed at Keble to read for a DPhil in organometallic chemistry, which encompassed synthetic, mechanistic and molecular modelling work. After post-doctoral spells in Strasbourg and Munich, David became a research fellow at the University of Reading funded by AWE. David joined the materials modelling team at AWE in 2002 after spending eight years at Reading University. He currently works on the ageing of materials, is a published author and is a full Member of the Royal Society of Chemistry and Chartered Chemist.
27
T
he reaction of hydrogen with uranium and plutonium is important to both the nuclear power and nuclear weapons communities1. Given this fact, it is surprising how poorly understood this reaction actually remains. Although experimental work can provide us with valuable insight into the reaction mechanism, it can never provide the full picture. To answer questions such as ‘which species of hydrogen diffuses to the metal surface’ requires the use of molecular modelling. Here, we will illustrate some of these molecular modelling techniques and their application to uranium hydriding.
the material towards further oxidation. Thus, metal oxides are not only important as the external surfaces of most metals, but also as materials in their own right. In order to quantify the hydrogen-uranium reactions, it is necessary to understand the mechanisms involved in gas transport through the surface oxide layers. Specifically, one requires knowledge of the species of the gas involved in diffusion through the oxide, the diffusion rate and its temperature dependence, the solubility of the gas species in the oxide and the mechanism by which it diffuses, e.g. via interstitial sites, vacancies or grain-boundaries within the oxide. Our modelling studies here at AWE, at Even ‘noble’ metals such as platinum Reading University and at Imperial develop some form of oxide coating College London are aimed at answering upon exposure to an oxygen-containing these questions using some of the atmosphere. The phase and composition methods explained below. of the surface oxide layer can greatly affect the reaction of the underlying Ab initio Calculations metal with gases, e.g. the oxide film on the surface of aluminium acts as a barrier The most desirable approach to to oxygen transport through to the modelling a system involves using a set metal surface and so rapidly passivates of physical constants, fundamental
18
equations and basic assumptions and predicting physical and chemical properties of the system in an ab initio manner without regard to experimental data. In this scenario, the rôle of experimental data is to provide a means by which modelling results may be validated, and not to predicate the modelling. This can be achieved if one uses methods based on quantum mechanics, which enable the calculation of structural, electronic and thermodynamic properties of materials and quantum molecular dynamics, which enable transport properties and other time and temperature-dependent properties to be calculated. Each of these techniques has proven to be a very powerful predictive tool when applied to small systems containing only a relatively low number of electrons. However the associated calculations rapidly become impossible when the number of electrons increases beyond approximately 500-1000. In fact, quantum molecular dynamics calculations can even become intractable problems when this number is much smaller.
Discovery • The Science & Technology Journal of AWE
In addition to the problems encountered with systems containing a relatively large number of electrons, actinides contain f electrons that are problematic in their own right for a number of reasons:
stoichiometric complexity. Thus, a calculation of the diffusion coefficient of a hydrogen atom in UO2 would
require a simulation cell that contained of the order of 300 uranium ions (27,000 electrons), which far exceeds f electrons are quite diffuse in what is technically possible on today’s nature, i.e. their associated super-computers. In fact, it is unlikely electron density is spread relatively that such calculations will be feasible evenly over a large radius from the within the next ten years, even if nucleus, making overlap and computer power increases at the rapid consequent interaction with other pace it has in the past 20 years. atoms and ions weak; partially filled f shells lead to low-lying unoccupied states and these cause convergence problems;
Molecular Modelling Calculations
when summed with the coulombic interactions between charged species the energy of a particular configuration of the ions is calculable. In systems which contain highly charged species, the van der Waals interactions represented by the pair-potentials detailed above, contribute approximately 10% of the interaction energy of the system. Thus, small errors in the pair-potentials in these systems can easily be accommodated and will not be significantly detrimental to the quality of the results2.
Atomistic simulations include interactions between the electrons and nuclei Fortunately, one can simplify the implicitly and not explicitly and are problem of calculating diffusion consequently thousands of times the f electrons can be localised on coefficients and solubility curves for quicker to perform on a system with an an individual ion or spread over atoms and molecules in materials by equivalent number of atoms. Key to several ions (itinerancy) which can setting up models of the system in the success of atomistic simulations change with charge state; which the atoms and ions are is the quality of the models for the represented by point charges which relativistic effects are very important interact with each other according to a individual ions (shell model or better, for example) and the pair-potentials set of functions or tabulated values. to actinide electronic properties. between the ions, though in highly This affects the energetic ordering This approach, sometimes called charged systems the pair-potentials molecular modelling or atomistic of the orbitals and means that simulation, requires parameters for the are of lesser importance. In systems spin-orbit coupling should be charge on each ion, or pair of charges containing uncharged atoms, such as included in the calculations since this will have a large effect on the for each ion in the shell model (Figure 1). hydrogen atoms, there are no permanent coulombic interactions, only instantaneous spin properties of the actinide atoms The charge in the rigid-ion model is ones from temporary polarisation of the or ions and magnetic properties of normally set to the formal charge on the ion, but this is not always the case, atoms by charged ions in the vicinity of the materials they constitute; the atom. In these cases the success of and using partial charges on ions can the simulations will rely heavily upon lead to unpredictable results and is many computer programs cannot the quality of the pair-potentials used another variable that needs to be handle f orbitals. defined in the model. The interactions to describe the van der Waals interactions between the atoms and ions and on the A purely ab initio or quantum mechanical between the ions and atoms in the treatment of these systems is impractical. system are defined by a set of functions shell-model parameters for the atoms themselves representing the potential energy for The accurate calculation of diffusion coefficients and solubility curves of gas each combination of ions and atoms Atomistic simulations can be at different separations; known as atoms and molecules in ceramics parameterised by a number of methods: pair-potentials since they act upon requires models for the ceramic that pairs of ions or atoms (Figure 2). These by fully quantum mechanical calculations contain of the order of 1000 atoms; of isolated systems such as pairs of pair-potentials represent the van der more if the oxide does not have a Waals interactions between species and ions or atoms; by semi-empirical simple crystal structure or has a
19
Figure 1 Core shell spring Shell model qcore q shell Rigid-Ion Model q
±
a) The rigid-ion model for ions and atoms, and b) the shell model for atoms and ions. In the rigid-ion model the charge on the point charge representing the ion is normally the formal charge on the ion, but some parameter sets do vary. b) In the shell model, the formal charge on the ion is usually the sum of the charge on the core and shell, qion = qcore + qshell
Equation 1
Etotal = E coulombic + E bond + E rotation + E torsion + E vdW + ( E H −bond )
Equation 2
qi q j e 2 r 4πε o ij j =1 i > j
N −1 N
Ecoulombic = where:
∑ ∑
calculations, such as electron-gas calculations, on isolated systems containing pairs of ions, where the electron density distribution of the ions is fixed at some unperturbed value and by fitting to bulk properties of the material to be simulated, such as bulk modulus elastic constants and lattice energy. In systems such as metal oxides the latter method is most desirable since it guarantees that the parameter set will reproduce the bulk properties of the material and thus its response to external stimuli such as isotropic and non-isotropic pressure changes and heating. In systems where no crystal structure data can exist, such as hydrogen interactions with an actinide oxide, the former two methods must be applied to the problem. Molecular modelling has been used to compute the energy of systems with the ions being described by one or other of the two models outlined above. The Force-Field method relies on the separation of energy terms into component parts, each of which can be equated to a physically meaningful quantity (equation 1). In the predominantly ionic systems under investigation in this work, the bond stretch term, Ebond, the bond
N=total number of ions qi, qj = formal charge on ion x
rotation term, Erotation, and the bond
rij = separation between ions i and j
torsion term, Etorsion, can be neglected and only the coulombic term, Ecoulombic, and van der Waals term, EvdW, need be
Equation 3
E vdW =
− r ij B A − exp ∑∑ ρ r ij6 j =1 i > j
N −1 N
20
considered. The coulombic term takes the form of the function that describes the potential energy between two charged particles, but summed over all possible pair interactions (equation 2). This coulombic term converges only very slowly as the separation between ions increases, so a mathematical
Discovery • The Science & Technology Journal of AWE
method known as the Ewald summation is used to model the interactions. This requires three terms: a real space cut-off distance; a reciprocal space cut-off and an Ewald sum constant. The mathematical formalism can be found in standard textbooks3.
Figure 2 5.0 Atomic Separation ( Å) R 4.0
3.0
The van der Waals term can take one of several forms, but the commonest of these is the Buckingham potential shown in equation 3.
D
2.0
1.0
These functions produce a potential energy curve somewhat like that plotted in Figure 2, having a potential minimum of depth, Dφ, at some separation of ions, Rφ. The repulsive part of the curve comes from the exponential term in the Buckingham function, and emanates from the repulsion of electron gas clouds around ions as they are brought together. The attractive part of the curve derives from the r-6 term and is a response to the dispersion forces sometimes called van der Waals attraction. Many properties that we wish to calculate depend up time and temperature, but molecular modelling as described above only deals with static systems. In order to calculate dynamic properties, molecular dynamics simulations must be carried out. These simulations involve solving the Newtonian equations of motion at very small time intervals (time steps) to calculate the positions of the atoms in a system as a series of ‘snapshots’ in time. Infinite systems such as crystals can be examined by using periodic boundary conditions, in which the simulation cell is repeated on all sides so that atoms leaving the cell on one side are never lost to the system because they re-enter the cell from the opposite side (Figure 3). The average
φ
0.0 φ -1.0 0.0
2.0
4.0
6.0
8.0
10.0
Potential Energy (kJmol -1 ) A plot of the Buckingham potential
Figure 3
The periodic boundary condition as applied to a simulation of UO2-x. The simulation cell is surrounded by replicas of itself such that an oxide ion leaving the cell on the left re-enters the central cell from the right hand cell
21
Figure 4
kinetic energy that the atoms in the system possess is determined by the temperature at which the simulation is carried out. Molecular dynamics simulations allow us to calculate properties of materials such as diffusion coefficients, melting temperatures and vibrational spectra.
X
Example of the Problems Addressed
The unit cell of UO2, with uranium ions represented by blue spheres, oxide ions represented by red spheres and the interstitial site indicated by a green cross
Figure 5
The reaction of hydrogen with uranium is a complex process that depends on many variables, one of which is the oxide over-layer1. The oxide that is formed on uranium adopts the fluorite structure which has an empty space at its centre, known as an interstitial site (marked by a green X in Figure 4). Though it is known that the thickness of the surface oxide layer influences the rate of reaction of hydrogen with uranium, it is not known what hydrogen species is transported through the oxide layer or the mechanism by which it diffuses. It is important to understand these factors in order to quantify the reaction of hydrogen with uranium and we are addressing these questions with our molecular modelling studies. We are investigating whether the hydrogen diffuses via interstitial sites, vacancies in the crystal or grain-boundaries of the oxide. We have carried out molecular dynamics simulations at a number of temperatures of perfect crystals of UO2 with different species of hydrogen in order to calculate the propensity of hydrogen to diffuse via the interstitial sites within the crystal (Figure 5).
Four unit cells of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. Atom diffusion via interstitial sites as indicated by curved arrow
22
We have built models of imperfect crystals of UO2 containing oxide vacancies
Discovery • The Science & Technology Journal of AWE
and hydrogen species and calculated the diffusion coefficients using molecular dynamics simulations at a range of temperatures (Figure 6). The lowest energy structures of 20 low Millar-index surfaces have been determined by molecular dynamics simulations at 300 K and this has revealed that some surfaces are very stable while some are very unstable as created. Those that are the most stable generally require little movement of the atoms in the top five surface layers and thus do not change much in appearance. One such example of a naturally stable surface is the [1 1 1] surface shown in Figure 7. Higher energy surfaces, such as the [4 2 0] shown in Figure 8, tend to reconstruct to minimise the surface energy. The reconstruction of the high energy surfaces can lead to the formation of microfacets, which are ‘mini-surfaces’ on top of the main surface. In Figure 8 the [4 2 0] surface has reconstructed to form [1 1 0] and [1 0 0] ‘mini-surfaces’ on top leading to the saw-tooth effect observed in this case.
Figure 6
Four unit cells of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. An arrow indicates the path of a diffusing atom through an oxide vacancy site (indicated by a green dot)
Figure 7
A number of models of surfaces have been combined to create models for grain-boundaries in UO2 and used in molecular dynamics simulations at 300K to determine which are energetically the most stable and might provide pathways for hydrogen diffusion (Figure 9). Recently we have carried out molecular dynamics simulations of helium atoms in our models for UO2 grain-boundaries and determined activation barriers from Arrhenius plots of the diffusion coefficients at different temperatures (Figure 10). These simulations allow us to place an upper limit of the rates of diffusion of hydrogen at different temperatures.
The calculated lowest energy structure of the [1 1 1] surface of UO2 at 300K with the surface running horizontally perpendicular to the page
23
Figure 8
The calculated lowest energy structure of the [4 2 0] surface of UO2 at 300K exhibiting the rumpling of microfaceting while the surface runs horizontally perpendicular to the page
Figure 9
Parameterisation of Materials The physical properties of actinide materials define their response to external stresses and strains imposed by the environment and intrinsic properties of the materials such as self-diffusion, defect energies and specific heat capacity. The oxides, carbides and hydrides of plutonium and the hydrides of uranium fall into this category. Little experimental data are available for these materials since they are difficult to prepare pure, either in a given stoichiometry or as a particular allotrope, are difficult to handle, and, of course, they are radioactive. Molecular modelling studies of these materials require physical properties to parameterise the models, but these data are unavailable for the reasons stated. Thus, they are ideal candidates for pioneering ab initio computational studies of their electronic, structural and mechanical properties. Workers at Reading University4 are using plane-wave density functional theory to calculate the ground-state structures of plutonium oxides such as PuO2, α-Pu2O3 and β-Pu2O3, and from
X
X
X
A model of the [1 1 0]-[2 2 1] grain-boundary of UO2, with uranium ions represented by blue spheres and oxide ions represented by red spheres. The grain-boundary runs horizontally across the page through the centre of the picture. Three green crosses indicate channels running into the page, down which diffusion can occur
24
these their mechanical properties such as bulk modulus and elastic constants, C11, C12 and C44. Very high cut-off energies are required in simulations of these systems, which are equivalent to extremely high quality basis sets in traditional ab initio calculations. This requirement for high cut-off energies means that these calculations are extraordinarily computer resource intensive and so Reading University makes use of the super-computer power at AWE. (See Box 1).
Discovery • The Science & Technology Journal of AWE
In Conclusion
they test current methodology to its limits and we require state-of-the-art Molecular modelling allows us to probe computer software to carry out our materials that are difficult or dangerous simulations. to handle and assess conditions that Our initial results suggest that both are impossible to achieve in the molecular and atomic hydrogen are laboratory. We routinely carry out immobilised in perfect crystals of UO2, simulations of radioactive materials
quantum mechanical methods and will use this parameter set in our future simulations of hydrogen solubility and diffusion in imperfect crystals and grain-boundaries of uranium and plutonium oxides. These simulations will contribute to both a better understanding of the mechanism of being trapped in the interstitial (vacant) hydrogen transport through actinide sites of the lattice. Our current studies, oxides and data for numerical models of uranium and plutonium hydriding. however, with imperfect crystals, With every simulation, we are gaining surfaces and grain-boundaries of UO2 more information about materials that indicate that helium atoms are much are technically very challenging to more mobile in these non-ideal experimentalists and theoreticians alike, environments, with activation barriers -1 to diffusion in the range 20-80 kJmol . actinide oxides.
such as plutonium and uranium oxides in order to determine diffusion coefficients of atoms, ions and molecules at temperatures up to 3000 K, which would be incredibly difficult to achieve in the laboratory. Our simulations are, however, very demanding of computational resources and require the huge memory capacity and processor We are currently defining a parameter speed of the AWE super-computers to complete successfully. Not only are our set for hydrogen interactions with actinide oxides using high level calculations computer-intensive, but
Figure 10 -21.0
He diffusion coefficients
ln D(D in m 2s -1 )1/T (10 -3 K)
-21.5 -22.0 -22.5 -23.0 -23.5 -24.0 -24.5 -25.0
0.6
0.8
1.0
1.2
Activation Energy = 25.9
1.4
kJmol -1+/-
2.5
1.6
1.8
kJmol -1
An Arrhenius plot of our helium diffusion data for the [3 2 1]-[2 2 1] grain-boundary in UO2
25
Box 1
Computer Resources at AWE
Figure 11
At AWE we have two super-computers which are constantly upgraded. A massively parallel IBM SP3 super computer5 shown in Figure 11 is equipped with 1920 Power3 II processors, 2 Terabytes (2 million Megabytes) of memory, and connected by Nighthawk II network switches provides 3 TeraFLOPS (3 x 1012 floating point operations per second) of computing power. This machine is used for calculations that can be run in parallel and require a lot of memory. Plane-wave DFT simulations of actinide oxides typically run on 40 to 250 processors over a period of 2 to 14 days and require 20 to 120 Gigabytes of memory! The 3 TeraFLOPS IBM super-computer at AWE
AWE’s beowolf super-computer shown in Figure 12 has recently been upgraded to 2.8 GHz Pentium IV Xeon CPUs delivering 1.2 TeraFLOPS with Gigabit and Myrinet networking, 552 Gigabytes of memory and 8.4 Terabytes of disk space. This machine is used for serial and parallel molecular dynamics calculations up to 32 processors requiring days or weeks of computational time.
Figure 12
The beowolf super-computer at AWE
26
Discovery • The Science & Technology Journal of AWE
Acknowledgements
References
Author Profile
The author wishes to thank Dr. Joseph 1 Glascott, AWE, Prof. Michael G.B. Drew and Dr. Mark T. Storr, University of Reading and the AWE super-computing staff of 2 Ash Vadgamma, John Dolphin, Mark Roberts, Paul Tomlinson and Pete Brown.
J. Glascott, “Hydrogen and Uranium”, Discovery, Issue 6, January 2003
3
M. P. Allen and D. J. Tildesley, “Computer Simulations of Liquids”, Clarendon Press, Oxford
4
M. G. B. Drew, M. T. Storr and D. W. Price, unpublished work.
5
R. W. Grimes, Imperial College, London University, personal communication
A. Baxter, “New Supercomputer”, Discovery, Issue 6, January 2003; J. Taylor, “History of Computing”, Discovery, Issue 3, July 2001
David can be contacted on e-mail: [email protected]
David Price David was educated at Kingshurst Comprehensive then Solihull Sixth Form College before going up to Keble College Oxford to read for a BA(Hons) in Chemistry. He stayed at Keble to read for a DPhil in organometallic chemistry, which encompassed synthetic, mechanistic and molecular modelling work. After post-doctoral spells in Strasbourg and Munich, David became a research fellow at the University of Reading funded by AWE. David joined the materials modelling team at AWE in 2002 after spending eight years at Reading University. He currently works on the ageing of materials, is a published author and is a full Member of the Royal Society of Chemistry and Chartered Chemist.
27
