Materials and Process Simulation Center,
Beckmann Institute, (139-74),
California Institute of Technology, Pasadena, CA 91125
Over the last decade increasing number of researchers from different fields of science attempted to investigate various aspects of silica, structure, melting, phase behavior, energetics using molecular mechanics, molecular dynamics and monte carlo techniques. If silica is exposed to shock wave, various low density silica polymorphs and its glass transform to high density stishovite along the Hugoniot isentrope. Several years ago Stolper and Ahrens pointed out a difference in the alpha-quartz to stishovite vs. glass to stishovite transformation, "this transformation is more displacive than reconsructive". In this study, we aim at shedding some light into this sound conjecture by analizing the microstructures of the model systems along the transformation path.
Central to all classical computational approaches is the representation of atomic interactions. We recently developed a simple 2-body interaction potential with long range electrostatic interactions are handled via variable charge model whereas a short range morse term supplies the combination of valence and van der Waals interactions. The partial charges are determined by the charge equilibration method of Rappe and Goddard. We have succesfully applied this potential to characterize silica polymorphs and especially the microstructure of glassy form. 
Our objective is to present a systematic study on crystal-crystal, glass- crystal transformations observed in silica and assess the role of each factor, such as temperature, pressure loading rate on the transition pressure. Also using the microscopic level information to investigate the details of transformation to relate to the experiments, which are mostly the shock wave experiments at these pressures, and loading rates. Furthermore, to develop an understanding of behavior of various silica polymorphs and glasses subjected to high pressures.
In the Theory and Method section, we briefly describe the interaction potentials used in this study, and procedures used in molecular mechanics and molecular dynamics simulations. In the Results section, we present and discuss the simulation results we obtained.
Interaction Potentials for Silica
Silica can exist in several crystalline polymorphs as well as in the amorphous state. A number of interatomic potential functions have been developed which can, with reasonable accuracy, reproduce the structural characteristics of one or more of the silica polymorphs. For studying silica most widely used forms of interaction potentials have a 2-body Born-Mayer-Huggins[4-8] term along with an attractive dispersion term, i.e. an exp-6 potential. However they may slightly differ in handling the electrostatic interaction, screened coulomb (erfc -complementary error function) to formal charge models. Some model potentials also include a Stillinger-Weber type 3-body term. 
We developed the new pure 2-body interaction potential for simulation of all phases of silica, glasses and melts and extended it to alimuna as well. This new MS-Q interatomic potential differs from the previous potentials by allowing the charges to adjust self-consistently to changes in the environment. These charges are determined via the Charge Equilibration (Qeq) method of Rappe and Goddard.
A Morse type function is assumed for short range bonding forces. Thus the MS-Q two body interatomic potential function has the form:
where qi are the charges and Rij is the distance between atoms i and j. Do, Ro and gamma are adjustable parameters.
The charges for the Si atoms in the optimized silica polymorphs varied between 1.22 and 1.4. This compares with potential derived charges and Mulliken charges obtained from a Si(OH)4 cluster ab initio Hartree Fock calculation using a 6-31G** basis set of 1.528 and 1.493 respectively. Thus, the Qeq charges appear reasonable. The Qeq charges for O vary from -0.61 to -0.70.
In our model glass simulations with this potentials, we found the pair distribution functions, g(r), and angular distribution functions for O-Si-O and Si-O-Si are in good agreement with experiment.
Model Systems and Simulation Procedure
In our studies we employed both molecular mechanics (for optimization of the structure and lattice parameters/volume, i.e. objective function is the enthalpy) and variable cell Parrinello-Rahman-Nos\' e molecular dynamics. Model systems are supercells of unit alpha-quartz, coesite and stishovite, and glasses which were obtained from melt through lengthy quench procedure with sizes ranging 576 atom to 640 atoms. Energy, force and stresses in the simulations described below all employed Ewald summation method.
In molecular mechanics simulations we did not impose any crystal symmetry constraints. We used an external hydrostatic pressure, and the pressure on the model systems were varied by 1 GPa intervals upto 10 GPa, 2 GPa intervals upto 30 GPa, and 5 GPa intervals beyond 30 GPa. In order to explore the release ``isentrope'' we also re-traced the loading curve backwards for stishovite and alpha-quartz cases only. Termination criteria for optimization was 0.1 kcal/mol/A. The results of these are presented in the next section.
Molecular dynamics simulations are performed under constant temperature condition, in contrast to shock wave experiments. Furthermore, the imposed pressure load or profile was isotropic rather than uniaxial. In these simulations 1 femtosecond is used as the integration step. In the course of simulations the pressure is incremented with a specific period ( step function compressions) to give approximate pressure loading rates 0.1, 0.25, 0.5, and 1.0 GPa/ps. The end points of each pressure segment were further continued for additional 25 ps subjected to the same pressure in order to study the micro structure and obtain averages.
Figure 1. Molecular Mechanics simulation of Equation of State for silica polymorphs and a model silica glass: a) Pressure vs. Density for alpha-quartz , coesite, stishovite and glass models. b) Loading and release isotherms (at 0 K) for quartz and stishovite
The pressure-density relationships reveal sharp transitions from alpha-quartz to stishovite (with a number of defects) at 24 GPa, and at 20 GPa for coesite, whereas transition in glass is gradual, as the pressure is increased the 5- 6- coordinated local structures slowly emerges. Even at a pressure as high as 100 GPa, the ratio of 6-coordinated atoms are less than that of structures obtained from quartz and coesite. When the loading curves retraced for stishovite and alpha-quartz and the transformation to stishovite from alpha-quartz is found to be not reversible. Resulting structure has a density of 3.7 g/cm3. From shock wave release isentropes similar observations are made.
We studied the rate and temperature dependence of the phase transformation using isothermal-isotension molecular dynamics, to mimic a non-equilibrium process where in the course of simulation the system is subjected to increasing isotropic pressure as step functions. The time interval between steps are varied to investigate the loading rate dependence of pressure loading. The rates are varied from 0.1 GPa to 1.0 GPa per picoseconds. The simulations for 1 GPa/ps rate were repeated at two higher temperatures 500 and 700 K, to investigate the temperature dependence of the transition preassure in the case of alpha-quartz. These temperatures are considerably lower than the temperatures observed in shock wave experiments. The results of pressure loading rate and temperature effects are displayed in Fig. 2.a and 2.b respectively.
Figure 2. Temperature and pressure loading rate dependence of transformation pressure for alpha-quartz to stishovite; a) Temperature dependence, pressure loading rate is 1.0 GPa/ps, b) Pressure loading rate dependence, T = 300 K
The temperature dependence of transformation pressure is stronger than the rate dependence, at least for the rates we have utilized in this study. The phase transition pressure changes by about 3 GPa for a temperature increase of 400 K. Experimental Hugoniot data show evidence of transition starting at about 14 GPa for 500 K , which is in agreement with our result, 15-16 GPa at 500 K.
Finally, to check the conjecture of Stolper and Ahrens, we analyzed the trajectories for microstructural evolution. alpha-quartz or glass to stishovite transition can easily be characterized by the variation of the distribution of O-Si-O angles in the structure. This angle is 100 % tetrahedral for quartz, dominantly tetrahedral, but not all, for glass structure, however in stishovite this angle is octahedral, i.e. 90 degrees. In Figure 3.a and 3.b we present the variation of the O-Si-O angles for quartz and silica glass. We see an abrupt change in the distribution for quartz whereas the change in glass is continuous even at pressures upto 120 GPa. The angle distribution for stishovite like structure obtained from the model glass is almost twice broader at the octahedral angle value. As a note, the peaks at 45 and 135 degrees are due to the silicon bonds to oxygens on the bisectors of 90 degree O-Si-O angles. We found that the transformations from alpha-quartz to stishovite are reconstructive whereas the transformation in the glass is displacive, with the peaks of the angle distributions progressively approaching to the 6-fold coordinated values which correspond to stishovite.
Figure 3. Microscopic analysis of pressure induced transformations using O-Si-O angle distribution as a precursor: a) alpha-quartz, b) model silica glass