# Densification and Devitrification of Fused Silica Induced by Ballistic Impact: A Computational Investigation

Authors: Mica Grujicic; Jennifer Snipes; S. Ramaswami; Ramin Yavari; Bryan Cheeseman
Journal: Journal of Nanomaterials
 (2015)

---

## Abstract

A molecular-level computational investigation is carried out to determine the dynamic response and material topology changes of fused silica subjected to ballistic impact by a hard projectile. The analysis was focused on the investigation of specific aspects of the dynamic response and of the topological changes such as the deformation of highly sheared and densified regions, and the conversion of amorphous fused silica to SiO2 crystalline polymorphs (in particular, α-quartz and stishovite). The topological changes in question were determined by carrying out a postprocessing atom-coordination procedure. This procedure suggested the formation of stishovite (and perhaps α-quartz) within fused silica during ballistic impact. To rationalize the findings obtained, the all-atom molecular-level computational analysis is complemented by a series of quantum-mechanics density functional theory (DFT) computations. The latter computations enable determination of the relative potential energies of the fused silica, α-quartz and stishovite, under ambient pressure (i.e., under their natural densities) as well as under imposed (as high as 50 GPa) pressures (i.e., under higher densities) and shear strains. In addition, the transition states associated with various fused-silica devitrification processes were identified. The results obtained are found to be in good agreement with their respective experimental counterparts.


---

## Body

## 1. Introduction

The present work deals with the problem of formation of crystalline phases (specificallyα-quartz and stishovite) within a fused-silica (ceramic glass containing a high purity SiO2) target during a ballistic impact. Hence, the main aspects of the present work include (a) short- and intermediate-range order topology of fused silica, (b) (crystalline) polymorphs of SiO2, and (c) devitrification (i.e., crystallization of glass taking place under the ballistic-impact loading conditions). A brief overview of these aspects of the problem at hand is presented in the remainder of this section.Short- and Intermediate-Range Order Topology of Fused Silica. Ceramic glasses (such as fused silica), like metallic glasses and glassy polymers, are amorphous materials. The molecular-level topology of these materials involves entities such as a random network of covalently bonded atoms, atomic free volumes, (network former) bridging and nonbridging oxygen atoms, and cations of glass-modifier species. However, despite the absence of a crystalline structure, the topology of ceramic glasses is not completely random. Rather, it involves different extents of short- and intermediate-range order spanning over a range of length-scales (from the quantum-mechanical to the continuum-level). To describe the structure of ceramic glasses as determined using various experimental techniques (e.g., neutron-diffraction, nuclear magnetic resonance, and small angle X-ray scattering (SAXS)), the so-called “random network model” [1] is typically employed. Such a model represents an amorphous material as a three-dimensional linked network of polyhedra. The character (number of facets) of the polyhedra is controlled by the species-specific coordination of the central (glass-forming) atom (cation). In the case of silicate-based glasses like fused silica and soda-lime glass, the polyhedron-center atoms are all silicon and each silicon atom is surrounded by four oxygen atoms (while each oxygen atom is connected to or bridges two silicon atoms) forming a SiO44- tetrahedron. In glasses of different formulations, other types of polyhedra may exist. Since silicon has a tendency to form a continuous network with (bridging) oxygen atoms, SiO2 is commonly referred to as a “network former.”Polymorphs of SiO2. Previous investigations (e.g., [2–14]) established that, under high-rate (shockwave or ballistic) loading conditions, fused-silica targets can, at least in the vicinity of the impact region, experience transformation of the amorphous structure into a crystalline one. Examination of the standard temperature-pressure phase diagram for SiO2 reveals that, at room temperature and slightly above it, SiO2 can appear in one of the following three polymorphs: α-quartz, coesite, and stishovite. Consequently, during a ballistic impact into a fused-silica target, these three SiO2 polymorphs are most likely to form (provided such impact produces a local crystallization of the fused silica). Since the main crystallographic features of these three polymorphs can be readily found in the public-domain literature, they will not be overviewed here but will be used in the computational analyses presented later in this paper. The atomic arrangements within the nonprimitive unit cells of (a) α-quartz and (b) stishovite are given in Figures 1(a) and 1(b), respectively. The same is not given for coesite since this SiO2 polymorph does not generally form under ballistic-impact conditions.Figure 1
The atomic arrangements within the nonprimitive unit cells of (a)α-quartz and (b) stishovite.
(a)
(b)Dynamic-Loading-Induced Crystallization of Glass. A detailed overview of the public-domain literature carried out as part of the present work revealed a number of experimental and computational investigations dealing with the mechanical response of fused silica to dynamic loading. Some of these studies revealed the formation of shear bands within otherwise amorphous glass and others established the formation of crystalline phases (α-quartz and stishovite, but not coesite), while still others demonstrated increased hardness of the material surrounding the impact region without establishing the topological cause for this property change.Chakraborty et al. [15] carried out a series of nanoindentation tests and showed that as the loading rate increases, the extent of shear band formation in the region surrounding the indentation decreases, while the hardness value increases. No crystal-structure analysis was carried out to determine potential formation of any of the crystalline phases as a result of loading or to provide a rationale for the observed effect of the loading rate.Tschauner et al. [16] investigated formation of stishovite in soda-lime glass during 57 GPa shock loading experiments and the reversion of this phase during subsequent release/unloading. They demonstrated that, upon loading, high-density non-fully-crystallized SiO2 phase was present in the “shocked” fused silica. Upon static loading to only 13 GPa, this phase was fully converted into the crystalline stishovite, suggesting that the shock loading was able to devitrify fused silica and form crystalline stishovite.Salleo et al. [17] demonstrated the formation of a defective form of stishovite in the surface region of fused-silica wafers through irradiation with a high-powered laser. The formation of such a phase and its continuous growth has been found to be the main cause of failure in optics used for high-power photonics.Mantisi et al. [18] carried out a comprehensive atomic-scale simulation of fused silica under combined pressure/shear loading conditions. The results obtained established permanent/irreversible densification of the fused-silica test sample and the change of the silicon and oxygen coordination relative to that present in the as-received fused silica. However, these topological changes did not appear to be the result of glass devitrification/crystallization processes; that is, the glass, while compacted, remained amorphous. One of the potential reasons for the differences in the findings reported in the present paper and in the work of Mantisi et al. regarding the response of glass to high pressure and shear is the use of different interatomic force-field potentials. It is well-established that relatively minor modifications in the force-field potentials may have profound effects on the outcome of the all-atom molecular-level simulations.Kubota et al. [19] used molecular dynamics simulations to infer the atomic-scale structural changes in fused silica induced by shock-compression loading. The results obtained revealed that shock-compressive loading involving stress levels exceeding the Hugoniot Elastic Limit gives rise to dramatic changes in the structure and topology of the fused-silica network and densifications in excess of 20%. Coordination analysis of the as-shocked fused silica revealed the formation of under- and overcoordinated Si atoms. While undercoordinated Si-atom regions could be interpreted as shock-induced fused-silica flaws, overcoordinated Si-atom regions showed some resemblance to the crystalline stishovite. In addition to the coordination changes just described, changes in glass topology (such as increases in the number of the threefold, fourfold, sevenfold, and larger rings) were observed.Main Objectives. The main objective of the present work is to carry out a series of all-atom molecular-level computational investigations of the ballistic impact by a hard projectile onto a fused-silica target-plate and to characterize the target-plate material response to such impact. In addition, in order to help rationalize some of the findings regarding the topology of the fused silica following the impact, a series of quantum-mechanical density functional theory (DFT) analyses will be carried out. These analyses will help reveal the relative potential energies of the SiO2 amorphous state and two SiO2 crystalline polymorphs, that is, α-quartz and stishovite, and the changes in these energies as a function of the extent of material compression and shear.Paper Organization. Details regarding the all-atom molecular-level computational procedure used to simulate the ballistic impact by a hard projectile onto a fused-silica target-plate are presented in Section 2. The quantum-mechanical DFT procedure used to determine the relative potential energies of the SiO2 amorphous state and the two SiO2 crystalline polymorphs, as well as the associated transition states (to be defined later), is provided in Section 3. Key results obtained in the present work are presented and discussed in Section 4, while a summary of the main findings and conclusions is provided in Section 5.


## 2. Molecular-Level Analysis of Ballistic Impact

