Open Access
Add to BibSonomy
Abstract
The laser damage induced by nano-absorbing centers generally results in a local high temperature and pressure environment, leading to denser phases and complex hydrodynamic processes. Here we parameterize the metal-organic framework force field to overcome the notorious unphysical agglomeration at small atomic distance in a Buckingham term. The structure and shock properties of amorphous silica are predicted well by the parameterized force field. By avoiding the Ewald summation of long-range coulomb interaction, the periodic boundary condition is not in such demand that the computational efficiency is greatly improved. The parameterized force field implicates a prospect for the atomic investigation of laser-induced hydrodynamic processes around the free surface or interface.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Silica is composed by the most abundant elements on the Earth and widely used from architectural glass to high-power laser components. The laser beam in the high-power laser system is usually manipulated by various coatings or glasses consisting of amorphous silica [1–5]. The energy deposition initiated by nano-absorbing centers generally results in a local high temperature and pressure environment, leading to the denser phases and complex hydrodynamic processes [6,7]. Many efforts have been devoted to clarifying the shock dynamics and topological properties of amorphous silica under the extreme conditions of laser irradiation in both experimentally and theoretically [8–13]. There is always a necessity to explore the topological and kinetic properties at large scale by the classical molecular dynamic (MD) simulation with force field that could well reproduce the experimental results. A state-of-art pair potential model proposed by van Beest, Kramer, and van Santen (BKS) [14] has been reported to well reproduce properties of various silica phase in both the standard and shock conditions [11,13,15–20]. Although various modifications have been tried to solve the notorious unphysical agglomeration at small atomic distance with an arbitrary spring function, the fundamental physical meaning is still lacking. Moreover, the long-range Coulombic interaction is often calculated via reciprocal space techniques (such as the Ewald method) with the periodic boundary conditions (PBCs), demanding enormous computational resource. PBCs means that a vacuum slab is required to simulate the laser-induced hydrodynamic process around the free surface or interface, which would bring in the extra computational burden. Recently, Bureekaew et al. [21] presented a metal-organic framework force field (MOF-FF) with consistent fashion from first principles data. The unphysical dispersion term could be corrected analog to the Grimme correction used in density functional theory (DFT) [22]. Moreover, the damped shifted force (DSF) model has been reported to be a promising method which considers the long-range interaction without the computational resource from the reciprocal space method [23]. Herein, a combination of MOF-FF and DSF model is parameterized to reproduce the structure and shock properties of amorphous silica, which might provide an efficient yet versatile solution to clarify the kinetic event during laser-induced damage.
2. Parameterization of MOF-FF
The BKS potential combines the Buckingham and Coulombic term. Firstly, the standard Buckingham term denoted as $\phi _{ij}^{buck}(r_{ij})$ is represented as a function of atomic distance $r_{ij}$ in a short-range [14]
where $A_{ij}$ and $C_{ij}$ are force constants, $b_{ij}$ is an ionic-pair dependent parameters. These parameters are summarized in Table 1 [14]. The Buckingham term vanishes when the atomic distance is beyond the cutoff distance $r_c$. In Fig. 1, the Buckingham term shows the notorious unphysical agglomeration at small atomic distance, although the atoms attempt to repel each other actually when they are close enough.
Fig. 1. Comparison of pair potentials for Si-O and O-O interactions. The dash and solid line depicts the original and modified Buckingham term, respectively. The Buckingham term in each version is truncated at the cutoff distance $r_c$ = 6 Å.
Table 1. Parameters for the standard and modified Buckingham term.
Then, the two-body interaction of the MOF-FF model denoted as $\phi _{ij}^{MOF-FF}$ modifies the dispersion term of Buckingham potential with a damped coefficient [21]
Finally, the Coulombic interaction is considered with atomic charges $q_{\textrm{Si}}$ = +2.4$e$ and $q_{\textrm{O}}$ = −1.2$e$ due to the contribution of ionic bond, in which $e$ is the elementary charge. Here the DSF model is used to accelerate the simulation without the computational burden from the reciprocal space method [23]. The MOF-FF model uses the spherical Gaussian type charge distributions to dampen the electrostatic potential for high charges at close distances [21]. The electrostatic potential $\phi _{ij}^{coul}$ in the MOF-FF model is given as [21]
In this work, a model combining the Coulombic term $\phi _{ij}^{coul}$ and modified Buckingham term $\phi _{ij}^{MOF-FF}$ is used to construct the amorphous silica and calculate the corresponding properties, which will be named as MOF-FF in the remained context. The smooth, damping, and cutoff parameters are parameterized as $\kappa$ = 0.5, $\alpha$ = 1.5 Å$^{-1}$, and $r_c$ = 6 Å, respectively. All the atomic simulations are implemented in the LAMMPS package (Version 3Mar20) [24].
3. Structure properties
To construct the amorphous silica with the new parameterized model via the standard melt-quenching schematic, the $\beta$-cristobalite structure suggested by Wright et al. is used as the raw materials [25]. The unit cell of $\beta$-cristobalite consists of 24 atoms and the duplication of 10 $\times$ 10 $\times$ 10 (24,000 atoms) is built as the simulation box with PBCs. The isothermal-isobaric (NPT) ensemble is used and the timestep is 1.6 fs in this section. Before initiating the melt-quenching schematic, the constructed simulation box is examined by relaxing for 200 ps at 300K. The simulated density is about 2.192 $\pm$ 0.002 g/cm$^3$, which agrees well with the experimental density (2.201 g/cm$^3$) of the $\beta$-cristobalite [25]. To calculate the structure properties, 100 independent configurations are generated with 1 ps interval under the equilibrium conditions for both the melt and solid silica.
At the first step of the melt-quenching schematic, the built simulation box is heated for 200 ps at 6060 K. The simulated melting point of silica is well known to be higher than the actual temperature for most potentials. The short-range order is dominated by the arrangement of atoms within the $\left [\rm {SiO_4}\right ]$ tetrahedron, which could be revealed from the partial radial distribution function (RDF) and the O-Si-O bond angle distribution (BAD). In Fig. 2(a) and (b), the melt silica still maintains the typical short-range order around the experimental values ( $r_\textrm {Si-O}$ = 1.62 Å, $r_\textrm {O-O }$ = 2.65 Å, $r_\textrm {Si-Si}$ = 3.12 Å [26], $\theta _\textrm {O-Si-O}$ = 109.7$^\circ$, and $\theta _\textrm {Si-O-Si}$ = 152.0$^\circ$ [27]) of fused silica at 300 K. The broad distribution of the structure properties is resulted from the large disorder of relative orientation among the $\left [\rm {SiO_4}\right ]$ tetrahedra in the diffusive melt.
Fig. 2. (a) Partial radial distribution function and (b) bond angle distribution for silica melt at 6060 K. (c) Silica density during the quenching process. (d) Total correlation function of amorphous silica. (e) Partial radial distribution function and (f) bond angle distribution for the quenched amorphous silica at 300 K. The dash lines in (a), (b), (e), and (f) label the corresponding experimental values of amorphous silica at 300 K.
The melt is then quenched to 300 K gradually with a cooling rate 0.6 K/ps, producing a solid density 2.222 $\pm$ 0.001 g/cm$^3$ within 1$\%$ of the experimental values (2.200 g/cm$^3$). In Fig. 2(c), the typical anomalous density of amorphous silica is also observed during quenching process. The simulated total correlation function in Fig. 2(d) agrees well with that from neutron scattering in both peak position and intensity [28]. The first peak of total correlation function from MOF-FF model is narrower and larger than the measured results, indicating a better short-range order of the $\left [\rm {SiO_4}\right ]$ tetrahedron. The excellent short-range order can be also depicted by the partial RDF of Si-O (Fig. 2(e)), in which the peak intensity at $r_\textrm {Si-O}$ agrees well with the previous results by other empirical potentials [18,20]. The cutoff distance for Si-O bond is chosen as 1.83 Å according to the partial RDF of Si-O, which is used to build the network consisted of Si and O atoms. Both the short-range order of the $\left [\rm {SiO_4}\right ]$ tetrahedron and the orientation order among the $\left [\rm {SiO_4}\right ]$ tetrahedra are significantly improved after quenching (Fig. 2(e) and (f)). The BAD of the quenched amorphous silica converges to the experimental values with a narrow distribution. The O-Si-O ranges from 90$^\circ$ to 130$^\circ$, while the Si-O-Si ranges from 120$^\circ$ to 180$^\circ$. The simulated BAD agrees well with that of other verified empirical potentials [18,20,29]. The excellent BAD indicates a ideal network linked by the corner-sharing tetrahedra with few structure defects.
To evaluate the topological properties of the Si–O network constructed from the parameterized MOF-FF model, the coordination and ring size distribution is analyzed. The Si and O atoms with coordination number $m$ are labeled as Si$^m$ and O$^m$, respectively. For an ideal $\left [\rm {SiO_4}\right ]$ tetrahedron, the network should be only consisted of Si$^4$ and O$^2$. Therefore, the atoms with other coordination numbers are considered as coordination defects. Due to the excellent short-range order and BAD, the proportion of Si$^4$ and O$^2$ is about 99.68$\%$ and 99.49$\%$, respectively, indicating pretty low density defects. The proportion of over-coordinated atoms is almost equal to that of low-coordinated atoms. In addition, the medium-range order is generally depicted by the distribution of ring size $n$. The $n$-ring means a ring consisted of $n$ Si atoms. A $n$-ring is identified with the King’s shortest-path criterion [30,31]. The accumulated number of $n$-ring is averaged by the total number of Si atoms. The distribution of $n$-ring is illustrated in Fig. 3. As expected, the 6-ring dominates the topological structure of the amorphous silica, whereas the number of the larger or smaller rings decreases rapidly. The simulated topological properties agree well with the previous reported results [18,20,29,30].
The density of amorphous silica is analyzed at different static pressures ($P$) to predict the equation of state for cold compression. As shown in Fig. 4(a), the density increases rapidly above 7 GPa and the growth rate slows down above 20 GPa. The pressure-dependent density agrees well with the experimental results in the whole range [32,33]. The minor deviation of density at the higher pressure might be related with the selection of cutoff distance in this work. The density of silica family predicted by the standard BKS potential has been reported to be sensitive with the cutoff distance [34]. The high pressure phase of silica generally requires a slightly larger cutoff distance to exactly reproduce the experimental density. Since the demonstrated MOF-FF here is modified from the BKS potential, the similar deviation could be expected. Nevertheless, the relative error of density at high pressure is small enough and shows little influence on the formation of high pressure phase.
Fig. 4. Pressure-dependent (a) density, (b) Si$^m$, (c) ring size ($n$ = 3 – 6), and (d) ring size ($n$ = 7 – 9). The experimental results are from Ref. [32,33]
The pressure-dependent coordination number of Si atom is shown in Fig. 4(b). The four-fold coordinated Si atoms dominates below 7 GPa. The proportion of Si$^4$ decreases rapidly and transforms into the Si$^5$ and Si$^6$ above 7 GPa, which is corresponding to the rapid increase of density. The proportion of Si$^5$ reaches the maximum at 20 GPa, while the proportion of Si$^4$ almost vanishes. Further increase the pressure will enable the Si$^5$ $\rightarrow$ Si$^6$ transformation to be the dominant process, leading to the decrease of Si$^5$. Since the available Si$^5$ is only a part of the original Si$^4$, the increase of density due to Si$^5$ $\rightarrow$ Si$^6$ transformation slows down. The proportion of Si$^6$ could surpass the Si$^5$ above 35 GPa. The existence of the Si$^5$ and Si$^6$ usually indicates the formation of high pressure phase such as coesite and stishovite. Since the Si$^5$ and Si$^6$ almost initiates simultaneously, the weak Si$^4$ $\rightarrow$ Si$^6$ transformation could also be predicted by the parameterized potential. Nevertheless, the Si$^5$ is the major intermediary to bridge the transformation from Si$^4$ to Si$^6$. In addition, the pressure-dependent density could also be revealed by the variation of ring size during the cold compression. As shown in Fig. 4(c) and (d), the small size ($n<6$) ring increases whereas the large size ($n>6$) ring decreases. The characteristic pressure follows those of Si$^m$, indicating an intimate relationship between small size ring and high pressure phase. The simulated pressure-dependent coordination number and ring size agrees well with those of the previous results [36].
In addition to the static structure properties, the dynamic structure properties of amorphous silica is also essential in reproducing the shock response during laser-material interaction. One of the typical dynamic structure properties is the atomic diffusion above the melting point. Here the mean-square displacement (MSD) of the Si and O atoms at 4000 K – 6000 K is used to extract the diffusion coefficients ($D$). For example, the MSD of atoms at 6000 K varies linearly with time in Fig. 5(a) and the slope indicates the diffusion coefficient. Similar phenomenon could be observed at other temperatures. The Si atom diffuses more slowly than the O atom due to the larger mass. The diffusion coefficients at various temperatures ($T$) are summarized in Fig. 5(b) together with the previous DFT results [35], which are linearly fitted according to the Arrhenius law. The slope of the fitting line indicates the diffusion activation energies ($Q$), where $Q_\textrm {Si}$ = 2.87 eV and $Q_\textrm {O}$ = 2.73 eV. The predicted diffusion activation energies are slightly smaller than the DFT results, which might be due to the non Arrhenius behavior at the higher temperature [11,35].
Fig. 5. (a) MSD of Si and O atoms at 6000K. (b) Diffusion coefficients (symbols) versus temperature at 300K. Lines represent the Arrhenius law. The DFT results are from Ref. [35].
4. Shock properties
During the explosive or laser-driven hydrodynamic processes, the expansion of the core region under extreme conditions leads to the movement of core/cold region boundary, known as particle or piston velocity. The compression of the core/cold region boundary drives the shockwave propagation in the cold materials, leading to a wide range of energy deposition [7]. The thermodynamic properties of a material under shock compression is usually described as Hugoniostat equation of state (HEOS), which is the relationship among shock pressure along the propagating direction ($\sigma _{zz}$), shock temperature ($T$), piston velocity ($u_p$), and shock velocity ($u_s$). Here the constant-stress Hugoniostat method is used to satisfy the Rankine-Hugoniot jump conditions for steady shocks [37]. The shock load is described as $\sigma _{zz}$ or $u_p$, while the shock response is described as $T$ or $u_s$. Generally, there is a Hugoniostat elastic limit (HEL) of shock load to separate the shock response into the elastic and shock region. The shock region could be further classified as elastic-plastic transition region and plastic region depending on the shock load. Compared with the traditional nonequilibrium molecular dynamics, reaching the steady shock state with the Hugoniostat method is far more time-saving. To get accurate shock properties, a timestep of 1 fs is used in this section with the MOF-FF model.
It takes a much longer time to reach the steady shock state from the initial amorphous silica under the shock load in the elastic-plastic transition region. The shock density and temperature has been examined under various shock pressure so that a shock duration of 60 ps is used to reach the steady state of the thermodynamic properties. The shock pressure ranges from 1 GPa to 100 GPa. The shock temperature stays at the room temperature below 12 GPa and increases almost linearly with the shock pressure above 12 GPa (Fig. 6(a)). The shock temperature surpasses the melting point of amorphous silica easily above 38 GPa. The pressure-dependent shock temperature agrees well with the available experimental results [38]. The liquid silica above the melting point has been reported to behave as the metal for the typical 351 nm laser, leading to the solid-absorptive front [41–43]. As expected, the piston velocity increases as the shock pressure monotonously (Fig. 6(b)). At the pressure less than 12 GPa, the shock velocity decreases gradually with the shock pressure. Considering negligible temperature variation, there is no plastic work during this stage. As a result, the shockwave is an elastic wave below 12 GPa. Above 12 GPa, the shock velocity increase gradually with the shock pressure. The rapid increase of temperature in this region indicates the plastic work during the compress process. Therefore, the HEL locates around 12 GPa, which is consistent with the reported values in experimentally [39,40]. Above 38 GPa, the piston velocity is greater than 2.82 km/s and the corresponding shock velocity is comparable or larger than the nominal acoustic speed ($\sim$ 6 km/s), indicating the plastic work. As shown in Fig. 6(b), the Hugoniostat is much more complex between HEL and 38GPa. The shock-induced plastic work becomes obvious due to the temperature rise and non-zero shock velocity. Moreover, the shock velocity is less than the acoustic velocity. Therefore, the Hugoniostat in this pressure interval comprises of the elastic and plastic waves.
Fig. 6. (a) Shock temperature versus shock pressure. (b) HEOS of $u_s$ and $u_p$. The experimental results from Lyzenga et al.[38], Renou et al. [39] and Sugiura et al. [40] are listed. The mapping relationship between shock pressure and piston velocity is also given.
The density or compressibility ($V/V_0$) depending on shock pressure is also investigated with the MOF-FF model in Fig. 7(a). Below the HEL, the volume decreases linearly with the shock pressure. Furthermore, the volume decreases rapidly as long as the shock pressure surpasses the HEL. Nevertheless, the volume is stable when the shock pressure is above 38 GPa, indicating a negligible compressible state. The pressure-dependent compressibility agrees well with the experimental results [44]. In Fig. 7(b), the compressibility below the HEL is also comparable with the laser-driven shock experimental and theoretical results by Renou et al. [39]. The shock pressure dependent coordination is also summarized in Fig. 7(c). Below the HEL, the fraction of Si$^4$ is stable, explaining the elastic interaction of shockwave. The fraction of Si$^4$ decreases rapidly and transforms into Si$^3$ and Si$^5$. The variation of Si$^m$ is stable above 38 GPa, which is responsible for the incompressible state. Similar phenomenon could be also observed on the shock pressure dependent distribution of ring size in Fig. 7(d). The increase of the smaller size ($n$ < 5) ring is generally related with the high-pressure phase. The increase of the larger size ($n$ > 8) ring is related with the fracture formation during the shock propagation [12]. Obviously, the elastic-plastic transition during shock compression is intimately related with the topological reconstruction in the window of HEL – 38 GPa. The negligible compressibility could be also revealed by the insensitive topological properties above 38 GPa.
Fig. 7. (a) Shock pressure dependent compressibility. The experimental result is also listed [44]. (b) Compressibility below the HEL. The experimental result from Renou et al. is also listed [39]. shock pressure dependent (c) coordination and (d) distribution of ring size.
5. Discussion
The larger-size ring decomposes into the smaller rings and the low-coordination Si atoms transforms into the over-coordination in the cold compression. The topological variation of shock compression is significantly different from that of cold compression. The shock temperature is sufficiently high due to the shock pressure, which exerts a resistance against the compression. The topological variation is caused by the slightly reconstruction among the standard $\left [\rm {SiO_4}\right ]$ tetrahedra. Both the number of the smaller and bigger size ring increases at sufficient high shock pressure. The Si$^3$ also increases with the shock pressure, which is not observed in the cold compression. The occurrence of the Si$^3$ might be related with the formation or the larger size ring during laser shock. The nominal Si$^6$ of the high-pressure phase is not observed during the shock process, which is similar with the reports in Ref. [12,39]. The formation of typical high-pressure phase during the shock propagation is a relatively slow kinetic process of nucleation, in which the timescale is one or two order of magnitudes longer than that of the current work [11]. Therefore, the parameterized MOF-FF could well reproduce the structure and shock properties of amorphous silica.
In addition to the excellent reproduction of the experimental results, the computation speed of the parameterized MOF-FF is much faster than that of the standard BKS potential. To compare the computation speed, the constructed amorphous silica (24,000 atoms) in Section 3. is relaxed for 100 ps with a timestep of 1.6 fs under NPT ensemble. The speed of the MOF-FF and BKS potential is about 33.5 ns/day and 5.2 ns/day under the same computational resource (Intel Xeon E5-2692 v2, 24 cores), respectively. Similar efficiency enhancement could be also observed for the quenching process and Hugoniostat simulation. The excellent reproduction and efficiency enhancement benefits from the truly consideration of long-range Coulombic interaction without the burden from the Ewald summation. By avoiding the Ewald summation, PBCs are not demanded any more, which provides the possibility to directly study the laser-induced hydrodynamic processes around the free surface or interface. Again, the computational burden would be further decreased significantly since the general vacuum slab in the PBCs is not required. Moreover, the Tersoff potential developed by Munetoh et al. has been reported as the fastest model for the Si–O system since the Coulomb interaction term is excluded and the charge-transfer between Si and O atoms is considered with a parameter $\chi _{Si-O}$[45]. With the same benchmark, the speed of the Tersoff potential is about 59.3 ns/day. However, the amorphous silica constructed by the Tersoff potential with the same cooling rate shows more coordination defects. The proportion of Si$^4$ and O$^2$ is about 87.10$\%$ and 96.46$\%$, respectively. Therefore, the parameterized MOF-FF in this work not only shows the high computation speed close to that of the Tersoff potential but also predicts the topological properties well as those of the BKS potential.
6. Conclusion
In this work, a combination of MOF-FF and DSF model is parameterized to describe the mixed ionic-valance bond and long-range Coulombic interaction interaction among atoms in amorphous silica. The structure and shock properties of amorphous silica are well reproduced with the new parameterized potential. As expected, the free of computation burden in reciprocal space results in significant enhancement of computation efficiency. Without the requirement of PBCs, the parameterized potential could be directly employed to study the large-scale hydrodynamic and kinetic processes along with the fracture, interface and free surface during the laser-induced damage of silica components. Since the hydrodynamic and kinetic process is critical to the formation of high-pressure phase and topological defects, our work would also be meaningful for the lifetime and mitigation of laser components.
Funding
Shanghai Science and Technology Commission (20511107400); National Key Research and Development Program of China (2018YFE0118000); Shanghai Sailing Program (18YF1426400); Youth Innovation Promotion Association of the Chinese Academy of Sciences (2020253); NSAF Joint Fund (U1930119); National Natural Science Foundation of China (11904376, 61805263, 61975220).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. S. Ly, N. Shen, R. A. Negres, C. W. Carr, D. A. Alessi, J. D. Bude, A. Rigatti, and T. A. Laurence, “The role of defects in laser-induced modifications of silica coatings and fused silica using picosecond pulses at 1053 nm: I Damage morphology,” Opt. Express 25(13), 15161 (2017). [CrossRef]
2. H. Wang, H. Qi, J. Zhao, Y. Chai, B. Wang, and J. Shao, “Intrafilm separation of solgel film under nanosecond irradiation,” Appl. Opt. 54(35), 10504 (2015). [CrossRef]
3. J. Zhao, H. Qi, H. Wang, H. He, and W. Zhang, “Scattering analysis for random antireflective structures on fused silica in the ultraviolet,” Opt. Lett. 40(22), 5168 (2015). [CrossRef]
4. H. Wang, H. Qi, W. Zhang, J. Sun, Y. Chai, F. Tu, J. Zhao, Z. Yu, B. Wang, M. Zhu, K. Yi, and J. Shao, “Suppression of nano-absorbing precursors and damage mechanism in optical coatings for 3ω mirrors,” Opt. Lett. 41(6), 1209 (2016). [CrossRef]
5. H. Wang, H. Qi, J. Zhao, B. Wang, Y. Chai, Z. Yu, and J. Shao, “Origin of the plasma scalds in dielectric coatings induced by 1ω laser,” Appl. Phys. Lett. 108(14), 141603 (2016). [CrossRef]
6. P. Grua, L. Lamaignère, M. Chambonneau, R. Courchinoux, and J. Néauport, “Nanosecond laser damage initiation at 0.35 μm in fused silica,” Opt. Lett. 43(11), 2692 (2018). [CrossRef]
7. P. DeMange, R. A. Negres, R. N. Raman, J. D. Colvin, and S. G. Demos, “Role of phase instabilities in the early response of bulk fused silica during laser-induced breakdown,” Phys. Rev. B 84(5), 054118 (2011). [CrossRef]
8. A. Salleo, S. T. Taylor, M. C. Martin, W. R. Panero, R. Jeanloz, T. Sands, and F. Y. Génin, “Laser-driven formation of a high-pressure phase in amorphous silica,” Nat. Mater. 2(12), 796–800 (2003). [CrossRef]
9. A. E. Gleason, C. A. Bolme, H. J. Lee, B. Nagler, E. Galtier, D. Milathianaki, J. Hawreliak, R. G. Kraus, J. H. Eggert, D. E. Fratanduono, G. W. Collins, R. Sandberg, W. Yang, and W. L. Mao, “Ultrafast visualization of crystallization and grain growth in shock-compressed SiO2,” Nat. Commun. 6(1), 8191 (2015). [CrossRef]
10. C. E. Wehrenberg, D. McGonegle, C. Bolme, A. Higginbotham, A. Lazicki, H. J. Lee, B. Nagler, H.-S. Park, B. A. Remington, R. E. Rudd, M. Sliwa, M. Suggit, D. Swift, F. Tavella, L. Zepeda-Ruiz, and J. S. Wark, “In situ X-ray diffraction measurement of shock-wave-driven twinning and lattice dynamics,” Nature 550(7677), 496–499 (2017). [CrossRef]
11. Y. Shen, S. B. Jester, T. Qi, and E. J. Reed, “Nanosecond homogeneous nucleation and crystal growth in shock-compressed SiO2,” Nat. Mater. 15(1), 60–65 (2016). [CrossRef]
12. A. Kubota, M. J. Caturla, J. S. Stölken, and M. D. Feit, “Densification of fused silica due to shock waves and its implications for 351 nm laser induced damage,” Opt. Express 8(11), 611 (2001). [CrossRef]
13. R. Su, M. Xiang, J. Chen, S. Jiang, and H. Wei, “Molecular dynamics simulation of shock induced ejection on fused silica surface,” J. Appl. Phys. 115(19), 193508 (2014). [CrossRef]
14. B. W. H. van Beest, G. J. Kramer, and R. A. van Santen, “Force fields for silicas and aluminophosphates based on ab initio calculations,” Phys. Rev. Lett. 64(16), 1955–1958 (1990). [CrossRef]
15. L. V. Woodcock, C. A. Angell, and P. Cheeseman, “Molecular dynamics studies of the vitreous state: Simple ionic systems and silica,” J. Chem. Phys. 65(4), 1565–1577 (1976). [CrossRef]
16. S. Tsuneyuki, M. Tsukada, H. Aoki, and Y. Matsui, “First-principles interatomic potential of silica applied to molecular dynamics,” Phys. Rev. Lett. 61(7), 869–872 (1988). [CrossRef]
17. E. Demiralp, T. Çağin, and W. A. Goddard, “Morse stretch potential charge equilibrium force field for ceramics: application to the quartz-stishovite phase transition and to silica glass,” Phys. Rev. Lett. 82(8), 1708–1711 (1999). [CrossRef]
18. T. F. Soules, G. H. Gilmer, M. J. Matthews, J. S. Stolken, and M. D. Feit, “Silica molecular dynamic force fields—A practical assessment,” J. Non-Cryst. Solids 357(6), 1564–1573 (2011). [CrossRef]
19. S. Izvekov and B. M. Rice, “Mechanism of densification in silica glass under pressure as revealed by a bottom-up pairwise effective interaction model,” J. Chem. Phys. 136(13), 134508 (2012). [CrossRef]
20. S. Sundararaman, W.-Y. Ching, and L. Huang, “Mechanical properties of silica glass predicted by a pair-wise potential in molecular dynamics simulations,” J. Non-Cryst. Solids 445-446, 102–109 (2016). [CrossRef]
21. S. Bureekaew, S. Amirjalayer, M. Tafipolsky, C. Spickermann, T. K. Roy, and R. Schmid, “MOF-FF - A flexible first-principles derived force field for metal-organic frameworks,” Phys. Status Solidi B 250(6), 1128–1141 (2013). [CrossRef]
22. S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” J. Chem. Phys. 132(15), 154104 (2010). [CrossRef]
23. C. J. Fennell and J. D. Gezelter, “Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics,” J. Chem. Phys. 124(23), 234104 (2006). [CrossRef]
24. S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117(1), 1–19 (1995). [CrossRef]
25. A. F. Wright and A. J. Leadbetter, “The structures of the β-cristobalite phases of SiO2 and AlPO4,” Philos. Mag. 31(6), 1391–1401 (1975). [CrossRef]
26. “SiO2,” in Handbook of Glass Data, vol. 15 of Physical Sciences DataO. V. Mazurin, M. V. Streltsina, and T. P. Shvaiko-Shvaikovskaya, eds. (Elsevier, 1983), pp. 1–185.
27. R. Zallen, The Physics of Amorphous Solids (John Wiley and Sons, New York, 1983).
28. P. A. Johnson, A. C. Wright, and R. N. Sinclair, “Neutron scattering from vitreous silica II. Twin-axis diffraction experiments,” J. Non-Cryst. Solids 58(1), 109–130 (1983). [CrossRef]
29. J. P. Rino, I. Ebbsjö, R. K. Kalia, A. Nakano, and P. Vashishta, “Structure of rings in vitreous SiO2,” Phys. Rev. B 47(6), 3053–3062 (1993). [CrossRef]
30. S. V. King, “Ring configurations in a random network model of vitreous silica,” Nature 213(5081), 1112–1113 (1967). [CrossRef]
31. D. S. Franzblau, “Computation of ring statistics for network models of solids,” Phys. Rev. B 44(10), 4925–4930 (1991). [CrossRef]
32. O. B. Tsiok, V. V. Brazhkin, A. G. Lyapin, and L. G. Khvostantsev, “Logarithmic kinetics of the amorphous-amorphous transformations in SiO2 and GeO2 glasses under high pressure,” Phys. Rev. Lett. 80(5), 999–1002 (1998). [CrossRef]
33. T. Sato and N. Funamori, “Sixfold-coordinated amorphous polymorph of SiO2 under high pressure,” Phys. Rev. Lett. 101(25), 255502 (2008). [CrossRef]
34. B. J. Cowen and M. S. El-Genk, “On force fields for molecular dynamics simulations of crystalline silica,” Comput. Mater. Sci. 107, 88–101 (2015). [CrossRef]
35. B. B. Karki, D. Bhattarai, and L. Stixrude, “First-principles simulations of liquid silica: Structural and dynamical behavior at high pressure,” Phys. Rev. B 76(10), 104205 (2007). [CrossRef]
36. A. Zeidler, K. Wezka, R. F. Rowlands, D. A. J. Whittaker, P. S. Salmon, A. Polidori, J. W. E. Drewitt, S. Klotz, H. E. Fischer, M. C. Wilding, C. L. Bull, M. G. Tucker, and M. Wilson, “High-pressure transformation of SiO2 glass from a tetrahedral to an octahedral network: a joint approach using neutron diffraction and molecular dynamics,” Phys. Rev. Lett. 113(13), 135501 (2014). [CrossRef]
37. R. Ravelo, B. L. Holian, T. C. Germann, and P. S. Lomdahl, “Constant-stress Hugoniostat method for following the dynamical evolution of shock matter,” Phys. Rev. B 70(1), 014103 (2004). [CrossRef]
38. G. A. Lyzenga, T. J. Ahrens, and A. C. Mitchell, “Shock temperatures of SiO2 and their geophysical implications,” J. Geophys. Res.: Solid Earth 88(B3), 2431–2444 (1983). [CrossRef]
39. R. Renou, L. Soulard, E. Lescoute, C. Dereure, D. Loison, and J.-P. Guin, “Silica glass structural properties under elastic shock compression: experiments and molecular simulations,” J. Phys. Chem. C 121(24), 13324–13334 (2017). [CrossRef]
40. H. Sugiura, K. Kondo, and A. Sawaoka, “Dynamic response of fused quartz in the permanent densification region,” J. Appl. Phys. 52(5), 3375–3382 (1981). [CrossRef]
41. K. Saito and A. J. Ikushima, “Absorption edge in silica glass,” Phys. Rev. B 62(13), 8584–8587 (2000). [CrossRef]
42. B. Sadigh, P. Erhart, D. Åberg, A. Trave, E. Schwegler, and J. Bude, “First-principles calculations of the urbach tail in the optical absorption spectra of silica glass,” Phys. Rev. Lett. 106(2), 027401 (2011). [CrossRef]
43. C. W. Carr, J. D. Bude, and P. DeMange, “Laser-supported solid-state absorption fronts in silica,” Phys. Rev. B 82(18), 184304 (2010). [CrossRef]
44. S. P. Marsh, LASL Shock Hugoniot Data, vol. 5 (University of California, 1980).
45. S. Munetoh, T. Motooka, K. Moriguchi, and A. Shintani, “Interatomic potential for Si–O systems using Tersoff parameterization,” Comput. Mater. Sci. 39(2), 334–339 (2007). [CrossRef]
