Development of a New Sr-O Parameterization to Describe the Influence of SrO on Iron-Phosphate Glass Structural Properties Using Molecular Dynamics Simulations

Iron-phosphate glasses, due to their properties, have many potential applications. One of the most promising seems to be nuclear waste immobilization. Radioactive 90Sr isotope is the main short-lived product of fission and, due to its high solubility, it can enter groundwater and pose a threat to the environment. On the other hand, Sr is an important element in hard tissue metabolic processes, and phosphate glasses containing Sr are considered bioactive. This study investigated the effect of SrO addition on a glass structure of nominal 30Fe2O3-70P2O5 chemical composition using classical molecular dynamics simulations. To describe the interaction between Sr-O ion pairs, new interatomic potential parameters of the Buckingham-type were developed and tested for crystalline compounds. The short-range structure of the simulated glasses is presented and is in agreement with previous experimental and theoretical studies. The simulations showed that an increase in SrO content in the glass led to phosphate network depolymerization. Analysis demonstrated that the non-network oxygen did not take part in the phosphate network depolymerization. Furthermore, strontium aggregation in the glass structure was observed to lead to the non-homogeneity of the glass network. It was demonstrated that Sr ions prefer to locate near to Fe(II), which may induce crystallization of strontium phosphates with divalent iron.


Introduction
Phosphate glasses are an important group of materials that possess many interesting features. Their properties mean that they can be applied in many different fields. Exemplary glasses containing calcium are biocompatible and can be used as bone or dental implants [1][2][3]. The addition of rare earths makes them interesting in the production of optoelectronic devices [4][5][6]. Phosphate glasses are also used as eco-fertilizers with a controlled dissolution rate [7][8][9]. However, the existence of easily hydrated P-O-P bridges leads to their corrosion, induced by water originating, for example, in a humid environment [10,11]. This effect may limit such glasses' potential applications. The problem can be solved by the addition of Al 2 O 3 or Fe 2 O 3 , which very strongly increase the glasses' chemical durability. As such, a superior water-resistant material can be obtained and the glasses can be considered as matrixes in the waste immobilization process [12,13]. One of the most promising materials for this purpose is glass of the composition 60P 2 O 5 -40Fe 2 O 3 [12][13][14]. Thus, by changing the iron-phosphate glass (IPG) composition, one can adjust the water resistance of the final material depending on the applicational purpose and obtain, e.g., a material of controlled dissolution rate.
In the case of radioactive waste immobilization, one of the major sources of nuclear waste activity is 90 Sr isotope, which is the source of β radiation. Due to its high solubility, it can enter the groundwater from waste. Such 90 Sr can be incorporated into bones and thus remain in the body. [12,13]. On the other hand, the similarity of strontium to calcium makes