As mentioned earlier, one of the objectives of the present work is to carry out a series of all-atom molecular-level computational analyses of the problem of a ballistic impact by a hard projectile onto a fused-silica target-plate. Within the all-atom molecular-level computational methods and tools, every atom and every bond is explicitly accounted for and molecular mechanics and dynamics algorithms are used to quantify the state and behavior of the material under investigation. All-atom molecular-level simulation problems typically require the specification of the following: (a) a molecular-level computational model consisting of atoms, ions, functional groups, and/or molecules; (b) a set of force-field functions (mathematical expressions describing bonding and nonbonding interactions between the model constituents, e.g., atoms and ions); and (c) computational method(s) to be used in the simulation. More details of these three aspects of the molecular-level modeling and simulation of fused silica are provided below.


### 2.1. Computational Model

The computational model used in this portion of the work consists of two distinct subdomains: (a) the projectile subdomain and (b) the target-plate subdomain. The two subdomains are shown and labeled in Figure2(a).Figure 2
(a) The computational model used in this portion of the work, consisting of two distinct subdomains for the projectile and the target-plate, and (b) close-up of the resulting fused-silica molecular-level random network topology. Larger ball sizes are used in this figure to highlight a pair ofSiO44- tetrahedra sharing a common oxygen atom.
(a)
(b)The projectile subdomain is of a right-circular solid cylindrical geometry (height over diameter ratio = 1.0, axis aligned with thez-direction). While, at least, the core of ballistic projectiles is typically made of hard and heavy metallic materials (e.g., tungsten), metallic materials could not be used in the construction of the projectile in the present work. The reason for this is the absence of metallic force-field functions (in the pure metallic environment) within the force-field function database used in the present work. Consequently, the projectile was made of another hard material, diamond. Typically, the projectile subdomain contained 150 covalently bonded carbon atoms forming a perfect single-crystalline diamond structure.As far as the target-plate subdomain is concerned, it is of a rectangular parallelepiped (plate-like) shape and is made of fused silica. At the molecular level, fused silica is modeled as a discrete-particle based material consisting of silicon (Si) and oxygen (O) atoms mutually bonded via a single covalent bond and forming a connected, nonstructured/amorphous network of silica (SiO44-) tetrahedra. While fused silica is an amorphous material and does not possess any long-range regularity in its atomic/molecular structure, modeling of bulk behavior of fused silica is typically done at the molecular level by assuming the existence of a larger (amorphous) unit cell. Repetition of this cell in the three orthogonal directions (the process also known as application of the “periodic boundary conditions”) results in the formation of an infinitely large bulk-type material. This procedure was adopted in the present work. The parallelepiped-shaped target-plate computational subdomain used in the present (ballistic-impact) analysis contained 9600 particles (3200 Si atoms and 6400 O atoms). The edge-lengths of the computational cell were initially set to a=b~7.9 nm, c~2.9 nm (approximately), yielding a fused-silica initial nominal density of 2.19 g/cm3. The three edges (a, b, and c) of the cell were aligned, respectively, with the three coordinate axes (x, y, and z) (with the target-plate thickness aligned with the z-direction).To create the ambient temperature/pressure “equilibrium” atomic configuration within the computational cell, the following procedure was implemented within the Visualizer [20] program from Accelrys.(a) A starting computational cell was first constructed by stacking the appropriate number ofα-quartz unit cells in the three orthogonal directions. This was followed by the affine distortion of the computational cell in order to obtain the correct mass density of the fused-silica amorphous state (obtained using the following procedure).(b) To convert the crystalline material into an amorphous one, a stochastic bond-switching algorithm was then implemented [21] using a Monte Carlo computational procedure. Within this algorithm, two neighboring Si–O pairs were randomly selected from the computational cell and the potential energy change ΔE resulting from the Si–O bond switching computed. In the ΔE<0 case, the bond switching in question was accepted without any additional conditions. On the contrary, in the ΔE>0 case, a Boltzmann probability factor PB=exp⁡-ΔE/3NkT/2 (where N is the number of atoms within the computational cell, k is the Boltzmann constant, and T is the absolute temperature) was first calculated and compared with a random number RN drawn from a (0,1) uniform distribution function. The bond switching in question was then adopted only if PB>RN.(c) The resulting structure was then subjected to a carefully devised set ofNVT (where N (=9600) is the (fixed) number of atoms within the computational cell, V is the computational cell volume (also fixed), and T is a fixed temperature) molecular dynamics simulations. Specifically, the NVT simulations were started at a temperature of 5300 K and carried out in such a way that the temperature was controlled using the following “simulated-annealing” scheme: (i) a particle-velocity scaling algorithm was applied every time step for the first 6000 steps. This enforced strict control of the temperature but produced particle velocities which were inconsistent with the target-temperature Maxwell-Boltzmann distribution function. (ii) Within the next 6000 NVT simulation steps, the frequency of particle-velocity scaling was decreased to every 40 time steps while a Nosé-Hoover [22] temperature-control algorithm (“thermostat”) was applied between the particle-velocity scaling steps. A brief description of the Nosé-Hoover thermostat could be found in our prior work [23]. (iii) During the final 8000 steps, the temperature was controlled using only the Nosé-Hoover thermostat. This procedure ensured an efficient temperature control while yielding an equilibrium state of the material (i.e., a particle-velocity distribution consistent with the target-temperature Maxwell-Boltzmann distribution function).Upon establishing the thermodynamic equilibrium at 5300 K, the target temperature was reduced by 500 K and then this procedure was reapplied at progressively (by 500 K) lower temperatures until the final temperature of 300 K was reached. The total system equilibration procedure typically involved simulation times on the order of 500 ps resulting in an average cooling-rate of ~10 K/ps. A close-up of the resulting fused-silica molecular-level random network is displayed in Figure2(b). Larger ball sizes are used in this figure to highlight a pair of SiO44- tetrahedra sharing a common oxygen atom.


### 2.2. Force Fields

The behavior of a material system at the molecular level is governed by the appropriate force fields which describe, in an approximate manner, the various interactions taking place between the constituent particles, atoms, ions, charge groups, and so forth. In other words, the knowledge of force fields enables determination of the potential energy of a system in a given configuration. In addition, gradients of the force-field functions quantify the net forces experienced by the particles, the information that is needed in the molecular dynamics simulations.In general, the potential energy of a system of interacting particles can be expressed as a sum of the valence (or bond),Evalence, cross-term, Ecross-term, and nonbond, Enon-bond, interaction energies as(1)Etotal=Evalence+Ecross-term+Enon-bond.The valence energy generally accounts for the contribution of valence electrons bonding and contains the following components: (a) a bond length/stretching term; (b) a two-bond angle term; (c) a three-bond dihedral/torsion angle term; (d) an inversion (or a four-atom out-of-plane interaction) term; and (e) the so-called “three-atom Urey-Bradley term” (i.e., the interaction of two atoms which are bonded to a common atom).The cross-term interacting energy accounts for the cross-interactions between the aforementioned valence-energy components and includes terms like (a) stretch-stretch interactions between two adjacent bonds; (b) stretch-bend interactions between a two-bond angle and one of its bonds; (c) bend-bend interactions between two valence angles associated with a common vertex atom; (d) stretch-torsion interactions between a dihedral angle and one of its end bonds; (e) stretch-torsion interactions between a dihedral angle and its middle bond; (f) bend-torsion interactions between a dihedral angle and one of its valence angles; and (g) bend-bend-torsion interactions between a dihedral angle and its two valence angles.The nonbond interaction term accounts for the interactions between nonbonded atoms and includes the van der Waals energy and the Coulomb electrostatic energy.In the present work, the so-called “COMPASS” (Condensed-Phase Optimized Molecular Potentials for Atomistic Simulation Studies) force field is used [24, 25]. This highly accurate force field is of anab initio type since most of its parameters were determined by matching the predictions made by theab initio quantum-mechanics calculations to the condensed-matter experimental data. A summary of the COMPASS force-field functions can be found in our previous work [26].


### 2.3. Computational Method(s)

