⇦Chinese Physics Letters, 2020, Vol. 37, No. 7, Article code 076201 Pressure Effects on the Transport and Structural Properties of Metallic Glass-Forming Liquid Qi-Long Cao (曹启龙)*, Duo-Hui Huang (黄多辉), Jun-Sheng Yang (杨俊升), and Fan-Hou Wang (王藩侯) Affiliations Key Laboratory of Computational Physics, Yibin University, Yibin 644007, China Received 9 April 2020; accepted 3 May 2020; published online 21 June 2020 Supported by the National Natural Science Foundation of China (Grant No. 11704329.)
*Corresponding author. Email: qlcao@mail.ustc.edu.cn Citation Text: Cao Q L, Huang D H, Yang J S and Wang F H 2020 Chin. Phys. Lett. 37 076201    Abstract Transport and structural properties of metallic glass-forming liquid Cu$_{50}$Zr$_{50}$ are investigated by molecular dynamics simulation, under high pressures from 1 bar to 70 GPa. The following results have been obtained: (i) reversals of component diffusion coefficients ($D_{\rm Cu}$ and $D_{\rm Zr}$) are observed at the reversion pressure. At low pressures below the reversion pressure, $D_{\rm Cu}/D_{\rm Zr}$ decreases from about 1.4 to 1.0. At high pressures above the reversion pressure, $D_{\rm Cu}/D_{\rm Zr}$ decreases more rapidly from 1.0 to about 0.7. (ii) Component diffusion coefficients decay exponentially with pressure up to reversion pressure, then the strength of the exponential dependence changes, while the pressure-dependent behavior of viscosity can be well described by a single exponential relation over the full range of pressure. (iii) The Stokes–Einstein relation (SER) works well at low pressures and starts to be violated at the breakdown pressure. For glass-forming liquid Cu$_{50}$Zr$_{50}$ along the 2000 K isotherm, the breakdown pressure equals the reversion pressure of component diffusion coefficients and is about 35 GPa. (iv) The pressure dependences of the ratio between component diffusion coefficients can be used to predict the breakdown pressure of SER along isotherm. The validity of SER and the reversals of component diffusion coefficients are found to be related to the pressure dependence of the relative total fractions of predominant Voronoi polyhedrons around individual components. DOI:10.1088/0256-307X/37/7/076201 PACS:62.50.-p, 66.20.-d, 61.20.-p © 2020 Chinese Physics Society Article Text Transport properties of metallic glass-forming liquids and the relationships between transport properties have attracted tremendous interest because of the scientific importance and technological relevance. Temperature and pressure are two important thermodynamic variables that could be used to tune transport properties of liquids.[1,2] Metallic glasses are conventionally obtained by cooling metallic liquid with a sufficiently high cooling rate to avoid crystallization. Consequently, numerous studies focus on temperature effects on the transport properties of liquids. For metallic liquids, under ambient pressure, it is widely (but not uniformly or completely) believed that: (i) the Arrhenius relation represents the temperature dependence of transport coefficients in the high-temperature region well, and it will be replaced by the super-Arrhenius behavior in the low-temperature region.[3–7] (ii) The dynamics crossover from high-temperature Arrhenius to low-temperature super-Arrhenius manner is located at about twice of glass transition temperature.[5,6,8] (iii) Below the dynamics crossover temperature, the familiar Stokes–Einstein relation (SER) is broken down and replaced by a fractional Stokes–Einstein relation.[5–11] On the other hand, although the study on the relaxation of liquids using elevated pressure has a long history, the pressure dependence of properties of metallic glass-forming liquid has drawn much less attention than their temperature dependence. Answers to many questions regarding pressure dependence are either not known or ill-understood. Previous studies on this topic show that pressure may have different effects on transport property and local structure for metallic liquids.[12,13] For instance, Arnab et al. found that the transport coefficients show an exponential dependence on pressure (which is weaker than temperature dependence) and there is a clear break in the strength of the exponential dependence at a certain pressure, in the Kob–Andersen binary mixture model.[12] Molecular dynamics simulations on the Zr$_{50}$Cu$_{44}$Al$_{6}$ shows that the temperature- and pressure-induced glass transitions differ in the formation of local coordination structures and the variation of fragility.[13] Hence, it is natural to expect that pressure may have different effects on the relationships between transport properties. However, the direct measure of transport property is still lacking and of challenge. The paucity of transport coefficients at high pressure has limited the study of this topic, and the effects of pressure on the validity of SER and the dynamics crossover are still unclear. Cu–Zr alloys are usually considered as ideal model systems to study the metallic glass-forming liquids, and a large number of articles have reported on the transport properties, local structure, and SER in liquid and amorphous Cu–Zr alloys with various compositions under ambient pressure. Based on an optimized embedded-atom-method interatomic potential provided by Sheng et al.,[14] molecular dynamics (MD) simulation has been proven to be an effective way of simulating transport and structure properties of Cu–Zr alloys at ambient-pressure as well as high-pressure conditions.[13,15] It motivates us to carry out the MD simulations for liquid Cu$_{50}$Zr$_{50}$ alloy in a pressure range from 1 bar to 70 GPa keeping the temperature constant at 2000 K, thereby performing a systematic study of the pressure effects on the transport properties and the relationships between transport properties (SER) for liquid metallic alloy. The choice of the temperature 2000 K permits the long pressure range needed in this study, at low pressures the system keeps in the stable liquid state and turns to undercooled liquid state at high pressures. In this work, all MD simulations were performed using the open source code-LAMMPS,[16] Initially, 3375 Cu atoms and 3375 Zr atoms were generated in a cubic box randomly under three-dimensional periodic boundary conditions, and equilibrated at high temperature under different external hydrostatic pressures. Samples are prepared by a quench from liquid at high temperature to 2000 K with a rate of $1\times 10^{12}$ K/s. The NPT (constant number of particles, pressure, temperature) ensemble was applied to reach constant pressure and temperature with a time step of 1 fs, and the ensemble was switched to the canonical ensemble (NVT) for production runs after equilibration with the same time step. The pair distribution functions and Voronoi polyhedrons are obtained by averaging 40 configurations.
cpl-37-7-076201-fig1.png
Fig. 1. (a) Self-diffusion coefficients of components and (b) the viscosity of Cu$_{50}$Zr$_{50}$ versus pressure at 2000 K. The solid and dashed lines are the best-fit lines. The inset in (a) shows the pressure dependences of the ratio $D_{\rm Cu}$/$D_{\rm Zr}$ along 2000 K isotherm, and the inset in (b) shows the stress autocorrelation functions (SACFs).
The self-diffusion coefficients of Cu and Zr ($D_{\rm Cu}$ and $D_{\rm Zr}$) in liquid Cu$_{50}$Zr$_{50}$, calculated from the long-time evolution of respective mean-square displacements are shown in Fig. 1(a), versus pressure. Due to the lack of available data about the self-diffusion coefficients at 2000 K for liquid Cu$_{50}$Zr$_{50}$, it is difficult to give a direct comparison of our simulation with previous studies. Han et al.[4] studied the self-diffusion coefficients of Cu in liquids Cu$_{8}$Zr$_{3}$ and CuZr$_{2}$ using a modified embedded atom method and $D_{\rm Cu}=4.7\times 10^{-9}$ m$^{2}$/s and $3.52\times 10^{-9}$ m$^{2}$/s at 2000 K, respectively, which agree quite well with the result in this work, $D_{\rm Cu}=4.02\times 10^{-9}$ m$^{2}$/s at 2000 K under ambient pressure condition. As is expected, both $D_{\rm Cu}$ and $D_{\rm Zr}$ decreases with increasing pressure. A striking note in this figure is the reversal of self-diffusion coefficients of Cu and Zr. At low pressures below 35 GPa, the self-diffusion coefficients of Cu are always larger than that of Zr, and the difference between them decreases with increasing pressure. At high pressures above 35 GPa, the self-diffusion coefficients of Cu is always smaller than that of Zr, and the difference between them increases with increasing pressure. To make it clearer, we plot the ratio between $D_{\rm Cu}$ and $D_{\rm Zr}$ versus pressure in the inset of Fig. 1(a). It shows that the reversion pressure is about 35 Gpa and there is a break in the slope of the pressure dependence of the ratio at 35 GPa. Another interesting feature is that the pressure dependence of self-diffusion coefficients cannot be described by a single Arrhenius or exponential relation over the full range of pressure. As suggested by Bagchi et al.,[12] the exponential dependence on pressure, i.e., $$ D_{\alpha }=\xi_{\alpha }\exp \left(\lambda_{\alpha }P \right) +D_{\alpha }^{0},~~ \tag {1} $$ was used to fit the self-diffusion coefficients. In the above equation, $\zeta_{\alpha}$, $\lambda_{\alpha}$ and $D_{\alpha }^{0}$ are fitting parameters. Figure 1(a) also depicts the fitted lines at low pressures (solid lines) and at high pressures (dashed lines). We can find that there is a change in the coefficient $\lambda_{\alpha}$ after about 35 GPa. This means that the strength of pressure dependence changes at the reversion pressure. This result is in agreement with the previous study on the Kob–Andersen binary mixture model for which there is a change in the pressure dependence of self-diffusion coefficients.[12] Figure 1(b) shows the viscosities of liquid Cu$_{50}$Zr$_{50}$ at 2000 K, calculated according to the Green–Kublo equation by a time integral over the stress autocorrelation functions (SAFs),[17] versus pressure. In order to determine the correlation time at different pressures, the normalized SAFs over the whole pressure range are given in the inset in Fig. 1(b). A glance at this figure shows that the SACFs decay monotonically toward zero as the correlation time increases, and the correlation time increases with increasing pressure. At 1 bar, our calculated viscosity is about 2.28 mPa$\cdot$s, which is in reasonable agreement with the previously reported data 2.72 mPa$\cdot$s given by Hu et al.[18] As is expected, viscosity increases with increasing pressure at constant temperature. Contrasting to the pressure-dependent behavior of self-diffusion coefficients, the pressure-dependent behavior of viscosity can be well described by a single exponential relation over the full range of pressure with $$ \eta =\sigma \exp \left(\kappa P \right)+\eta_{0},~~ \tag {2} $$ where $\sigma$, $k$ and $\eta_{0}$ are fit parameters. This may imply the existence of decoupling of self-diffusion coefficient and viscosity behind the reversion pressure. This result is different from the result on the Kob–Andersen binary mixture model for which there is a change in the pressure dependence of viscosity, but in agreement with the result on the liquid Zr$_{50}$Cu$_{44}$Al$_{6}$ for which the relaxation time can be well described by the pressure counterpart of the Vogel–Fulcher–Tamman equation in the whole pressure range.[13] SER is often used to describe the relationship between the solute diffusion coefficients $D$ and size $R$ and the solvent viscosity $\eta$ as follows: $$ D\eta =\frac{k_{_{\rm B}}T}{4\pi R}, $$ where $T$ is the temperature of the system. Figure 2 presents the pressure dependence of the product $D_{\rm Zr}\eta R_{\rm Zr}$ keeping temperature at 2000 K, and the effective diameter $R_{\rm Zr}$ is selected by the first peak position $r_{\rm ZrZr}$ of the partial pair distribution function $g_{\rm ZrZr}(r)$. Two pressure regimes can be distinguished: the product $D_{\rm Zr}\eta R_{\rm Zr}$ is independent of pressure at low pressures below 35 GPa, followed by a repaid increase with increasing pressure at high pressures up to 70 GPa. As predicted by the SER, as long as $D_{\rm Zr}\eta R_{\rm Zr}$ is independent of pressure along the isotherm, the SER is deemed to hold, otherwise the SER has been broken down. This means that the SER holds well at low pressures below 35 GPa, is broken down at high pressures above 35 GPa, when keeping temperature at 2000 K. It should be noted that the breakdown pressure equals the reversion pressure of component self-diffusion coefficients. Therefore, we can presume that the pressure dependences of the ratio between component self-diffusion coefficients can be used to predict the breakdown pressure of SER along isotherm. This feature is similar to the case under ambient pressure that the breakdown temperature can be predicted quantitatively by the divergence behavior of component self-diffusion.[5,19,20] As mentioned above, at ambient pressure, numerous studies have demonstrated that the breakdown temperature of SER is about two times the glass-transition temperature $T_{\rm g}$ for metallic liquids. It is interesting to detect the proposition under high-pressure condition. Our calculated $T_{\rm g}$, determined by intersecting two linear parts of the potential energy curve, is about 1165.2 K at 35 GPa for Cu$_{50}$Zr$_{50}$. In other words, the reduced breakdown temperature of the SER is about 1.72 (2000 K/1165.2 K) at 35 GPa, which is slightly smaller than 2.0. This is unsurprising since the Cu$_{50}$Zr$_{50}$ liquids will become more fragile with higher pressure,[15] and the reduced breakdown temperature will decrease with increasing fragility as previously reported.[8]
cpl-37-7-076201-fig2.png
Fig. 2. The product $D_{\rm Zr}\eta R_{\rm Zr}$ versus pressure at 2000 K for Cu$_{50}$Zr$_{50}$ melt. Lines are guides for the eyes.
cpl-37-7-076201-fig3.png
Fig. 3. Behavior of the first peak position $r_{\alpha \alpha}$ in the partial pair distribution function $g_{\alpha \alpha }(r)$ versus pressure along isotherm (a), and pressure dependence of the first peak amplitude $g(r_{1})_{\alpha \alpha}$ in the partial pair distribution function $g_{\alpha \alpha }(r)$ (b).
Following the conventional explains, the abnormal phenomena of self-diffusion coefficients can be traced back to the structure anomalies.[20–26] As the most commonly used parameter for describing the structure of melts and the only one can be measured from experiments, the pair distribution function (or structure factor) was used to analyze the pressure dependences of atomic structure in Cu$_{50}$Zr$_{50}$ melts. For a binary alloy of components $\alpha$ and $\beta$, the first peak position $r_{\alpha \alpha}$ of the partial pair distribution function $g_{\alpha \alpha }(r)$ is the most frequently used representation of the effective diameter of $\alpha$ components. Figure 3(a) shows the pressure dependence of the effective diameters of Cu and Zr along 2000 K isotherm. We can find that, at given temperature, the effective diameters of both Zr and Cu decrease as pressure increases. This means that atomic pairs become more densely packed when the pressure increases. This is the typical characteristic of liquids along compressing. In addition, the effective diameter of Zr is much larger than that of Cu. In Fig. 3(b), we also show the first peak intensity $g(r_{1})_{\rm CuCu}$ and $g(r_{1})_{\rm ZrZr}$ against pressure. We note that, as is expected, $g(r_{1})_{\rm CuCu}$ increases as pressure increases. However, the pressure-dependent behavior of $g(r_{1})_{\rm ZrZr}$ is an unusual, it increases (not decreases, but small in spite of the increment) at low pressures, and decreases at high pressure as pressure increases. This result can be partially verified by the previous study.[15] This is an indication that the pressure dependences of the topological order surrounding Cu and Zr are indeed different. Therefore, a careful consideration needs to be made to determine the structure difference surrounding Cu and Zr. The structure orders around Cu and Zr atoms were further characterized by employing the Voronoi tessellation. In the Voronoi tessellation analysis, a Voronoi polyhedron (VP) is denoted by a set of indices ($n_{1}$, $n_{2}$, $n_{3}$,$\ldots$, $n_{i}$,$\ldots$), where $n_{i}$ is the number of faces with $i$ vertices, and their sum $\sum n_{i}$ corresponds to the coordination number (CN) of the center atom. Four indices ($n_{3}$, $n_{4}$, $n_{5}$, $n_{6}$) are sufficient for our purpose in this work. Previous works have shown that icosahedral clusters (including ideal icosahedron and icosahedral-like polyhedrons) and the mixed polyhedrons are more abundant in the stable and supercooled liquids.[24–30] Figure 4 shows the fractions of the most popular VPs around Cu (a) and Zr (b) atoms in the Cu$_{50}$Zr$_{50}$ as functions of pressure. As is expected, the most popular VPs in liquid make up a smaller fraction of the total.[29,30] Predominant types of Cu-centered VPs are ideal icosahedron, icosahedral-like polyhedrons and the mixed polyhedrons, of which the CN are in the range of 11–14 such as $\langle 0, 0, 12,0\rangle $, $\langle 0, 2, 8,1\rangle $, $\langle 0, 2, 8,2\rangle $, $\langle 0, 2, 8,4\rangle $, $\langle 0, 1, 10,2\rangle $, $\langle 0, 3, 6,3\rangle $, and $\langle 0, 3, 6,4\rangle $.[27] Predominant types of Zr-centered VPs are icosahedral-like polyhedrons and mixed polyhedrons, of which the CN are in the range of 13–16 such as $\langle 0, 2, 8,4\rangle $, $\langle 0, 2, 8,5\rangle $, $\langle 0, 2, 8,6\rangle $, $\langle 0, 1, 10, 3\rangle $, $\langle 0, 1, 10, 4\rangle $, $\langle 0, 1, 10, 5\rangle $ and $\langle 0, 3, 6,4\rangle $. The fraction of ideal icosahedron around the Zr atom with $\langle 0, 0, 12, 0\rangle $ index is very small. Around Cu atoms, the pressure dependences of the fractions of VPs with small CN ($\langle 0, 2, 8,1\rangle $, $\langle 0, 2, 8,2\rangle $, and $\langle 0, 3, 6,3\rangle $) can be separated into three pressure regimes. In regime I with pressures lower than 20 GPa, the fractions increase rapidly. In pressure regime II between 20 and 45 GPa, the increase slows down; In regime III with pressures higher than 45 GPa, there is a dramatic decrease with increasing pressure. As shown in Fig. 4(b), similar to the behavior of VPs with small CN around the Cu atom, the pressure dependences of the VPs around Zr atoms can also be separated into three pressure regimes at the same turning points. Moreover, the fractions of VPs with larger CN ($\langle 0, 2, 8, 4\rangle $, $\langle 0, 1, 10, 2\rangle $, and $\langle 0, 3, 6, 4\rangle $) around Cu show different pressure dependences. They increase with pressure over the whole pressure range. In Fig. 4(c), we plot the ratio between the total fractions of the seven predominant types of Zr-centered VPs and that of the seven predominant types of Cu-centered VPs (VP$_{\rm Zr}$/VP$_{\rm Cu}$) versus pressure. Similar to the pressure dependence of $D_{\rm Cu}/D_{\rm Zr}$ shown in the inset of Fig. 1(a), the pressure-dependent behavior of the ratio VP$_{\rm Zr}$/VP$_{\rm Cu}$ can be separated into two pressure regimes, a low decrease at low pressures below 35 GPa and a more rapid one at high pressures above 35 GPa. Considering the fact that the component diffusion coefficients will increased with the decreasing VPs,[24–26] we suggest that the reversal of the component diffusion coefficients seems to be in accordance with the different pressure-dependent behavior of VPs around individual components. Moreover, the relative total fractions of predominant VPs around components can also be used to explain and predict the validity of SER along the isotherm.
cpl-37-7-076201-fig4.png
Fig. 4. The pressure dependences of the fractions of predominant VPs in Cu$_{50}$Zr$_{50}$ around Cu (a) and Zr (b), and the ratio between the total fractions of the predominant types of Zr-centered VPs and that of the predominant types of Cu-centered VPs (VP$_{\rm Zr}$/VP$_{\rm Cu}$) (c). Lines are guides for the eyes.
In conclusion, using MD simulations we have studied the pressure effect on transport coefficients in the glass-forming liquid Cu$_{50}$Zr$_{50}$, and the validity of SER along the isotherm is also detected up to 70 GPa. It is found that the pressure-dependent behavior of viscosity can be described well by a single exponential relation over the full range of pressure, while there is a change in the pressure dependence of component diffusion at reversion pressure. At low pressures below the reversion pressure, the self-diffusion coefficients of Cu are always larger than those of Zr, and the difference between them becomes smaller with increasing pressure. At high pressures above the reversion pressure, the feature becomes polar opposite, the self-diffusion coefficients of Cu are always smaller than those of Zr, and the difference between them becomes larger with increasing pressure. The SER relation between self-diffusion coefficients, effective diameter, and the viscosity works well at low pressures and starts to be violated at the breakdown pressure. For glass-forming liquid Cu$_{50}$Zr$_{50}$ along the 2000 K isotherm, the reverse pressure equals the breakdown pressure of SER and is about 35 GPa. At the same pressure, there is also a change in the slope of the ratio between the relative total fractions of the seven predominant VPs around components to pressure. Therefore, the validity of SER and the reverse behavior of component diffusion coefficients can be predicted by the pressure dependence of the relative total fractions of predominant VPs around individual components. References Supercooled dynamics of glass-forming liquids and polymers under hydrostatic pressureImpact of the application of pressure on the fundamental understanding of glass transitionA molecular dynamics examination of the relationship between self-diffusion and viscosity in liquid metalsHigh temperature breakdown of the Stokes-Einstein relation in a computer simulated Cu-Zr meltRevisiting the Stokes–Einstein relation for glass-forming meltsAtomic-scale dynamics of a model glass-forming metallic liquid: Dynamical crossover, dynamical decoupling, and dynamical clusteringStructural signatures evidenced in dynamic crossover phenomena in metallic glass-forming liquidsCorrelation between Fragility and the Arrhenius Crossover Phenomenon in Metallic, Molecular, and Network LiquidsAppearance of a fractional Stokes–Einstein relation in water and a structural interpretation of its onsetStructural origin of fractional Stokes-Einstein relation in glass-forming liquidsSignatures of fragile-to-strong transition in a binary metallic glass-forming liquidPressure and temperature dependence of viscosity and diffusion coefficients of a glassy binary mixturePressure effects on structure and dynamics of metallic glass-forming liquidHydrodynamic Relaxation of an Electron Plasma to a Near-Maximum Entropy StateAnomalous structure-property relationships in metallic glasses through pressure-mediated glass formationFast Parallel Algorithms for Short-Range Molecular DynamicsFive-fold symmetry as indicator of dynamic arrest in metallic glass-forming liquidsTransport properties and Stokes-Einstein relation in a computer-simulated glass-forming Cu 33 . 3 Zr 66 . 7 meltTransport properties and abnormal breakdown of the Stokes-Einstein relation in computer simulated Al72Ni16Co12 metallic meltCorrelation between Dynamic Heterogeneity and Medium-Range Order in Two-Dimensional Glass-Forming LiquidsLocal order and dynamic properties of liquid Au x Si 1− x alloys by molecular dynamics simulationsInterplay between the structure and dynamics in liquid and undercooled boron: An ab initio molecular dynamics simulation studyRelationship between structure, dynamics, and mechanical properties in metallic glass-forming alloysLocal order and dynamic properties of liquid and undercooled Cu 55 Hf 45 and Cu 62 Hf 38 alloys by ab initio molecular dynamicsStructural aspects of the Stokes-Einstein relation breakdown in high temperature meltsAtomic packing and medium-range order in Ni 3 Al metallic glassAb initio molecular dynamics simulation of the atom packing and density of Al-Ni amorphous alloysEvaluation of the local environment for nanoscale quasicrystal formation in Zr 80 Pt 20 glassy alloy using Voronoi analysisExperimental and computer simulation determination of the structural changes occurring through the liquid–glass transition in Cu–Zr alloys
[1] Roland C M, Hensel-Bielowka S, Paluch M and Casalini R 2005 Rep. Prog. Phys. 68 1405
[2] Ngai K L and Capaccioli S 2008 J. Phys.: Condens. Matter 20 244101
[3] Lü Y, Cheng H and Chen M 2012 J. Chem. Phys. 136 214505
[4] Han X J, Li J G and Schober H R 2016 J. Chem. Phys. 144 124505
[5] Cao Q, Wang P and Huang D 2020 Phys. Chem. Chem. Phys. 22 2557
[6] Jaiswal A, Egami T and Zhang Y 2015 Phys. Rev. B 91 134204
[7] Hu Y C, Li F X, Li M Z, Bai H Y and Wang W H 2016 J. Appl. Phys. 119 205108
[8] Jaiswal A, Egami T, Kelton K F, Schweizer K S and Zhang Y 2016 Phys. Rev. Lett. 117 205701
[9] Xu L, Mallamace F, Yan Z, Starr F W, Buldyrev S V and Eugene Stanley H 2009 Nat. Phys. 5 565
[10] Pan S, Wu Z W, Wang W H, Li M Z and Xu L 2017 Sci. Rep. 7 39938
[11] Lad K N, Jakse N and Pasturel A 2012 J. Chem. Phys. 136 104509
[12] Mukherjee A, Bhattacharyya S and Bagchi B 2002 J. Chem. Phys. 116 4577
[13] Hu Y, Guan P, Wang Q, Yang Y and Bai H 2017 J. Chem. Phys. 146 024507
[14] Cheng Y Q, Ma E and Sheng H W 2009 Phys. Rev. Lett. 102 244501
[15] Ding J, Asta M and Ritchie R O 2016 Phys. Rev. B 93 140204(R)
[16] Plimpton S 1995 J. Comput. Phys. 117 1
[17]Allen M P and Tildesley D J 1987 Computer Simulation of Liquids (Oxford: Clarendon Press)
[18] Hu Y C, Li F X, Wang W H, Li M Z, Bai H Y and Wang W H 2015 Nat. Commun. 6 8310
[19] Han X J and Schober H R 2011 Phys. Rev. B 83 224201
[20] Zhou Y H, Han X J and Li J G 2019 J. Non-Cryst. Solids 517 83
[21] Kawasaki T, Araki T and Tanaka H 2007 Phys. Rev. Lett. 99 215701
[22] Jakse N, Nguyen T L T and Pasturel A 2012 J. Chem. Phys. 137 204504
[23] Jakse N and Pasturel A 2014 J. Chem. Phys. 141 234504
[24] Cheng Y Q, Sheng H W and Ma E 2008 Phys. Rev. B 78 014207
[25] Jakse N, Nguyen T L T, Pasturel A, Jakse N, Nguyen T L T and Pasturel A 2013 J. Appl. Phys. 114 063514
[26] Li C H, Luan Y W, Han X J and Li J G 2017 J. Non-Cryst. Solids 458 107
[27] Trady S, Hasnaoui A and Mazroui M 2017 J. Non-Cryst. Solids 468 27
[28] Yu C, Hui X, Chen X, Liu X, Lin D, Liu Z and Chen G 2010 Sci. Chin. Technol. Sci. 53 3175
[29] Saida J, Itoh K, Sato S and Imafuku M 2009 J. Phys.: Condens. Matter 21 375104
[30] Mendelev M I, Kramer M J, Ott R T, Sordelet D J, Besser M F, Kreyssig A, Goldman A I, Wessels V, Sahu K K, Kelton K F, Hyers R W, Canepari S and Rogers J R 2010 Philos. Mag. 90 3795