Simulation Methods
This study examined glasses with the formal composition (100-x)(0.3Fe 2 O 3 -0.7P 2 O 5 )-xSrO mol. %, where x = 0, 5, 10, 15, . . . , 50. During the synthesis of iron-phosphate glasses, some of the iron atoms are reduced from Fe(III) to Fe(II). Since the share of Fe(II) atoms is usually in the range of 20-40% of the total iron in simulations, we used a constant value for the share equal to 30% [28,30]. It should be pointed out that the iron oxidation state depends on the glass synthesis conditions (melting temperature, atmosphere, etc.). According to this fact, glass composition may influence the glass melting conditions, and thus the redox state may change. It has been observed that an increase in modifier content the IPG glasses can lead to a decrease in the Fe(II) share [31,32]. Thus, we can assume that the Fe(II) content will decrease with the increase of the SrO in the glass.
The simulations were conducted using classical molecular dynamics with Lammps software [33]. In the simulation, approximately 50,000 ions were placed randomly in a cubic simulation box. The number of ions was set to reflect the tested glass chemical composition. The box lengths were set according to the experimental glass density. To remove the surface effects, periodic boundary conditions were applied.
The pair interactions of i-j ions were described using the Buckingham potential of the equation: where → r ij is the distance between ions i-j, ε 0 is the permittivity of free space, q i q j is the product of effective charges, and A ij , ρ ij , and C ij are parameters, given in Table 1, that were taken from the literature [19,34], except for the Sr-O potential parameters, which were fitted as described below. A short-range cut-off of 14 Å was set for the calculations of the non-Coulomb part of the potential. The Coulombic long-range forces were calculated by using the particle-mesh Ewald method [35] with a precision of 10 −6 . Table 1. Values of the Buckingham potential parameters of the i-j pair [19,34].
Pair (i-j) q i (e) A ij (eV) q ij (Å) C ij (eV·Å 6  To properly describe the covalent character of P-O-P and O-P-O bonds, an additional three-body term was added to the potential energy in the form proposed by Stillinger and Weber. The term parameters were taken directly from [21,36]. The simulations were conducted with the 1 fs timestep in the NVT ensemble. At the first stage, the starting temperature of 10,000 K was gradually lowered to 1600 K with relaxation steps at 5000 K and 1600 K. After 100,000 timesteps, the temperature was reduced to 300 K, at which point the system was relaxed through an additional 100,000 timesteps. A detailed description of the simulation protocol can be found in [20,21].
As an example, the glass network fragments were visualized using VESTA 3 software for three-dimensional visualization of the crystal, volumetric, and morphology data [37].

Development of Sr-O Interatomic Potential Parameters
According to our knowledge, there are no studies in the literature that deal with the Sr-O interaction parameters for the partial charges used in the simulations. That is why we had to determine the A Sr-O , ρ Sr-O , and C Sr-O parameters of the interatomic potential. The parameters were fitted to reflect as closely as possible the crystal structure parameters of several crystalline compounds containing Sr through the use of the GULP program [38]. The set of compounds consisted of simple oxide (SrO), strontium phosphates (Sr 2 P 2 O 7 , Sr(PO 3 ) 2 ), and iron phosphates (SrFe(III) 3 (PO 4 ) 3 O, SrFe(III)Fe(II) 2 (PO 4 ) 3 , SrFe(III) 2 (P 2 O 7 ) 2 ). The crystalline phases were previously detected as the crystallizing compounds in the strontium-containing iron-phosphate glasses [39]. The compounds' crystal structure parameters were taken from the Crystallography Open Database and the appropriate COD card numbers are given in Table 2. All the structures were fitted simultaneously and the relaxed fitting procedure was employed [38]. The determined parameters are collected in Table 1, whereas the reference and the simulated crystal structure parameters are summarized in Table 2.
It should be pointed out here that the fitted parameters could reproduce the crystal structure parameters of the compounds very well. The differences between the experimental and the simulated results were about a few percent or less. The Sr-O potential parameters were fitted based on the crystal structures, which were detected as crystallizing compounds in iron-phosphate glasses of similar compositions [39], and similar compounds were expected during the crystallization of the studied system. Taking into account the fact that the glass network structure should somehow reflect the possible crystallizing compounds, we assumed that the parameters would also properly describe the structural features of the simulated systems.

Glass Synthesis
An important point of the glass network simulations was the application of proper glass density, which in the proposed simulation protocol was maintained during the simu-lations. In the study, we used an experimental glass density. Therefore, at the beginning, we had to synthesize the simulated glasses. The glass synthesis was conducted with a conventional melting and quenching technique. The appropriate amounts of chemicalgrade purity NH 4 H 2 PO 4 , Fe 2 O 3 , and SrCO 3 were carefully homogenized in the planetary ball mill. To compensate for P 2 O 5 losses during melting, approximately 20% overweight NH 4 H 2 PO 4 was used. The batches were melted in an electric laboratory furnace in Al 2 O 3 crucibles in an air atmosphere. The melting temperature was 1473 K and glasses were kept at this temperature for 2 h. The melts were vitrified by casting onto steel plates. The synthesis procedure is known to give glass compositions close to the limit of about one percent, and Al 2 O 3 contamination is below 1% with an Fe(II) content of about 30% of the total iron. More details can be found in [22,28,40,41].
The glasses density was measured using the Archimedes method and deionized water was used as the immersed liquid. The density was measured for five pieces of glass from every composition and mean values and mean standard deviations were obtained.

Results and Discussion
The experimentally determined glass densities and corresponding molar volumes are presented in Figure 1. The molar volume was calculated according to the formula V = M/ρ, where ρ is the glass density and M is the molar mass of the glass.
into account the fact that the glass network structure should somehow reflect the possible crystallizing compounds, we assumed that the parameters would also properly describe the structural features of the simulated systems.

Glass Synthesis
An important point of the glass network simulations was the application of proper glass density, which in the proposed simulation protocol was maintained during the simulations. In the study, we used an experimental glass density. Therefore, at the beginning, we had to synthesize the simulated glasses. The glass synthesis was conducted with a conventional melting and quenching technique. The appropriate amounts of chemical-grade purity NH4H2PO4, Fe2O3, and SrCO3 were carefully homogenized in the planetary ball mill. To compensate for P2O5 losses during melting, approximately 20% overweight NH4H2PO4 was used. The batches were melted in an electric laboratory furnace in Al2O3 crucibles in an air atmosphere. The melting temperature was 1473 K and glasses were kept at this temperature for 2 h. The melts were vitrified by casting onto steel plates. The synthesis procedure is known to give glass compositions close to the limit of about one percent, and Al2O3 contamination is below 1% with an Fe(II) content of about 30% of the total iron. More details can be found in [22,28,40,41].
The glasses density was measured using the Archimedes method and deionized water was used as the immersed liquid. The density was measured for five pieces of glass from every composition and mean values and mean standard deviations were obtained.

Results and Discussion
The experimentally determined glass densities and corresponding molar volumes are presented in Figure 1. The molar volume was calculated according to the formula V = M/ρ, where ρ is the glass density and M is the molar mass of the glass. The glass density increased according to the formula ρ(x) = 1.5·10 −4 x 2 + 0.0856x + 2.7472 (g/cm 3 ). More information concerning the glass network changes can help determine the molar volume. The molar volume decreased according to the relation V(x) The glass density increased according to the formula ρ(x) = 1.5·10 −4 x 2 + 0.0856x + 2.7472 (g/cm 3 ). More information concerning the glass network changes can help determine the molar volume. The molar volume decreased according to the relation V(x) = −0.0012x 2 −0.2747x + 51.0002 (cm 3 /mol). The decrease indicated that the glass network had become more compacted. That may have been related to the placing of the Sr atoms in the open voids of the glass network; in this way, the network free space would have decreased. This strongly suggests the modifier role of strontium in the glass network.
The short-range order in the glasses can be described by the partial distribution function (PDF). Exemplary = −0.0012x 2 −0.2747x + 51.0002 (cm 3 /mol). The decrease indicated that the glass network had become more compacted. That may have been related to the placing of the Sr atoms in the open voids of the glass network; in this way, the network free space would have decreased. This strongly suggests the modifier role of strontium in the glass network.
The short-range order in the glasses can be described by the partial distribution function (PDF). Exemplary PDF curves for P-O, Fe(III)-O, Fe(II)-O, and Sr-O pairs for x = 5 glass are presented in Figure 2a. As would be expected for amorphous systems, the curves were characterized by a single well-developed peak which as the effect of the short-range order in the closest coordination shell. The further peaks were considerably less intense, which proved the glassy character of the simulated systems. For the rest of the tested materials, the obtained curves were very similar. It should be noted that the first PDF peak for the P-O pair was composed of two peaks (as shown in Figure 2a), which were attributed to two main P-O bond lengths, depending on the role of oxygen in the [PO4] tetrahedron. The first peak of the lower P-O length was attributed to shorter P-Ot or P-ONB bonds with the maximum at ca. 1.475 Å. This corresponds fairly well with the literature data, where the bond length is in the range of 1.43-1.49 Å [20][21][22][41][42][43][44][45]. The distance increased linearly with the SrO content, as presented in Figure 3 (P-ONB curve). This increase may be attributed to the increased content of P-ONB bonds compared to P-Ot bonds. P-Ot bonds are shorter than P-ONB bonds [22,41,45]. Thus, with the increase of the P-ONB content, the position of the PDF curve component moves toward the higher lengths. The effect may also indicate that Sr ions prefer to break double P-Ot bonds and transform them into P-ONB-Sr joinings, instead of disrupting P-OB-P bridges. As would be expected for amorphous systems, the curves were characterized by a single well-developed peak which as the effect of the short-range order in the closest coordination shell. The further peaks were considerably less intense, which proved the glassy character of the simulated systems. For the rest of the tested materials, the obtained curves were very similar. It should be noted that the first PDF peak for the P-O pair was composed of two peaks (as shown in Figure 2a The P-OB bond length was greater and the corresponding component of the first PDF peak was located for the longer distances (Figure 2.). The position of the peak was between 1.52 and 1.53 Å and was almost constant (Figure 3). The P-OB distances followed previous MD results indicating ca.  The P-O B bond length was greater and the corresponding component of the first PDF peak was located for the longer distances (Figure 2.). The position of the peak was between 1.52 and 1.53 Å and was almost constant (Figure 3). The P-O B distances followed previous MD results indicating ca. 1.51-1.64 Å [17,19,21,29,[41][42][43][45][46][47][48]. Nevertheless, it should be noted that the MD values are lower than the values from ab initio simulations, which are closer to 1.6 Å [22,29,41,45]. This may be an effect of the simulation model and the P-O interatomic potential parameters used.
The integration of the PDF curve gives the dependence of the coordination number of the given central atom on the distance. Therefore, the integration gives the average coordination number (CN) of P to oxygen as a function of the distance from the central P atom. The obtained CN functions are presented in Figure 2b. The CN P appeared as a step-like curve that is characteristic of a glass network species, and the P atoms took, in a step-like manner, a constant coordination number of 4 in relation to oxygen; further, the value was independent of the assumed cut-off radius in a wide region from above 1.6 Å up to ca. 2.6 Å. This confirmed that all the phosphorous atoms were placed in the middle of the oxygen tetrahedrons [PO 4 ]. Above 2.6 Å the second coordination sphere started.
The longer bond length created an Fe(III)-O pair. Therefore, the maximum of the PDF curve lay at the higher lengths ( Figure 2). The position of the maximum was about 1.9 Å and slightly increased with the SrO content, as presented in Figure 4. The value was in accordance with earlier studies of IPGs where the distance was in the range of 1.87-1.91 Å [19,21,26,29,45,49].  The CN curve for the Fe(III)-O (Figure 2b) also showed step-like behavior. Nevertheless, in this case, the step was not flat and increased with the distance. Therefore, the coordination number depended on the value of the cut-off radius taken. In the further analysis, the cut-off radius was defined as the distance which had the minimum value between the first and the second peaks on the corresponding PDF curve. For Fe(III)-O the minimum value was at 2.4 Å, and this value was used as the cut-off radius. Thus, the average coordination numbers for iron in relation to oxygen could be obtained and are shown in Figure 5. The CN curve for the Fe(III)-O (Figure 2b) also showed step-like behavior. Nevertheless, in this case, the step was not flat and increased with the distance. Therefore, the coordination number depended on the value of the cut-off radius taken. In the further analysis, the cut-off radius was defined as the distance which had the minimum value between the first and the second peaks on the corresponding PDF curve. For Fe(III)-O the minimum value was at 2.4 Å, and this value was used as the cut-off radius. Thus, the average coordination numbers for iron in relation to oxygen could be obtained and are shown in Figure 5.
It can be seen that the CN number for Fe(III) linearly increased with the SrO content in the glass and was about 4.5. This is in good agreement with earlier MD simulations [19][20][21]. In iron-phosphate glasses, iron atoms have a dual role depending on the crystal-chemical surroundings. They may substitute phosphorous atoms in their network positions in a tetrahedral coordination or may function as glass network modifiers with the higher coordination numbers. In the simulated glasses, the mean coordination of ca. 4.5 was realized by about 52% of the Fe(III) in tetrahedral positions, 45% in distorted bi-pyramidal sites (CN = 5), and only 3% in octahedral sites. Exemplary Fe(III) surroundings are presented in Figure 6.
Nevertheless, in this case, the step was not flat and increased with the distance. Therefore, the coordination number depended on the value of the cut-off radius taken. In the further analysis, the cut-off radius was defined as the distance which had the minimum value between the first and the second peaks on the corresponding PDF curve. For Fe(III)-O the minimum value was at 2.4 Å, and this value was used as the cut-off radius. Thus, the average coordination numbers for iron in relation to oxygen could be obtained and are shown in Figure 5. It can be seen that the CN number for Fe(III) linearly increased with the SrO content in the glass and was about 4.5. This is in good agreement with earlier MD simulations [19][20][21]. In iron-phosphate glasses, iron atoms have a dual role depending on the crystalchemical surroundings. They may substitute phosphorous atoms in their network positions in a tetrahedral coordination or may function as glass network modifiers with the higher coordination numbers. In the simulated glasses, the mean coordination of ca. 4.5 was realized by about 52% of the Fe(III) in tetrahedral positions, 45% in distorted bi-  Thus, we assumed that about half of the iron Fe(III) atoms occupied glass network positions. Nevertheless, it should be noted that this assumption was based on the value of the cut-off radius used and the geometry of the sites. A much more precise determination of the coordination number could be obtained based on bond order analysis, in which only interacting atoms are taken into account. It has been shown that in iron-aluminum-phosphate glasses, some oxygen atoms lying closer, e.g., as in aluminum, do not create a chemical bond, whereas those lying further away do create a bond [22,41]. Such an analysis is beyond the capacities of classical methods and can only be done by utilizing ab initio simulations.
Similar behavior to that of Fe(III) was evidenced for Fe(II). However, the CN function was less step-like, which was the reason for the lower short-range ordering. In this case, the distances were longer and the coordination number was also higher. This was a reason for the Fe(II) atoms' greater size and lower contribution to the directional, covalentcharacter interaction in the Fe(II)-O chemical bond compared to Fe(III)-O [22,29].
The longest distance was the Sr-O distance, which increased with SrO content in the glass (Figure 7a), and the mean coordination number for oxygen also increased from ca. 6.4 Å to ca. 7 Å (Figure 7b). It should be noted that in the fitted crystal structures strontium atoms were in coordination 6 or 7 in relation to oxygen, with the mean bond length being between ca. 2.56 and ca. 2.69 Å, which is in agreement with previous results. The Sr-O bond distance was comparable with but higher than that of Ca-O, which is about 2.43 Å in similar IPGs [20]. The difference seems to be reasonable taking into account the fact that strontium is larger than calcium. Thus, we assumed that about half of the iron Fe(III) atoms occupied glass network positions. Nevertheless, it should be noted that this assumption was based on the value of the cut-off radius used and the geometry of the sites. A much more precise determination of the coordination number could be obtained based on bond order analysis, in which only interacting atoms are taken into account. It has been shown that in iron-aluminum-phosphate glasses, some oxygen atoms lying closer, e.g., as in aluminum, do not create a chemical bond, whereas those lying further away do create a bond [22,41]. Such an analysis is beyond the capacities of classical methods and can only be done by utilizing ab initio simulations.
Similar behavior to that of Fe(III) was evidenced for Fe(II). However, the CN function was less step-like, which was the reason for the lower short-range ordering. In this case, the distances were longer and the coordination number was also higher. This was a reason for the Fe(II) atoms' greater size and lower contribution to the directional, covalent-character interaction in the Fe(II)-O chemical bond compared to Fe(III)-O [22,29].
The longest distance was the Sr-O distance, which increased with SrO content in the glass (Figure 7a), and the mean coordination number for oxygen also increased from ca. 6.4 Å to ca. 7 Å (Figure 7b). It should be noted that in the fitted crystal structures strontium atoms were in coordination 6 or 7 in relation to oxygen, with the mean bond length being between ca. 2.56 and ca. 2.69 Å, which is in agreement with previous results. The Sr-O Materials 2021, 14, 4326 9 of 16 bond distance was comparable with but higher than that of Ca-O, which is about 2.43 Å in similar IPGs [20]. The difference seems to be reasonable taking into account the fact that strontium is larger than calcium. We can make the general conclusion that the increase of SrO content increased the mean coordination numbers for oxygen for both forms of iron and strontium. The increase was related to the extension of the mean distances from oxygen. In all the cases, this was a linear growth that confirmed the correlation between the distance and the coordination number. Additionally, we may assume that the chemical bond lengths followed the trend. Similar behavior has been observed for calcium atoms in IPGs, although the calcium CN was higher by about 0.3 [20].
In the simulated glasses, P2O5 was the only pure network oxide. Therefore, in the further description of the glass network, we focus mostly on this cation. The important parameters describing the network are the O-P-O and P-O-P angle distribution functions (ADFs). The ADFs for the simulated glasses are presented in Figure 8.  We can make the general conclusion that the increase of SrO content increased the mean coordination numbers for oxygen for both forms of iron and strontium. The increase was related to the extension of the mean distances from oxygen. In all the cases, this was a linear growth that confirmed the correlation between the distance and the coordination number. Additionally, we may assume that the chemical bond lengths followed the trend. Similar behavior has been observed for calcium atoms in IPGs, although the calcium CN was higher by about 0.3 [20].
In the simulated glasses, P 2 O 5 was the only pure network oxide. Therefore, in the further description of the glass network, we focus mostly on this cation. The important parameters describing the network are the O-P-O and P-O-P angle distribution functions (ADFs). The ADFs for the simulated glasses are presented in Figure 8. We can make the general conclusion that the increase of SrO content increased the mean coordination numbers for oxygen for both forms of iron and strontium. The increase was related to the extension of the mean distances from oxygen. In all the cases, this was a linear growth that confirmed the correlation between the distance and the coordination number. Additionally, we may assume that the chemical bond lengths followed the trend. Similar behavior has been observed for calcium atoms in IPGs, although the calcium CN was higher by about 0.3 [20].
In the simulated glasses, P2O5 was the only pure network oxide. Therefore, in the further description of the glass network, we focus mostly on this cation. The important parameters describing the network are the O-P-O and P-O-P angle distribution functions (ADFs). The ADFs for the simulated glasses are presented in Figure 8.   The O-P-O distribution describes the bond angles in [PO 4 ] tetrahedrons. In an ideal tetrahedron, the distribution would have a maximum for the angle 109.47 0 , and a very similar maximum was observed in this study. The position of the maximum was independent of the glass composition. The intensity of the peak slightly increased with the SrO content and followed a very slight decrease of the half-width. This was an effect of the decrease in the number of O t oxygen atoms and their transformation to O nb . The difference in bond length between P-O NB and P-O B was lower compared to P-O t . Thus, the [PO 4 ] tetrahedron was less distorted and the corresponding ADF was narrower. The P-O-P angle is the angle at which [PO 4 ] tetrahedrons are joined by each other. In this case, the distribution was wider due to a much higher level of randomness. The position of the maximum linearly decreased with the SrO content in the glass (Figure 8-inset), which may have been an effect of the phosphate chains' shortening.
Glass network connectivity (NC) is another important parameter that can be used in the description of a glass network. The NC is defined as the number of bridging oxygen atoms per glass network-forming cation [24]. To initially determine the, one must establish the distribution of different oxygen types in the glass. In the studied system, there were four different oxygen atom types. There were terminal oxygen atoms (O T ), which create double P=O bonds; non-bridging atoms (O NB ), which are the atoms that create linkages for the P-O-Fe/Sr bonds; bridging oxygen atoms (O B ), which form P-O-P bridges; and non-network oxygen atoms (O NN ), which are not connected with phosphorous. Examples of atom clusters where the oxygen atom types are present are shown in Figure 9. The O-P-O distribution describes the bond angles in [PO4] tetrahedrons. In an ideal tetrahedron, the distribution would have a maximum for the angle 109.47 0 , and a very similar maximum was observed in this study. The position of the maximum was independent of the glass composition. The intensity of the peak slightly increased with the SrO content and followed a very slight decrease of the half-width. This was an effect of the decrease in the number of Ot oxygen atoms and their transformation to Onb. The difference in bond length between P-ONB and P-OB was lower compared to P-Ot. Thus, the [PO4] tetrahedron was less distorted and the corresponding ADF was narrower. The P-O-P angle is the angle at which [PO4] tetrahedrons are joined by each other. In this case, the distribution was wider due to a much higher level of randomness. The position of the maximum linearly decreased with the SrO content in the glass (Figure 8-inset), which may have been an effect of the phosphate chains' shortening.
Glass network connectivity (NC) is another important parameter that can be used in the description of a glass network. The NC is defined as the number of bridging oxygen atoms per glass network-forming cation [24]. To initially determine the, one must establish the distribution of different oxygen types in the glass. In the studied system, there were four different oxygen atom types. There were terminal oxygen atoms (OT), which create double P=O bonds; non-bridging atoms (ONB), which are the atoms that create linkages for the P-O-Fe/Sr bonds; bridging oxygen atoms (OB), which form P-O-P bridges; and non-network oxygen atoms (ONN), which are not connected with phosphorous. Examples of atom clusters where the oxygen atom types are present are shown in Figure 9. The distribution of the oxygen atom types as a function of the glass composition is shown in Figure 10a. The distribution of the oxygen atom types as a function of the glass composition is shown in Figure 10a. It can be seen that the number of OB oxygen atoms decreased with the SrO content, and the decrease accompanied an increase in the numbers of ONB and ONN oxygen atoms. The OB decrease together with the increase of ONB demonstrated the glass network's  (Figure 10b). Thus, the phosphate network was almost continuous only for the base glass where the NC was close to 1.
The distribution of different Q i structural units in the simulated glasses and the dimensionality of the glass network, defined as the mean i parameter of the Q i structural units, are given in Figure 11. It can be seen that the number of OB oxygen atoms decreased with the SrO content, and the decrease accompanied an increase in the numbers of ONB and ONN oxygen atoms. The OB decrease together with the increase of ONB demonstrated the glass network's gradual depolymerization due to the increase of the SrO content in the glass. This was confirmed by the decrease of the NC parameter, which followed the linear trend of OB changes. Based on the calculated OB values, the NC parameter decreased from 0.97 for x = 0 to 0.31 for x = 50 glass (Figure 10b). Thus, the phosphate network was almost continuous only for the base glass where the NC was close to 1.
The distribution of different Q i structural units in the simulated glasses and the dimensionality of the glass network, defined as the mean i parameter of the Q i structural units, are given in Figure 11. Exemplary phosphate network elements for the selected glasses are presented in Figure 12. Exemplary phosphate network elements for the selected glasses are presented in Figure 12. The base glass was built mostly of Q 2 with some comparable additions of Q 1 and Q 3 structural units. Thus, the phosphate network was composed of long chains that were joined to each other by Q 3 units and terminated by Q 1 units (Figure 12a). The addition of SrO led at first to a decrease in the number of Q 3 units. Therefore, the connected chains The base glass was built mostly of Q 2 with some comparable additions of Q 1 and Q 3 structural units. Thus, the phosphate network was composed of long chains that were joined to each other by Q 3 units and terminated by Q 1 units (Figure 12a). The addition of SrO led at first to a decrease in the number of Q 3 units. Therefore, the connected chains were being separated. Then, for the higher SrO content, the decrease in Q 2 units became more rapid and, in its place, formation of Q 1 and Q 0 occurred. Thus, in the middle range of SrO in the glass, the long chains were being shortened and even separated Q 0 units were distinguished (Figure 12b). The high modifier-content glass was built of very short chains, like Q 1 dimers and separated Q 0 units (Figure 12c). This followed the decrease in the glass network dimensionality (Figure 11b), which was the average <i> value of the Q i unit distribution [24,41]. Initially, the glass created a 2D phosphate network, in which the average dimensionality was decreased to one-dimensional chains for x = 40.
The glass network dimensionality could also be simply predicted based on the O/P ratio. In the case of the pure P 2 O 5, the O/P = 2.5 and <i> = 3. The parameter decreased linearly with the ratio through O/P = 3, <i> = 2 (metaphosphate) and O/P = 3.5, <i> = 1 (pyrophosphate) until O/P = 4, <i> = 0 (orthophosphate) [17]. According to the ratio, the base glass (x = 0) had O/P = 3.14, which resulted in <i> = 1.71. The value was lower compared to the simulated <i> of 2.02 (Figure 11b). For the glasses, we assumed that 30% of the Fe(III) atoms were reduced to Fe(II), which would have been connected with some oxygen loss. After considering this effect, the O/P ratio was judged to equal 3.08 and the corresponding <i> 1.84. Additionally, in the system, there were also O NN oxygen atoms that did not take part in the network depolymerization. Taking into account the O NN oxygen atoms led to an <i> value of 1.92, which was much closer to the simulated <i> value of 2.02. A similar result was found for all the simulated glasses. For example, for x = 25, considering only the iron reduction led to <i> = 1.37, whereas including O NN atoms led to <i> = 1.48. The value was close to the simulated value of <i> = 1.45. Thus, comparing the experimentally obtained dimensionality with the value predicted based on the glass stoichiometry provided the possibility of estimating the quantity of O NN oxygen, as proposed in [28,41].
Interesting information was gained from the analysis of the average coordination number of P for Fe(III) and Fe(II). The dependence on the SrO content in the glass is presented in Figure 13. The number of Fe(III) around P initially increased and achieved a maximum for about x = 0.25 before decreasing. The curve for Fe(II) demonstrates the opposite tendency. It suggests that, for low Sr contents, Fe(II) atoms could be removed from the phosphorous surroundings. This shows that strontium may prefer to be located near Fe(II) sites in a The number of Fe(III) around P initially increased and achieved a maximum for about x = 0.25 before decreasing. The curve for Fe(II) demonstrates the opposite tendency. It suggests that, for low Sr contents, Fe(II) atoms could be removed from the phosphorous surroundings. This shows that strontium may prefer to be located near Fe(II) sites in a network. On the other hand, the withdrawing of Fe(II) from P was compensated for by the affinity of Fe(III). This was confirmed by analyzing coordination numbers for Sr for metal cations. The compositional dependence is shown in Figure 14. The number of Fe(III) around P initially increased and achieved a maximum for about x = 0.25 before decreasing. The curve for Fe(II) demonstrates the opposite tendency. It suggests that, for low Sr contents, Fe(II) atoms could be removed from the phosphorous surroundings. This shows that strontium may prefer to be located near Fe(II) sites in a network. On the other hand, the withdrawing of Fe(II) from P was compensated for by the affinity of Fe(III). This was confirmed by analyzing coordination numbers for Sr for metal cations. The compositional dependence is shown in Figure 14. It can be clearly seen that the coordination numbers for all the cations increased with increasing SrO content. The fastest increase was for Sr around P, which showed a preference for strontium to be located around [PO4] tetrahedra. This is quite natural and has been It can be clearly seen that the coordination numbers for all the cations increased with increasing SrO content. The fastest increase was for Sr around P, which showed a preference for strontium to be located around [PO 4 ] tetrahedra. This is quite natural and has been previously detected for sodium and calcium [20,21]. Despite the fact that Fe(II) only comprised 30% of the total iron in the glass structure, strontium had a higher preference for locations near Fe(II) than Fe (III).
This supports the idea of Sr's affinity for Fe(II). It should be pointed out that such behavior may lead to a preference for the crystallization of compounds containing Sr and Fe(II). This was observed experimentally, for which in the comparable Sr-containing IPGs the major crystalline compound was SrFe(II) 2 P 2 O 7 [39].
The calculated values of the coordination numbers of Sr for Sr were used to obtain the aggregation parameter (R Sr-Sr ), as defined in [50][51][52]. The parameter was a simple ratio of the coordination number to a random homogenous distribution. A value over 1 indicated ions clustering. The estimated values of the aggregation parameter are shown in Figure 15.
The values of the parameter show that, with low Sr content, very strong aggregation of the atoms could be observed. This led to the formation of Sr-rich clusters and the inhomogeneity of the glass structure. The increase of SrO content led to an increase of the regions and, in this way, the whole glass network became more homogenous. A similar effect has been observed for Ca [20]. However, in the case of Ca, the parameter value was lower and, for the low modifier, the contents were about twice the Ca value. Thus, we judged that the aggregation of Sr was stronger than that of Ca.
In light of the obtained results, it can be observed that the glass network was inhomogeneous, with regions enriched and depleted in specific species alongside more ordered regions, like in Lebedev crystalline theory, which interlaces with more disordered regions of Zachariasen-Waren theory. Thus, the network was close to the model domain proposed by Görlich [53,54].
Fe(II). This was observed experimentally, for which in the comparable Sr-containing IPGs the major crystalline compound was SrFe(II)2P2O7 [39].
The calculated values of the coordination numbers of Sr for Sr were used to obtain the aggregation parameter (RSr-Sr), as defined in [50][51][52]. The parameter was a simple ratio of the coordination number to a random homogenous distribution. A value over 1 indicated ions clustering. The estimated values of the aggregation parameter are shown in Figure 15. The values of the parameter show that, with low Sr content, very strong aggregation of the atoms could be observed. This led to the formation of Sr-rich clusters and the inhomogeneity of the glass structure. The increase of SrO content led to an increase of the regions and, in this way, the whole glass network became more homogenous. A similar effect has been observed for Ca [20]. However, in the case of Ca, the parameter value was lower and, for the low modifier, the contents were about twice the Ca value. Thus, we judged that the aggregation of Sr was stronger than that of Ca.
In light of the obtained results, it can be observed that the glass network was inhomogeneous, with regions enriched and depleted in specific species alongside more ordered regions, like in Lebedev crystalline theory, which interlaces with more disordered regions of Zachariasen-Waren theory. Thus, the network was close to the model domain proposed by Görlich [53,54].

Conclusions
New interatomic potential parameters for Sr-O interaction were developed and tested for crystalline compounds.

Conclusions
New interatomic potential parameters for Sr-O interaction were developed and tested for crystalline compounds.
The influence of increasing SrO content on the structural properties of 30Fe 2 O 3 -70P 2 O 5 glasses was tested using classical molecular dynamics. The short-range structure results were in good agreement with previous experimental and theoretical studies.
The glasses with low SrO contents were built of chains that can have branches and form three-dimensional structures. The increase in SrO content in the glass led to phosphate network depolymerization, which was characterized by decreases in glass connectivity and dimensionality parameters.
It was observed that non-network oxygen atoms existed in the studied system that did not take part in the phosphate network depolymerization. The quantity of the atoms can be estimated through a comparison of the experimental dimensionality parameter with the value predicted based on the glass stoichiometry.
It was detected that Sr ions prefer to be located near to Fe(II), which can induce crystallization of strontium phosphates with divalent iron.
Evidence of strontium aggregation was found, and it was the most intense for low SrO contents. This can lead to the inhomogeneity of the glass network with the formation of Sr-rich and Sr-depleted regions.