All the all-atom molecular-level calculations were carried out within the present work using Discover [27] (an atomic simulation program from Accelrys). Both equilibrium and nonequilibrium molecular-dynamics (MD) analyses were employed in the present work. Within the molecular dynamics method, negative gradient of the potential energy evaluated at the location of each atom is first used to compute forces acting on each atom. Then, the associated Newton’s equations of motion (three equations for each atom) are integrated numerically over a femtosecond-long time interval. This procedure is repeated over a pico-to-nanosecond-long simulation time in order to determine the temporal evolution of the material molecular-level configuration. Equilibrium MD was used in the aforementioned Monte Carlo based bond-switching procedure to generate amorphous state from the crystalline state in SiO2. Nonequilibrium MD was used during the simulation of the ballistic impact by a hard projectile onto a fused-silica target-plate.Within the equilibrium MD method, the system under consideration is coupled to an (external) environment (e.g., constant-pressure piston and constant-temperature reservoir) which ensures that the system remains in equilibrium (i.e., the system is not subjected to any thermodynamic fluxes). In the present work,NVT, NPT (P is pressure), and NVE (E is the total energy) equilibrium MD simulations were used. Equilibrium MD calculations enable determination of the (equilibrium) thermodynamic properties of a material system through the use of time averages of the state variables sampled along the calculated system trajectories.Within nonequilibrium MD, the system is subjected to large mechanical and/or thermal perturbations (momentum transfer from the moving projectile to the initially stationary target-plate, in the present case). As a consequence, the system experiences large fluxes of its thermodynamic quantities (mass, momentum, and energy, in the present case). Since Discover was initially designed to carry out equilibrium MD simulations, a procedure had to be devised to deactivate “equilibration” portions of this algorithm so that nonequilibrium MD calculations can be carried out. This procedure was implemented using a Discover input file [25]. This file is written using the Basic Tool Command Language (BTCL) which enabled the use of a scripting engine that provides very precise control of simulation tasks (e.g., specification of the projectile constant incident velocity and deactivation of the thermal-equilibration algorithm).


### 2.4. Problem Formulation

The problem addressed in the present work involves all-atom molecular-level computational modeling of the ballistic impact by a hard projectile onto a fused-silica target-plate, and the coordination/topology analysis of the as-impacted fused silica. The projectile is of a solid right-circular cylindrical shape and impacts the target-plate normally (i.e., at a 0° obliquity angle through-the-thickness direction). The projectile is driven at a constant velocity of 1500 m/s and the target-plate is confined only along its bottom rim.


## 2.1. Computational Model

The computational model used in this portion of the work consists of two distinct subdomains: (a) the projectile subdomain and (b) the target-plate subdomain. The two subdomains are shown and labeled in Figure2(a).Figure 2
(a) The computational model used in this portion of the work, consisting of two distinct subdomains for the projectile and the target-plate, and (b) close-up of the resulting fused-silica molecular-level random network topology. Larger ball sizes are used in this figure to highlight a pair ofSiO44- tetrahedra sharing a common oxygen atom.
(a)
(b)The projectile subdomain is of a right-circular solid cylindrical geometry (height over diameter ratio = 1.0, axis aligned with thez-direction). While, at least, the core of ballistic projectiles is typically made of hard and heavy metallic materials (e.g., tungsten), metallic materials could not be used in the construction of the projectile in the present work. The reason for this is the absence of metallic force-field functions (in the pure metallic environment) within the force-field function database used in the present work. Consequently, the projectile was made of another hard material, diamond. Typically, the projectile subdomain contained 150 covalently bonded carbon atoms forming a perfect single-crystalline diamond structure.As far as the target-plate subdomain is concerned, it is of a rectangular parallelepiped (plate-like) shape and is made of fused silica. At the molecular level, fused silica is modeled as a discrete-particle based material consisting of silicon (Si) and oxygen (O) atoms mutually bonded via a single covalent bond and forming a connected, nonstructured/amorphous network of silica (SiO44-) tetrahedra. While fused silica is an amorphous material and does not possess any long-range regularity in its atomic/molecular structure, modeling of bulk behavior of fused silica is typically done at the molecular level by assuming the existence of a larger (amorphous) unit cell. Repetition of this cell in the three orthogonal directions (the process also known as application of the “periodic boundary conditions”) results in the formation of an infinitely large bulk-type material. This procedure was adopted in the present work. The parallelepiped-shaped target-plate computational subdomain used in the present (ballistic-impact) analysis contained 9600 particles (3200 Si atoms and 6400 O atoms). The edge-lengths of the computational cell were initially set to a=b~7.9 nm, c~2.9 nm (approximately), yielding a fused-silica initial nominal density of 2.19 g/cm3. The three edges (a, b, and c) of the cell were aligned, respectively, with the three coordinate axes (x, y, and z) (with the target-plate thickness aligned with the z-direction).To create the ambient temperature/pressure “equilibrium” atomic configuration within the computational cell, the following procedure was implemented within the Visualizer [20] program from Accelrys.(a) A starting computational cell was first constructed by stacking the appropriate number ofα-quartz unit cells in the three orthogonal directions. This was followed by the affine distortion of the computational cell in order to obtain the correct mass density of the fused-silica amorphous state (obtained using the following procedure).(b) To convert the crystalline material into an amorphous one, a stochastic bond-switching algorithm was then implemented [21] using a Monte Carlo computational procedure. Within this algorithm, two neighboring Si–O pairs were randomly selected from the computational cell and the potential energy change ΔE resulting from the Si–O bond switching computed. In the ΔE<0 case, the bond switching in question was accepted without any additional conditions. On the contrary, in the ΔE>0 case, a Boltzmann probability factor PB=exp⁡-ΔE/3NkT/2 (where N is the number of atoms within the computational cell, k is the Boltzmann constant, and T is the absolute temperature) was first calculated and compared with a random number RN drawn from a (0,1) uniform distribution function. The bond switching in question was then adopted only if PB>RN.(c) The resulting structure was then subjected to a carefully devised set ofNVT (where N (=9600) is the (fixed) number of atoms within the computational cell, V is the computational cell volume (also fixed), and T is a fixed temperature) molecular dynamics simulations. Specifically, the NVT simulations were started at a temperature of 5300 K and carried out in such a way that the temperature was controlled using the following “simulated-annealing” scheme: (i) a particle-velocity scaling algorithm was applied every time step for the first 6000 steps. This enforced strict control of the temperature but produced particle velocities which were inconsistent with the target-temperature Maxwell-Boltzmann distribution function. (ii) Within the next 6000 NVT simulation steps, the frequency of particle-velocity scaling was decreased to every 40 time steps while a Nosé-Hoover [22] temperature-control algorithm (“thermostat”) was applied between the particle-velocity scaling steps. A brief description of the Nosé-Hoover thermostat could be found in our prior work [23]. (iii) During the final 8000 steps, the temperature was controlled using only the Nosé-Hoover thermostat. This procedure ensured an efficient temperature control while yielding an equilibrium state of the material (i.e., a particle-velocity distribution consistent with the target-temperature Maxwell-Boltzmann distribution function).Upon establishing the thermodynamic equilibrium at 5300 K, the target temperature was reduced by 500 K and then this procedure was reapplied at progressively (by 500 K) lower temperatures until the final temperature of 300 K was reached. The total system equilibration procedure typically involved simulation times on the order of 500 ps resulting in an average cooling-rate of ~10 K/ps. A close-up of the resulting fused-silica molecular-level random network is displayed in Figure2(b). Larger ball sizes are used in this figure to highlight a pair of SiO44- tetrahedra sharing a common oxygen atom.


## 2.2. Force Fields

