Iron phosphate glasses: structure determination and displacement energy thresholds, using a fixed charge potential model

Abstract Iron phosphate glass is a versatile matrix for the immobilisation of various radioactive elements found in high-level nuclear waste (HLW). Quenched glass structures of iron phosphate glasses with Fe/P ratios of 0.33, 0.67 and 0.75 and with a composition of 40 mol% Fe 2 O 3 and 60 mol% P 2 O 5 , with 4% and 17% Fe 2 + ion concentrations were generated using molecular dynamics and the threshold displacement energies calculated. In the minimum energy structures, we found that in nearly all cases the P atoms were 4-fold coordinated. The potential energy per atom increased with increasing concentration of Fe 2 + ions with similar Fe/P ratio, suggesting that decreasing the Fe 2 + content is a stabilising factor. The average bond distances between Fe 2 + –O, Fe 3 + –O, P–O and O–O were calculated as 2.12, 1.88, 1.5 and 2.5 A respectively. The threshold displacement energy (E d ) was found to be dependent upon the ion specie, less for Fe 2 + ions compared to Fe 3 + ions, and was overall slightly lower than that determined for borosilicate glass.


Introduction
Phosphate glasses, due to their favourable properties such as: reasonably low liquid and glass transition temperatures, low viscosity, high thermal expansion coefficient, high electrical conductivity and high ultraviolet transmission, find application in a wide range of fields. For example, phosphate glasses are used in lasers [1], solid electrolytes [2], bio-medical devices [3] and nuclear waste immobilisation [4]. A good review up to the year 2000 is given in [5]. According to the literature [6], in spite of the good glass forming characteristics of phosphate glasses, their relatively poor chemical durability limits their application, especially in the field of nuclear waste immobilisation. However, a new group of phosphate glasses, iron phosphate glass, is being considered as a promising matrix for the immobilisation of high-level waste rich in alkali oxides, sulphates and chrome oxides [7][8][9][10]. Iron phosphate glass and its waste form containing simulated fast reactor waste were synthesised, characterised and reported by us earlier [11,12]. Higher waste loading, better chemical durability and better corrosion resistance [13,14] are certain promising features of iron phosphate glass compared to other phosphate glasses. Among the various compositions of iron phosphate glass, the one with 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 (referred as IPG in the present paper) is found to be more chemically durable [15][16][17]. It also has the ability to accommodate large amounts of certain nuclear wastes, especially those that are not well suited for borosilicate glasses [16]. The better chemical durability of iron phosphate glass is attributed to the presence of more hydration resistant Fe-O-P bonds compared to P-O-P bonds available in other phosphate glasses [16]. The synthesis, characterisation and experimental determination to elucidate the structure of a wide variety of iron phosphate glasses are available in the literature [7,[17][18][19][20]. However, the available literature on the structural modelling of iron phosphate glasses is limited [21,22]. This is because, the structure of iron phosphate glasses, not only depends on composition, quenching temperature of the melt and quenching technique [23], it depends also on the concentration of Fe 2+ and Fe 3+ ions in the glass.
The iron phosphate glass with the composition 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 , has been reported by various researchers to have different ratios of Fe 2+ /Fe [7,12,24,25]. The concentration of Fe 2+ /Fe in glass with the same atomic composition varies from 4 to 20%, while the density varies between 2.9 and 3.0 g cm −3 with the uncertainty in density ranging from ±0.005 to 0.02 g cm −3 [7,24]. The concentrations of Fe 2 + /Fe and density of IPG reported by us [12] were 4% and 2.9 g cm − 3 respectively. The promising composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 (IPG) is found with the varying density depending on the concentration of Fe 2+ in the glass. It is evident from the literature [7,24,25] that as the concentration of Fe 2+ in IPG increases, the density also increases for the same atomic composition. It has been shown experimentally by Mössbauer spectroscopy [24], that as the melting temperature of IPG increased from 1423 to 1673 K, the density of the glass increased along with the increase in concentration of Fe 2+ . Thus, it becomes essential to model the structure of iron phosphate glasses specific to the composition and Fe 2+ /Fe ratio.
It is well known that amorphous structures are more difficult to handle computationally than crystalline lattices since many different atomic configurations are possible. Therefore averaging over many different structures is very important. It is also essential to have good inter-atomic potentials that describe the amorphous systems. A recent paper [21] reported one of the first studies of iron phosphate glasses using fixed charge potentials.
Other authors have also developed potentials of a form that could be used to model phosphate glasses. However potentials that involve shell models such as that by Ainsworth et al. [26] are not really suitable for radiation cascade studies and a previous potential formulation by Pedone et al. [27], which also included the capability to model Fe-P-O systems had not been explicitly tested on these glasses. Since a key aim of our research programme is to investigate radiation effects in phosphate glasses, the potential developed in [21] was chosen as the underlying model for the work.
Thus we modify the potential in [21] so that it is suitable for radiation studies, compare the results with the previous work, and then test other compositional structures with various Fe/P ratio and Fe 2+ / Fe ratios that we have produced experimentally. Furthermore, we determine the threshold displacement energies as the first step towards the investigation of collision cascades in these systems.

Methodology
Molecular dynamics simulation studies were carried out to model the structure of iron phosphate glasses. The code (LBOMD) has been used over many years for radiation damage studies, see e.g. [28] for work involving radiation damage in spinels. The two body interactions, i.e. Fe-O, P-O and O-O were modelled using a Buckingham rigid ion potential (Eq. (1)), together with a Coulomb term to model the long-range interactions between ionic charges.
The subscripts, i and j, refer to each ion; q i and q j are the ion charges, r ij is the inter-ionic distance and ε 0 is the permittivity of free space. The parameters: A ij , ρ ij and C ij for each ionic bond are given in reference [21].
In addition to the two-body terms, three body terms were also used to control the local bond angles. This three-body potential is especially very important for P-O bonds due to their ionic-covalent nature similar to that of Si-O bonds found in silicate glasses [29].
The original authors [21] used a harmonic three-body potential, V(θ ijk ) = 1/2 k(θ ijk − θ 0 ) 2 , where k = 3.5 eV and θ 0 = 109.47°for O-P-O, and k = 3.0 eV and θ 0 = 135.5°for P-O-P. This function presents some computational difficulties for radiation damage studies, as it requires the computation of the derivative of arccos(θ) to obtain the forces. This derivative is infinite when θ = 180°. Furthermore, there is no smooth cut-off to zero as the atomic separation increases.
Among the various analytical forms of three-body potentials listed in the literature [30][31][32], the three-body Stillinger-Weber potential [31] (Eq. (2)) was chosen instead to model the three body terms. This is more suitable to use in radiation damage simulations and the construction of a potential that can be used in such simulations is also a key aim of the work.
For an atom triplet (j-i-k), r ij and r ik , are the two internal atomic separations, and θ jik , is the bond angle at the central atom 'i'. θ 0 represents the angle towards which the angle θ jik is constrained and r c is the cut-off radius beyond which the three-body terms do not apply. λ, γ and θ 0 are adjustable parameters.
The Stillinger-Weber potential was chosen because the cos(θ) term is obtainable directly from the atomic separations via the dot product (also, the derivatives are a function of cos(θ)). The Stillinger-Weber potential has continuous derivatives and a built-in smooth cut-off to zero at r = r c . We fit the Stillinger-Weber potential to the original authors' function by computing the parameter λ. We do this by computing the Taylor series expansion of Eq. (2) at θ = θ 0 to the second order. For this calculation we set γ = 0.5 eV and r ij = r ik = 1.5 Å (the equilibrium bond length).
For both triplets (P-O-P and O-P-O), we set this equal to the original authors' function and rearrange to find λ. These parameter values (λ, γ and θ 0 ) for both triplets are summarised in Table 1.
For radiation damage studies, the two-body potential cannot model the repulsion between nuclei when the inter-particle separation is small. A screened Coulomb potential is normally used to model such interactions and the ZBL model [33] is the commonly used model when the separation is small. As a result we have joined the two-body potential given in Eq. (1) to the ZBL potential using a splining function. The details and parameters for the splining function are given in Table 2.
Each of the glass structures was prepared by distributing the required number of atoms of each species randomly within a cubic box. We then use a simple temperature-rescaling algorithm [34] to quench the system. This works by measuring the temperature at each timestep. If the temperature exceeds the desired value by 7%, then, the velocities of all the atoms are rescaled, such that, the temperature is the desired value. We quench the system from 6000 K to 10 K at a rate of Table 1 Three-body potential terms used in the present study. Table 2 The splining parameters for joining the Buckingham + Coulomb potential to the ZBL potential. The splining function is of the form exp(f 0 + f 1 x + f 2 x 2 + f 3 x 3 + f 4 x 4 + f 5 x 5 ) joined to the ZBL potential at x = a and to the Buckingham + Coulomb at x = b so that the function and its first 2 derivatives are continuous. The units of f i are Å −i . In the case of Fe 3+ -O and O-P, there were offsets of 13 and 50 eV respectively added to the potential (and later subtracted) to make the potential positive at x = b. 5 × 10 12 Ks −1 using this temperature-rescaling algorithm with constant volume molecular dynamics. The total simulation time required was 1.2 ns. Finally, the structures were minimised to 0 K using the conjugate gradient method. This cooling rate is slower than that given in [21] by a factor of 2. This was found to be necessary to generate the lowest energy structures. Molecular dynamics computation was carried out on these iron phosphate glasses contained in a cubic box with a small (~150) number of atoms with a varying Fe/P atomic ratio. The total number of atoms varied depending on the Fe/P atomic ratio. The potential energy of iron phosphate glasses with a varying Fe/P atomic ratio was analysed as a function of density to obtain the structure with the lowest potential energy, and hence the optimal density for each case. We then computed quenched structures for each system with a larger number of atoms (1100-1550) at the optimal density. Having determined the optimised structure of IPG, the mechanical properties (bulk modulus, shear modulus and Young's modulus) were calculated using GULP [30] and compared with that of experimental values available in the literature for IPG with a composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 .
The threshold displacement energies of each species: iron (Fe 2+ and Fe 3+ ), phosphorous and oxygen were computed. Since only iron ions exist in ferrous and ferric states, the oxidation state is specified as Fe 2+ and Fe 3+ and the ionic nature of ions (P 5+ and O 2− ) with single oxidation state are not specified throughout the paper. We consider each atomic species, oxidation state and local coordination separately for the calculation of threshold energies.
To find the threshold energies for each case, firstly, using the quenched lattices, a list of candidate atoms to test is generated. We then randomly select a target atom from this list and pick a random direction. To find the threshold energy, we begin by performing a trial cascade by injecting 100 eV of kinetic energy into the target atom. After 2 ps of MD simulation time we compare the initial and final lattices. If the target atom has been displaced, we know that the threshold energy is lower than 100 eV, otherwise, it must be higher. We then perform additional cascades to zero-in on the threshold energy using the binary search algorithm. The cut-off tolerance used was 0.5 eV, meaning that the results obtained are within 0.5 eV of the actual threshold energy. We then select a new random target atom from our list, and then pick a new random direction; we then find the threshold energy using the method above. We repeat this process at least 900 times to obtain good statistics.

Results and discussion
The details of the composition of the smaller systems of iron phosphate glass, Fe/P atomic ratio and the number of atoms quenched are given in Table 3. Fig. 1 shows the average potential energy obtained as a function of glass density for systems Fe 3 (P 2 O 7 ) 2 , Fe(PO 3 ) 3 and Fe 4 (P 2 O 7 ) 3 . The composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 is the same as that of the Fe 4 (P 2 O 7 ) 3 glass and can also be referred as IPG. However, Fe 4 (P 2 O 7 ) 3 does not have any Fe 2+ ions in the glass. For the small systems, the lowest energy was obtained with densities of 3.04 and 3.2 g cm −3 for Fe 4 (P 2 O 7 ) 3 and Fe 3 (P 2 O 7 ) 2 glasses respectively. The lowest energy structure of Fe 3 (P 2 O 7 ) 2 glass (containing both Fe 2+ and Fe 3+ ions) had a higher density compared to the other systems containing only Fe 3+ ions. Larger systems of Fe 3 (P 2 O 7 ) 2 and Fe 4 (P 2 O 7 ) 3 were generated with the density that showed the lowest energy in the smaller atom systems. The computation was done a minimum of three times for each large system at the fixed density. Whenever the difference in energy between the structures was small (b1%), the structure with the lowest potential energy was chosen for further analysis based on certain logical guidelines.
The logical guidelines used to choose the best structure when the difference in potential energy was small are given below: 1. The coordination of P should be 4 for a maximum number of P ions in the glass. For Fe(PO 3 ) 3 glass (130 atoms), the difference in potential energy between the structures with densities of 2.9 and 3.04 g cm −3 was very small, i.e., 3.5 × 10 −3 %. Thus larger systems were modelled using both the densities (2.9 and 3.04 g cm −3 ) for Fe(PO 3 ) 3 glass (1300 atoms). For Fe(PO 3 ) 3 glass (1300 atoms), based on the guidelines above, the structure with a density of 3.04 g cm − 3 was chosen for Table 3 The composition (in atomic %), of each atomic species, in the three small system glass structures studied.    3 and Fe(PO 3 ) 3 , are presented in Table 4. Table 4 shows that all the phosphorous ions have 4-fold coordination in Fe 3 (P 2 O 7 ) 2 and Fe 4 (P 2 O 7 ) 3 glass. However, two five-fold coordinated phosphorous ions were observed in the Fe(PO 3 ) 3 glass. This is significantly lower than the 1.9% that has been reported previously [21]. Fe 2+ is present only in Fe 3 (P 2 O 7 ) 2 glass and it is found with coordination numbers 4 to 6. Generally, the Fe 3+ ion is found with 4 to 6 coordination in Fe 3 (P 2 O 7 ) 2 , Fe 4 (P 2 O 7 ) 3 and Fe(PO 3 ) 3 glasses. The oxygen ion coordination was found to vary between 1 and 3 in Fe 3 (P 2 O 7 ) 2 , Fe 4 (P 2 O 7 ) 3 and Fe(PO 3 ) 3 glasses. Fe 3 (P 2 O 7 ) 2 glass also showed 0.5% of four coordinated oxygen ions. The snapshots, shown in Fig. 2, show the coordination of ions within the Fe 3 (P 2 O 7 ) 2 glass structure.
In addition to these three iron phosphate glasses of varying Fe/P atomic ratio, it is useful to model the structure of the chemically durable iron phosphate glass with a composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 (IPG) [15]. Before describing the results of the simulations, it is useful to understand the experimental observations on crystallisation of iron phosphate glasses. Iron phosphate glass crystallises to FePO 4 , Fe 3 (P 2 O 7 ) 2 , Fe(PO 3 ) 3 , Fe(PO 3 ) 2 , Fe 7 (PO 4 ) 6 , and Fe 2 P 2 O 7 depending on the composition, atmosphere of crystallisation and ratio of Fe 2+ /Fe [24,25,[35][36][37][38][39]. It was reported by us [40] that, on heating, IPG containing 4% Fe 2+ , crystallises to a mixture of Fe 3 (P 2 O 7 ) 2 , Fe(PO 3 ) 3 and Fe 4 (P 2 O 7 ) 3 under flowing argon atmosphere. C. S. Ray et al. [37] reported a fraction of 17% of Fe 2+ ions within the IPG glass with each phase also of the form Fe 3 (P 2 O 7 ) 2 , Fe(PO 3 ) 3 and Fe 4 (P 2 O 7 ) 3 . Thus these are the data that are presented in Table 5.
The structure of IPG was computed using the combination of atoms of Fe, P and O to provide concentrations of (a) 4% Fe 2+ and (b) 17% Fe 2+ . The atomic percentage of each ion used for constructing the structure of IPG involving 4% Fe 2+ and 17% Fe 2+ , is also presented in Table 5. The Fe/P atomic ratio of the final glass (Table 5), differs slightly from the starting batch composition of Fe/P atomic ratio = 0.67 in both the cases. However, it is well within the reported experimental uncertainty [37,41]. Computation was carried out to obtain structures with the lowest potential energy using densities 2.9 and 3.04 g cm −3 for IPG containing 4% Fe 2+ ions. These densities are chosen based on the experimental result [11] and the value obtained from ref. [21]. For IPG containing 17% Fe 2+ , the densities of 2.9, 2.95 and 3.04 g cm −3 were used to compute the structure of glass to obtain the glass with the lowest energy. A density of 2.95 g cm −3 was additionally used in IPG containing 17% Fe 2+ ions since the experimental density reported by C. S. Ray et al. [37] was 2.95 g cm −3 .
As explained earlier, the lowest energy structure of IPGs containing 4% Fe 2+ and 17% Fe 2+ was computed using similar MD methods. The potential energy obtained for densities of 2.9 and 3.04 g cm −3 is presented in Table 6. The contribution of various ions and their coordination numbers is also given in Table 6 for both the densities (2.9 and 3.04 g cm −3 ). As per the guidelines, the structure with the lowest number of Fe 2+ ions with 3-fold coordination has a density of 2.9 g cm −3 . In the structure with a density of 3.04 g cm −3 , there are a greater number of 3-fold coordinated Fe 2+ ions. This reduces the average coordination number of Fe 2+ below 4.0 and is therefore not considered for computation of the threshold displacement energy. Under a similar analogy, the structure of IPG with 17% Fe 2+ was the structure with a density of 2.95 g cm − 3 and that the optimised structure was used for further threshold energy displacement calculations. The details of potential energy, coordination number of various ions in IPG with 17% Fe 2+ , are also presented in Table 6. The computed density of the structures containing 4% and 17% Fe 2+ is in good agreement with the reported experimental values [11,37].
The IPG containing 4% Fe 2+ showed Fe 2+ and Fe 3+ with coordination number 3, whereas, the IPG containing 17% Fe 2+ showed no 3 coordinated iron ions. However, the IPG containing 17% Fe 2+ was found to contain one five-fold coordinated phosphorous ion.
In order to compare the structures, the potential energy per atom was determined and this is presented in Fig. 3. The potential energy per atom increases linearly with the Fe/P ratio. This plot also clearly shows that as the concentration of Fe 2+ ions increases in the glass, the potential energy per atom also increases.
In the experiment, the starting material contains only Fe 3+ ions but during the preparation process, some Fe 2+ ions form. Ammonium dihydrogen phosphate is used as the source of phosphorous during the preparation of the iron phosphate glass. Evolution of ammonia, a reducing agent, during the preparation of the iron phosphate glass (due to the decomposition of ammonium di-hydrogen phosphate), reduces some Fe 3+ to Fe 2+ . However, the reducing effect can be minimised by a low temperature pre-calcination, as indicated in ref. [40]. Hence the most stable structure would be expected to be the one with the lowest concentration of Fe 2+ ions, that is the IPG containing 4% Fe 2+ concentration.
Modelling studies on the structure of iron phosphate glasses are given in refs. [21,22]. During the structural evolution of the IPG, Stoch et al. reported [22] the formation of crystallised products in air as Eq. (4) also justifies the formation of Fe 2+ ions (Fe 2 P 2 O 7 ) from the starting composition containing only Fe 3+ ions with the loss of oxygen. Thus some Fe 2+ are to be expected despite the fact that the model gives a higher cohesive energy when the amount of Fe 2+ is minimised.
The radial distribution functions of Fe 3 (P 2 O 7 ) 2 glass and IPG containing 4% Fe 2+ concentration are shown in Fig. 4. The bond distance between Fe 3+ -O in all the five iron phosphate glasses did not vary significantly. This agrees with the bond distance determined in [21]. As expected, the Fe 2+ -O distance was longer than that of Fe 3+ -O. These values are in good agreement with the experimental values reported in the literature [19,23,[42][43][44] for iron phosphate glasses. Table 7 presents details of the mean bond distance between Fe 2+ -O, Fe 3+ -O and P-O in the iron phosphate glasses that we have modelled.
The calculated elastic properties (bulk modulus and shear modulus) of IPG containing 4% and 17% concentrations of Fe 2+ ions are presented in Table 8. The values of these properties are obtained as the average of the three methods employed in GULP [30]. Experimentally determined bulk modulus and Young's modulus of IPG with a composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 were 47 GPa and 70-72 GPa respectively as reported in the literature [45,46]. The calculated properties are in good agreement with the experimentally reported values.
The threshold displacement energy (E d ) of various glasses Fe 3 (P 2 O 7 ) 2 (1260 atoms), Fe 4 (P 2 O 7 ) 3 (1550 atoms) and Fe(PO 3 ) 3 (1300 atoms) was computed. E d was found to change with the coordination number and species of the particular ion in Fe 3 (P 2 O 7 ) 2 glass. The effect of the coordination of Fe 2+ ions on E d is shown in Fig. 5. As would be expected, the values of E d are generally lower for the lower coordination numbers. This was found to be true for other two systems (Fe 4 (P 2 O 7 ) 3 and Fe(PO 3 ) 3 ) also. E d also varied as a function of the coordination number and nature of     Fig. 6, Fe 2+ has a lower E d than Fe 3+ . Since phosphorous is 4-fold coordinated in all the iron phosphate glasses, the values of E d for phosphorous ions in all these glasses are compared and shown in Fig. 7. It is apparent from the plot, that E d remains similar for phosphorous ions in all the iron phosphate glasses.
The peak values of E d of the different ions as a function of their coordination number are presented in Table 9. Although a large number of calculations were performed to determine these values, there is a wide statistical variation, and the values for a given species are very similar for the different coordination numbers.
The values of E d for the O ions given in Table 9 are very similar to those calculated for the borosilicate glasses [34] but the thresholds for P and Fe 2+ are lower than those calculated for boron and silicon suggesting that the phosphate glasses might not be as radiation resistant as the borosilicate.

Conclusion
The structures of iron phosphate glasses with varying Fe/P atomic ratios were modelled as a function of density using the molecular dynamics simulation method. Furthermore,the structures of iron phosphate glass with a composition of 40 mol% Fe 2 O 3 -60 mol% P 2 O 5 (IPG) containing the concentrations of 4 and 17% Fe 2+ ions, were also computed and compared in order to understand the effect of Fe 2+ concentration. The potential energy per atom of all these glasses indicated that the system energy increases with both increasing Fe/P atomic ratio and Fe 2+ concentration. The computed elastic properties of   40 mol% Fe 2 O 3 -60 mol% P 2 O 5 glass were found to be similar to that of the experimental values reported in the literature. The threshold displacement energy varies with the nature and as well as with the coordination of ions in these iron phosphate glasses. However, similar threshold energies were obtained for IPG containing 4 and 17% Fe 2+ ions. The threshold energy of phosphorous ions remains the same for all the iron phosphate glasses. The results show that the peak displacement energy threshold of the Fe 2+ ions is lower than that for Fe 3+ and O ions. The cohesive energy of the glasses is also reduced as the amount of Fe 2+ increases. This suggests that for the purposes of nuclear waste immobilisation, it would be preferable to produce the glasses with as low a Fe 2+ content as possible. This is not the only consideration since the glasses must also be suitable to incorporate heavy atoms such as caesium and the associated structure is also important. However, the bond lengths and densities have been shown to vary only slightly as a function of composition.
Future work will investigate the structural changes induced by collision cascades in the glasses. It was for this purpose that the modified potential was developed since it provides a better description of the collision dynamics than the simple Teter model developed only for structural studies.