Molecular Dynamic Simulation of Large Model of Silica Liquid - Nguyen Thi Thanh Ha

VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 19 Molecular Dynamic Simulation of Large Model of Silica Liquid Nguyen Thi Thanh Ha1,*, Phan Quan1, Tran Van Hong2, Le Van Vinh1 1Department of Computational Physics, Hanoi University of Science and Technology, Vietnam 2Department of Physics, Thai Nguyen University of Education, Thai Nguyen, Vietnam Received 01 January 2018 Revised 30 February 2018; Accepted 20 March 2018 Abstract: We perform a molecular dynamics simulation to study the microstructure and dynamical properties in large silica model at liquid state. The models consisting of 19998 atoms were constructed under a wide range of pressure (0-20 GPa) and at 3500K temperature. Structural characteristics were clarified through the pair radial distribution function (PRDF), the distribution of SiOx coordination units and network structure. The result shows that these liquids consist of identical units SiO4, SiO5 and SiO6 and have common partial O―Si―O angle distribution. Furthermore, the major change in the diffusion mechanism under pressure is also considered and discussed. Keywords: Molecular dynamics, structure, coordination units, diffusion, network structure. 1. Introduction Silica and silicate minerals (mixture of SiO2 and other metal oxides) play an important role in geosciences and technology. So, the complete knowledge of the structure and dynamical properties in silica and silicate under conditions of high pressure and temperature is quite necessary. The results of research studies will optimize the process of manufacturing new materials for intended use and controlling geological activities [1-4]. As we have known, the silicate glass- and melt- structures consists of SiO4 tetrahedra linking to each other form continuous random tetrahedral-network in three-dimensional space [5-7]. In pure silica glass and melt, each SiO4 tetrahedron connects to four other adjacent SiO4 tetrahedra via bridging-oxygen (BO). The addition of other oxides (such as Na2O, CaO, MgO or PbO) into pure silica (SiO2) disrupts the basic silica network by breaking part of the Si-O bonds, generating non-bridging oxygen (NBO). The percentage of NBOs in the system increases with the alkali content [8-9]. The detailed degree of mixing between Si-O network and metal ions as well as the distribution of BO and NBO species yields insights into the atomic scale structure that govern the mechanical-physical-chemical properties of silicate glass system. The degree of polymerization  Corresponding author: Tel.: 84-983012387. Email: ha.nguyenthithanh1@hust.edu.vn https//doi.org/ 10.25073/2588-1124/vnumap.4231 N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 20 depends on the ratio between the number of bridging oxygen (BO) and SiO4 tetrahedral units (the abundance of different Q(n) species, where Q represents the SiO4 tetrahedron; n is the number of BO). The response of viscosity to pressure is found to be strikingly distinct between polymerized and simple liquid [10-11]. The degree of polymerization is also closely related to the diffusion of atoms in silicate. The neutron, X-ray diffraction (XRD) techniques, magic-angle spinning (MAS) nuclear magnetic resonance (NMR) [12-13] and simulation methods [14-15] have observed these phenomena. Liquid silica is a typical network forming system. Its structure consists of basic structural units SiOx (x=4, 5, 6) and the SiO4 units is dominant at low pressure. When the liquid is compressed to smaller volume, the units SiO5 and SiO6 become prevalent. Due to its network-forming ability, the liquid silica exhibits a number of peculiar properties, which have been observed in both simulation and experiment. Namely, the compressibility of liquid silica in the interval of pressure up to 10 GPa is substantially higher than that of quartz [16]. Moreover, anomalous behavior diffusion and spatially heterogeneous dynamics are also observed [17-18]. In our previous works, the micro-structure and dynamical properties in silica liquid also have been investigated via through the pair radial distribution function (PRDF), the distribution of SiOx coordination units, the bond angle distribution and structure network. However, the research model has only 2000 atoms. To check the accuracy of previous research results and evaluate the effect of model size, we have conducted a SiO2 model consisting of 19998 atoms by mean of molecular dynamic simulation. This is a large model and quite difficult to build it. We will investigate the structural characteristic and coefficient diffusion in silica liquid under different pressure, compare with the experimental results and simulations. Moreover, the major change in the diffusion mechanism under pressure is also considered and discussed. 2. Calculation We have prepared the model which consists of 6666 Si and 13332 O by means of MD simulation. In this paper, we simulate silica liquid with the BKS potential. The integration is performed using a velocity-Verlet algorithm with time step of 1.0 fs. Initial configuration of MD silica liquid model is constructed by randomly in simulation box and heated up to 6000K to remove possible memory effects. Then the obtained model is cooled down to 3500K. Next, a long relaxation has been done in NPT ensemble (moles (N), pressure (P) and temperature (T) are constant) to produce a model at 3500K and upon ambient pressure. We obtain a sample at ambient pressure which is denoted to model M1. The model M1 has been compressed to different pressure (5, 10, 15, 20 GPa). In order to improve statistics, all quantities of considered structural data were calculated by averaging over the 5.000 configurations during the last simulation (105 MD steps). To observe the dynamical processes, two above models are relaxed in NVE ensemble (the system is isolated from changes in moles (N), volume (V) and energy (E)) for a long time (106 MD steps). 3. Results and discussion 3.1. Radial distribution function Firstly, we examine the structural characteristics of constructed models. The Fig.1 presents the PRDF of Si–Si, Si–O and O-O at different compressions. It can be seen that the first peak all atomic pairs decreases in amplitude and becomes broader under compression. The position of the first peak of Si–Si, and O–O pairs decreases meanwhile Si–O pairs, the position of the first peak of Si- O increases. N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 21 Moreover, the shift of the first peak of gSi–O ( r ) is the least and this indicates that the short-range order of silica liquid is not sensitive to the compression at pressures ranging from 0 to 20 GPa. The detail result is showed in Table 1. The characteristic of the PRDF is in good agreement with the reported data in refs. [19,20] Fig 1. Radial distribution functions of Si−Si, Si−O, and O−O pairs at different pressures. 3.2. Coordination units The structure of SiO2 liquid consists of structural units SiOx (x = 4, 5, 6) and OSiy (y = 2,3). Figure 2 (a) distribution of coordination units SiO4 in the model. In which, only coordination units SiO4 are drawn, while the others are removed. Similarly, Figure 2 (b)–(c) shows the distributions of coordination units SiO5, SiO6. The detail dependence of the fraction of coordination units SiOx and OSiy on pressure was presented in Table 2. 0 1 2 3 Si-O 0 3 6 9 g (r ) O-O 0 GPa 5 GPa 10 GPa 15 GPa 20 GPa 0 2 4 6 8 10 0 1 2 3 r, Å Si-Si N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 22 Table 1. Structural characteristics of SiO2 liquid, rlk is positions of first peak of PRDF, glk is high of first peak of PRDF Model 0GPa 5 GPa 10 GPa 15 GPa 20 GPa Ref [19,20] rSi-Si, [Å] 3.16 3.12 3.14 3.12 3.12 3.12 rSi-O, [Å] 1.62 1.62 1.62 1.64 1.62 1.62 rO-O, [Å] 2.64 2.52 2.56 2.5 2.48 2.65 gSi-Si 2.86 2.67 2.45 2.37 2.35 - gSi-O 8.62 7.14 6.19 5.54 5.26 - gO-O 2.72 2.44 2.37 2.39 2.37 - It has been seen that the number of SiO4 and OSi2 unit is domain at ambient. When increasing pressure, the fraction of SiO5 and SiO6 increases meanwhile fraction of OSi3 and OSi2 also increases. This demonstrates that a transition from SiO4 to SiO5 or SiO6 should be accompanied by a transition Fig 2. Snapshots of coordination unit distributions SiOx in model at ambient pressure: (a) distribution of coordination units SiO4, (b) distribution of coordination units SiO5, (c) distribution of coordination units SiO6. (a) (b) (c) N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 23 from OSi2 to OSi3. It means that increasing pressure, there is a transformation from four-fold coordination (SiO4) to five- and six-fold coordination (SiO5 and SiO6). Table 2. The percentage fraction of the coordination units in silica liquid at different pressure Fig.3. Partial bond angle distributions for structural units SiOx (x = 4, 5, 6) at different pressures. Model 0GPa 5 GPa 10 GPa 15 GPa 20 GPa SiO4 93.94 77.26 59.89 44.89 35.69 SiO5 3.90 19.34 34.65 45.66 50.13 SiO6 0.03 0.79 3.51 8.14 12.77 OSi2 95.72 86.26 75.53 66.01 59.95 OSi3 4.28 13.74 24.47 33.99 40.05 0 3 6 SiO 5 0GPa 5GPa 10GPa 15GPa 20GPa 0 3 6 SiO 6 F ra ct io n 40 60 80 100 120 140 160 180 0 3 6 Angle ( degree) SiO 4 N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 24 To investigate the short-range order, the O−Si−O bond angle and Si−O bond length distribution in the coordination units SiOx (main coordination units) have been calculated. Figure 3 presents the partial bond angle distributions for structural units SiOx (x = 4, 5, 6) at different pressures. The results show that the partial O-Si−O bond angle distribution in each kind of coordination unit SiOx is almost the same for different pressure. This means that the distribution of partial bond angle in SiO4, SiO5, and SiO6 units is not dependent on pressure. Here angle distribution in SiO4 units has a form of Gauss function and a pronounced peak at 105° and 90° with SiO5 unit. This result is similar to experimental and other simulated data reported in [19, 21] and indicates a slightly distorted tetrahedron with a Si atom at the center and four O atoms at the vertices. In the case of angle distribution in SiO6, there are two peaks: a main peak locates at 90° and small one at 165°. The partial Si−O bond length distribution in coordination units SiO4, SiO5, and SiO6 is shown in Fig. 4. It shows that The Si – O bond length distributions in SiO4, SiO5, and SiO6 units have peaks at 1.62, 1.64, and 1.74 Å, respectively. This structural characteristic has been explained via the Coulomb repulsive force between anion and anion (O-2 and O-2 ions).The force in SiO5 and SiO6 units is much stronger than the one in SiO4. This leads an increase of Si–O bond length in SiO5 and SiO6. Fig 4. The Si−O bond length distributions for structural units SiOx (x = 4, 5, 6) at different pressures. 0 3 6 9 0GPa 5GPa 10GPa 15GPa 20GPa SiO 5 0 3 6 SiO 6 F ra c ti o n 1.0 1.5 2.0 2.5 0 3 6 Distance (Å ) SiO 4 N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 25 The diffusion process in silica under compression The diffusion coefficient of particles is determined from the mean square displacement (MSD) of atom via Einstein equation 2( ) lim 6t R t D t    Where t=N.TMD; N is number of MD steps; TMD is MD steps and equal 1.0 fs, The Fig 5 describes the time dependence of Silicon (or Oxygen) MSD under different pressure. One can see that the data fall on straight- line plot. We have calculated coefficient diffusion from the slope of these lines. Fig 5. The step dependence of mean square displacement in silica liquid at different pressure. The anomalous diffusion is observed via Fig.6. It can be seen that there is a monotonous growth of coefficient diffusion with upon compression and a maximum point near P=10 GPa that consistent with previous simulation [22-23]. The self-diffusion coefficient increases with increasing pressure (0÷ 10 GPa) meanwhile the self-diffusion coefficient decreases with pressure at 10÷20 GPa. We find that the change in diffusion mechanism between low and high-pressure in silica liquid. The change in diffusion mechanism can be explained by that upon compression the SiO4 unit transforms to SiO5 and SiO6 where the Si–O bond is much weaker. This leads an increase in mobility of Si and O atoms. On the other hand the liquid becomes much denser under pressure and diffusion process more difficult. One can see that the bond weakening and the liquid densification affect the diffusion coefficient in opposite directions. At the low-pressure configuration, the influence of the bond weakening of Si, O in the newly formed structure units SiO5 and SiO6 is larger than that of liquid densification. Therefore, the coefficient diffusion increases. However, at the high-pressure configuration, rate of transition from SiO4 to SiO5 or SiO6 decreases and effect of liquid densification predominate. The coefficient diffusion decreases and Fig.6 has a maximum point. 0 100000 200000 300000 400000 0 100 200 300 400 M S D (Å 2 ) n (step) 0 GPa 5 GPa 10 GPa 15GPa 20GPa Silicon 0 100000 200000 300000 400000 0 100 200 300 400 500 600 M S D (Å 2 ) n (step) 0 GPa 5 GPa 10 GPa 15GPa 20GPa Oxygen N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 26 Fig 6. The pressure dependence of diffusion coefficient of silicon and oxygen in silica liquid. 3. Conclusion In this work, microstructure and dynamic properties of large model of silica liquid at pressures ranging from 0 to 20 GPa were investigated in detail. Results reveal that the structure of SiO2 liquid consists of structural units SiOx (x = 4, 5, 6) and OSiy (y = 2,3). At ambient, the number of SiO4 and OSi2 unit is domain. Upon compression, the fraction of SiO5 and SiO6 increases meanwhile fraction of OSi3 and OSi2 also increases. The distribution of partial bond angle and bond length in SiO4, SiO5, and SiO6 units is not dependent on pressure. The cation–cation, anion–anion Coulomb repulsion is the origin of the increase of the mean Si-O bond length under compression. Furthermore, the diffusion in silica liquid shows anomalous behavior and the origin of anomalous diffusivity is due to the change of SiO5 and SiO6 concentration under compression. These results are similar to our previous simulation. These demonstrate that the model silica liquid consisting of 2000 atoms was good enough for us to study structural properties and diffusion. However, the large model will help us further study the original of dynamical heterogeneity that is difficult with the small model (2000 atoms). Acknowledgement The authors are grateful for support by the Ministry-level projects (grant No B2018-BKA-57) References [1] Gergely Molnár, Patrick Ganster, Anne Tanguy, Physical review E 95, 043001 (2017) [2] M. M. Smedskjaer,Frontier Mater, 1(23),1,(2014) [3] B. Hehlen and D. R. Neuville. J Phys Chem B. 119 (10), 4093,( 2015) [4] T. Kawasaki, H. Tanaka, J. Phys.: Condens. Matter 22, 232102 (2010). [5] G. Calas, L. Galoisy, L. Cormier, G. Ferlat, G. Lelong,Procedia Materials Science 7, 23 (2014) 0 5 10 15 20 0.04 0.08 0.12 0.16 0.20 D [ 1 0 5 c m 2 /s ] Pressure (GPa) Silicon Oxygen N.T.T. Ha et al. / VNU Journal of Science: Mathematics – Physics, Vol. 34, No. 4 (2018) 19-27 27 [6] J. Badro, D. M. Teter, R. T. Downs, P. Gillet, R. J. Hemley, and J.L. Barrat, Phys. Rev. B 56, 5797 (1997) [7] C. Weigel, L. Cormier, G. Calas, L. Galoisy, D.T. Bowron, Phys. Rev. B 78, 064202 (2008) [8] H. Jabraoui, Y.Vaills, A. Hasnaoui, M. Badawi and S. Ouaskit, J. Phys. Chem. B 120, 13193 (2016). [9] T. K. Bechgaard eltal., J. Non-Cryst. Solids 441, 49 (2016) [10] S. K. Baggain, D. B. Ghosh, B. B. Karki, Phys. Chem. Min. 42, 393 (2015). [11] A. W.Cooper, P. Harrowell, and H. Fynewever, Phys. Rev. Lett 93, 135701 (2004). [12] J.R. Allwardt, J.F. Stebbins, B.C. Schmidt, D.J. Frost, A.C. Withers, M.M. Hirschmann, Am. Mineral. 90, 1218 (2005) [13] M. Bauchy, M. Micoulaut, Physical review B 83, 184118 (2011) [14] H. Jabraoui, E.M. Achhal, A. Hasnaoui, J.-L. Garden,Y.Vaills, S. Ouaskit, J.Non-Cryst-Solid 448,16(2016) [15] S.K. Lee, G.D. Cody, Y. Fei, B.O. Mysen, Chem. Geol. 229,162 (2006). [16] I. Jackson, Phys. Earth Planet. Inter. 1, 218 (1976) [17] B.T. Poe etal., Science 276, 1245 (1997) [18] M. Scott Shell, G.D.Pablo , Z.P. Athanassios, Phys. Rev.E 66, 011202 (2002) [19] D.I. Grimley, A.C. Wright, R.N. Sinclair, J. Non-Cryst. Solids 119, 49 (1990). [20] Mozzi R L and Warren B E, J. Appl. Crystallogr. 2 164 (1969) [21] Bauchy M, J Chem Phys. 141, 024507 (2014) [22] T. Morishita, Phys. Rev. E 72, 021201 (2005) [23] P.K. Hung, N.T.T. Ha, N.V. Hong, J. Non-Cryst. Solids 358, 1649 (2012)