The behavior of a material system at the molecular level is governed by the appropriate force fields which describe, in an approximate manner, the various interactions taking place between the constituent particles, atoms, ions, charge groups, and so forth. In other words, the knowledge of force fields enables determination of the potential energy of a system in a given configuration. In addition, gradients of the force-field functions quantify the net forces experienced by the particles, the information that is needed in the molecular dynamics simulations.In general, the potential energy of a system of interacting particles can be expressed as a sum of the valence (or bond),Evalence, cross-term, Ecross-term, and nonbond, Enon-bond, interaction energies as(1)Etotal=Evalence+Ecross-term+Enon-bond.The valence energy generally accounts for the contribution of valence electrons bonding and contains the following components: (a) a bond length/stretching term; (b) a two-bond angle term; (c) a three-bond dihedral/torsion angle term; (d) an inversion (or a four-atom out-of-plane interaction) term; and (e) the so-called “three-atom Urey-Bradley term” (i.e., the interaction of two atoms which are bonded to a common atom).The cross-term interacting energy accounts for the cross-interactions between the aforementioned valence-energy components and includes terms like (a) stretch-stretch interactions between two adjacent bonds; (b) stretch-bend interactions between a two-bond angle and one of its bonds; (c) bend-bend interactions between two valence angles associated with a common vertex atom; (d) stretch-torsion interactions between a dihedral angle and one of its end bonds; (e) stretch-torsion interactions between a dihedral angle and its middle bond; (f) bend-torsion interactions between a dihedral angle and one of its valence angles; and (g) bend-bend-torsion interactions between a dihedral angle and its two valence angles.The nonbond interaction term accounts for the interactions between nonbonded atoms and includes the van der Waals energy and the Coulomb electrostatic energy.In the present work, the so-called “COMPASS” (Condensed-Phase Optimized Molecular Potentials for Atomistic Simulation Studies) force field is used [24, 25]. This highly accurate force field is of anab initio type since most of its parameters were determined by matching the predictions made by theab initio quantum-mechanics calculations to the condensed-matter experimental data. A summary of the COMPASS force-field functions can be found in our previous work [26].


## 2.3. Computational Method(s)

All the all-atom molecular-level calculations were carried out within the present work using Discover [27] (an atomic simulation program from Accelrys). Both equilibrium and nonequilibrium molecular-dynamics (MD) analyses were employed in the present work. Within the molecular dynamics method, negative gradient of the potential energy evaluated at the location of each atom is first used to compute forces acting on each atom. Then, the associated Newton’s equations of motion (three equations for each atom) are integrated numerically over a femtosecond-long time interval. This procedure is repeated over a pico-to-nanosecond-long simulation time in order to determine the temporal evolution of the material molecular-level configuration. Equilibrium MD was used in the aforementioned Monte Carlo based bond-switching procedure to generate amorphous state from the crystalline state in SiO2. Nonequilibrium MD was used during the simulation of the ballistic impact by a hard projectile onto a fused-silica target-plate.Within the equilibrium MD method, the system under consideration is coupled to an (external) environment (e.g., constant-pressure piston and constant-temperature reservoir) which ensures that the system remains in equilibrium (i.e., the system is not subjected to any thermodynamic fluxes). In the present work,NVT, NPT (P is pressure), and NVE (E is the total energy) equilibrium MD simulations were used. Equilibrium MD calculations enable determination of the (equilibrium) thermodynamic properties of a material system through the use of time averages of the state variables sampled along the calculated system trajectories.Within nonequilibrium MD, the system is subjected to large mechanical and/or thermal perturbations (momentum transfer from the moving projectile to the initially stationary target-plate, in the present case). As a consequence, the system experiences large fluxes of its thermodynamic quantities (mass, momentum, and energy, in the present case). Since Discover was initially designed to carry out equilibrium MD simulations, a procedure had to be devised to deactivate “equilibration” portions of this algorithm so that nonequilibrium MD calculations can be carried out. This procedure was implemented using a Discover input file [25]. This file is written using the Basic Tool Command Language (BTCL) which enabled the use of a scripting engine that provides very precise control of simulation tasks (e.g., specification of the projectile constant incident velocity and deactivation of the thermal-equilibration algorithm).


## 2.4. Problem Formulation

The problem addressed in the present work involves all-atom molecular-level computational modeling of the ballistic impact by a hard projectile onto a fused-silica target-plate, and the coordination/topology analysis of the as-impacted fused silica. The projectile is of a solid right-circular cylindrical shape and impacts the target-plate normally (i.e., at a 0° obliquity angle through-the-thickness direction). The projectile is driven at a constant velocity of 1500 m/s and the target-plate is confined only along its bottom rim.


## 3. Quantum-Mechanical Analysis of Glass Devitrification

In order to rationalize the molecular-level computational results pertaining to the ballistic impact of a projectile onto the fused-silica target-plate, and the potentially accompanying changes in the glass topology, a quantum-mechanical DFT analysis of the fused-silica devitrification process is carried out in the present work. The main purpose of this analysis was (a) the establishment of the energy difference between the as-received glass material state and the crystallineα-quartz and stishovite SiO2 states; (b) determination of the transition state energy barriers associated with the two (glass → α-quartz or glass → stishovite) fused-silica devitrification processes. The transition state is the material state along the devitrification pathway which is associated with the maximum potential energy; and (c) the establishment of the effect of high pressure and shear stresses on the energetics of the initial, transition, and final states.


### 3.1. Computational Models

As mentioned earlier, coesite does not typically form under ballistic-loading conditions and, hence, will not be the subject of the present DFT investigation. Consequently, only two glass devitrification product states,α-quartz and stishovite, will be investigated. The respective computational cells are shown in Figures 1(a) and 1(b). Since fused silica possesses an amorphous structure, without long-range order, a substantially larger computational cell had to be used for this (initial) state of SiO2, in order to more accurately determine its average potential energy. An example of such a computational cell is given in Figure 2(a), the target-plate computational subdomain.In the portion of the analysis dealing with determination of the transition state, the initial and final states of the material have to contain the same number of atoms and species. For that reason, fused-silica computational cells containing 48 Si and 96 O atoms (the numbers correspond to a2×2×2=8 stack of the unit cells depicted in Figure 1(a)) are used to investigate fused-silica → α-quartz transition state. Likewise, fused-silica computational cells containing 16 Si and 32 O atoms (the numbers correspond to a 2×2×2=8 stack of the unit cells depicted in Figure 1(b)) are used to investigate fused-silica → stishovite transition state. In order to account for the statistical effects associated with the extraction of such small fused-silica computational cells from a larger computational cell, ten small fused-silica computational cells (not shown for brevity) are generated for each of the two transition state analyses and the results obtained are averaged out. It should be noted, however, that due to the fact that the initial state of the material within the computational cells is amorphous, the cell sizes used are relatively small, negatively affecting the accuracy of the results obtained.


### 3.2. Computational Method

All calculations in the present work are carried out using theabinitio density-functional theory code DMol3 developed by Accelrys Inc. [28]. In this code each electronic wave function is expanded in a localized atom-centered basis set with each basis function defined numerically on a dense radial grid. No pseudopotential approximation is used for the near-core electrons. Instead, all-electron calculations are performed with a double numerical polarized (DNP) basis set, the most complete basis set available in the DMol3 code. This basis set is equivalent to the commonly used analytical 6-31g∗∗ basis set, a split-valence basis set with polarization functions p to H and d to C, and the halogens [29].Within the exchange-correlation potential of the density functional theory (DFT), the variation of the electron-charge density at the atomic scale is represented using the so-called Generalized Gradient Approximation (GGA). Following Grujicic et al. [30, 31], the Perdew-Burke-Ernzerhof (PBE) gradient-corrected functional [32] is used in the present work. A standard value of 5.5 Å is assigned to the finite basis-set cutoff radius.


### 3.3. Determination of the Transition States

In general, the transition state lies on a Minimum Energy Pathway connecting the initial and final states’ potential energy minima, where the pathway and the two minima all reside on the associated potential energy hypersurface. The defining feature of the transition state is that it is associated with a minimum potential energy in all the directions but one (in which the potential energy experiences a maximum). In the case of a system with two degrees of freedom, the transition state corresponds to a saddle point, as depicted in Figure3. Moving the system under consideration from the transition state point (in either direction) along the steepest-descent path leads to the initial/final states.Figure 3
The initial-to-final state transition in a two-degree-of-freedom system and identification of the corresponding transition state (the saddle point).There are several algorithms which are commonly used for determination of the transition state. The three most frequently used include (a) linear synchronous transit (LST) [34], used in the present work; (b) quadratic synchronous transit (QST) [34]; and (c) nudged elastic band [33].The LST method constructs the initial-state → final-state transition pathway by connecting each atom in its initial and final states using a straight-line pathway and by constructing the intermediate configurations by linearly and synchronously interpolating the atomic positions along the pathway of each atom. In other words, the structure of the intermediate states is a rule of mixtures of the initial structure and the final structure, with a weighting factor for the product statef=0,1. The potential energy is next determined for all the intermediate states, and the one associated with the largest potential energy is taken as a first approximation of the system transition state. The location of the transition state is further improved by carrying out a constrained optimization of its first approximation.


## 3.1. Computational Models

