Article pubs.acs.org/JPCC
Quantum Nuclear Effects in Stishovite Crystallization in ShockCompressed Fused Silica Yuan Shen† and Evan J. Reed*,‡ †
Department of Physics and ‡Department of Materials Science and Engineering, Stanford University, 496 Lomita Mall, Stanford, California 94305, United States ABSTRACT: Quantum nuclear effects include zero-point energy for each vibrational mode and potentially significant deviations of the heat capacity from classical values. While these effects play important roles in shock-induced chemical reactions and phase transitions, they are absent in classical atomistic shock simulations. To address this shortcoming, we have studied the quantum nuclear effects for stishovite crystallization in shock-compressed fused silica by employing the quantum bath multiscale shock technique, which couples the shock simulation with a colored-noise Langevin thermostat. We find that this semiclassical approach gives shock temperatures as much as 7% higher than classical simulations near the onset of stishovite crystallization in silica. We have also studied the impact of this approach on the kinetics of stishovite crystallization and the position of high-pressure silica melt line. We further describe a systematic way of setting up the parameters for the quantum thermal bath and quantum bath multiscale shock techniques. colored noise Langevin thermostat17 has been developed by Qi et al. (QBMSST) and applied to shocked methane.18 The onset of chemistry is reported to be at a pressure on the shock Hugoniot that is 40% lower than observed with classical molecular dynamics and temperature corrections are approximately a factor of 2 (∼1000 K) for some density range on the Hugoniot. These are sizable effects that may exceed typical errors associated with representation of the interatomic potential energy surface. These works have shown significant quantum nuclear effects in shock induced chemical reactions for materials with a large population of hydrogen atoms and thus more likely to exhibit quantum effects. Here we study the role of quantum nuclear effects in shock-induced solid−solid phase transitions in fused silica by employing the quantum bath multiscale shock technique (QBMSST) to study the quantum nuclear effects in shock induced stishovite crystallization.19 The multiscale shock technique employs approximate hydrodynamics to reproduce the kinetics of shocked materials and enables the use of expensive energy models including density functional theory.1−6 Large scale MSST simulations of stishovite formation in shocked fused silica and quartz have been recently reported by the present authors.19 These simulations employed classical dynamics of the atoms. The mechanism of stishovite formation for both fused silica and quartz is found to be a homogeneous nucleation and grain growth process for a chemically pure material. Recent time-resolved X-ray diffraction measurements20 of the shock-induced fused silica to stishovite transition have been reported to suggest a homogeneous nucleation and growth mechanism on nanosecond time scales at pressures lower than
1. INTRODUCTION Shock compression of materials enables the study of materials at pressures and temperatures higher than can be obtained under static conditions. It has historically been one of the primary tools used to study materials at high pressures. Shock experiments are often complicated by chemical reactions, phase transitions, and other kinetic processes that are intrinsic to the relatively short time scales of these experiments. Such kinetic processes can complicate deduction of the thermodynamically stable states. Atomistic simulations play a complementary role to these experiments, supplying microscopic information that is difficult to measure experimentally. The past decades have witnessed atomistic modeling methods for shock compression push to increasingly accurate approaches, including electronic structure methods and a variety of methods aimed at speeding the simulations up using specialized thermostats and boundary conditions.1−9 As these methods become more accurate and introduce less error in the results, the error associated with the usual approximation of the nuclei as classical particles warrants consideration. Classical molecular dynamics (MD) simulations do not include any quantum nuclear effects. Quantum treatment of the vibrational modes will introduce zero-point energy into the system, alter the energy power spectrum, and change the heat capacity from the classical value. Path integral simulations can accurately treat nuclear vibrational effects,10,11 but the computational effort required for these methods can be prohibitive in many cases. Theoretical postprocessing quantum corrections of DFT-based simulations of shock compressed water and methane have been conducted by Goldman et al.,12 where Hugoniot temperature corrections as much as 30% are reported, providing better agreement with experimental measurements. Subsequently, a self-consistent method that couples a multiscale shock simulation technique13−16 (MSST) to a quantum thermal bath described by a © 2016 American Chemical Society
Received: May 19, 2016 Revised: July 12, 2016 Published: July 13, 2016 17759
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C
temperature changes must be accounted for. In classical simulations, the MSST scheme has been used to describe the thermodynamics and kinetics of shock processes out to nanosecond and longer time scales.6,13,19 The thermodynamic constraints of a steady shock in a continuum are applied to MD simulations, restraining the system energy to the shock Hugoniot energy and the system stress to Rayleigh line stress.6 The computational cell can be thought of as a Lagrangian particle embedded in material that the shock is passing by. Longer simulation time corresponds to states further behind the shock front. In the QBMSST, coupling between the quantum thermal bath temperature and the shocked system is described by the equation of motion for the time-dependent target quantum temperature Tqm through a unitless coupling constant η,18
the molecular dynamics. This phase transition and its kinetics have fundamental impact for laser damage processes and meteorite impact events.
2. COMPUTATIONAL METHODS 2.1. Simulation Details. Methods that incorporate some aspects of quantum nuclear effects using different Langevin-type colored thermostats have been introduced by Dammak et al.17 and Ceriotti et al.21 to circumvent the high computational cost of more accurate path integral approaches. These methods make the quasiharmonic assumption that each vibrational mode of the system is a harmonic oscillator and that each vibrational mode should exhibit the free energy of a quantum harmonic oscillator. The thermostat developed by Dammak et al.17 has been used to compute temperature dependence of the lattice constant and heat capacity of MgO crystal in a temperature regime where quantum statistical effects cannot be neglected.17 Barrat et al.22 have described an efficient implementation of the quantum thermal bath due to Dammak et al.17 that is utilized by the QBMSST algorithm. We have also implemented this method as a user package quantum thermal bath (QTB) in the LAMMPS molecular dynamics simulator.23 We briefly summarize the algorithm of ref 22 as follows. The atoms obey an equation of motion of the Langevin form subject to random force Riλ(t) and a dissipative force −γmiẋiλ in addition to the interatomic force f iλ mixï λ = fiλ − γmixi̇ λ + R iλ(t )
qm
Ṫ (t ) = γη
Nf − 1
∑ k =−Nf
( ) sin( )
α sin
kπ 2αNf
kπ 2Nf
(1)
E(|ωk|) cos(ωktn) (2)
⎡ ⎤ 1 1 ⎥ is the energy expectation where E(ω) = ℏω⎢ 2 + ℏω ⎢⎣ exp( k T qm ) − 1 ⎥ ⎦ B of a quantum harmonic oscillator with kB being the Boltzmann kπ constant and ωk = αN δt provides a grid in the frequency domain f
of size 2Nf. The random force is given by a convolution with H(tn) R iλ(tn) =
2γmi αδt
Nf − 1
∑ n ′=−Nf
(4)
Here M is the total mass of atoms in the computational cell, N is the total number of atoms, and ẽ is the instantaneous augmented total energy (shock Hugoniot energy)13,16,18 per unit mass. Time propagation of the QBMSST is carried out according to eqs 19 and 22 in ref 18 with the random and dissipative force terms given in eq 3. Both MSST and QBMSST are implemented in the publically available distribution of the LAMMPS molecular dynamics code.23 We find that the computational expense of this implementation of QBMSST is comparable to MSST, with wall clock run times for the simulations in this work (144 000 atoms parallelized over 16 cores) differing by approximately only 40%. 2.2. Modeling Amorphous Silica. The analytical pair potential we use is proposed by van Beest, Kramer, and van Santen24 (BKS) to model atomic interactions in silicon dioxide. This potential has been employed to describe properties of silica melt,25 quartz,24,26 stishovite,26 and amorphous states from low27 to high28 pressures. It has also been used to study silica phase transitions, including the amorphization of α-quartz,29 the α,βquartz phase transition,30 and shock induced stishovite formation.19 The BKS potential was fitted to experimentally measured lattice and elastic properties of α-quartz at ambient pressure and thus incorporates some indirect aspect of the quantum nature of silica in the parameters. Quantum nuclear effects most relevant here are derived from higher frequency vibrational frequencies, which are not explicitly influenced by any quantum nature of experimental fitting data. Furthermore, we are starting simulations from fused silica where a cutoff potential is employed. The use of classical MD with an interatomic potential fit to experimental data does not describe the quantum thermodynamics of the vibrational modes, including zero point energy, power spectrum and quantum heat capacity, that are here brought into the simulations by the integration schemes, i.e., QTB and QBMSST. The BKS potential, i.e., eq 5, consists of a Coulomb term and a Buckingham term.24 This potential has been reported to provide microscopic quantities in reasonable agreement with experimental measurements of radial distribution function, bond− bond angle distribution function, coordination numbers, ring sizes, and vibrational spectrum.27 Force field parameters of eq 5 are listed in Table 1.
for the i-th atom of mass mi and direction λ (1, 2, or 3). Here γ is the dissipation constant. In a molecular dynamics simulation with step length δt, random forces on the atoms can be updated on a coarser time grid every α steps, at times tn = nαδt. The random forces are generated though convolution with a noise filter H(tn) 1 H(tn) = 2Nf
M[e (̃ t ) − e (̃ t = 0)] 3NkB
riλ(tn − tn ′)H(tn ′) (3)
where riλ is a pseudorandom number series with zero mean, unit variance, and uncorrelated for different iλ pairs. For a single harmonic oscillator of a particular frequency, eqs 1 and 3 can be shown to provide the quantum energy expectation for the same frequency,17,22 i.e., E(ω), if γ is chosen to be sufficiently small. This method is memory usage efficient because only 2Nf steps of the noise filter H(tn) and the random number sequence riλ(tn) will be stored in the memory, independent of the total simulation length. While the quantum thermal baths of Dammak et al.17 and Ceriotti et al.21 apply to systems at a fixed temperature, the temperature varies dramatically during shock loading and 17760
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C Table 1. Force Field Parameters i-j
Aij (eV)
ρij (Å)
Cij (eV·Å6)
rmij (Å)
wij (eV)
O−O Si−Si Si−O
1388.773
0.36232
175
1.439
20.782
18003.7572
0.20521
133.5381
1.1936
−27.3050
⎧ κ(r − r m)2 + w ; r < r m ij ij ij ij ⎪ ij ⎪qq ⎛ ⎞ ⎪ i j e 2 + A e−rij / ρij − Cij − ⎜A e−rc / ρij − Cij ⎟ ; r m ij ij ij 6 6 ⎪ r rij rc ⎠ ⎝ ϕij(rij) = ⎨ ij ⎪ < rij < rc ⎪ ⎪ qiqj 2 e ; rc < rij ⎪ ⎩ rij
the number of time steps between updates of the random force, is chosen to ensure that the random noise spectrum frequencies exceed the Debye frequency or the highest vibrational frequency spectral content of the system. π max 0 ≤ k ≤ Nf ωk = > fD (7) αδt Here we take α = 4 and δt = 1 fs which provides π = 785.4 αδt THz, greater than the 50 THz maximum spectral content that we find for silica under initial and shock conditions. The friction coefficient γ is chosen so that the thermostat can establish a vibrational energy spectrum that closely resembles the spectrum expected if all of the vibrational modes are quantum harmonic oscillators. Reference 17 studied the infrared absorption spectrum of a MgO crystal and reported convergence for γ below 0.3 THz. For liquid methane, ref 18 showed that the classical temperature converges when γ is smaller than 5 THz for liquid methane. For any value of γ, we can obtain the vibrational energy or power spectrum by computing the velocity autocorrelation function ⟨miẋiλ(0)ẋiλ(t)⟩iλ where brackets denote averaging over atom indices i and directions λ, and then performing a Fourier transformation.
(5)
rmij
Within small distances as listed in Table 1, we replace the BKS form with a harmonic term to prevent unphysical particle fusion reported in ref 27. For temperatures below 6000 K and pressures below 130 GPa, we have not observed particle fusion. The Buckingham part is truncated and shifted at a cutoff distance rc = 6 Å to improve agreement with the experimentally measured density of the ambient condition amorphous phase.27 Fused silica is obtained using an annealing process. First we melt αquartz with 144 000 atoms at 6000 K, and then we quench it to 300 K with a cooling rate γ = 6 × 1011 K·s−1. The macroscopic and microscopic properties of the glass at this cooling rating has been reported to yield microscopic quantities of fused glass in ref 27 in closest agreement with experimental results, including radial distribution function, bond−bond angle distribution function, coordination numbers, ring sizes, and vibrational spectrum. This construction of the silica glass yields an initial density of 2260 kg·m−3, 2.7% higher than reported experimentally measured values.31 The elastic properties, melting temperature, activation energies for self-diffusion, and high pressure phase boundaries of this potential are quantified in ref 19 and its Supporting Information. 2.3. Simulation Parameters. Setting up the QBMSST simulation starts with setting the parameters for a classical MSST simulation. References 13 and 19 contain discussions about how to set up six key parameters. These are the cell mass-like parameter Q = 40.0 amu2·Å6 which is chosen to enable the initial elastic (fastest portion) compression on a time scale ∼1 ps, the optional artificial viscosity μ which is set to be zero here, the LAMMPS parameter tscale = 0.01 that will assist the compression but have no statistical impact on the subsequent simulation, p0 = 0 GPa as the initial pressure, v0 = 0.4425 cm3·g−1 as the initial volume. Values of these parameters for the shocked fused silica simulations performed here are the same as those used for the simulations in ref 19. The next step is to determine the parameters associated with the quantum bath and QBMSST. We re-equilibrate the fused silica of 144 000 atoms obtained through the classical annealing process of ref 27 with the quantum thermal bath at a specified initial temperature Tinit = 300 K. We choose Nf = 100 which is shown in refs 22 and 18 to be large enough to converge the classical temperature Tcl T cl =
1 3NkB
∑ mixi̇ 2λ iλ
κ = 100 eV·Å2 qO = −1.2 qSi = +2.4
+∞
P(ω) =
∫−∞
⟨mixi̇ λ(0)xi̇ λ(t )⟩iλ eiωt
dt 2π
(8)
For a classical MD simulation of a harmonic system, all vibration modes share the same energy, namely kBTcl. Therefore, the classical power spectrum Pcl(ω) is proportional to the phonon density of states (PDoS), i.e., g(ω). P cl(ω) = g (ω)kBT cl
(9)
Under the quantum quasi-harmonic assumption, each vibrational mode has on average an energy that depends on its ⎡ ⎤ 1 1 ⎥. Therefore, the frequency ω, E(ω) = ℏω⎢ 2 + ℏω ⎢⎣ exp( k T qm ) − 1 ⎥ ⎦ B vibrational energy spectrum of a harmonic system equilibrated with a quantum thermal bath is thus targeted to be eq 10.
P
⎡ ⎢1 (ω) = g (ω)ℏω⎢ + ⎢2 exp ⎣
PDoS + QHO
1
( ) ℏω kBT qm
⎤ ⎥ ⎥ − 1⎥ ⎦
(10)
The choice of γ impacts the degree to which the Langevin thermostat can achieve the spectrum PPDoS+QHO(ω). When γ is too large, the equation of motion eq 1 is sufficiently altered to provided unphysical dynamics and an energy spectrum that deviates from the target spectrum.22 Setting γ = 0 THz reduces eq 1 to the unthermostated case since the term −γmiẋiλ + Riλ becomes zero because the random forces Riλ are proportional to √γ. Another consideration for choice of γ is anharmonicity. While the random force terms pump energy into the system toward the spectrum PPDoS+QHO(ω), the anharmonic terms that couple modes of different frequencies re-equilibrate the energy
(6)
Minimization of the value of Nf subject to an accurate spectrum minimizes the memory requirement. The parameter α, 17761
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C spectrum to its classical form Pcl(ω) at a higher temperature. The choice of γ must be sufficiently large that rate of achieving the quantum spectrum on each mode exceeds the rate at which anharmonic couplings force the system toward the classical distribution Pcl(ω). The effect of γ on the energy spectrum is illustrated in Figure 1 for the initial state of fused silica. Here we have set the target
Figure 1. Power spectrum for different values of parameter γ. The power spectra are computed by making the quantum harmonic oscillator correction to the phonon density of states (PDoS+QHO), and by applying the quantum bath thermostat (γ values). PDoS+QHO is supposed to be the target power intensity. The quantum bath thermostat spectrum is closest to the approximate target spectrum (PDoS+QHO) when γ is neither too small (γ = 0.01 THz) nor too large (γ = 100 THz). Units of the power intensity are given as K·ps per atom so that the integral of the power spectrum yields the classical temperature.
Figure 2. Temperature profiles during shock compression for different values of parameter η, showing (a) quantum temperatures and (b) classical temperatures, i.e., eq 6, during the first 30 ps of shock compression. Values of η near 0.2 minimize over and under damping caused by rapid initial compression. Quantum and classical temperatures exhibit no oscillations at later times where volume changes are much slower as the phase change occurs.
3. RESULTS AND DISCUSSION 3.1. Power Spectrum of Fused Silica and Stishovite. With a choice of friction coefficient γ, we test this value to provide adequate description of the stishovite crystal after phase transition and is thus valid throughout the shock loading. Figure 3 shows the power density spectrum P(ω) of silicon (blue) and oxygen (red) atoms for both the initial and final phases. The
quantum bath temperature to be 300 K and plotted 2πP(ω)/kB such that the areas under curves correspond to classical temperatures. PPDoS+QHO(ω) is denoted by the green curve and shows the ideal target spectrum of the quantum thermal bath. The choice of γ = 1 THz (blue curve) shows a good agreement with the target PPDoS+QHO(ω) spectrum. When γ = 0.001 THz (orange curve), anharmonic terms overwhelm the thermostat and the power density becomes similar to a classical spectrum at Tcl = 530 K. On the large side, when γ = 100 THz (red curve), all vibrations are overdamped and the spectrum is washed out. For this work, we choose γ = 1 THz. After equilibrating the initial state with the quantum thermal bath at 300 K, the QBMSST initial energy parameter ẽ0 = −9.143 × 105 kJ·g−1 is set as the initial total energy. To run the QBMSST, the coupling constant η of eq 4 must be chosen. Equation 4 can be modified to smooth thermal fluctuations of the quantum bath temperature and minimize computational expense by replacing ẽ(t) by its time average over the past β MD steps.18 β
δT qm = γη
M ∑l = 1 [e (̃ t − lδt ) − e0̃ ] 3NkB
δt
(11)
The value of β should be chosen to be small compared to the 1 time scale for thermostat coupling, i.e., β ≪ γδt to allow Tqm to Figure 3. Power spectrum of Si and O atoms in (a) fused silica and (b) high pressure stishovite (60 GPa) at a quantum temperature of 300 K. The PDoS+QHO power spectrum is computed by making the quantum harmonic oscillator correction to the classical phonon density of states and the QTB curves are computed using the quantum bath thermostat. Units of the power intensity are given as K·ps per atom so that the integral of the power spectrum yields the classical temperature. In these simulations, γ is set to be 1 THz.
track changes in thermodynamic state, which happen at a time scale longer than 1/γ. For γ = 1 THz and δt = 1 fs, we have used β = 50 in our simulations. Figure 2 shows the first 20 ps of shock compressed fused silica for different values of η. In our silica simulations, we have chosen η = 0.2, similar to the value η = 0.33 employed in earlier work on liquid methane.18 17762
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C
in these two cases can be described by the shock Hugoniot energy equation
initial amorphous state in Figure 3a is at 1 atm, and the stishovite crystal phase in Figure 3b is at 60 GPa. Figure 3 also shows the target quasi-harmonic spectrum PPDoS+QHO(ω) of silicon and oxygen atoms for comparison. 3.2. Shift of the Hugoniot Induced by Quantum Nuclear Effects. Figure 4 shows Hugoniots for a range of
1 (p + pi )(Vf − Vi ) ≈ Efqm − E iqm 2 f
Efcl − E icl =
(12)
Here are energies corresponding to the initial and final state in classical MSST and QBMSST. Let us assume that the nonthermal portion of the energy (i.e., the T = 0 portion) change during the shock loading is ΔEcl,qm and we further neglect any heat exchange with the environment. Ecl,qm i,f
Efcl − E icl = ΔE cl +
∫T
Tfcl
C cl dT
i
≈ ΔE
qm
+
∫T
Tfqm
C qm dT
i
= Efqm − E iqm
(13)
Here ΔE − ΔE is the nonthermal portion of the energy change, which includes the zero point energy difference between the initial state and the final state. The heat capacities Ccl and Cqm describe the temperature dependence of the classical and quantum heat capacities, which can vary considerably with T because the initial temperature is well below the Debye cl
Figure 4. MD shock Hugoniot for fused silica using MSST (red, classical dynamics) and QBMSST (blue, with quantum nuclear effects). Open symbols show the temperature and pressure 10 ps (immediately) after shock pressure is achieved while solid symbols are at 1.5 ns. There are two states along the Hugoniot: crystallized (cyan block) and disordered (green block). For a fixed pressure or shock speed, the QBMSST result is at a higher temperature than the MSST result for both the disordered and crystallized states. Experimental data sets are plotted for comparison: up pointing triangles34 and the gray band35 are from laser shock experiments, down pointing triangles38,41 are gas-gun data. Kinetic effects likely play a role in the comparison to experimental data. See text for details.
temperature
qm
hfD kBTi
≈ 7 ≫ 1. Using the mean value theorem, the
quantum nuclear correction of the final state temperature can be qm cl cl qm expressed as follows with Tqm f > T* > Ti, Tf > T* > Ti, and C* ≡ cl cl cl Cqm(Tqm ), C ≡ C (T ) * * * ⎛ C cl − C qm ⎞ ⎛ ΔE qm − ΔE cl ⎞ Tfqm − Tfcl = ⎜⎜ * qm * ⎟⎟(Tfcl − Ti) − ⎜ ⎟ C qm ⎝ ⎠ ⎝ C* ⎠ *
shock speeds computed after 1.5 ns with QBMSST (solid blue) and classical MSST (solid red). Figure 4 also shows comparison of the simulations to available experimental data, including gasgun32,33 (purple symbols) and laser-driven34,35 (brown symbols and gray band) methods. Two distinct states characterize the simulated Hugoniots: a crystallized state (cyan box) and a disordered state (green box). Shortly after fused silica is initially compressed by the shock (10 ps), the temperature of the system is relatively low and the shock Hugoniot lies along the disordered state line. Subsequent stishovite nucleation and grain growth are observed only for simulations within a specific temperature range. In these cases, the state of the system after 10 ps is marked by the open symbols in Figure 4. Incomplete crystallization is observed in the lowest pressure simulation. Quantum nuclear effects are likely to be significantly less important at the temperatures of Figure 4 than they are at room temperature. This is particularly true for this system lacking hydrogen atoms. However, Figure 4 indicates that quantum nuclear corrections can be as great as a few hundred degrees Kelvin. For the disordered state, QBMSST shock temperatures are shifted by 150 to 250 K (4.5% to 7%) from the classical MSST shock simulations. In the crystallized stishovite phase, QBMSST temperatures are shifted by up to 300 K (8%) from the MSST temperatures. The quantum nuclear correction brings the shock Hugoniot closer to the experimental reports of the disordered phase Hugoniot. To obtain some qualitative insight into the factors that determine the sign and magnitude of the quantum nuclear shifts, we perform an approximate analysis here. Let us denote the temperature of the initial state as Ti, which is the same for both shock simulations with and without quantum nuclear correction. We further take the pressures and volumes of the classical and quantum simulations to be the same, so that the energy increase
(14)
The two terms on the right of eq 14 each independently have the potential to contribute to the temperature shift. The term −
(
ΔE qm − ΔE cl C qm *
) contains the difference between zero-point
vibrational energies of the initial and final states. In cases where the final state phase has a higher zero-point energy than cl the initial state, this term is negative, giving Tqm f − Tf < 0. This might be the case for covalently or ionically bonded solids where the Debye frequency is expected to increase. To quantify eq 15, we have computed the respective zeropoint energies for fused silica and high pressure stishovite. This can be done by integrating the zero temperature power spectrum ℏω g(ω) 2 . Alternatively, evaluating the classical temperature Tcl at extremely low target quantum temperature can also give the zero point energy under the quasiharmonic assumption. limT qm → 0K kBT cl =
∫ g(ω) ℏ2ω dω
(15)
Both methods indicate that fused silica has a zero point energy of 450 K/atom at 1 atm and stishovite has a zero-point energy of 520 K/atom at 60 GPa. In the case of shocked fused silica transforming to stishovite, the term in eq 14 −
(
ΔE qm − ΔE cl C qm *
)≈
−70 K leads to a lower temperature. The other term
(
C cl − C qm * qm * C *
)(T
cl f
− Ti) in eq 14 reflects the difference in heat
capacity between the quantum and classical systems. It is expected to be positive because Ccl* − Cqm * > 0 due to the monotonic increase of Cqm(T) with T to the classical limit. The 17763
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C QBMSST temperature shifts in Figure 4 are positive, indicating that this term dominates the overall correction. Moreover, this term indicates that higher final state temperatures Tclf will result cl in a bigger value of Tqm f − Tf , consistent with the QBMSST data in Figure 4. 3.3. Position of the Melting Line. The kinetics observed in these simulations shed some new light on the position of the melting line relative to the earlier obtained experimental data. The laser shock data reported in ref 35 (gray band in Figure 4) shows a temperature decrease with increasing pressure between 55 and 65 GPa. This temperature drop is also observed in earlier laser34 and gas gun32,33 experiments. If the crystalline phase Hugoniot points are superheated and metastable, the melt line will intersect the lower part of the cusp, as it has been located in earlier work.35,36 However, the present MD simulations and recent time-resolved X-ray diffraction experiments20 indicate that the amorphous phase in this region is metastable and will crystallize given sufficient time, placing the melt line at the top of the cusp. The latter interpretation is in closer agreement with earlier reported DFT calculations of the melt line in ref 37. One would not expect to observe a superheated crystalline phase in this case because the material is initially amorphous, is still amorphous immediately after the shock front, and will crystallize only if the amorphous state is metastable. A similar temperature drop has been reported for shocked quartz.36,38 Our earlier work indicates that quartz amorphizes immediately behind the shock front and thus follows the same dynamics as shocked fused silica, with the melt line intersecting the top of the cusp.19 3.4. Kinetics in the QBMSST Scheme. The colored thermostat has the potential to maintain the system out of classical equilibrium, potentially leading to kinetic rates for barrier hopping processes that do not correspond with physically relevant rates. For sufficiently high temperatures, the colored noise spectrum becomes flat like a classical Langevin thermostat, and one might expect the thermostat to do a better job with kinetics in this regime. Earlier work on stishovite crystallization19 indicates that the kinetics are controlled by a diffusive process. To study the capacity of the thermostat to describe these diffusion effects, Figure 5 shows the computed self-diffusion coefficients in the temperature range where we observe crystallization. Here we
compare the self-diffusion coefficients of silicon (solid) and oxygen (open) for classical MD (red) and QTB (blue) simulations within a wide range of temperatures. They are both computed by evaluating the time derivatives of mean squared displacement at constant pressure 40 GPa for 144 000 atoms and become identical as the temperature increases. The effective activation energy of the modified BKS potential is found through linear fitting to be QSi = 2.36 ± 0.1 eV, QO = 2.43 ± 0.1 eV for classical MD and QSi = 1.96 ± 0.1 eV, QO = 2.03 ± 0.1 eV for QTB. Indirect comparisons with the DFT results (green) can also be made from Figure 5. Self-diffusion coefficients for silicon (solid symbols) and oxygen (open symbols) are evaluated at different temperature at pressures of 25 to 60 GPa. These pressures differ from the 40 GPa pressure of the BKS results. Classical self-diffusion coefficients for the BKS potential (red) at 40 GPa fall within the range of the DFT results at lower and higher pressures for the range of temperatures relevant to this work. As a further comparison of the kinetics, we compare the crystallization rates for MSST and QBMSST simulations with similar temperature profiles. We consider a MSST shock with speed 6.8 km/s and a QBMSST shock with speed 6.6 km/s. They exhibit comparable temperatures of both the disordered and crystallized states. These correspond to the solid and open symbols at 47 GPa (blue) and 51 GPa (red) in Figure 4. We have performed the shape matching analysis19,40 on these two simulations that was employed in ref 19. A crystal phase has the neighboring atomic arrangement for individual atoms resembling each other. However, a disordered phase may only have short-range order but lacks longer range resemblance. The shape matching analysis explores this distinction and allows the separation of crystallized silica from the disordered state. Grain size, density, and the fraction of crystallized atoms were computed as described in ref 19. In Figure 6, we have plotted
Figure 5. Self-diffusion coefficients of fused silica vs temperature at 40 GPa. Si (solid) and O (open) self-diffusion coefficients are computed by evaluating the time derivatives of mean squared displacement at constant pressure 40 GPa and listed temperatures. Red symbols are results from classical MD within the NPT ensemble, while blue symbols are results from the quantum thermostat MD. The green symbols are diffusion coefficients from reported DFT results within the pressure range of the shock Hugoniot.39 The pressures (GPa) are also marked close to the DFT data points.
Figure 6. Temperature, grain density, maximum grain size, and crystallization fraction of shocked fused silica using MSST (red) and QBMSST (blue). Shock speeds of the MSST and QBMSST simulations are 6.8 km/s (51 GPa) and 6.6 km/s (47 GPa). Similar temperature profiles lead to similar rates of nucleation and grain growth in these two simulations. White, green, and cyan regions show different stages of stishovite formation for shocked fused silica: the nucleation stage (white), the explosive grain growth stage (green), and the coalescence stage (cyan) as described in ref 19. 17764
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C
simulations. For linear explosive growth, classical simulations yield ln(Ke0/Å·ns−1) = 16.9(±0.6) and Qe = 3.59 eV(±0.2 eV) while the colored thermostat gives ln(Ke0/Å·ns−1) = 16.3(±1.0) and Qe = 3.42 eV(±0.3 eV). For the t1/3 coarsening stage, classical simulations yield ln(Kc0/nm3·ns−1) = 15.2(±0.5) and Qe = 3.85 eV(±0.2 eV) while the colored thermostat gives ln(Kc0/ nm3·ns−1) = 15.5(±0.4) and Qe = 3.89 eV(±0.2 eV). Taken together, these comparisons of kinetics indicate that the shock induced nucleation and grain growth are similar in the QBMSST scheme at the same quantum temperature as classical simulations. The effect of quantum nuclear effects on the kinetics in these simulations can be primarily understood to be a temperature shift of a few hundred Kelvins, which can lead to factor of 2 or 3 differences in kinetic rates due to the exponential forms of both nucleation and grain growth.
temperature profiles, grain densities, maximum grain sizes, and crystallization fractions for these two simulations. Comparable behavior exhibited by these quantities indicates that the QBMSST is likely to provide an acceptable description of the kinetics in this case. The kinetics for QBMSST and MSST are similar for similar pressures and temperatures. However, the Hugoniots for QBMSST and MSST differ, leading to apparently different kinetics for the same shock speed. Figure 6 shows a MSST shock speed 6.6 km/s and a QBMSST shock speed of 6.8 km/s. Figure 7 further validates this observation by studying the temperature dependence of the grain growth rate in both the
4. CONCLUSIONS In this work, we find that quantum nuclear effects play a role in the Hugoniot and kinetics of shock-induced phase transition in silica. While it is not surprising that such effects occur in hydrogen-containing materials, it might be surprising that they play any role in silica with its heavier atoms. Though the shockcompressed state ends up at a temperature where it is well described by classical mechanics, significant quantum nuclear effects in the preshocked state still play a role in determining the final temperature. We find that quantum nuclear corrections shift the Hugoniot curve by a few hundred Kelvins at the pressures studied here, which may hasten the phase transformation to a modest degree. The kinetics observed in these simulations suggests that the melt line is located at the high temperature side of the cusps in the Hugoniot. We provide details and figures describing a systematic way of setting up parameters for the quantum thermal bath and quantum bath multiscale shock technique schemes implemented in the publicly available LAMMPS code. We have also shown that, in the classical regime, the colored-noise thermostat reproduces the classical kinetics of nucleation and grain growth for shocked fused silica.
Figure 7. Homogeneous nucleation and grain growth model obtained through QBMSST MD simulations. Different colors indicate different shock speeds: 6.4 km/s (red), 6.6 km/s (cyan), 6.8 km/s (blue), 7.0 km/ s (purple), 7.2 km/s (green), 7.4 km/s (pink), and 7.6 km/s (red). Orange curves are nucleation and grain growth models fitted to QBMSST simulations. Dashed back curves are those obtained in classical MSST simulations plotted for comparison.19 (a and b) Grain growth rate vs temperature for explosive growth and deep coalescence stages. Here the explosive growth rate has a unit of Å/s while the coalescence rate has a unit of nm3/s. Error bars indicate temperature variations within the corresponding stage.
■
Corresponding Author
*E-mail:
[email protected]. Tel: 650-723-2971. Fax: 650725-4034.
explosive growth and coalescence stages. Identifications of these two stages through the time evolution of the grain density are given in ref 19. Since it is unlikely that the Coulombic form of the BKS can provide a reliable description of phase separation into silicon-rich and oxygen rich phases, we do not expect to find any stoichiometric phase separation in the calculations performed here. In this case, a general grain growth law for the average grain size can be written as n −n = K Δt
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Our work was supported in part by the Army High Performance Computing Research Center. Y.S. is supported by a William R. Hewlett Stanford Graduate Fellowship.
■
(16)
where we take n = 1 for the interface controlled explosive growth stage and n = 3 for the coarsening stage based on MSST results.19 Temperature dependence of the grain growth rate K is hypothesized to follow the Arrhenius form. ⎛ Q ⎞ K = K 0 exp⎜ − ⎟ ⎝ kBT ⎠
AUTHOR INFORMATION
REFERENCES
(1) Goldman, N.; Reed, E. J.; Kuo, I.-F. W.; Fried, L. E.; Mundy, C. J.; Curioni, A. Ab Initio Simulation of the Equation of State and Kinetics of Shocked Water. J. Chem. Phys. 2009, 130 (12), 124517. (2) Reed, E. J.; Riad Manaa, M.; Fried, L. E.; Glaesemann, K. R.; Joannopoulos, J. D. A Transient Semimetallic Layer in Detonating Nitromethane. Nat. Phys. 2008, 4 (1), 72−76. (3) Mundy, C. J.; Curioni, A.; Goldman, N.; Kuo, I.-F. W.; Reed, E. J.; Fried, L. E.; Ianuzzi, M. Ultrafast Transformation of Graphite to Diamond: An Ab Initio Study of Graphite under Shock Compression. J. Chem. Phys. 2008, 128 (18), 184701.
(17)
Figure 7 compares the temperature dependence of the growth rate for MSST (black dashed) and QBMSST (orange solid) 17765
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
Article
The Journal of Physical Chemistry C (4) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. Nitrogen-Rich Heterocycles as Reactivity Retardants in Shocked Insensitive Explosives. J. Am. Chem. Soc. 2009, 131 (15), 5483−5487. (5) Goldman, N.; Reed, E. J.; Fried, L. E.; Kuo, I.-F. W.; Maiti, A. Synthesis of Glycine-Containing Complexes in Impacts of Comets on Early Earth. Nat. Chem. 2010, 2 (11), 949−954. (6) Reed, E. J. Electron-Ion Coupling in Shocked Energetic Materials. J. Phys. Chem. C 2012, 116 (3), 2205−2211. (7) Bolesta, A. V.; Zheng, L.; Thompson, D. L.; Sewell, T. D. Molecular Dynamics Simulations of Shock Waves Using the Absorbing Boundary Condition: A Case Study of Methane. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76 (22), 224108. (8) Maillet, J.-B.; Mareschal, M.; Soulard, L.; Ravelo, R.; Lomdahl, P. S.; Germann, T. C.; Holian, B. L. Uniaxial Hugoniostat: A Method for Atomistic Simulations of Shocked Materials. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 63 (1), 016121. (9) Zhakhovsky, V. V.; Budzevich, M. M.; Inogamov, N. A.; Oleynik, I. I.; White, C. T. Two-Zone Elastic-Plastic Single Shock Waves in Solids. Phys. Rev. Lett. 2011, 107 (13), 135502. (10) Chakravarty, C. Path Integral Simulations of Quantum LennardJones Solids. J. Chem. Phys. 2002, 116 (20), 8938−8947. (11) Militzer, B.; Ceperley, D. M. Path Integral Monte Carlo Calculation of the Deuterium Hugoniot. Phys. Rev. Lett. 2000, 85 (9), 1890−1893. (12) Goldman, N.; Reed, E. J.; Fried, L. E. Quantum Mechanical Corrections to Simulated Shock Hugoniot Temperatures. J. Chem. Phys. 2009, 131 (20), 204103. (13) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. A Method for Tractable Dynamical Studies of Single and Double Shock Compression. Phys. Rev. Lett. 2003, 90 (23), 235503. (14) Reed, E. J.; Fried, L. E.; Henshaw, W. D.; Tarver, C. M. Analysis of Simulation Technique for Steady Shock Waves in Materials with Analytical Equations of State. Phys. Rev. E 2006, 74 (5), 056706. (15) Reed, E. J.; Maiti, A.; Fried, L. E. Anomalous Sound Propagation and Slow Kinetics in Dynamically Compressed Amorphous Carbon. Phys. Rev. E 2010, 81 (1), 016607. (16) Manaa, M. R. Chemistry at Extreme Conditions; Elsevier: Amsterdam, 2005. (17) Dammak, H.; Chalopin, Y.; Laroche, M.; Hayoun, M.; Greffet, J.-J. Quantum Thermal Bath for Molecular Dynamics Simulation. Phys. Rev. Lett. 2009, 103 (19), 190601. (18) Qi, T.; Reed, E. J. Simulations of Shocked Methane Including SelfConsistent Semiclassical Quantum Nuclear Effects. J. Phys. Chem. A 2012, 116 (42), 10451−10459. (19) Shen, Y.; Jester, S. B.; Qi, T.; Reed, E. J. Nanosecond Homogeneous Nucleation and Crystal Growth in Shock-Compressed SiO2. Nat. Mater. 2015, 15 (1), 60−65. (20) Gleason, A. E.; Bolme, C. A.; Lee, H. J.; Nagler, B.; Galtier, E.; Milathianaki, D.; Hawreliak, J.; Kraus, R. G.; Eggert, J. H.; Fratanduono, D. E.; et al. Ultrafast Visualization of Crystallization and Grain Growth in Shock-Compressed SiO2. Nat. Commun. 2015, 6, 8191. (21) Ceriotti, M.; Bussi, G.; Parrinello, M. Nuclear Quantum Effects in Solids Using a Colored-Noise Thermostat. Phys. Rev. Lett. 2009, 103 (3), 030603. (22) Barrat, J.-L.; Rodney, D. Portable Implementation of a Quantum Thermal Bath for Molecular Dynamics Simulations. J. Stat. Phys. 2011, 144 (3), 679−689. (23) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (24) Van Beest, B. W. H.; Kramer, G. J.; Van Santen, R. A. Force Fields for Silicas and Aluminophosphates Based on Ab Initio Calculations. Phys. Rev. Lett. 1990, 64 (16), 1955. (25) Saika-Voivod, I.; Poole, P. H.; Sciortino, F. Fragile-to-Strong Transition and Polyamorphism in the Energy Landscape of Liquid Silica. Nature 2001, 412 (6846), 514−517. (26) Saika-Voivod, I.; Sciortino, F.; Grande, T.; Poole, P. H. Phase Diagram of Silica from Computer Simulation. Phys. Rev. E 2004, 70 (6), 061507.
(27) Vollmayr, K.; Kob, W.; Binder, K. Cooling-Rate Effects in Amorphous Silica: A Computer-Simulation Study. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54 (22), 15808. (28) Tse, J. S.; Klug, D. D.; Le Page, Y. High-Pressure Densification of Amorphous Silica. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46 (10), 5933. (29) Tse, J. S.; Klug, D. D. Mechanical Instability of α-Quartz: A Molecular Dynamics Study. Phys. Rev. Lett. 1991, 67 (25), 3559. (30) Müser, M. H.; Binder, K. Molecular Dynamics Study of the α−β Transition in Quartz: Elastic Properties, Finite Size Effects, and Hysteresis in the Local Structure. Phys. Chem. Miner. 2001, 28 (10), 746−755. (31) Mozzi, R. L.; Warren, B. E. The Structure of Vitreous Silica. J. Appl. Crystallogr. 1969, 2 (4), 164−172. (32) Sugiura, H.; Kondo, K.; Sawaoka, A. Shock Temperatures in Fused Silica Measured by Optical Technique. J. Appl. Phys. 1982, 53 (6), 4512−4514. (33) Lyzenga, G. A.; Ahrens, T. J.; Mitchell, A. C. Shock Temperatures of SiO2 and Their Geophysical Implications. J. Geophys. Res. 1983, 88 (B3), 2431−2444. (34) Hicks, D. G.; Boehly, T. R.; Eggert, J. H.; Miller, J. E.; Celliers, P. M.; Collins, G. W. Dissociation of Liquid Silica at High Pressures and Temperatures. Phys. Rev. Lett. 2006, 97 (2), 025502. (35) Millot, M.; Dubrovinskaia, N.; Č ernok, A.; Blaha, S.; Dubrovinsky, L.; Braun, D. G.; Celliers, P. M.; Collins, G. W.; Eggert, J. H.; Jeanloz, R. Shock Compression of Stishovite and Melting of Silica at Planetary Interior Conditions. Science 2015, 347 (6220), 418−420. (36) Akins, J. A.; Ahrens, T. J. Dynamic Compression of SiO2: A New Interpretation. Geophys. Res. Lett. 2002, 29 (10), 31−1. (37) Usui, Y.; Tsuchiya, T. Ab Initio Two-Phase Molecular Dynamics on the Melting Curve of SiO2. J. Earth Sci. 2010, 21 (5), 801−810. (38) Lyzenga, G. A.; Ahrens, T. J.; Mitchell, A. C. Shock Temperatures of SiO2 and Their Geophysical Implications. J. Geophys. Res. 1983, 88 (B3), 2431−2444. (39) Karki, B. B.; Bhattarai, D.; Stixrude, L. First-Principles Simulations of Liquid Silica: Structural and Dynamical Behavior at High Pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76 (10), 104205. (40) Keys, A. S.; Iacovella, C. R.; Glotzer, S. C. Characterizing Structure Through Shape Matching and Applications to Self-Assembly. Annu. Rev. Condens. Matter Phys. 2011, 2 (1), 263−285. (41) Sugiura, H.; Kondo, K.; Sawaoka, A. Shock Temperatures in Fused Silica Measured by Optical Technique. J. Appl. Phys. 1982, 53 (6), 4512−4514.
17766
DOI: 10.1021/acs.jpcc.6b05083 J. Phys. Chem. C 2016, 120, 17759−17766
