Study of the structure of MgSiO3 system under compression by using ring statistics and voronoi analysis

The structural change of MgSiO3 liquid under compression is still one of the most interesting challenges. In this paper, we perform the molecular dynamics simulation to study the structural change of MgSiO3 liquid from 0 to 200 GPa. Ring statistics are analyzed to clarify the intermediate-range order, to explain why the second peak of Si–Si PRDFs splits into 2 subpeaks at 200 GPa, and to show the heterogeneity of MgSiO3. Large rings which form at high pressures would capture the oxygen atoms. Oxygen atoms which have negative charge attract Mg2+ ions, creating magnesium-rich regions. Besides, the Voronoi and Q n distribution changes on the ring with pressure are clarified to give more information about the rings.


Introduction
The magnesium silicate system is one of the main components of the Earth's mantle, chondrite meteorites and interplanetary dust [1][2][3]. It also plays an important application in technological and scientific fields such as materials science or geoscience. Therefore, it is necessary to study the structural properties of this material at different physical conditions. The continuous silicate network is mainly constructed by SiO 4 units through corner linkages. Adding MgO leads to the reduction of the polymerization of -Si-O-network and the formation of non-bridging oxygens (NBOs) and free oxygens (FOs). In MgSiO 3 , Mg can act as both network former and network modifier [4][5][6]. Previous works state that there are structural changes in the MgO-SiO 2 system by compression [7][8][9], by the molar ratio of MgO/SiO 2 content [10] and by temperature changing [11][12][13].
The radial distribution function and coordination number are among typical criteria for evaluating the local structure of materials. The average Si-O bond length in both MgSiO 3 glass and liquid is about 1.66 Å and does not depend much on temperature [14,15]. Under mild compression, the mean bond lengths of Si-O and Mg-O increase slightly [13,16]. The first peak of Si-O pair radial distribution functions (PRDF) shifts from 1.64 Å at 0 GPa to 1.68 Å at 30 GPa and the one of Mg-O PRDF shifts from 1.98 Å to 2.02 Å [17]. From neutron, x-ray diffraction data and a reserve Monte Carlo (RMC) simulation, authors in previous work suggested that Mg-O bond length has two values of 1.99 ± 0.01 and 2.21 ± 0.01 Å and Mg-O correlation is skewed to the high-r side [18]. Also by x-ray diffraction, M. Wilding et al indicate that the coordination number of Si-O is 4.0 while the one of Mg-O is 4.5 ± 0.2 due to the mixture of MgO x with x = 4, 5, 6 [10]. In particular, the fraction of MgO 4 , MgO 5 and MgO 6 is 68.8%, 27.8% and 3.4% respectively [18,19].
The intermediate-range order (IRO) is related to the connection of polyhedrons. There is no clear boundary between short and intermediate-range order. The IRO is usually considered relevant from 3 to 10 Å, that is ∼2π/Q First Sharp Diffraction Peak [20]. The SiO x polyhedrons bond to others by corner-, edge-and face-sharing linkages. The number of these linkages strongly change under compression [21][22][23][24]. At low pressure, the polyhedrons are connected mainly through corner-sharing linkages. When the pressure increases, edge-and face-sharing linkages are created, participating in the network formation. The same goes for Mg-Si pairs. The formation of face-sharing bonds has resulted in unequal density distribution and heterogeneity in the material [14,17,25]. The existence of Mg-rich regions under pressure is analyzed through the number of FO atoms and the formation of OMg x clusters [26]. In this paper, we carry out a different approach to the analysis of Mg-rich regions, which is to use ring statistics. Although ring statistics has been studied extensively in SiO 2 and MgSiO 3 glass [24,[27][28][29][30][31][32], studies on rings in MgSiO 3 liquid have been still in question.
The network information such as IRO [18] or the degree of polymerization [33] can be examined by the Q n species, in which, n is the number of bridging oxygens (BOs). In MgSiO 3 glass, the distribution of Q n units changes as a function of pressure [26]. In particular, when pressure rises up, the number of Q n , with 0 n 2, decreases while the number of Q n , with 3 n 6, increases. In MgSiO 3 liquid, however, the Q n distribution under compression, especially the Q n distribution on rings, has not been studied in detail. Therefore, in order to see the connectivity characteristic of rings into the network, we analyze the Q n distribution of SiO x units on different types of rings (here, n is still the number of BOs of SiO x units).
In terms of Voronoi, by RMC modelling and density function theory (DFT) simulation, previous work shows that the Voronoi volumes of Mg are greater than Bader volumes [34]. In MgSiO 3 glass, the Mg-Voronoi volume is about 13.36 Å for the RMC model and 12.69 Å for the DFT model. In Mg 2 SiO 4 , these values decrease to 11.83 Å and 12.52 Å, respectively. There have been some works on Voronoi in MgSiO 3 , however, the Voronoi analysis on rings in this material has not been clarified yet.
In this work, 5000-atom MgSiO 3 liquid models from 0 to 200 GPa are built by molecular dynamics simulation. We analyze the local structural change under pressure of the material through radial distribution function and coordination number. The Si-O ring statistics is clarified to explain the second peak splitting of the Si-Si radial distribution function at 200 GPa and the formation of Mg-rich regions under pressure. Q n and Voronoi are analyzed to provide additional information on rings. The article is arranged according to the layout: section 2 is simulation method, section 3 is the result and discussion, section 4 is the conclusion.

Methodology
Pressure-induced MgSiO 3 liquid models are erected by molecular dynamics simulation. Each model contains 5000 atoms (1000 Si, 3000 O and 1000 Mg atoms) with the periodic boundary condition. The temperature of the model is 3000 K and pressure is in the range from 0 to 200 GPa. The Oganov potential which is applied is presented as: in which, q , i q j are effective charges of atom i and atom j; r ij is the distance between i and j atoms; e is the charge of the electron; A ij and C ij are parameters attributing van de Waals repulsion and attraction, respectively; B ij is an e-folding length characterizing the radially-symmetric decay of electron repulsive Born energy between species i and j [35][36][37][38]. We use the Verlet algorithm with the time step of 0.47 fs to run the model. At first, atoms are randomly scattered inside a box and heated up to 6000 K to lose the memory of the initial configuration. The model then is quenched to 3000 K with the cooling rate of 2.5 K/ps. During that process, the model is relaxed at 5000 and 4000 K. At 3000 K, a long relaxation of 10 7 MD time steps is done in the NPT ensemble. The MgSiO 3 liquid model at 0 GPa is obtained. Finally, the model is compressed to desired pressures. The physical quantities of each model are analyzed based on averaging data from 1000 last MD steps.
To determine the correlations among atoms, the r cut is chosen based on the position of the first notch after the first peak of the PRDFs and adjusted to specific pressures. In this work, r cut for Si-O pair is 2.  Table 1. The first peak position of radial distribution functions (PRDFs) in MgSiO 3 liquid, at 0 GPa.
To analyze ring statistics, the -Si-O-network is considered as the undirected graph, whereas each atom is a node. Within the undirected graph, the path is the series of nodes and linkages that are connected consequently without overlap. There are a number of ring definitions. Some typical ring definitions that could be referred to are classic King's criterion [42], the irreducible [43,44] /primitivity [45,46]/no shortcut path [47] criterion, the shortest paths criterion [48,49]. Using different ring definitions could result in different number of rings. In this work, we use the no shortcut path criterion and algorithm which are introduced in the study of M Matsumoto et al [47]. n-fold ring consists of n T atoms (T is Mg or Si) and n O atoms. Typical Si rings are displayed in figure 2.

Structural change under compression
Under compression, the atoms are packed into small domains [20]. Therefore, both the coordination numbers of Si and Mg strongly depend on pressure (see ). In addition, there is a tight relation between PRDF and bond angle distribution [26]. Hence, the mean O-O bond length also decreases significantly under compression. Figure 5 shows the change in the fraction of SiO x (x 4) units. At 0 GPa, SiO 4 accounts for about 75%, SiO 5 is about 20%, and the rest is mostly SiO 3 . SiO x units with x 6 are then almost non-existent in the system. The formation of SiO 3 and SiO 5 in MgSiO 3 at 0 GPa is due to doping MgO [41]. The ratio of SiO x units varies drastically with pressure. Specifically, the proportion of SiO 4 , SiO 5 , SiO 6 , and SiO 7 units alternately account for the highest proportion at pressures from 0 to 6 GPa, 6 to 27 GPa, 27 to 125 GPa and from 125 onwards. The

Ring analysis
The ring analysis is carried out to somewhat clarify the IRO in MgSiO 3 . Figure 7    pressure. The Mg species would break the short paths connecting two Si atoms. Therefore, the ring statistics are not in the gaussian shape at pressures above 0 GPa.
In the case of thering, at ambient pressure, the 4-and 5-fold rings account for around 24% and these results. The 6-fold, 7-fold and 8-fold rings are around 15.2%, 6.2% and 2.2%, respectively. At higher pressures, the ring size distribution shifts to the left. From 50 GPa, the 3-fold ring accounts for the highest proportion, it is even more than 50% at pressure higher than or equal 100 GPa. Besides, the large rings (size 6) almost vanish at high pressure. At 200 GPa, the fraction of the 5-fold ring is about 3%. Under compression, the mean coordination number of Si and Mg increased (see figure 3). Atoms are more tightly packed. Therefore, the percentages of small T-O rings increasingly rise, and the percentages of largerings go down. For Si-O rings, there is a massive change in the ring statistics under compression. At 0 GPa, the Si-O ring statistic has a gaussian form, with a peak at 5-fold. Paradoxically, when the pressure increases and reaches 20 GPa, the 5-fold ring has the lowest proportion. The 11-fold ring which has the lowest proportion at 0 GPa now is dominant. As the   pressure increases further, the ring statistics continue changing significantly but 4-and 5-fold rings still account for the smallest ratio. This shows the strong influence of pressure on the structure of the Si-O network. At 200 GPa, the 9-fold ring fraction is surprisingly higher than others. This will affect the intermediate range order of the -Si-O-network. In particular, we use that phenomenon to explain the splitting peak of the second peak in Si-Si PRDF.
For convenience, we consider the distribution of bond length from a Si atom on a ring to the nearest second one on the same ring as the Si-Si-Si distance distribution. The Si-Si-Si distance distribution depends very strongly on certain rings (see figure 8). The 4-fold, 5-fold and 6-fold rings have a gaussian form distribution, with peaks at 4.1, 4.8 and 4.7 Å, respectively. The Si-Si-Si distance distribution has a narrower full width at half maximum (FWHM) on 4-fold rings than the ones on 5-fold rings and 6-fold rings. The Si-Si-Si distance distribution of 10 and 11-fold rings varies significantly compared to the small rings. Therefore, large rings are very topologically flexible. Large rings are constructed by many Si-O paths. The SiO x units on large rings have more ability to arrange than on small rings with a constraint minimizing the cohesive energy of the whole configuration. Conversely, low-fold rings, typically 2-fold and 3-fold rings, are formed by stretching Si-O bond length, bending of Si-O-Si and O-Si-O bond angle to optimize the raised energy [27]. The smaller the rings are, the more the rings are affected by the energy constraint. Compared to small-size rings, higher-fold ones have more varied shapes. Hence, Si-Si-Si distance distribution in large rings does not have the Gaussian shape. Although the distance distribution of Si-Si-Si of the 9-fold ring is asymmetric, it has a sharp peak at about 5.3 Å. Ring size distribution, as shown in figure 7, indicates that at 200 GPa, 9-fold rings are enormous in proportion compared to others. The 9-fold ring proportion is almost 24% while the second most dominant one is the 11fold ring which accounts for only about 15%. The combination of the dominance in the ratio and the high peak position of Si-Si-Si distance distribution of the 9-fold ring regulates Si-O PRDF. The second peak of Si-Si PRDF tends to shift to the left under compression (see figure 4), and at 200 GPa, the second peak splits into two subpeaks. The first subpeak is at about 4.6 Å while the second one is at about 5.3 Å. Based on the above analysis, the phenomenon that the second peak splits into two small subpeaks at 200 GPa is attributed to the massive formation of 9-fold rings. The small subpeak corresponds to the reduction of averaged bond length among interactive atoms (similar to the cases in pressure range beneath 200 GPa). The high subpeak is ascribed to the position of the sharp peak in Si-Si-Si distance distribution. Figure 9 demonstrates the Q n distributions of rings under compression to have a better understanding of rings. At 0 GPa, Q n distribution of 2-, and 3-fold rings have peaks at n = 4 (see figure 9(a)). On 2-fold rings, the value of Q 3 is close to the value of Q 5 , about 27%. The fraction of Q 4 in 2-fold rings is higher than the one in other rings. For rings of size 4 or higher, Q 3 accounts for the highest proportion. Q n has almost the same distribution across the different rings at the same pressure, except for the distribution of Q n on the 11-fold ring at 100 GPa (see figure 9(c)). The peak shift of Q n distribution from 3 or 4 at ambient pressure to 6 or 7 at pressure of 50 GPa or higher indicates that rings tend to merge under compression.

Mg-rich region
The ring statistical analysis can also explain the heterogeneity of MgSiO 3 . When large rings are formed on the -Si-O-network, inside those large rings there is no silicon atom because there is no shortcut path [47]. Therefore, the domain inside the large ring will have a high density of oxygen atoms creating a negative charge region (see figure 10(a)). The Mg 2+ ions will therefore be attracted inward to neutralize the charge (see figure 10(b)) forming an Mg-rich region. The spatial distribution of Mg 2+ ions is not uniform, leading to the heterogeneity of the material. As the pressure increases, the greater the proportion of large Si-O rings appears (see ring statistics in figure 7(a)), the higher the ability to capture Mg atoms. In other words, the heterogeneity of MgSiO 3 increases with pressure.

Voronoi
The Si-and O-Voronoi volume distributions on the rings are also analyzed to explain the arrangement of atoms on the -Si-O-network. Figure 11 shows the change in the volume of the Voronoi on particular rings with pressure. Figure 12 is an example of Si and O-Voronoi on 6-fold ring at pressures of 50 and 200 GPa, and the spatial distribution of Voronoi on the 6-fold is shown in figure 13. At 0 GPa, the Si-voronoi volume ranges from 7.50 to 7.76 Å 3 . In contrast, at ambient pressure, O-Voronoi are strongly dependent on ring type and have a much larger volume than Si-voronoi. At 0 GPa, the volume of O increases from about 13 Å 3 to almost 16 Å 3 , as the ring size increases from 2 to 6. However, with rings having size greater than 6, the larger the ring size, the  For the larger rings, the atoms are more flexible to arrange to form highly symmetric SiO x units, which helps to reduce the energy of the system [27]. This means that the SiO x units on the large rings (size > 5) have less distortion. There is a relation between Si-O-Si angle and O-Voronoi volume (see figure 11(b) and figure 14). It is noticed that the mean O-Voronoi volume is proportional to the mean angle of Si-O-Si. The fact that mean Si-O-Si angle reaches a value of 148.2°then gradually decreases to 137.5°could correspond to the volume of O-Voronoi decreasing when the ring size is greater than 6. Hence, the mean O-Voronoi volume in 6-fold rings is greater than the ones in other rings.
At pressures of 20 GPa or more, the mean Voronoi volumes of atoms of the same type on different rings have almost the same value. The Si atoms are surrounded by O atoms, while the O atoms are bonded to both Mg and Si atoms. Besides, the mean Si-O bond length is shorter than the mean Mg-O bond length at the same pressure condition [17]. The average distance from the center point of Si-Voronoi is then shorter than the one of O-voronoi. Thus, the O-Voronoi would have a larger volume than the Si-Voronoi. Specifically, at pressures at  50, 100 and 200 GPa, the Si-Voronoi volume is 6.7, 6.1, 4.8 Å 3 while the O-Voronoi volume is 9.0, 7.8 and 6.0 Å 3 , respectively. As the pressure increases, the number of edge-sharing and face-sharing bonds increases (see appendix -table A1), the size of the model decreases due to compression, and the polymerization of the -Si-O-  network increases. The atoms will be tightly packed together, and the rings connect to others or even overlap (for example, 2-fold rings overlap a 10-fold ring in figure 10(a) and 6-fold rings merge under compression in figure 13). Therefore, Voronoi of the same type on different rings have the same average volume.

Conclusion
We conduct a study on the local structural change and intermediate-range order of MgSiO 3 liquid at 3000 K under compression. The Si-Si and O-O PRDFs rapidly change by pressure while the Si-O PRDF is quite stable, because of the formation of edge-and face-sharing linkages as well as the increase in the mean coordination number of Si. At 200 GPa, the 9-fold rings dominate the system, Si-Si-Si distance distribution on 9-fold rings has a sharp peak at around 5.3 Å and Si-O linkages shrink due to high pressure. This infers the splitting peak of Si-Si PRDF into 2 subpeaks at 4.6 Å and 5.3 Å. The formation of large rings, at high pressures, sheds light on the formation of Mg-rich regions. In particular, the O atoms inside the large rings and there is no shortcut path which contains Si species. They would generate a high negative charge domain, attracting the Mg 2+ ions and forming Mg-rich regions. The non-uniform distribution of Mg in space leads to the heterogeneity of the MgSiO 3 . In terms of Voronoi analysis, our results indicate that at ambient pressure, the mean Si-Voronoi volume is in the range of 7.5-7.8 Å 3 on different types of rings. The mean O-Voronoi volume increases from about 13.0 Å 3 on 2-fold ring to nearly 16.0 Å 3 on 6-fold ring, then decreases to 14.3 Å 3 on 11-fold ring. This phenomenon is attributed to the topological distortion of SiO 4 units of small rings, the mean Si-O-Si angle and the way SiO x units arrange on large rings. Table A1. The number of corner-(N c ), edge-(N e ), face-sharing (N f ) network per the number of Si atoms at different pressure, corresponding with the mean bond length of them, corner-(CSBL), edge-(ESBL), face-sharing bond length (FSBL).