As mentioned earlier, coesite does not typically form under ballistic-loading conditions and, hence, will not be the subject of the present DFT investigation. Consequently, only two glass devitrification product states,α-quartz and stishovite, will be investigated. The respective computational cells are shown in Figures 1(a) and 1(b). Since fused silica possesses an amorphous structure, without long-range order, a substantially larger computational cell had to be used for this (initial) state of SiO2, in order to more accurately determine its average potential energy. An example of such a computational cell is given in Figure 2(a), the target-plate computational subdomain.In the portion of the analysis dealing with determination of the transition state, the initial and final states of the material have to contain the same number of atoms and species. For that reason, fused-silica computational cells containing 48 Si and 96 O atoms (the numbers correspond to a2×2×2=8 stack of the unit cells depicted in Figure 1(a)) are used to investigate fused-silica → α-quartz transition state. Likewise, fused-silica computational cells containing 16 Si and 32 O atoms (the numbers correspond to a 2×2×2=8 stack of the unit cells depicted in Figure 1(b)) are used to investigate fused-silica → stishovite transition state. In order to account for the statistical effects associated with the extraction of such small fused-silica computational cells from a larger computational cell, ten small fused-silica computational cells (not shown for brevity) are generated for each of the two transition state analyses and the results obtained are averaged out. It should be noted, however, that due to the fact that the initial state of the material within the computational cells is amorphous, the cell sizes used are relatively small, negatively affecting the accuracy of the results obtained.


## 3.2. Computational Method

All calculations in the present work are carried out using theabinitio density-functional theory code DMol3 developed by Accelrys Inc. [28]. In this code each electronic wave function is expanded in a localized atom-centered basis set with each basis function defined numerically on a dense radial grid. No pseudopotential approximation is used for the near-core electrons. Instead, all-electron calculations are performed with a double numerical polarized (DNP) basis set, the most complete basis set available in the DMol3 code. This basis set is equivalent to the commonly used analytical 6-31g∗∗ basis set, a split-valence basis set with polarization functions p to H and d to C, and the halogens [29].Within the exchange-correlation potential of the density functional theory (DFT), the variation of the electron-charge density at the atomic scale is represented using the so-called Generalized Gradient Approximation (GGA). Following Grujicic et al. [30, 31], the Perdew-Burke-Ernzerhof (PBE) gradient-corrected functional [32] is used in the present work. A standard value of 5.5 Å is assigned to the finite basis-set cutoff radius.


## 3.3. Determination of the Transition States

In general, the transition state lies on a Minimum Energy Pathway connecting the initial and final states’ potential energy minima, where the pathway and the two minima all reside on the associated potential energy hypersurface. The defining feature of the transition state is that it is associated with a minimum potential energy in all the directions but one (in which the potential energy experiences a maximum). In the case of a system with two degrees of freedom, the transition state corresponds to a saddle point, as depicted in Figure3. Moving the system under consideration from the transition state point (in either direction) along the steepest-descent path leads to the initial/final states.Figure 3
The initial-to-final state transition in a two-degree-of-freedom system and identification of the corresponding transition state (the saddle point).There are several algorithms which are commonly used for determination of the transition state. The three most frequently used include (a) linear synchronous transit (LST) [34], used in the present work; (b) quadratic synchronous transit (QST) [34]; and (c) nudged elastic band [33].The LST method constructs the initial-state → final-state transition pathway by connecting each atom in its initial and final states using a straight-line pathway and by constructing the intermediate configurations by linearly and synchronously interpolating the atomic positions along the pathway of each atom. In other words, the structure of the intermediate states is a rule of mixtures of the initial structure and the final structure, with a weighting factor for the product statef=0,1. The potential energy is next determined for all the intermediate states, and the one associated with the largest potential energy is taken as a first approximation of the system transition state. The location of the transition state is further improved by carrying out a constrained optimization of its first approximation.


## 4. Results and Discussion

### 4.1. Validation of the As-Received Fused-Silica Material State

In this section, an attempt is made to validate the fused-silica room-temperature/ambient-pressure structure obtained through application of the previously described bond-switching and simulated-annealing computational procedures. In particular, the material mass density, the partial radial distribution functions for the three (Si–Si, Si–O, and O–O) atomic pairs, and distributions of the Si-atom and O-atom coordination numbers are calculated and compared with their experimental counterparts.To compute the material mass density, the partial radial distribution functions, and atomic coordination combinations,NPT equilibrium molecular dynamics simulations were run at the room temperature and the ambient pressure. Simulations were carried out for over 10000 (0.1 fs-long) time steps while maintaining the temperature using the Nosé-Hoover thermostat and scaling the particle velocities every ten time steps. Pressure was controlled using a Berendsen barostat [35]. A brief overview of this barostat is provided in our prior work [36].


#### 4.1.1. Mass Density

The average mass density of fused silica computed from the computational-cell average volume and the cell mass has been found to be about 1.0% greater than its commonly cited experimental counterpart of 2.19 g/cm3. This computation/experiment agreement has been deemed to be reasonably good.


#### 4.1.2. Radial Distribution Functions

The partial radial distribution (often also referred to as the partial pair correlation) function provides a measure of the probability that, given the presence of an atom of typeα at the origin of an arbitrary reference frame, there will be an atom of type β within a spherical shell of infinitesimal thickness dr at a distance r from the reference atom. Alternatively, this function can be considered as a function which defines a ratio of the probability of finding an α-β atomic pair with the separation distance r, and the average probability of finding an α-β atomic pair at the same distance. In amorphous materials like fused silica, the partial pair correlation functions are quite important since they (a) provide an insight into the short-range order of the system, (b) can be used in the assessment of continuum-level thermodynamic material properties, and (c) provide a way of validating the molecular-level calculations since these quantities can also be determined experimentally using X-ray diffraction.The computed partial radial distribution functions for the as-received/initial state of fused silica, Figure4(a), are compared with their counterparts based on the shell model molecular dynamics calculations [37], Figure 4(b). This comparison reveals that the present computational results are qualitatively similar to the ones reported in [38–40]. As far as the quantitative agreement between the two sets of results is concerned, it could be characterized as being fair to good.Figure 4
A comparison between the fused-silica Si–Si, O–O, and Si–O partial radial distribution functions computed in (a) the as-received/initial state, determined in the present work, and (b) the equilibrium initial state as reported in Henkelman and Jónsson [33].
(a)
(b)


#### 4.1.3. Si- and O-Atom Coordinations

Since fused silica does not contain any network modifiers and Si is the only network-forming element, each Si atom is expected to be bonded to four O atoms while each O atom is expected to be bonded to two Si atoms. This is confirmed through postprocessing of the molecular dynamics results, as shown in Figures5(a)-5(b) (the data labeled “initial state”). The results shown in Figure 5(a) pertain to the Si-atom coordination while those displayed in Figure 5(b) refer to the O-atom coordination.Figure 5
Fractional distribution of the (a) Si-atom and (b) O-atom coordination numbers in the initial and the as-impacted states of fused silica.
(a)
(b)


### 4.2. Analysis of the Ballistic Impact

