- Email: [email protected]
Borosilicate glass potentials for radiation damage simulations
Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
Contents lists available at ScienceDirect
Nuclear Instruments and Methods in Physics Research B journal homepage: www.elsevier.com/locate/nimb
Borosilicate glass potentials for radiation damage simulations Kenny Jolley a,⇑, Roger Smith a, Kitheri Joseph b a b
Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, UK Materials Chemistry Division, Indira Gandhi Centre for Atomic Research, Kalpakkam 603 102, India
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 25 June 2014 Received in revised form 8 December 2014 Accepted 9 December 2014 Available online 26 December 2014
Three borosilicate glass (SiO2–B2O3) fixed charge potentials from the literature are compared (Delaye and Ghaleb, 1996; Kieu et al., 2011; Rushton, 2006) and their suitability for use in simulations of radiation damage is assessed. For a range of densities, we generate glass structures by quenching at 5 1012 K/s using constant volume Molecular Dynamics. In each case, the bond lengths, mean bond angles, bulk modulus, melting point and displacement energy thresholds are calculated, and where possible compared to experimental data. Whereas the bond lengths and mean bond angles are reasonably well predicted, we find that the potentials predict melting temperatures, bulk moduli and densities that are higher than experimental data. The displacement energy thresholds are generally lower than those for ionic crystalline materials, but show a wider spread of values. However, the barriers for atomic rearrangements, after atoms have been displaced in the equilibrium structures, are very high. This indicates, that the radiation damage produced in the ballistic phase of a collision cascade, is likely to persist for extended time scales. This is in contrast to crystals, where interstitials and vacancies can diffuse rapidly between successive radiation events. Ó 2014 Elsevier B.V. All rights reserved.
Keywords: Molecular Dynamics Borosilicate glass Inter-atomic potentials Radiation damage
1. Introduction High-Level nuclear Wastes (HLW) from spent nuclear power station fuels and decommissioned nuclear weapons must be safely removed from the environment. One leading method of achieving this, is long term storage by immobilisation in glass waste forms prior to permanent disposal in a geologically stable repository [4,5]. A primary concern with this method, is that over long time-scales, the radioactive material could leach out of the glass and contaminate the environment. It is for this reason we must understand how glass responds to radiation damage over long time scales. In the early stages of waste encapsulation the material is heated due to b-decay, however, in this work it is the damage caused by a-decay that will be the main focus of the study. Experiments on physical samples are limited to relatively short times (10 years), compared to geological time scales. Experiments are also carried out at higher dose rates, such as by using ion beams to simulate the collision cascade process. As a result, atomistic simulation methods can play an important role in the understanding of how these glasses behave after a radiation event. Although Molecular Dynamics (MD) simulations can only access very short ⇑ Corresponding author. E-mail addresses: [email protected] (R. Smith), [email protected] (K. Joseph).
(K.
http://dx.doi.org/10.1016/j.nimb.2014.12.024 0168-583X/Ó 2014 Elsevier B.V. All rights reserved.
Jolley),
[email protected]
times ( 1 ns), kinetic Monte Carlo (kMC) simulation techniques can be used to access longer time scales. In order to simulate the radiation event, it is necessary to have a good model of the amorphous glass structure. Amorphous materials are more difficult to handle computationally than crystals, since there is no equivalent concept of a vacancy or interstitial. In addition, although the structures have short range order, there are many different atomic configurations that are possible, so statistical averaging over different structures is very important. In addition, good inter-atomic potentials to describe these systems is essential. In this work, we compare three inter-atomic potentials from the literature, and test their suitability for Molecular Dynamics (MD) simulations of radiation effects in borosilicate glass. We investigated two potentials from Delaye’s group (Delaye1996 [1] and Kieu2011 [2]), and one potential from Rushton (Rushton2006 [3]). Only the potentials of Kieu were developed specifically for the simulation of large sets of compositions around the SiO2–B2O3–Na2O glass. The Delaye1996 potential was developed for simulating R7T7 glass, and the Rushton2006 potential was developed for Magnox glass. The ‘best’ potential model for simulating borosilicate glass was then used to determine threshold displacement energies, and some typical energy barriers for atomic movement after such atomic displacements.
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
For the boron–oxygen interaction, Aij ¼ 180390:53 eV, qij ¼ 0:124 Å and C ij ¼ 35:0019 eV Å6.
2. Methodology Computer simulations of borosilicate glass and sodium borosilicate glass were performed using MD. For borosilicate glass, the structures consisted of 70 mol% SiO2 and 30 mol% B2O3. We refer to this glass as, BSi, in this paper. For sodium borosilicate glass, we used the SBN12 blend which is described in Ref. [2]. This glass consists of 59.66 mol% SiO2, 28.14 mol% B2O3 and 12.20 mol% Na2O. We refer to this glass as, NaBSi, in this paper. A brief overview of each inter-atomic potential that we have investigated, is given in Sections 2.1–2.3 below. These potentials are suitable for near equilibrium conditions and must be joined to a screened Coulomb potential for close particle separation. Our splining method is described in Section 2.4, and the method for generating glass structures is described in Section 2.5. 2.1. Delaye1996 The Delaye1996 potential [1] is used to model BSi glass. This potential contains both 2 and 3 body terms. The 2-body interaction is modelled by the Born–Mayer–Huggins potential. The potential energy, Uðrij Þ, of two ions separated by r ij is given by Eq. (1) below:
Uðrij Þ ¼ Aij exp
r ij
! þ
qij
1 4p0
zi zj r ij
2.2.2. Kieu2011_v2 For the sodium borosilicate glass that we studied (NaBSi), R = 0.4335 and K = 2.12. This is described as the SBN12 blend in the original paper. The ion charges were: zSi ¼ 1:869128; zO ¼ 0:965872; zB ¼ 1:528756 and zNa ¼ 0:451628. The parameter, Aij , for the boron–oxygen bond, is 205751.1464 eV. The remaining parameters (qij and C ij ) are the same as before. 2.3. Rushton2006 The Rushton2006 potential [3,9] is used for modelling BSi glass, and also contains only 2-body potentials. In this case, a LennardJones pair potential plus the Coulomb term (Eq. 3) is used to model the atomic interactions.
pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi Bi Bj Ai Aj 1 zi zj Uðr ij Þ ¼ 12 6 þ 4p0 r ij r ij r ij
ð3Þ
The parameters A and B for each ion are given in the original paper. zi and zj are the full formal charges for each ion. 2.4. Application to collision cascades
ð1Þ
where Aij , and qij are parameters of the model and can be found in the original paper [1]. zi and zj are the ion charges. This potential uses full formal charges for each ion. A Stillinger–Weber [6] 3-body term is applied to the O–Si–O and Si–O–Si bond angles. We choose not to put any restriction on the O–B–O bond angle as was done by Connelly [7]. The adjustable parameters are listed in the original paper [1] in Table 4. Due to the inclusion of the 3-body terms, this potential is the most computationally intensive. The three body term also induces an extremely large penalty when the bond angles move away from their equilibrium values. 2.2. Kieu2011
Some of these potentials tend to 1 at short range due to the electrostatic coulomb term (zi zj =rij ) and the Buckingham term (C ij =r 6ij ). Therefore, the potential energy reaches a maximum value at short range, which can be as low as 150 eV (for a Si–O bond). Hence, during a radiation damage cascade, it becomes possible for some atoms to unnaturally fuse together. The short range interactions between atomic nuclei are more accurately modelled by the ZBL potential [10]. We therefore join these potentials to the ZBL potential at short range with a spline function given by Eq. (4).
FðrÞ ¼ exp a0 þ a1 r þ a2 r 2 þ a3 r3 þ a4 r 4 þ a5 r5
ð4Þ
The constants a0 – a5 are set such that the potential energy and its first and second derivatives are continuous. 2.5. Modelling of glass formation
The original paper by Kieu et al. [2], describes potentials for modelling both borosilicate (BSi) and sodium borosilicate (NaBSi) glass. This potential contains only 2-body terms. A Buckingham pair potential (given by Eq. (2) below), is used to model these interactions.
Uðrij Þ ¼ Aij exp
141
r ij
qij
! þ
1 zi zj C ij 6 4p0 r ij r ij
ð2Þ
The adjustable parameters Aij ; qij and C ij for each ion interaction are listed in the original paper [2]. The boron–oxygen interaction is a special case, where the parameter, Aij , depends upon the composition of the glass. The ion charges (zi and zj ), also depend upon the glass composition. This composition is characterised by the two molar ratios, R = [Na2O]/[B2O3] and K = [SiO2]/[B2O3]. Since the parameter values and charges depend upon the composition of the glass, we denote the potential used for BSi glass as Kieu2011_v1, and potential used for NaBSi glass as Kieu2011_v2.
The glass structures are amorphous and have no long range order, therefore, we cannot simply generate a structure as we would for a crystal lattice. The procedure for generating a glass structure is to perform an MD simulation of a quench from a molten state. Atoms are placed at random positions inside a cubic box, ensuring no two atoms are closer than 1 Å. Such a configuration has a large potential energy, and therefore a high initial temperature when simulated with MD, typically exceeding 10,000 K. At every time-step the temperature is controlled by rescaling the velocities of all the atoms. If the temperature exceeds the desired value by 7%, then the velocities are rescaled to the desired temperature. We equilibrate at 6000 K for 10 ps then quench at 5 1012 K/s to 50 K. The final minimisation to 0 K is obtained via conjugate gradient minimisation. A schematic of this process is shown in Fig. 1. 3. Results 3.1. Density
2.2.1. Kieu2011_v1 For borosilicate glass (BSi), R = 0 and K = 2.33, and therefore, the charges obtained are close to the values used in the Guillot–Sator potential model [8]: zSi ¼ 1:89; zO ¼ 0:945 and zB ¼ 1:4175. We used these charges in the Kieu2011_v1 potential.
For each inter-atomic potential we performed a set of quenches using constant volume MD. We vary the volume to obtain sets of quenched glass structures at a range of different densities. Fig. 2 shows the plots of the average potential energy per atom as a
142
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
The original authors generated their glass structures by quenching at a density 5–10% below the experimental value, and then relaxed the system at constant pressure. This enabled them to find a stable local minimum close to the experimental density. In our simulations, we used GULP [12] to perform constant pressure relaxations of the glass structures for a range of densities. The glass structure with the optimum density did not change in volume significantly. The low density structures, increased in density on relaxation, but, still had a lattice energy that was higher than the minimum found at the optimum density. Typically, relaxations at constant pressure using GULP caused the density to change by less than 10%. Low density structures relaxed by increasing their density, while the high density structures relaxed by lowering their density. Constant pressure minimisation did not yield a better optimal structure, therefore, we believe a large sample of constant volume quenches is sufficient to obtain the optimal structure. The value of the average potential energy per atom of the optimum glass structures also varies between each potential.
Fig. 1. A schematic of the quenching procedure.
3.2. Structural properties
Fig. 2. Plots of the average total potential energy per atom of quenched lattices vs density for each inter-atomic potential. Each curve has been offset on the vertical axis so that they appear within one plot.
function of glass density for each inter-atomic potential that we studied. Each plot has been offset on the vertical axis such that all curves appear within the same figure. In each case, the data about the minimum has been fitted with a 2nd order polynomial, from which we estimate the global minimum lattice energy. We find that borosilicate glass modelled with the Delaye1996 potential, has a global minimum at a density of 1.78 g/cm3. This is lower than the experimental value of 2.04 g/cm3 [11]. All the other potentials have global minima that are higher than the experimental values. The Kieu2011_v1 and Rushton2006 potentials have minima at 2.76 and 3.67 g/cm3, respectively, for borosilicate glass. The global minimum lattice energy found for NaBSi glass, modelled with the Kieu2011_v2 potential was, 2.728 g/cm3. The Rushton2006 potential produced a glass structure that contained many over co-ordinated silicon and boron atoms. This is the reason for the significantly higher optimum density found. The higher densities that we found in our constant volume simulations using the Kieu2011 potentials, could be a high density phase, which, during an actual experiment would not be realised. This is because the cooling process takes place over longer timescales at constant pressure, and could require large energy barriers to be overcome.
We measured some structural properties of each glass. This data is summarised in Table 1. The Delaye1996 and Kieu2011 potentials give bond length and bond angle distributions that are in good agreement with experimental data. The Rushton2006 potential, however, contains Si–O and B–O bond lengths that are too long. The bond angle distributions in this case also differ due to the over-coordinated structure. The Bulk moduli, K B , of the optimum quenched glass structures for each potential is found using GULP [12], and are reported in Table 1. We see that in all cases the bulk moduli for borosilicate glass exceeds the experimental value of 23.7 GPa [2]. The sodium borosilicate glass (modelled using Kieu2011_v2) structure also has a bulk modulus that exceeds the experimental value of 42 GPa. We found that the bulk modulus decreases with decreasing density, and that when quenching a glass using the original authors method, we were able to replicate the results in the original paper, which were much closer to the experimental values. 3.3. Melting points We also measured the melting temperature predicted by each inter-atomic potential. The method for estimating the melting point is to look at the peak of a plot of the specific heat capacity as a function of temperature. Fig. 3 shows this for each interatomic potential that we investigated. To generate this plot: 1. For a range of temperatures (0–14,000 K), using the optimum quenched lattice for each potential, we thermalise the system using the Berendsen thermostat for 10 ps. 2. Then we run a constant energy MD (NVE) simulation for each case, and compute the time average of the temperature and total lattice energy. 3. The specific heat capacity is then the gradient of a plot of the total lattice energy vs temperature.
Table 1 Summary of simulated glass structural properties for each inter-atomic potential. Glass
BSi BSi NaBSi BSi
Potential
Delaye1996 Kieu2011_v1 Kieu2011_v2 Rushton2006
Density (g/cm3)
1.78 2.76 2.728 3.67
K B (GPa)
90 79 79 –
Bond length (Å)
Bond angle (°)
Si–O
B–O
O–Si–O
O–B–O
1.595 1.607 1.608 1.93
1.355 1.365 1.428 1.468
109.43 109.36 109.34 100.89
119.26 119.10 112.93 103.10
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
143
3. After 1 ps we compare the initial and final lattices to see if the atom has become permanently displaced. The cut-off for displacement is 1 Å. 4. If it has not been displaced, we double the trial energy and rerun the cascade. 5. We find the minimum threshold energy required to permanently displace the atom using the binary search algorithm. We compute the threshold to an accuracy of 1 eV in each case.
Fig. 3. A Plot of the temperature dependence of the specific heat capacity. The peak occurs at 8000 K for Delaye1996, 6000 K for Rushton2006 and 2700 K for the Kieu2011 potentials.
4. To obtain a smooth curve, we use a least squares fit on several adjacent points to obtain a moving average of the gradient. We see that the Delaye1996 potential predicts the highest melting point at approximately 8000 K. The Rushton2006 potential predicts the melting point at about 6000 K, and both the Kieu2011_v1 and Kieu2011_v2 potentials give 2700 K. The experimental value of the melting temperature of sodium borosilicate glass ranges from 1310 to 1870 K depending upon the composition of the glass [11,13]. 3.4. Displacement energy thresholds
The distributions of threshold energies are shown in Fig. 4. Here we plot data for borosilicate glass at two densities (2.0 and 2.76 g/ cm3) modelled with the Kieu2011_v1 potential and sodium borosilicate glass (at the optimum density 2.74 g/cm3) modelled by Kieu2011_v2 potential. The average threshold energy does not have any directional dependence upon cascade direction, since the glass structures are amorphous. We also find that the lower density borosilicate glass has lower threshold energies for all atomic species, compared to the higher density case. Compared to ionic crystalline systems, the thresholds are lower [14] with a more continuous spectrum. The structures were also examined after a displacement event. It was found that there were no debonded atoms or voids formed in the system, a result confirmed also by running cascades up to 1 keV. Defective regions resulting from these low energy displacements, were best identified through the co-ordination number of each atom in the system. 3.5. Transition barriers In an ionic crystal one can determine the activation energies of defects such as vacancies and interstitials in a straightforward way.
We calculate the displacement threshold energies for each atom species in BSi and NaBSi quenched glass structures using the Kieu2011_v1 and Kieu2011_v2 potentials, respectively. The procedure we used is: 1. For each atom species, we choose an atom from the quenched lattice at random and pick a random direction in 3D space. 2. We run an initial trial cascade by adding 100 eV of kinetic energy to the selected atom in the random direction we picked earlier.
Fig. 4. The plot shows the displacement energy thresholds for each atomic species (B, Si and O) for 3 quenched glass structures. The average threshold energy for each species for each glass is summarised in the above table.
Fig. 5. Histogram of transition barriers for (A) borosilicate glass modelled by the Kieu2011_v1 potential and (B) sodium borosilicate glass modelled by Kieu2011_v2 potential. The insets in each figure show a zoomed in section of the data in the range 0–2 eV.
144
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
However, here there is no equivalent concept. Because rebonding occurred after displacement, we determined the longer time scale transition in two ways. In the first way, under co-ordinated atoms were identified and their motion to a higher co-ordinated state was examined. Such motion generally resulted in large energy barriers of 3 eV or more. At 520 K, which is the initial storage temperature of a waste glass due to self-heating from b-decay [5], a 3 eV barrier is equivalent to one event every 356 million years. So we adopted a different approach. Here we use the data we obtained from threshold calculations above to identify initial and final states of transitions within the glass structures. We use the nudged elastic band method [15] to compute the minimum energy pathway between these states. This gives an estimate of the barrier height for these transitions. Fig. 5 shows histograms of barriers found for two glass structures. For borosilicate glass, the minimum barrier found was 0.9 eV (52 ls per event at 520 K). For sodium borosilicate glass, there were more small barriers, however, these only involved small displacements of single Na atoms. The minimum barrier found for Si movements in this case was 0.76 eV (2.3 ls per event), and B movements, 0.67 eV (0.3 ls per event). 4. Conclusion The most recent potentials of Kieu: Kieu2011_v1 and Kieu2011_v2, predict the correct bond lengths and bond angles and has the lowest melting point. We therefore conclude that the best existing inter-atomic potential for borosilicate glass is the Kieu2011_v1 potential. However, the minimum energy structure has a density that is too high at q=2.76 g/cm3. It is possible to quench a glass at the experimental density using the original authors method. However, in the thermal spike induced by a collision cascade, there will be local melting and therefore, a likely transition to a higher density in some regions and a corresponding lower density (voids) in others. Our original intention in performing this comparison of potentials was to test their suitability for radiation damage simulations, especially long time scale recovery of the material after collision cascades. The fact that they all give melting points that exceed the experimental value means that the transition barriers determined in the previous section are also likely to be a little too high compared to experimental values. This implies a need further to develop inter-atomic potentials for these glass materials, but to obtain the correct structure, density and melting point of such complex multicomponent materials is an extremely challenging task. It is easier for atoms to be displaced by a radiation event in the amorphous glass structures, compared to other ionic systems, such as, spinel [16], but, the resulting damaged system is quickly rebonded, and due to the high transition barriers, does not result in fast diffusion of atoms. This indicates the possibility of modelling multiple radiation events within a system, without the necessity to perform kinetic Monte Carlo simulations of diffusion and
recombination between radiation events, as is required for crystalline materials [17]. Acknowledgements We thank Dr. J.G. Shah of the Bhabia Atomic Research Centre (BARC) and Dr. Karl Travis of Sheffield University for useful discussions. Computing resources were provided by Loughborough University HPC service. The work was funded as part of a joint UK-India Nuclear Collaboration through EPSRC Grant No. EP/ K007882/1. References [1] J.-M. Delaye, D. Ghaleb, Molecular dynamics simulation of SiO2 + B2O3 + Na2O + ZrO2 glass, J. Non-Cryst. Solids 195 (3) (1996) 239–248, http://dx.doi.org/ 10.1016/0022-3093(95)00527-7. [2] L.-H. Kieu, J.-M. Delaye, L. Cormier, C. Stolz, Development of empirical potentials for sodium borosilicate glass systems, J. Non-Cryst. Solids 357 (18) (2011) 3313–3321, http://dx.doi.org/10.1016/j.jnoncrysol.2011.05.024. [3] M. Rushton, Simulations of glass and ceramic systems for nuclear waste applications (Ph.D. thesis), Imperial College, University of London, 2006. [4] J.A.C. Marples, The preparation, properties, and disposal of vitrified high level waste from nuclear fuel reprocessing, Glass Technol. 29 (6) (1988). [5] W.J. Weber, R.C. Ewing, Radiation effects in glasses used for immobilization of high-level waste and plutonium disposition, J. Mater. Res. 12 (8) (1997), http:// dx.doi.org/10.1557/JMR.1997.0266. [6] F.H. Stillinger, T.A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B 31 (8) (1985) 5262–5271, http://dx.doi.org/ 10.1103/PhysRevB.31.5262. [7] A.J. Connelly, K.P. Travis, R.J. Hand, N.C. Hyatt, E. Maddrell, Compositionstructure relationships in simplified nuclear waste glasses: 1. Mixed alkali borosilicate glasses, J. Am. Ceram. Soc. 94 (1) (2011) 151–159, http:// dx.doi.org/10.1111/j.1551-2916.2010.04066.x. [8] B. Guillot, N. Sator, A computer simulation study of natural silicate melts. Part I: low pressure properties, Geochim. Cosmochim. Acta 71 (5) (2007) 1249– 1265, http://dx.doi.org/10.1016/j.gca.2006.11.015. [9] M.J.D. Rushton, R.W. Grimes, S.L. Owens, Predicted changes to alkali concentration adjacent to glasscrystal interfaces, J. Am. Ceram. Soc. 91 (5) (2008) 1659–1664, http://dx.doi.org/10.1111/j.1551-2916.2008.02352.x. [10] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping and Range of Ions in Matter, vol. 1, Pergamon, New York, 1985. [11] C. Maynell, G. Saunders, S. Scholes, Ultrasound propagation in glasses in the metastable immiscibility region of the sodium borosilicate system, J. NonCryst. Solids 12 (3) (1973) 271–294, http://dx.doi.org/10.1016/00223093(73)90001-X. [12] J.D. Gale, A.L. Rohl, The general utility lattice program (GULP), Mol. Simul. 29 (5) (2003) 291–341, http://dx.doi.org/10.1080/0892702031000104887. [13] R.L. Hervig, A. Navrotsky, Thermochemistry of sodium borosilicate glasses, J. Am. Ceram. Soc. 68 (6) (1985) 314–319, http://dx.doi.org/10.1111/j.11512916.1985.tb15232.x. [14] L. Kittiratanawasin, R. Smith, B. Uberuaga, K. Sickafus, Displacement threshold and Frenkel pair formation energy in ionic systems, Nucl. Instrum. Methods Phys. Res. Sect. B 268 (19) (2010) 2901–2906, http://dx.doi.org/10.1016/ j.nimb.2010.04.024. [15] G. Henkelman, B.P. Uberuaga, H. Jonsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (22) (2000) 9901, http://dx.doi.org/10.1063/1.1329672. [16] R. Smith, D. Bacorisen, B.P. Uberuaga, K.E. Sickafus, J.a. Ball, R.W. Grimes, Dynamical simulations of radiation damage in magnesium aluminate spinel, MgAl2O4, J. Phys. Condens. Matter 17 (6) (2005) 875–891, http://dx.doi.org/ 10.1088/0953-8984/17/6/008. [17] C. Scott, R. Smith, Modelling the sputtering of Au surfaces using a multi timescale technique, Proc. R. Soc. A: Math. Phys. Eng. Sci. 469 (2150) (2012) 20120480, http://dx.doi.org/10.1098/rspa.2012.0480.
Contents lists available at ScienceDirect
Nuclear Instruments and Methods in Physics Research B journal homepage: www.elsevier.com/locate/nimb
Borosilicate glass potentials for radiation damage simulations Kenny Jolley a,⇑, Roger Smith a, Kitheri Joseph b a b
Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, UK Materials Chemistry Division, Indira Gandhi Centre for Atomic Research, Kalpakkam 603 102, India
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 25 June 2014 Received in revised form 8 December 2014 Accepted 9 December 2014 Available online 26 December 2014
Three borosilicate glass (SiO2–B2O3) fixed charge potentials from the literature are compared (Delaye and Ghaleb, 1996; Kieu et al., 2011; Rushton, 2006) and their suitability for use in simulations of radiation damage is assessed. For a range of densities, we generate glass structures by quenching at 5 1012 K/s using constant volume Molecular Dynamics. In each case, the bond lengths, mean bond angles, bulk modulus, melting point and displacement energy thresholds are calculated, and where possible compared to experimental data. Whereas the bond lengths and mean bond angles are reasonably well predicted, we find that the potentials predict melting temperatures, bulk moduli and densities that are higher than experimental data. The displacement energy thresholds are generally lower than those for ionic crystalline materials, but show a wider spread of values. However, the barriers for atomic rearrangements, after atoms have been displaced in the equilibrium structures, are very high. This indicates, that the radiation damage produced in the ballistic phase of a collision cascade, is likely to persist for extended time scales. This is in contrast to crystals, where interstitials and vacancies can diffuse rapidly between successive radiation events. Ó 2014 Elsevier B.V. All rights reserved.
Keywords: Molecular Dynamics Borosilicate glass Inter-atomic potentials Radiation damage
1. Introduction High-Level nuclear Wastes (HLW) from spent nuclear power station fuels and decommissioned nuclear weapons must be safely removed from the environment. One leading method of achieving this, is long term storage by immobilisation in glass waste forms prior to permanent disposal in a geologically stable repository [4,5]. A primary concern with this method, is that over long time-scales, the radioactive material could leach out of the glass and contaminate the environment. It is for this reason we must understand how glass responds to radiation damage over long time scales. In the early stages of waste encapsulation the material is heated due to b-decay, however, in this work it is the damage caused by a-decay that will be the main focus of the study. Experiments on physical samples are limited to relatively short times (10 years), compared to geological time scales. Experiments are also carried out at higher dose rates, such as by using ion beams to simulate the collision cascade process. As a result, atomistic simulation methods can play an important role in the understanding of how these glasses behave after a radiation event. Although Molecular Dynamics (MD) simulations can only access very short ⇑ Corresponding author. E-mail addresses: [email protected] (R. Smith), [email protected] (K. Joseph).
(K.
http://dx.doi.org/10.1016/j.nimb.2014.12.024 0168-583X/Ó 2014 Elsevier B.V. All rights reserved.
Jolley),
[email protected]
times ( 1 ns), kinetic Monte Carlo (kMC) simulation techniques can be used to access longer time scales. In order to simulate the radiation event, it is necessary to have a good model of the amorphous glass structure. Amorphous materials are more difficult to handle computationally than crystals, since there is no equivalent concept of a vacancy or interstitial. In addition, although the structures have short range order, there are many different atomic configurations that are possible, so statistical averaging over different structures is very important. In addition, good inter-atomic potentials to describe these systems is essential. In this work, we compare three inter-atomic potentials from the literature, and test their suitability for Molecular Dynamics (MD) simulations of radiation effects in borosilicate glass. We investigated two potentials from Delaye’s group (Delaye1996 [1] and Kieu2011 [2]), and one potential from Rushton (Rushton2006 [3]). Only the potentials of Kieu were developed specifically for the simulation of large sets of compositions around the SiO2–B2O3–Na2O glass. The Delaye1996 potential was developed for simulating R7T7 glass, and the Rushton2006 potential was developed for Magnox glass. The ‘best’ potential model for simulating borosilicate glass was then used to determine threshold displacement energies, and some typical energy barriers for atomic movement after such atomic displacements.
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
For the boron–oxygen interaction, Aij ¼ 180390:53 eV, qij ¼ 0:124 Å and C ij ¼ 35:0019 eV Å6.
2. Methodology Computer simulations of borosilicate glass and sodium borosilicate glass were performed using MD. For borosilicate glass, the structures consisted of 70 mol% SiO2 and 30 mol% B2O3. We refer to this glass as, BSi, in this paper. For sodium borosilicate glass, we used the SBN12 blend which is described in Ref. [2]. This glass consists of 59.66 mol% SiO2, 28.14 mol% B2O3 and 12.20 mol% Na2O. We refer to this glass as, NaBSi, in this paper. A brief overview of each inter-atomic potential that we have investigated, is given in Sections 2.1–2.3 below. These potentials are suitable for near equilibrium conditions and must be joined to a screened Coulomb potential for close particle separation. Our splining method is described in Section 2.4, and the method for generating glass structures is described in Section 2.5. 2.1. Delaye1996 The Delaye1996 potential [1] is used to model BSi glass. This potential contains both 2 and 3 body terms. The 2-body interaction is modelled by the Born–Mayer–Huggins potential. The potential energy, Uðrij Þ, of two ions separated by r ij is given by Eq. (1) below:
Uðrij Þ ¼ Aij exp
r ij
! þ
qij
1 4p0
zi zj r ij
2.2.2. Kieu2011_v2 For the sodium borosilicate glass that we studied (NaBSi), R = 0.4335 and K = 2.12. This is described as the SBN12 blend in the original paper. The ion charges were: zSi ¼ 1:869128; zO ¼ 0:965872; zB ¼ 1:528756 and zNa ¼ 0:451628. The parameter, Aij , for the boron–oxygen bond, is 205751.1464 eV. The remaining parameters (qij and C ij ) are the same as before. 2.3. Rushton2006 The Rushton2006 potential [3,9] is used for modelling BSi glass, and also contains only 2-body potentials. In this case, a LennardJones pair potential plus the Coulomb term (Eq. 3) is used to model the atomic interactions.
pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi Bi Bj Ai Aj 1 zi zj Uðr ij Þ ¼ 12 6 þ 4p0 r ij r ij r ij
ð3Þ
The parameters A and B for each ion are given in the original paper. zi and zj are the full formal charges for each ion. 2.4. Application to collision cascades
ð1Þ
where Aij , and qij are parameters of the model and can be found in the original paper [1]. zi and zj are the ion charges. This potential uses full formal charges for each ion. A Stillinger–Weber [6] 3-body term is applied to the O–Si–O and Si–O–Si bond angles. We choose not to put any restriction on the O–B–O bond angle as was done by Connelly [7]. The adjustable parameters are listed in the original paper [1] in Table 4. Due to the inclusion of the 3-body terms, this potential is the most computationally intensive. The three body term also induces an extremely large penalty when the bond angles move away from their equilibrium values. 2.2. Kieu2011
Some of these potentials tend to 1 at short range due to the electrostatic coulomb term (zi zj =rij ) and the Buckingham term (C ij =r 6ij ). Therefore, the potential energy reaches a maximum value at short range, which can be as low as 150 eV (for a Si–O bond). Hence, during a radiation damage cascade, it becomes possible for some atoms to unnaturally fuse together. The short range interactions between atomic nuclei are more accurately modelled by the ZBL potential [10]. We therefore join these potentials to the ZBL potential at short range with a spline function given by Eq. (4).
FðrÞ ¼ exp a0 þ a1 r þ a2 r 2 þ a3 r3 þ a4 r 4 þ a5 r5
ð4Þ
The constants a0 – a5 are set such that the potential energy and its first and second derivatives are continuous. 2.5. Modelling of glass formation
The original paper by Kieu et al. [2], describes potentials for modelling both borosilicate (BSi) and sodium borosilicate (NaBSi) glass. This potential contains only 2-body terms. A Buckingham pair potential (given by Eq. (2) below), is used to model these interactions.
Uðrij Þ ¼ Aij exp
141
r ij
qij
! þ
1 zi zj C ij 6 4p0 r ij r ij
ð2Þ
The adjustable parameters Aij ; qij and C ij for each ion interaction are listed in the original paper [2]. The boron–oxygen interaction is a special case, where the parameter, Aij , depends upon the composition of the glass. The ion charges (zi and zj ), also depend upon the glass composition. This composition is characterised by the two molar ratios, R = [Na2O]/[B2O3] and K = [SiO2]/[B2O3]. Since the parameter values and charges depend upon the composition of the glass, we denote the potential used for BSi glass as Kieu2011_v1, and potential used for NaBSi glass as Kieu2011_v2.
The glass structures are amorphous and have no long range order, therefore, we cannot simply generate a structure as we would for a crystal lattice. The procedure for generating a glass structure is to perform an MD simulation of a quench from a molten state. Atoms are placed at random positions inside a cubic box, ensuring no two atoms are closer than 1 Å. Such a configuration has a large potential energy, and therefore a high initial temperature when simulated with MD, typically exceeding 10,000 K. At every time-step the temperature is controlled by rescaling the velocities of all the atoms. If the temperature exceeds the desired value by 7%, then the velocities are rescaled to the desired temperature. We equilibrate at 6000 K for 10 ps then quench at 5 1012 K/s to 50 K. The final minimisation to 0 K is obtained via conjugate gradient minimisation. A schematic of this process is shown in Fig. 1. 3. Results 3.1. Density
2.2.1. Kieu2011_v1 For borosilicate glass (BSi), R = 0 and K = 2.33, and therefore, the charges obtained are close to the values used in the Guillot–Sator potential model [8]: zSi ¼ 1:89; zO ¼ 0:945 and zB ¼ 1:4175. We used these charges in the Kieu2011_v1 potential.
For each inter-atomic potential we performed a set of quenches using constant volume MD. We vary the volume to obtain sets of quenched glass structures at a range of different densities. Fig. 2 shows the plots of the average potential energy per atom as a
142
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
The original authors generated their glass structures by quenching at a density 5–10% below the experimental value, and then relaxed the system at constant pressure. This enabled them to find a stable local minimum close to the experimental density. In our simulations, we used GULP [12] to perform constant pressure relaxations of the glass structures for a range of densities. The glass structure with the optimum density did not change in volume significantly. The low density structures, increased in density on relaxation, but, still had a lattice energy that was higher than the minimum found at the optimum density. Typically, relaxations at constant pressure using GULP caused the density to change by less than 10%. Low density structures relaxed by increasing their density, while the high density structures relaxed by lowering their density. Constant pressure minimisation did not yield a better optimal structure, therefore, we believe a large sample of constant volume quenches is sufficient to obtain the optimal structure. The value of the average potential energy per atom of the optimum glass structures also varies between each potential.
Fig. 1. A schematic of the quenching procedure.
3.2. Structural properties
Fig. 2. Plots of the average total potential energy per atom of quenched lattices vs density for each inter-atomic potential. Each curve has been offset on the vertical axis so that they appear within one plot.
function of glass density for each inter-atomic potential that we studied. Each plot has been offset on the vertical axis such that all curves appear within the same figure. In each case, the data about the minimum has been fitted with a 2nd order polynomial, from which we estimate the global minimum lattice energy. We find that borosilicate glass modelled with the Delaye1996 potential, has a global minimum at a density of 1.78 g/cm3. This is lower than the experimental value of 2.04 g/cm3 [11]. All the other potentials have global minima that are higher than the experimental values. The Kieu2011_v1 and Rushton2006 potentials have minima at 2.76 and 3.67 g/cm3, respectively, for borosilicate glass. The global minimum lattice energy found for NaBSi glass, modelled with the Kieu2011_v2 potential was, 2.728 g/cm3. The Rushton2006 potential produced a glass structure that contained many over co-ordinated silicon and boron atoms. This is the reason for the significantly higher optimum density found. The higher densities that we found in our constant volume simulations using the Kieu2011 potentials, could be a high density phase, which, during an actual experiment would not be realised. This is because the cooling process takes place over longer timescales at constant pressure, and could require large energy barriers to be overcome.
We measured some structural properties of each glass. This data is summarised in Table 1. The Delaye1996 and Kieu2011 potentials give bond length and bond angle distributions that are in good agreement with experimental data. The Rushton2006 potential, however, contains Si–O and B–O bond lengths that are too long. The bond angle distributions in this case also differ due to the over-coordinated structure. The Bulk moduli, K B , of the optimum quenched glass structures for each potential is found using GULP [12], and are reported in Table 1. We see that in all cases the bulk moduli for borosilicate glass exceeds the experimental value of 23.7 GPa [2]. The sodium borosilicate glass (modelled using Kieu2011_v2) structure also has a bulk modulus that exceeds the experimental value of 42 GPa. We found that the bulk modulus decreases with decreasing density, and that when quenching a glass using the original authors method, we were able to replicate the results in the original paper, which were much closer to the experimental values. 3.3. Melting points We also measured the melting temperature predicted by each inter-atomic potential. The method for estimating the melting point is to look at the peak of a plot of the specific heat capacity as a function of temperature. Fig. 3 shows this for each interatomic potential that we investigated. To generate this plot: 1. For a range of temperatures (0–14,000 K), using the optimum quenched lattice for each potential, we thermalise the system using the Berendsen thermostat for 10 ps. 2. Then we run a constant energy MD (NVE) simulation for each case, and compute the time average of the temperature and total lattice energy. 3. The specific heat capacity is then the gradient of a plot of the total lattice energy vs temperature.
Table 1 Summary of simulated glass structural properties for each inter-atomic potential. Glass
BSi BSi NaBSi BSi
Potential
Delaye1996 Kieu2011_v1 Kieu2011_v2 Rushton2006
Density (g/cm3)
1.78 2.76 2.728 3.67
K B (GPa)
90 79 79 –
Bond length (Å)
Bond angle (°)
Si–O
B–O
O–Si–O
O–B–O
1.595 1.607 1.608 1.93
1.355 1.365 1.428 1.468
109.43 109.36 109.34 100.89
119.26 119.10 112.93 103.10
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
143
3. After 1 ps we compare the initial and final lattices to see if the atom has become permanently displaced. The cut-off for displacement is 1 Å. 4. If it has not been displaced, we double the trial energy and rerun the cascade. 5. We find the minimum threshold energy required to permanently displace the atom using the binary search algorithm. We compute the threshold to an accuracy of 1 eV in each case.
Fig. 3. A Plot of the temperature dependence of the specific heat capacity. The peak occurs at 8000 K for Delaye1996, 6000 K for Rushton2006 and 2700 K for the Kieu2011 potentials.
4. To obtain a smooth curve, we use a least squares fit on several adjacent points to obtain a moving average of the gradient. We see that the Delaye1996 potential predicts the highest melting point at approximately 8000 K. The Rushton2006 potential predicts the melting point at about 6000 K, and both the Kieu2011_v1 and Kieu2011_v2 potentials give 2700 K. The experimental value of the melting temperature of sodium borosilicate glass ranges from 1310 to 1870 K depending upon the composition of the glass [11,13]. 3.4. Displacement energy thresholds
The distributions of threshold energies are shown in Fig. 4. Here we plot data for borosilicate glass at two densities (2.0 and 2.76 g/ cm3) modelled with the Kieu2011_v1 potential and sodium borosilicate glass (at the optimum density 2.74 g/cm3) modelled by Kieu2011_v2 potential. The average threshold energy does not have any directional dependence upon cascade direction, since the glass structures are amorphous. We also find that the lower density borosilicate glass has lower threshold energies for all atomic species, compared to the higher density case. Compared to ionic crystalline systems, the thresholds are lower [14] with a more continuous spectrum. The structures were also examined after a displacement event. It was found that there were no debonded atoms or voids formed in the system, a result confirmed also by running cascades up to 1 keV. Defective regions resulting from these low energy displacements, were best identified through the co-ordination number of each atom in the system. 3.5. Transition barriers In an ionic crystal one can determine the activation energies of defects such as vacancies and interstitials in a straightforward way.
We calculate the displacement threshold energies for each atom species in BSi and NaBSi quenched glass structures using the Kieu2011_v1 and Kieu2011_v2 potentials, respectively. The procedure we used is: 1. For each atom species, we choose an atom from the quenched lattice at random and pick a random direction in 3D space. 2. We run an initial trial cascade by adding 100 eV of kinetic energy to the selected atom in the random direction we picked earlier.
Fig. 4. The plot shows the displacement energy thresholds for each atomic species (B, Si and O) for 3 quenched glass structures. The average threshold energy for each species for each glass is summarised in the above table.
Fig. 5. Histogram of transition barriers for (A) borosilicate glass modelled by the Kieu2011_v1 potential and (B) sodium borosilicate glass modelled by Kieu2011_v2 potential. The insets in each figure show a zoomed in section of the data in the range 0–2 eV.
144
K. Jolley et al. / Nuclear Instruments and Methods in Physics Research B 352 (2015) 140–144
However, here there is no equivalent concept. Because rebonding occurred after displacement, we determined the longer time scale transition in two ways. In the first way, under co-ordinated atoms were identified and their motion to a higher co-ordinated state was examined. Such motion generally resulted in large energy barriers of 3 eV or more. At 520 K, which is the initial storage temperature of a waste glass due to self-heating from b-decay [5], a 3 eV barrier is equivalent to one event every 356 million years. So we adopted a different approach. Here we use the data we obtained from threshold calculations above to identify initial and final states of transitions within the glass structures. We use the nudged elastic band method [15] to compute the minimum energy pathway between these states. This gives an estimate of the barrier height for these transitions. Fig. 5 shows histograms of barriers found for two glass structures. For borosilicate glass, the minimum barrier found was 0.9 eV (52 ls per event at 520 K). For sodium borosilicate glass, there were more small barriers, however, these only involved small displacements of single Na atoms. The minimum barrier found for Si movements in this case was 0.76 eV (2.3 ls per event), and B movements, 0.67 eV (0.3 ls per event). 4. Conclusion The most recent potentials of Kieu: Kieu2011_v1 and Kieu2011_v2, predict the correct bond lengths and bond angles and has the lowest melting point. We therefore conclude that the best existing inter-atomic potential for borosilicate glass is the Kieu2011_v1 potential. However, the minimum energy structure has a density that is too high at q=2.76 g/cm3. It is possible to quench a glass at the experimental density using the original authors method. However, in the thermal spike induced by a collision cascade, there will be local melting and therefore, a likely transition to a higher density in some regions and a corresponding lower density (voids) in others. Our original intention in performing this comparison of potentials was to test their suitability for radiation damage simulations, especially long time scale recovery of the material after collision cascades. The fact that they all give melting points that exceed the experimental value means that the transition barriers determined in the previous section are also likely to be a little too high compared to experimental values. This implies a need further to develop inter-atomic potentials for these glass materials, but to obtain the correct structure, density and melting point of such complex multicomponent materials is an extremely challenging task. It is easier for atoms to be displaced by a radiation event in the amorphous glass structures, compared to other ionic systems, such as, spinel [16], but, the resulting damaged system is quickly rebonded, and due to the high transition barriers, does not result in fast diffusion of atoms. This indicates the possibility of modelling multiple radiation events within a system, without the necessity to perform kinetic Monte Carlo simulations of diffusion and
recombination between radiation events, as is required for crystalline materials [17]. Acknowledgements We thank Dr. J.G. Shah of the Bhabia Atomic Research Centre (BARC) and Dr. Karl Travis of Sheffield University for useful discussions. Computing resources were provided by Loughborough University HPC service. The work was funded as part of a joint UK-India Nuclear Collaboration through EPSRC Grant No. EP/ K007882/1. References [1] J.-M. Delaye, D. Ghaleb, Molecular dynamics simulation of SiO2 + B2O3 + Na2O + ZrO2 glass, J. Non-Cryst. Solids 195 (3) (1996) 239–248, http://dx.doi.org/ 10.1016/0022-3093(95)00527-7. [2] L.-H. Kieu, J.-M. Delaye, L. Cormier, C. Stolz, Development of empirical potentials for sodium borosilicate glass systems, J. Non-Cryst. Solids 357 (18) (2011) 3313–3321, http://dx.doi.org/10.1016/j.jnoncrysol.2011.05.024. [3] M. Rushton, Simulations of glass and ceramic systems for nuclear waste applications (Ph.D. thesis), Imperial College, University of London, 2006. [4] J.A.C. Marples, The preparation, properties, and disposal of vitrified high level waste from nuclear fuel reprocessing, Glass Technol. 29 (6) (1988). [5] W.J. Weber, R.C. Ewing, Radiation effects in glasses used for immobilization of high-level waste and plutonium disposition, J. Mater. Res. 12 (8) (1997), http:// dx.doi.org/10.1557/JMR.1997.0266. [6] F.H. Stillinger, T.A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B 31 (8) (1985) 5262–5271, http://dx.doi.org/ 10.1103/PhysRevB.31.5262. [7] A.J. Connelly, K.P. Travis, R.J. Hand, N.C. Hyatt, E. Maddrell, Compositionstructure relationships in simplified nuclear waste glasses: 1. Mixed alkali borosilicate glasses, J. Am. Ceram. Soc. 94 (1) (2011) 151–159, http:// dx.doi.org/10.1111/j.1551-2916.2010.04066.x. [8] B. Guillot, N. Sator, A computer simulation study of natural silicate melts. Part I: low pressure properties, Geochim. Cosmochim. Acta 71 (5) (2007) 1249– 1265, http://dx.doi.org/10.1016/j.gca.2006.11.015. [9] M.J.D. Rushton, R.W. Grimes, S.L. Owens, Predicted changes to alkali concentration adjacent to glasscrystal interfaces, J. Am. Ceram. Soc. 91 (5) (2008) 1659–1664, http://dx.doi.org/10.1111/j.1551-2916.2008.02352.x. [10] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping and Range of Ions in Matter, vol. 1, Pergamon, New York, 1985. [11] C. Maynell, G. Saunders, S. Scholes, Ultrasound propagation in glasses in the metastable immiscibility region of the sodium borosilicate system, J. NonCryst. Solids 12 (3) (1973) 271–294, http://dx.doi.org/10.1016/00223093(73)90001-X. [12] J.D. Gale, A.L. Rohl, The general utility lattice program (GULP), Mol. Simul. 29 (5) (2003) 291–341, http://dx.doi.org/10.1080/0892702031000104887. [13] R.L. Hervig, A. Navrotsky, Thermochemistry of sodium borosilicate glasses, J. Am. Ceram. Soc. 68 (6) (1985) 314–319, http://dx.doi.org/10.1111/j.11512916.1985.tb15232.x. [14] L. Kittiratanawasin, R. Smith, B. Uberuaga, K. Sickafus, Displacement threshold and Frenkel pair formation energy in ionic systems, Nucl. Instrum. Methods Phys. Res. Sect. B 268 (19) (2010) 2901–2906, http://dx.doi.org/10.1016/ j.nimb.2010.04.024. [15] G. Henkelman, B.P. Uberuaga, H. Jonsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (22) (2000) 9901, http://dx.doi.org/10.1063/1.1329672. [16] R. Smith, D. Bacorisen, B.P. Uberuaga, K.E. Sickafus, J.a. Ball, R.W. Grimes, Dynamical simulations of radiation damage in magnesium aluminate spinel, MgAl2O4, J. Phys. Condens. Matter 17 (6) (2005) 875–891, http://dx.doi.org/ 10.1088/0953-8984/17/6/008. [17] C. Scott, R. Smith, Modelling the sputtering of Au surfaces using a multi timescale technique, Proc. R. Soc. A: Math. Phys. Eng. Sci. 469 (2150) (2012) 20120480, http://dx.doi.org/10.1098/rspa.2012.0480.