Temporal evolution of the computational domain at four (0.5 ps, 1.5 ps, 2.5 ps, and 3.5 ps) times following the initial contact between the diamond solid right-circular cylindrical projectile moving at a high velocity and the fused-silica target-plate is depicted in Figures6(a)–6(d), respectively. For improved clarity, the symbols for the target-plate atoms are made smaller. Examination of the results displayed in Figures 6(a)–6(d) reveals that the target-plate material in the close vicinity of the projectile experiences major changes. However, the nature of these changes, that is, the accompanying alterations in the material topology, is difficult to infer from the results displayed in Figures 6(a)–6(d). To overcome this problem, a few selected results revealing the as-impacted fused-silica local topology in the region adjacent to the projectile are presented and discussed below.Figure 6
Temporal evolution of the computational domain at four times: (a) 0.5 ps, (b) 1.5 ps, (c) 2.5 ps, and (d) 3.5 ps, following the initial contact between the diamond solid right-circular cylindrical projectile moving at a high velocity and the fused-silica target-plate.
(a)
(b)
(c)
(d)Figure7 shows a thin slice of the fused-silica target-plate, having the faces parallel to the x-z plane, centered on the hole created by the projectile after penetrating the target by about half of its thickness. Examination of the fused-silica topology within this slice revealed the presence of numerous six-folded Si and three-folded O atoms (the Si and O coordinations characteristic of stishovite and not of fused silica). In Figure 7, two six-coordinated Si and two three-coordinated O atoms are highlighted by assigning a larger sphere radius to the atoms involved. In addition, six-folded Si atoms and three-folded O atoms are tagged with circular symbols. In addition to the highlighted atomic configurations, many more instances of six-coordinated Si and three-coordinated O are found in the region surrounding the penetration hole. Since, as mentioned earlier, the initial state of fused-silica only contained four-folded Si and two-folded O atoms, this finding suggests that ballistic impact can lead, at least locally, to the conversion of amorphous SiO2 into the crystalline (but highly deformed) stishovite structure.Figure 7
Local stishovite-like topology showing two six-folded Si and two three-folded O atoms. For clarity, a larger sphere radius is assigned to the atoms involved. In addition, six-folded Si atoms and three-folded O atoms are tagged with circular symbols.Additional changes observed in the as-impacted fused silica pertain to the distribution of the smallest Si–O rings. To quantify the size-distribution of the smallest Si–O rings in both the initial and the as-impacted fused-silica states, a computational method was developed. This method solves the class of so-called “shortest path problems” and is a simple modification of Dijkstra’s algorithm [41, 42]. The main modifications in this algorithm are associated with the fact that, in the present case, the starting point and the destination point of the path are identical. The smallest-ring size-distribution results for the initial and as-impacted fused-silica states are displayed in Figures 8(a)-8(b), respectively. Examination of these results reveals that ballistic impact alters the ring structure within fused silica. Specifically, while no five-membered rings were present in the initial state of fused silica, such rings are found in the as-impacted state of the same material. Additional changes observed pertain to the topology of the six-membered rings. That is, while in the initial state of fused silica these rings resemble the corresponding rings found in cristobalite (another polymorph of SiO2), Figure 9(a), in the as-impacted fused silica the topology of the six-membered rings was found to resemble more that found in α-quartz, Figure 9(b).Figure 8
Size-distribution function for the smallest Si–O rings in the fused-silica: (a) initial state and (b) as-impacted state.
(a)
(b)Figure 9
Two views of the atomic structure of six-membered Si–O rings in (a) cristobalite and (b) quartz.
(a)
(b)Three partial radial distribution functions for the fused-silica target-plate after the diamond impactor has penetrated halfway through the target-plate thickness are depicted in Figure10(a). A comparison of these results with their as-received fused-silica counterparts, Figure 4(a), reveals that the ballistic impact causes distinct changes in the short-range order and atomic coordination within this material. To help rationalize these changes, the same partial radial distribution functions are calculated for α-quartz, Figure 10(b), and stishovite, Figure 10(c). A comparison of the results displayed in Figure 4(a) and Figures 10(a)–10(c) further confirms that the ballistic impact-induced changes in the pair correlation functions are associated with the devitrification of fused silica, and the formation of α-quartz and stishovite. In other words, differences in the partial pair correlation functions in the as-impacted fused silica relative to those in the as-received fused silica appear to be caused by the presence of Si–Si, Si–O, and O–O atomic pairs with an atomic environment similar to those found in stishovite (and α-quartz).Figure 10
Three partial radial distribution functions for (a) the fused-silica region adjacent to the projectile, after the diamond impactor has penetrated approximately halfway through the target-plate thickness, (b)α-quartz and (c) stishovite.
(a)
(b)
(c)To further demonstrate the conversion of fused silica to stishovite as a result of the ballistic impact, Si- and O-atom coordinations in the fused-silica region surrounding the penetration hole are determined. The results of this procedure, labeled “as-impacted state,” are shown in Figures5(a)-5(b). Examination of the results displayed in these figures reveals the presence of six-coordinated Si and three-coordinated O atoms, the atomic coordination which characterizes the stishovite crystal structure.It should be noted that all the results reported in this subsection were obtained under the conditions of the projectile being driven at a constant velocity of 1500 m/s during the MD simulations. Under such simulation conditions, one cannot determine the extent of momentum transfer from the projectile to the target. In a separate analysis, the projectile was assigned an initial velocity of 1500 m/s rather than being driven at this velocity. In this simulation, it was found that about 60% of the incident linear momentum of the projectile is transferred to the target. In other words, the residual velocity of the projectile after passing through the target is approximately 600 m/s.


### 4.3. Relative Stability of Fused Silica,α-Quartz and Stishovite

To further help rationalize the changes in the partial pair correlation functions induced by the ballistic impact, the quantum-mechanical DFT calculation results pertaining to the relative stability of fused silica,α-quartz and stishovite, as well as of the energy barriers associated with the fused-silica → α-quartz and fused-silica → stishovite transition states, are presented and discussed in this section. The relative room-temperature potential energies of the fused silica, α-quartz and stishovite as a function of pressure are presented in Figure 11(a). Examination of the results displayed in this figure reveals that, at the ambient pressure, α-quartz is the most stable form of SiO2, followed by fused silica and then by stishovite. At a pressure of 50 GPa, on the other hand, the relative positions of α-quartz and fused silica, with respect to thermodynamic stability, have been exchanged.Figure 11
Relative room-temperature potential energies of fused silica,α-quartz and stishovite, as a function of pressure, in the (a) absence and (b) presence of 5% shear.
(a)
(b)In Figure10(b), the relative room-temperature potential energies of the fused silica, α-quartz and stishovite as a function of pressure, in the presence of 5% shear, are presented. A comparison of the results displayed in Figures 11(a)-11(b) reveals that the application of 5% shear strain did not change the relative stability of fused silica, α-quartz and stishovite at the ambient pressure. In sharp contrast, at 50 GPa pressure, the application of 5% shear has been found to make the stishovite the most stable phase, followed by α-quartz and then fused silica. This finding suggests that, in the presence of shear, fused silica is more likely to undergo a devitrification conversion into stishovite and, to a lower extent, into α-quartz. This finding is important since the fused-silica region surrounding the projectile is generally subjected to high shear, in addition to high pressures. This explains why, in Figure 7, this region was found to undergo an extensive fused-silica → stishovite conversion.The likelihood for the aforementioned devitrification processes is affected not only by the relative stabilities of fused silica,α-quartz and stishovite, but also by the size of the energy barrier associated with the respective transition state. Variations in the fused-silica → α-quartz and fused-silica → stishovite transition state energy barriers with pressure, in the absence and the presence of 5% shear, are shown, respectively, in Figures 12(a)-12(b). Examination of the results displayed in these figures, at a pressure of 50 GPa, reveals that the energy barrier for fused-silica devitrification, in the absence of shear, is slightly higher for the fused-silica → stishovite conversion. On the other hand, the energy barrier for fused-silica devitrification, in the presence of shear, is slightly higher for the fused-silica → α-quartz conversion. These findings suggest that, in the presence of shear, and at pressures as high as 50 GPa, fused silica is more likely to convert to stishovite than α-quartz.Figure 12
Variations in the fused-silica →α-quartz and fused-silica → stishovite transition state energy barriers with pressure, in the (a) absence and (b) presence of 5% shear.
(a)
(b)Figure13 shows the conversion of an initially amorphous SiO2 structure into a stishovite-like structure under the influence of high pressure and shear. The transition state associated with this conversion is also shown in this figure. For clarity, only a small SiO2 region consisting of three Si and fifteen O atoms is displayed in Figure 13. Due to such a small size of the region, some O atoms appear to be nonbonded. These O atoms are bonded to Si atoms, but the Si atoms that they are bonded to are not shown in this figure. The Si atoms in question are also bonded to some of the bonded O atoms displayed in Figure 13. Examination of this figure reveals that, as expected, the fused-silica state contains only fourfold coordinated Si and twofold coordinated O atoms, while the stishovite state contains sixfold coordinated Si and threefold coordinated O atoms. The Si- and O-atom coordination in the transition state involves fourfold and fivefold coordinated Si and twofold and threefold coordinated O atoms. To help with the interpretation of the results displayed in Figure 13, the same three Si atoms appearing in the three configurations are denoted using labels “Si 1,” “Si 2,” and “Si 3.” It should be recalled that the results presented in Figures 12(a)-12(b) reveal that the maximum energy associated with the transition state is lowered and, thus, the conversion of the fused-silica → stishovite becomes feasible, under high pressure and shear.Figure 13
Conversion of an initially amorphous SiO2 structure into a stishovite-like structure under the influence of high pressure and shear. The associated transition state is also shown.


## 4.1. Validation of the As-Received Fused-Silica Material State

In this section, an attempt is made to validate the fused-silica room-temperature/ambient-pressure structure obtained through application of the previously described bond-switching and simulated-annealing computational procedures. In particular, the material mass density, the partial radial distribution functions for the three (Si–Si, Si–O, and O–O) atomic pairs, and distributions of the Si-atom and O-atom coordination numbers are calculated and compared with their experimental counterparts.To compute the material mass density, the partial radial distribution functions, and atomic coordination combinations,NPT equilibrium molecular dynamics simulations were run at the room temperature and the ambient pressure. Simulations were carried out for over 10000 (0.1 fs-long) time steps while maintaining the temperature using the Nosé-Hoover thermostat and scaling the particle velocities every ten time steps. Pressure was controlled using a Berendsen barostat [35]. A brief overview of this barostat is provided in our prior work [36].


### 4.1.1. Mass Density

The average mass density of fused silica computed from the computational-cell average volume and the cell mass has been found to be about 1.0% greater than its commonly cited experimental counterpart of 2.19 g/cm3. This computation/experiment agreement has been deemed to be reasonably good.


### 4.1.2. Radial Distribution Functions

The partial radial distribution (often also referred to as the partial pair correlation) function provides a measure of the probability that, given the presence of an atom of typeα at the origin of an arbitrary reference frame, there will be an atom of type β within a spherical shell of infinitesimal thickness dr at a distance r from the reference atom. Alternatively, this function can be considered as a function which defines a ratio of the probability of finding an α-β atomic pair with the separation distance r, and the average probability of finding an α-β atomic pair at the same distance. In amorphous materials like fused silica, the partial pair correlation functions are quite important since they (a) provide an insight into the short-range order of the system, (b) can be used in the assessment of continuum-level thermodynamic material properties, and (c) provide a way of validating the molecular-level calculations since these quantities can also be determined experimentally using X-ray diffraction.The computed partial radial distribution functions for the as-received/initial state of fused silica, Figure4(a), are compared with their counterparts based on the shell model molecular dynamics calculations [37], Figure 4(b). This comparison reveals that the present computational results are qualitatively similar to the ones reported in [38–40]. As far as the quantitative agreement between the two sets of results is concerned, it could be characterized as being fair to good.Figure 4
A comparison between the fused-silica Si–Si, O–O, and Si–O partial radial distribution functions computed in (a) the as-received/initial state, determined in the present work, and (b) the equilibrium initial state as reported in Henkelman and Jónsson [33].
(a)
(b)


### 4.1.3. Si- and O-Atom Coordinations

Since fused silica does not contain any network modifiers and Si is the only network-forming element, each Si atom is expected to be bonded to four O atoms while each O atom is expected to be bonded to two Si atoms. This is confirmed through postprocessing of the molecular dynamics results, as shown in Figures5(a)-5(b) (the data labeled “initial state”). The results shown in Figure 5(a) pertain to the Si-atom coordination while those displayed in Figure 5(b) refer to the O-atom coordination.Figure 5
Fractional distribution of the (a) Si-atom and (b) O-atom coordination numbers in the initial and the as-impacted states of fused silica.
(a)
(b)


## 4.1.1. Mass Density

The average mass density of fused silica computed from the computational-cell average volume and the cell mass has been found to be about 1.0% greater than its commonly cited experimental counterpart of 2.19 g/cm3. This computation/experiment agreement has been deemed to be reasonably good.


## 4.1.2. Radial Distribution Functions

The partial radial distribution (often also referred to as the partial pair correlation) function provides a measure of the probability that, given the presence of an atom of typeα at the origin of an arbitrary reference frame, there will be an atom of type β within a spherical shell of infinitesimal thickness dr at a distance r from the reference atom. Alternatively, this function can be considered as a function which defines a ratio of the probability of finding an α-β atomic pair with the separation distance r, and the average probability of finding an α-β atomic pair at the same distance. In amorphous materials like fused silica, the partial pair correlation functions are quite important since they (a) provide an insight into the short-range order of the system, (b) can be used in the assessment of continuum-level thermodynamic material properties, and (c) provide a way of validating the molecular-level calculations since these quantities can also be determined experimentally using X-ray diffraction.The computed partial radial distribution functions for the as-received/initial state of fused silica, Figure4(a), are compared with their counterparts based on the shell model molecular dynamics calculations [37], Figure 4(b). This comparison reveals that the present computational results are qualitatively similar to the ones reported in [38–40]. As far as the quantitative agreement between the two sets of results is concerned, it could be characterized as being fair to good.Figure 4
A comparison between the fused-silica Si–Si, O–O, and Si–O partial radial distribution functions computed in (a) the as-received/initial state, determined in the present work, and (b) the equilibrium initial state as reported in Henkelman and Jónsson [33].
(a)
(b)


## 4.1.3. Si- and O-Atom Coordinations

Since fused silica does not contain any network modifiers and Si is the only network-forming element, each Si atom is expected to be bonded to four O atoms while each O atom is expected to be bonded to two Si atoms. This is confirmed through postprocessing of the molecular dynamics results, as shown in Figures5(a)-5(b) (the data labeled “initial state”). The results shown in Figure 5(a) pertain to the Si-atom coordination while those displayed in Figure 5(b) refer to the O-atom coordination.Figure 5
Fractional distribution of the (a) Si-atom and (b) O-atom coordination numbers in the initial and the as-impacted states of fused silica.
(a)
(b)


## 4.2. Analysis of the Ballistic Impact

Temporal evolution of the computational domain at four (0.5 ps, 1.5 ps, 2.5 ps, and 3.5 ps) times following the initial contact between the diamond solid right-circular cylindrical projectile moving at a high velocity and the fused-silica target-plate is depicted in Figures6(a)–6(d), respectively. For improved clarity, the symbols for the target-plate atoms are made smaller. Examination of the results displayed in Figures 6(a)–6(d) reveals that the target-plate material in the close vicinity of the projectile experiences major changes. However, the nature of these changes, that is, the accompanying alterations in the material topology, is difficult to infer from the results displayed in Figures 6(a)–6(d). To overcome this problem, a few selected results revealing the as-impacted fused-silica local topology in the region adjacent to the projectile are presented and discussed below.Figure 6
Temporal evolution of the computational domain at four times: (a) 0.5 ps, (b) 1.5 ps, (c) 2.5 ps, and (d) 3.5 ps, following the initial contact between the diamond solid right-circular cylindrical projectile moving at a high velocity and the fused-silica target-plate.
(a)
(b)
(c)
(d)Figure7 shows a thin slice of the fused-silica target-plate, having the faces parallel to the x-z plane, centered on the hole created by the projectile after penetrating the target by about half of its thickness. Examination of the fused-silica topology within this slice revealed the presence of numerous six-folded Si and three-folded O atoms (the Si and O coordinations characteristic of stishovite and not of fused silica). In Figure 7, two six-coordinated Si and two three-coordinated O atoms are highlighted by assigning a larger sphere radius to the atoms involved. In addition, six-folded Si atoms and three-folded O atoms are tagged with circular symbols. In addition to the highlighted atomic configurations, many more instances of six-coordinated Si and three-coordinated O are found in the region surrounding the penetration hole. Since, as mentioned earlier, the initial state of fused-silica only contained four-folded Si and two-folded O atoms, this finding suggests that ballistic impact can lead, at least locally, to the conversion of amorphous SiO2 into the crystalline (but highly deformed) stishovite structure.Figure 7
Local stishovite-like topology showing two six-folded Si and two three-folded O atoms. For clarity, a larger sphere radius is assigned to the atoms involved. In addition, six-folded Si atoms and three-folded O atoms are tagged with circular symbols.Additional changes observed in the as-impacted fused silica pertain to the distribution of the smallest Si–O rings. To quantify the size-distribution of the smallest Si–O rings in both the initial and the as-impacted fused-silica states, a computational method was developed. This method solves the class of so-called “shortest path problems” and is a simple modification of Dijkstra’s algorithm [41, 42]. The main modifications in this algorithm are associated with the fact that, in the present case, the starting point and the destination point of the path are identical. The smallest-ring size-distribution results for the initial and as-impacted fused-silica states are displayed in Figures 8(a)-8(b), respectively. Examination of these results reveals that ballistic impact alters the ring structure within fused silica. Specifically, while no five-membered rings were present in the initial state of fused silica, such rings are found in the as-impacted state of the same material. Additional changes observed pertain to the topology of the six-membered rings. That is, while in the initial state of fused silica these rings resemble the corresponding rings found in cristobalite (another polymorph of SiO2), Figure 9(a), in the as-impacted fused silica the topology of the six-membered rings was found to resemble more that found in α-quartz, Figure 9(b).Figure 8
Size-distribution function for the smallest Si–O rings in the fused-silica: (a) initial state and (b) as-impacted state.
(a)
(b)Figure 9
Two views of the atomic structure of six-membered Si–O rings in (a) cristobalite and (b) quartz.
(a)
(b)Three partial radial distribution functions for the fused-silica target-plate after the diamond impactor has penetrated halfway through the target-plate thickness are depicted in Figure10(a). A comparison of these results with their as-received fused-silica counterparts, Figure 4(a), reveals that the ballistic impact causes distinct changes in the short-range order and atomic coordination within this material. To help rationalize these changes, the same partial radial distribution functions are calculated for α-quartz, Figure 10(b), and stishovite, Figure 10(c). A comparison of the results displayed in Figure 4(a) and Figures 10(a)–10(c) further confirms that the ballistic impact-induced changes in the pair correlation functions are associated with the devitrification of fused silica, and the formation of α-quartz and stishovite. In other words, differences in the partial pair correlation functions in the as-impacted fused silica relative to those in the as-received fused silica appear to be caused by the presence of Si–Si, Si–O, and O–O atomic pairs with an atomic environment similar to those found in stishovite (and α-quartz).Figure 10
Three partial radial distribution functions for (a) the fused-silica region adjacent to the projectile, after the diamond impactor has penetrated approximately halfway through the target-plate thickness, (b)α-quartz and (c) stishovite.
(a)
(b)
(c)To further demonstrate the conversion of fused silica to stishovite as a result of the ballistic impact, Si- and O-atom coordinations in the fused-silica region surrounding the penetration hole are determined. The results of this procedure, labeled “as-impacted state,” are shown in Figures5(a)-5(b). Examination of the results displayed in these figures reveals the presence of six-coordinated Si and three-coordinated O atoms, the atomic coordination which characterizes the stishovite crystal structure.It should be noted that all the results reported in this subsection were obtained under the conditions of the projectile being driven at a constant velocity of 1500 m/s during the MD simulations. Under such simulation conditions, one cannot determine the extent of momentum transfer from the projectile to the target. In a separate analysis, the projectile was assigned an initial velocity of 1500 m/s rather than being driven at this velocity. In this simulation, it was found that about 60% of the incident linear momentum of the projectile is transferred to the target. In other words, the residual velocity of the projectile after passing through the target is approximately 600 m/s.


## 4.3. Relative Stability of Fused Silica,α-Quartz and Stishovite

To further help rationalize the changes in the partial pair correlation functions induced by the ballistic impact, the quantum-mechanical DFT calculation results pertaining to the relative stability of fused silica,α-quartz and stishovite, as well as of the energy barriers associated with the fused-silica → α-quartz and fused-silica → stishovite transition states, are presented and discussed in this section. The relative room-temperature potential energies of the fused silica, α-quartz and stishovite as a function of pressure are presented in Figure 11(a). Examination of the results displayed in this figure reveals that, at the ambient pressure, α-quartz is the most stable form of SiO2, followed by fused silica and then by stishovite. At a pressure of 50 GPa, on the other hand, the relative positions of α-quartz and fused silica, with respect to thermodynamic stability, have been exchanged.Figure 11
Relative room-temperature potential energies of fused silica,α-quartz and stishovite, as a function of pressure, in the (a) absence and (b) presence of 5% shear.
(a)
(b)In Figure10(b), the relative room-temperature potential energies of the fused silica, α-quartz and stishovite as a function of pressure, in the presence of 5% shear, are presented. A comparison of the results displayed in Figures 11(a)-11(b) reveals that the application of 5% shear strain did not change the relative stability of fused silica, α-quartz and stishovite at the ambient pressure. In sharp contrast, at 50 GPa pressure, the application of 5% shear has been found to make the stishovite the most stable phase, followed by α-quartz and then fused silica. This finding suggests that, in the presence of shear, fused silica is more likely to undergo a devitrification conversion into stishovite and, to a lower extent, into α-quartz. This finding is important since the fused-silica region surrounding the projectile is generally subjected to high shear, in addition to high pressures. This explains why, in Figure 7, this region was found to undergo an extensive fused-silica → stishovite conversion.The likelihood for the aforementioned devitrification processes is affected not only by the relative stabilities of fused silica,α-quartz and stishovite, but also by the size of the energy barrier associated with the respective transition state. Variations in the fused-silica → α-quartz and fused-silica → stishovite transition state energy barriers with pressure, in the absence and the presence of 5% shear, are shown, respectively, in Figures 12(a)-12(b). Examination of the results displayed in these figures, at a pressure of 50 GPa, reveals that the energy barrier for fused-silica devitrification, in the absence of shear, is slightly higher for the fused-silica → stishovite conversion. On the other hand, the energy barrier for fused-silica devitrification, in the presence of shear, is slightly higher for the fused-silica → α-quartz conversion. These findings suggest that, in the presence of shear, and at pressures as high as 50 GPa, fused silica is more likely to convert to stishovite than α-quartz.Figure 12
Variations in the fused-silica →α-quartz and fused-silica → stishovite transition state energy barriers with pressure, in the (a) absence and (b) presence of 5% shear.
(a)
(b)Figure13 shows the conversion of an initially amorphous SiO2 structure into a stishovite-like structure under the influence of high pressure and shear. The transition state associated with this conversion is also shown in this figure. For clarity, only a small SiO2 region consisting of three Si and fifteen O atoms is displayed in Figure 13. Due to such a small size of the region, some O atoms appear to be nonbonded. These O atoms are bonded to Si atoms, but the Si atoms that they are bonded to are not shown in this figure. The Si atoms in question are also bonded to some of the bonded O atoms displayed in Figure 13. Examination of this figure reveals that, as expected, the fused-silica state contains only fourfold coordinated Si and twofold coordinated O atoms, while the stishovite state contains sixfold coordinated Si and threefold coordinated O atoms. The Si- and O-atom coordination in the transition state involves fourfold and fivefold coordinated Si and twofold and threefold coordinated O atoms. To help with the interpretation of the results displayed in Figure 13, the same three Si atoms appearing in the three configurations are denoted using labels “Si 1,” “Si 2,” and “Si 3.” It should be recalled that the results presented in Figures 12(a)-12(b) reveal that the maximum energy associated with the transition state is lowered and, thus, the conversion of the fused-silica → stishovite becomes feasible, under high pressure and shear.Figure 13
Conversion of an initially amorphous SiO2 structure into a stishovite-like structure under the influence of high pressure and shear. The associated transition state is also shown.


## 5. Summary and Conclusions

Based on the results obtained in the present work, the following summary remarks and main conclusions can be drawn.(1) By determining the Si–Si, Si–O, and O–O partial radial distribution functions, and the coordination numbers for Si and O atoms, before and after a ballistic impact experienced by a fused-silica target-plate, potential devitrification of the initial material and its conversion intoα-quartz and stishovite are investigated.(2) The investigation is carried out computationally through the use of all-atom nonequilibrium molecular-dynamics simulations.(3) To rationalize the results obtained, quantum-mechanical density functional theory computational analyses are employed. These analyses helped to establish the effect of high pressure and shear on the relative stability of the amorphous fused silica, crystallineα-quartz and crystalline stishovite, and on the energy barriers associated with the fused-silica → α-quartz and fused-silica → stishovite conversion reactions.(4) The results obtained demonstrated that under pressures on the order of 50 GPa, and in the presence of shear, stishovite becomes the most stable form of SiO2 and the energy barrier associated with the fused-silica → stishovite conversion is the smallest. Since the fused-silica target-plate regions beneath and surrounding the hard projectile experience high pressures and shear, formation of stishovite and, to a lower extent, α-quartz from fused silica, as observed in the all-atom molecular-level computational analyses of the ballistic-impact problem, is justified.



---

*Source: 650625-2015-04-29.xml*