Simulations of Shocked Methane Including Self-Consistent

Article pubs.acs.org/JPCA

Simulations of Shocked Methane Including Self-Consistent Semiclassical Quantum Nuclear Effects Tingting Qi and Evan J. Reed* Department of Materials Science and Engineering, Stanford University, Stanford, California 94305, United States ABSTRACT: A methodology is described for atomistic simulations of shock-compressed materials that incorporates quantum nuclear effects on the fly. We introduce a modification of the multiscale shock technique (MSST) that couples to a quantum thermal bath described by a colored noise Langevin thermostat. The new approach, which we call QB-MSST, is of comparable computational cost to MSST and self-consistently incorporates quantum heat capacities and Bose−Einstein harmonic vibrational distributions. As a first test, we study shock-compressed methane using the ReaxFF potential. The Hugoniot curves predicted from the new approach are found comparable with existing experimental data. We find that the self-consistent nature of the method results in the onset of chemistry at 40% lower pressure on the shock Hugoniot than observed with classical molecular dynamics. The temperature shift associated with quantum heat capacity is determined to be the primary factor in this shift.

■

INTRODUCTION Shock compression techniques have traditionally been the primary approach for studying materials at pressures higher than those that can be achieved under static conditions. Shock waves also arise in fundamental materials phenomena like the detonation of an energetic material. Shock-induced chemistry results from a complicated interplay of mechanical and thermochemical processes. In recent years, atomistic simulations of shock compression have arisen as a promising tool for providing microscopic insight into the detailed chemistry and kinetics in these processes (see, e.g., ref 1 for a review). A challenge has been balancing the need to accurately describe chemical reactions with the nanosecond and longer time scales over which shockinduced chemistry can occur. Properly describing chemistry requires an accurate potential energy surface, which is usually computationally expensive. The multiscale shock technique (MSST)2−10 provides a partial solution to this problem by employing approximate hydrodynamics to enable the use of computationally expensive potential energy surfaces like those provided by density functional theory.11 MSST has been shown to approximately reproduce the kinetics of shocked materials5 and is among a number of recent methods aimed at improving the tractability of atomistic shock simulation.12−14 Much of the theoretical research on shock compression utilizes either analytical interatomic energy models or DFT. Great efforts have been made to improve the fit of interatomic potentials to DFT and other more accurate results (see, e.g., refs 15−17). While this is very important work, it is important to note that no matter how good the fit is, classical molecular dynamics lacks the quantum vibrational effects that can play an important role in shock compression. Quantum nuclear effects play several roles in shock compression. First, they provide a © 2012 American Chemical Society

temperature- and pressure-dependent heat capacity that impacts the temperature and pressure of the shocked material. Second, they provide a frequency-dependent energy distribution of the vibrational modes, changing the relative free energies of potential decomposition products. Other potential effects include quantum tunneling of nuclei. Classical MD simulations are only approximately correct for temperatures that are higher than the Debye temperature TD (TD = ℏΩD/kB, where ΩD is the Debye frequency), which is often at least on the order of a few hundred Kelvins. When T < TD, quantum effects cannot be neglected. For instance, in a water molecule at room temperature, the zero-point energy (ZPE) 1/2ℏω in each of the two O−H stretching modes and the single bending mode is more than 20 times kBTroom.18 Figure 1 schematically illustrates the classical and quantum energy vibrational spectrum at two different temperatures, showing the energy difference at typical organic molecule vibrational frequencies. Quantum nuclear effects in shocked matter have been studied by Goldman, Reed, and Fried by performing classical MSST calculations using a postcalculation temperature correction.19 They found that quantum temperature corrections to DFT calculations of the water and methane Hugoniots are as much as 20−30% of the temperature. The quantum corrections were found to bring the DFT-calculated temperatures closer to experimental measurements. In this work, we seek to employ a method that incorporates quantum nuclear effects in a self-consistent fashion, i.e., going beyond a postcalculation temperature correction. This is Received: August 14, 2012 Revised: September 24, 2012 Published: September 26, 2012 10451

dx.doi.org/10.1021/jp308068c | J. Phys. Chem. A 2012, 116, 10451−10459

The Journal of Physical Chemistry A

Article

linearly scales with atom velocity. We define the power spectral density (PSD) of the random force ∞

I R iκR jζ(ω) ≡

∫−∞

R iκ(t )R jζ(t + τ ) eiωτ dτ

(2)

where i and j are atom indices, κ and ζ are the labels for Cartesian directions. IRiκRjζ(ω) can be related to γ by the fluctuation−dissipation theorem22,24 I R iκR jζ(ω) = 2miγδijδκζ Θ(ω , T )

(3)

where ⎞ ⎛1 1 Θ(ω , T ) = ℏ|ω|⎜ + ⎟ exp(ℏ|ω| /kBT ) − 1 ⎠ ⎝2 Figure 1. Average energies of a harmonic oscillator shown at two different temperatures, T1 = 200 K and T2 = 1000 K. The black and red correspond to the classical and quantum mechanical expectation values, respectively. According to the equipartition theorem in classical statistics, the average energy of a harmonic oscillator in thermal equilibrium always equals kBT. When the system is at low temperature compared with TD, the quantum energy differs from the classical energy. The vibrational frequencies of some typical chemical bonds are marked approximately with dashed lines.

Equation 3 can be derived by considering the special case of a harmonic oscillator. We consider a 1-D harmonic oscillator along x with characteristic frequency ω0. The equation of motion, derived from eq 1, is mx(̈ t ) = −mω02x(t ) + R(t ) − mγx(̇ t )

(5)

Fourier transforming and solving for position and velocity give

motivated by a desire to enable description of the shifts in equilibrium chemical composition that quantum effects can potentially bring. We also desire an approach that is ideally no more complicated or expensive to use than classical MD. Among the most accurate self-consistent approaches to quantum nuclear effects is the path integral method. The required computational effort can be prohibitive, especially for low-temperature simulations due to the large number of required replicas to converge to the quantum limit.20 Another approach is to use a semiclassical theory, which takes numerically computed classical trajectories as the input and then postprocesses the data to recover quantum effects.18,21 For example, the classical MD temperature can be rescaled to an effective quantum temperature. However, this method lacks the self-consistency we seek. In this work, we adopt a so-called colored thermostat for quantum nuclear effects. Dammak et al.22 and Ceriotti et al.23 describe Langevin-type equations of motion that force the vibrational energies in the system to obey a Bose−Einstein spectrum. Instead of the white-noise spectrum random force employed in the usual Langevin thermostat, the random force obeys a “colored” spectral distribution. In this paper, we develop a new version of MSST that is combined with a quantum thermal bath (QTB) scheme similar to ref 22. This new approach enables the self-consistent inclusion of some approximate quantum nuclear effects “on the fly” for shock compression.

x(̃ ω) =

v (̃ ω) =

m(ω02

1 R̃(ω) − ω 2) + imγω

(6)

m(ω02

iω R̃(ω) − ω 2) + imγω

(7)

R̃ (ω) is the Fourier transform of time-domain random force R(t). The total energy is given by

∫−∞ d2ωπ ⎛⎝ 12 mω02|x(̃ ω)|2 + 12 m|v(̃ ω)|2 ⎞⎠ ∞

E=

⎜

⎟

ω2 + ω 2

∞

=

∫−∞ d2ωπ 2m(ω 2 − ω2)2 +0 2mγ 2ω2 |R̃(ω)|2 0

(8)

In the limit of γ ≪ ω0, the above integral done on a contour yields E=

|R̃(ω0)|2 2mγ

(9)

which is equivalent to eq 3, where |R̃ (ω)|2 = IRR satisfying the Wiener−Khinchin theorem and E(ω0) is equivalent to Θ(ω) in eq 3. In a common application of the Langevin equation, the random force is taken to be a white noise, which implies that Θ(ω, T) equals kBT. This is the classical limit where all harmonic modes have the same energy, regardless of frequency. In the QTB thermostat, Θ(ω, T) is chosen to be the Bose− Einstein spectrum for a quantum harmonic oscillator. The resulting average vibrational energy includes zero-point energy and a spectral dependence. By applying this new technique, Dammak et al. found that the thermal expansion and heat capacity of solid MgO can be successfully predicted at low temperatures, and it describes the liquid behavior of 4He above the λ point.25 The so-called quasiharmonic approximation employed here, where each mode is assumed to be a quantum harmonic oscillator, has been found to provide reasonable agreement with some experimental

■

METHODOLOGY Quantum Thermal Bath. In the quantum thermal bath scheme proposed by Dammak et al.,22 the equation of motion follows a Langevin form mi ri ⃗ ̈ = fi ⃗ + R⃗ i − miγ ri ⃗ ̇

(4)

(1)

where mi, ri⃗ , fi⃗ , and R⃗ i are the mass, position, force exerted by all other atoms, and random force of atom i. The parameter γ is the frictional coefficient, and the dissipative force (−miγr ⃗ ̇i) 10452



Article

transformation of H̃ (ω) gives H(t). Then we can calculate the random signal θ(t) and random force Riκ(t) by eqs 10 and 15. In practice, discrete Fourier transforms and convolutions are required instead of the continuous versions. The filter H̃ (ω) is an even function discretized on a grid with 2Nf points with δω spacing over the frequency range [−Ωmax, Ωmax]

results for other materials, such as liquid water,26 and approximate results for anharmonic potentials.25 The parameter γ controls the coupling strength of the system’s phonon modes to the bath. A suitable choice of γ is critical for the QTB thermostat, as is the case for any Langevin thermostat. There are two criteria to be considered. First, when γ is too big, eq 3 is no longer exact. The resulting spectrum deviates from the prescribed Θ(ω), potentially leading to inaccurate results. Large values of γ also impact the dynamics in eq 1. Second, when γ is too small, eq 1 reduces to the classical limit, where there is no coupling to the bath. Anharmonic coupling between phonon modes provides a finite lifetime, characterized by an effective damping Γanh. If the coupling strength γ is much stronger than the anharmonic coupling Γanh, then QTB thermostat should be able to impose the desired energy spectrum. Anharmonicity sets a limit on the minimum value of γ. To generate the random forces, we employ the method of Barrat et al.,27 which circumvents the need to generate all random forces for all times before the simulation.22 This algorithm employs a method for noise synthesis utilized by the signal-processing community. The memory requirement of this approach is less demanding and independent of the simulation duration. The goal is to calculate time-dependent random forces Riκ(t). Here, we define R (t ) θ(t ) ≡ iκ 2miγ

H̃k = H̃ −k = H̃ (kδω),

Hn = H(nδh) =

Θ(ω)

R iκ(nδh) =

(11)

(12)

(13)

(14)

However, eq 14 requires the inverse Fourier transform be performed at the start of the calculation and θ(t) stored in memory for the duration of the simulation. Instead, Barrat et al. take ∞

θ (t ) =

∫−∞ H(τ)r(t − τ)dτ

H̃k cos(πkn/Nf ) (17)

2miγ

∑

Hmrn − m (18)

where rn−m is a white noise array drawn from a uniform distribution with a variance 1/h. The random number array rn−m is first initialized by drawing 2Nf initial values. At each subsequent time step n, with time step size δh, rn−Nf+1 is discarded and a new random number rn+Nf is generated and stored. As a result, the storage size is therefore independent of the total number of time steps. In principle, the molecular dynamics integration time step δt does not have to equal the noise generation time step δh. In fact, we take a δh that is an integer multiple larger than δt, δh = αδt. The same random force can be exerted α times during the MD integration. The integration step δt primarily depends on the Debye frequency ΩD and is typically on the order of 0.01 × 2π/ΩD. For a higher Debye frequency, a finer δt is necessary to capture the fast motion of atoms and molecules. However, the choice of δh is determined by the Nyquist−Shannon sampling theorem, δh ≤ 0.5 × 2π/Ωmax. The time step h is inversely proportional to Ωmax, which is the sampling upper bound of spectrum Θ(ω). MSST with Quantum Thermal Bath (QB-MSST). The direct method of simulating a shock wave propagating through a material involves propagating the shock down a long tube of atoms. Achieving shock propagation time scales >10 ps requires use of a sufficiently large number of atoms to preclude the use of the more expensive methods like density functional theory.28,29 Instead of simulating a shock wave within a large computational cell with many atoms, the computational cell of MSST follows a Lagrangian material element through a shock wave at a specified shock speed, enabling simulation of the shock wave with significantly fewer atoms and lower computational cost. MSST simulates steady shock waves (i.e., constant shock speed) by time-evolving equations of motion for the atoms and volume of the computational cell to constrain the shock propagation direction stress to the Rayleigh line and energy to the Hugoniot energy conditions.2−4 For a specified shock speed, the latter two relations describe a steady planar

(10)

∞

∫−∞ d2ωπ H̃ (ω)r(̃ ω)e−iωt

∑ k =−Nf

m =−Nf

Taking the inverse fast Fourier transform of eq 11, θ(t) is computed as θ (t ) =

Nf − 1

Nf − 1

where the white noise property of r̃(ω) makes |r̃(ω)|2 a constant, which we take to be unity. H̃ (ω) can be chosen H(̃ ω) =

1 2Nf

The random force for atom i along Cartesian direction κ at the nth time step can be computed by performing a discrete convolution in the time range [−π/δω, π/δω]

where the frequency-dependent filter H̃ (ω) is to be determined and r̃(ω) is a white noise which is readily generated using standard random number generators. From eq 11, we get Θ(ω) = |θ (̃ ω)|2 = |H̃ (ω)|2 |r (̃ ω)|2

(16)

Ωmax is chosen such that Ωmax > ΩD. More details concerning the choice of Ωmax will be discussed in subsequent sections. The frequency sampling step size is δω = 2Ωmax/2Nf. The noise filter in the time domain H(t) can then be computed via discrete inverse Fourier transform of H̃ (ω), over the time range [−π/δω, π/δω] with a time step size δh = π/Ωmax. At the nth time step

where the power spectral density of θ(t) is Θ(ω). The Fourier transform of θ(t) is denoted as θ̃(ω). We want to express θ̃(ω) as the following product θ (̃ ω) = H̃ (ω)r (̃ ω)

k = −Nf ...Nf − 1

(15)

by using the convolution theorem. This approach is faster due to the relatively small spectral content of H̃ (ω) compared with r̃(ω) for practical simulations. The spectral content of r̃(ω) increases with simulation duration while H̃ (ω) is independent of duration. To summarize, in order to get the random force in the time domain, H̃ (ω) is obtained first according to eq 13. A Fourier 10453



Article

Figure 2. Flowchart for QB-MSST algorithm implemented in LAMMPS.

and v is the computational cell volume. A scaled coordinate scheme is utilized, where ri ⃗ ̇ ≡ Asi⃗ ̇. A contains the computational cell lattice vectors in columns, and si⃗ is the scaled position of atom. Uniaxial strain of the computational cell is used since the shock is planar, and dA/dv is a matrix corresponding to the direction to the uniaxial strain; i.e., A is only allowed to change with one degree of freedom. Equations 20 and 21 are the equations of motion derived from the Hamiltonian based on the conserved quantity ẽ (energy per unit mass)

shock wave within continuum theory. Uniaxial strain appropriate for a planar shock wave is imposed on the computational cell. The MSST has been demonstrated to reproduce accurately the sequence of thermodynamic states throughout the reaction zone of explosives with analytical equations of state, shock waves in amorphous Lennard-Jones, and amorphous Tersoff carbon.3,5 It has been used to simulate shocks in water, 6 nitromethane, 7 and other materials represented by DFT and tight-binding quantum energy models.8−10,30 MSST was derived from Navier−Stokes equations. Mass, momentum, and energy are conserved everywhere in the wave. Positions and velocities are updated at each time step to mimic a compressive shock wave passing over the system. It varies the cell volume and temperature in such a way as to restrain the system’s thermodynamic states to the shock Hugoniot and the Rayleigh line. These restraints correspond to the macroscopic conservation laws dictated by a shock.2,5 The MD equations of motion for volume of the computational cell and the atoms are Qv∼ ̈ =

∑ miAsi⃗̇· i

ẽ =

i

To couple the quantum thermal bath to MSST, the goal is to let the QTB thermostat regulate the velocities so that the velocity spectrum obeys the one for a quantum harmonic oscillator. The atomic equation of motion for MSST is updated with the random and dissipative force terms in eq 1, −1/ miA−1R⃗ i and −γs ⃗ ̇i

v2 dA ̇ dϕ v∼ ̇ si⃗ − − s2 (v0̃ − v )̃ − p0 − μ dv dv ṽ v0̃

si⃗ ̈ =

∼2

si⃗ ̇ μv ̇ −1 −1 ∂ϕ A − G−1Gṡ i⃗ ̇ + M mi miv ̃ ∑j |Asj⃗|̇ 2 ∂r ⃗ i

2 vs2 ⎛ 1 1 ṽ ⎞ 1 2 mi|Asi⃗|̇ 2 + ϕ({ ri }) ⃗ − p0 (v0̃ − v )̃ − ⎜1 − ⎟ + Qv∼ ̇ M 2M 2⎝ v0̃ ⎠ 2

(21)

(19)

si⃗ ̈ =

∑

2 si⃗ ̇ μv∼ ̇ − 1 −1 ⃗ − 1 −1 ∂ϕ + − G−1Gṡ i⃗ ̇ + M A A R i − γ si⃗ ̇ ∂ ri ⃗ mi miv ̃ ∑j |Asj⃗|̇ 2 mi

(22)

One challenge that must be addressed here is the choice of temperature in the QTB. The quantum thermal bath described in the previous section is a constant-temperature thermostat. However, the temperature varies, sometimes drastically, upon shock compression. We introduce a time-dependent temperature TQM(t) which is the target temperature for the bath. Θ(ω, TQM(t)) generated from the bath regulates atom

(20)

The shock wave moves at a speed vs. ṽ = 1/ρ is the specific volume. Q is a masslike parameter for the volume, and μ is an empirical macroscopic viscosity term. M is the total system mass. G is defined as ATA, ϕ({ri⃗ }) is the total potential energy, 10454



Article

Figure 3. Convergence tests with respect to Nf, γ, Ωmax, and η. Unless otherwise depicted, these plots take Nf = 1000, γ = 5 THz, and Ωmax/2π = 333 THz. In QTB, the classical temperature TC varies as a function of quantum bath temperature TQM and (a) frequency grid Nf, (b) frictional coefficient γ, and (c) upper frequency noise limit Ωmax. (d) In QB-MSST, the temperature trajectory upon shock (6 km/s) is set by the value of η. (e) The computational cost of QB-MSST is most strongly influenced by Nf. When Nf = 100, the cost of QB-MSST is only 6% more than classical MSST with no quantum mechanical nuclear effects. The parameters α = 6 and β = 600. When δt is 0.25 fs, β = 600 gives a 150 fs time step for integration of eq 23.

(23)

of TQM(t). Generally, the faster the volume change during shock compression, the more frequently TQM(t) needs to be adjusted. Updating the bath temperature around 10−15 times during the temperature ramp is found to enable TQM(t) to track thermodynamic changes.

where M(ẽ − ẽ0) defines the constraint energy quantity change. The term 3NkB in the denominator is a rough measure of heat capacity, converting the energy deviation to temperature. The parameter γ defines the typical time scale (1/γ) to reach thermal equilibrium at the temperature TQM(t). The rate for the bath to reach equilibrium and the maximum rate of bath temperature change should be comparable. The factor η is a new parameter introduced in this equation; it defines how fast the target quantum mechanical bath temperature responds to energy changes of the system. We find that it is unnecessary to generate a new TQM(t) at every integration step δt. TQM(t) is updated every β time steps (typically 50−150 fs), because the temperature is ill-defined on time scales shorter than the vibrational time scale determined by 2π/ΩD. Therefore, we take an averaged instead of an instantaneous ẽ (shown in Figure 2), as we propagate TQM(t). Taking a time-averaged ẽ reduces the magnitude of fluctuations

RESULTS AND DISCUSSION We perform our MD simulations using the LAMMPS code,31 in which MSST is already implemented. We modified the code to incorporate the QB-MSST scheme. We chose to study liquid methane because of its high-energy phonon modes and significant quantum effects. Also, methane is the most abundant organic molecule in the universe, and it is one of the major constituents of giant planets Uranus and Neptune. Hence, studying methane subjected to high pressure and high temperature is of great importance in planetary physics. Methane’s four harmonic normal modes are computed to be 41.0, 47.4, 90.7, and 94.6 THz by using hybrid density functional theory.32 We use a version of the ReaxFF potential,15 and the computational cell 23.87 × 23.87 × 23.87 Å3 consists of 216 CH4 molecules with 0.423 g/cm3 initial density and 111 K initial temperature, the same as gas-gun experimental

velocities of the system. Here, we utilize the equation of motion for updating bath temperature Ṫ

QM

(t ) = γη

M(e (̃ t ) − e0̃ ) 3NkB

■

10455



Article

Figure 4. (a) Comparison of methane heat capacities (ρ = 0.423 g/cm3) between QTB-MD (red dots) and classical MD (black solid line) to experimental data (blue solid line).36 (b) Temperature and the conserved energy quantity (in units of temperature, ΔT̃ = M(ẽ(t) − ẽ0)/3NkB) for methane shocked with 6 km/s shock speed.

the low-temperature range. A value of γ ≤ 5 THz gives reasonably converged results in this case. Because the highest vibrational frequency of methane is around 95 THz, we tested four different random force sampling upper bounds Ωmax/2π that are higher than 95 THz. In Figure 3c, we vary Ωmax with fixed Nf and γ. Ωmax is found sufficient at and above 250 THz, corresponding to about 2.5ΩD. In our case, a value of ΩD/2π = 250 THz converges the results. One word of caution is to be aware of the fact that the chemical bonds can get stiffened upon compression, particularly in solids. Figure 3 shows a variety of Ωmax values at ambient conditions. However, Ωmax appropriate for the high-pressure and hightemperature conditions of QB-MSST might be a larger value. In our case, we find Ωmax/2π = 333 THz good for converging methane QB-MSST simulations at shock pressures. Finally, we run our new approach QB-MSST to test η by performing a shock at 6 km/s. In Figure 3d, small η values cannot adjust the quantum mechanical target temperature fast enough during the fast temperature-ramping period; while large η leads to big temperature oscillations. We choose η = 0.33 for the calculations here. This value of η is likely to also be appropriate for simulating other systems under shock compression. We observe that different values of η lead to almost the same final thermodynamic state behind the shock, as expected. Furthermore, we note that the eight simulation parameters associated with the QB-MSST method that were determined in this work are expected to be appropriate for similar hydrocarbon systems involving the same chemical bonds. For systems of different chemical makeup, convergence tests are required for the three parameters Q, μ, and γ, using the following set of values for remaining five, relatively systemindependent, parameters as an initial guess: Nf ≈ 1000, Ωmax ≈ 3.5ΩD, η ≈ 1/3, α ≤ π/Ωmax/δt, and β ≈ 2π/ΩD·10/δt. Application to Methane. The temperature of a material determines the number of phonons occupying each phonon mode. The phonon occupancy has important implications for physical properties including specific heat and thermal conductivity. Classically, the heat capacity per atom for a 3D harmonic crystal is equal to 3kB, independent of the temperature. The heat capacity of a quantum harmonic solid converges to the classical limit 3kB only when temperature T ≫ TD and goes to zero as temperature approaches zero. Using the parameters determined from the convergence tests, we compute the heat capacity by differentiation of the mean total energy from simulations equilibrated with the quantum thermal bath. Figure 4a shows good agreement between our results and experimental measurements for methane at a

conditions.33 The integration time step we use varies from 0.25 to 0.12 fs, as the shock speed increases from 6 to 12.2 km/s. The total simulation time is approximately 200−300 ps. This time might not be sufficient for the system to reach chemical equilibrium. Some of the reactions, such as polymerization,34,35 may evolve on a much longer time scale. Convergence Tests. In this section we discuss the choice of simulation parameters. Convergence tests are necessary and therefore highly recommended before the actual QB-MSST simulations. The parameters that we will establish default values for are Q, μ, Nf, Ωmax, γ, η, α, and β. Q and μ can be determined by running MSST without the quantum bath. Nf, Ωmax, and γ can be determined by running QTB without MSST. η is determined by running QB-MSST directly. First, we run standard MSST simulations to determine Q and μ. While setting μ = 0, a value of Q is chosen to obtain an initial volume compression time of approximately 1 ps, characteristic of shocked liquids. Then if necessary, a nonzero μ can be set to damp any observed excessive volume oscillation at the start of the simulation. In our case, Q = 6.9 × 10−13 kg2/m4 and μ = 1.5 × 10−2 kg/m/s. Both Q and μ are transferable from MSST to QB-MSST. We note that Q and μ impact the initial (unreacted) compression but have little impact on subsequent slower changes associated with chemistry.5 Second, we run quantum thermal bath calculations at a series of temperatures distributed between 0 and 1500 K to determine the number of frequency grid points Nf, frictional coefficient γ, and upper cutoff frequency Ωmax. The frequency grid size Nf is the number of points from 0 to Ωmax over which the Bose−Einstein spectrum Θ(ω) is sampled. A large Nf provides a more accurate sampling of the spectrum. In Figure 3a, with fixed Ωmax and γ, we vary Nf and plot the classical temperature as a function of target quantum bath temperature. The instantaneous classical temperature TC is derived from the atom velocities at thermal equilibrium N

C

T =

∑i = 1 mivi2 3NkB

(24)

Figure 3a shows nonzero TC as the quantum temperature T approaches zero. This is associated with the zero-point energy at low temperature. A value of Nf = 100 converges the results well. In this paper, we choose to use Nf = 1000. Then, we run QTB to test γ with fixed Ωmax and Nf. As we discussed earlier, a γ chosen too large leads to an energy spectrum that differs from the desired Θ(ω). Figure 3b shows that large γ values lead to overestimation of TC, especially for QM

10456



Article

Figure 5. Comparison of our QB-MSST and MSST ReaxFF (a) temperature vs density Hugoniot and (b) pressure vs density Hugoniot with other theories19,37 and experiments.33,38 (b) Temperature, density, and pressure are taken to be the averaged values for the last 3 ps at the end of each simulation. Quantum corrections (QB-MSST) to ReaxFF potential systematically shift T−ρ and p−ρ Hugoniot curves to higher-temperature and higher-pressure phase region. Moreover, QB-MSST is found to improve agreement with the experimental measurements and quantum-corrected DFT results of refs 19 and 37 for most of the thermodynamic range considered. In (b), as pressure increases, both MSST and QB-MSST results start to deviate from the experimental measurement which suggests that the p−ρ quality of the ReaxFF potential could be further refined in the highpressure regime.

density of 0.423 g/cm3. It is known that the heat capacity plays a crucial role in shock compression. During dynamical compression, mechanical energy is converted into heat, yielding a temperature increase. The heat capacity impacts the magnitude of this temperature change. Figure 4b shows shock-induced temperature change of MSST and QB-MSST, for a 6 km/s shock speed. At this relatively low shock speed, no appreciable reactivity is observed in our simulations. The elevated final temperature for QBMSST, compared with MSST, results from the smaller specific heat of the quantum bath. Also, the conserved energy ẽ(t) (converted to an approximate temperature) from QB-MSST fluctuates around zero, with largest variation during the temperature ramp. Figures 5a and 5b show simulations with shock wave speeds from 6 to 12.2 km/s. Depicted are the thermodynamic states of the system after a simulation duration of 200−300 ps. We observe substantial temperature and pressure difference between MSST and QB-MSST. While our aim here is to study the role of quantum effects in the context of a particular energy model (ReaxFF) rather than benchmark the model against other available data, Figures 5a and 5b also include the principal Hugoniot curves from two other density functional theory-based theoretical works19,37 and gas-gun experiments.33,38 The uncorrected T−ρ Hugoniot (Figure 5a, data labeled by MSST) is consistent with the uncorrected DFT calculations from ref 19. The corrected T−ρ Hugoniot (data labeled by QB-MSST) exhibits reasonable agreement with other correction theories and available experimental measurements. A similar trend is observed in the p−ρ Hugoniot (Figure 5b). The p−ρ Hugoniot deviation between the ReaxFF QBMSST results and the experimental measurements increases as pressure increases. The volume equation of motion (eq 19) employed here utilizes a classical form of the thermal stress in the computational cell, given in part by the first term on the right-hand side of eq 19. The quantum nature of nuclear vibrations should impose some correction to this stress which we neglect. To obtain a rough estimate of the potential significance of this quantum correction, we consider the thermal term of the virial stress equation NkB(TQM − TC)/V where TQM and TC are the final quantum and classical temperatures of QB-MSST calculations. The value of this stress is approximately 0.7 GPa for the simulations performed here.

Chemical reactions are observed at sufficiently large shock wave speeds. The self-consistent nature of QB-MSST enables quantum corrections for reacting mixtures. Although coarsely determined, we find different reaction thresholds for QB-MSST (1.05% CH4 decomposed, shock speed 9 km/s, thermodynamic state 2120 K, 0.986 g/cm3, and 19.7 GPa) and MSST (2.03% CH4 decomposed, shock speed 11 km/s, thermodynamic state 2630 K, 1.139 g/cm3, and 32.2 GPa), where around 1−2% methane is decomposed. We note that at the end of the whole simulation duration (200−300 ps), the chemical compositions appear not to have reached steady values for all shocks where chemistry is observed. The major reaction products are H, H2, and a variety of hydrocarbons. CH4 → H + CH3 is the first reaction to take place in both the MSST and QB-MSST calculations, with and without quantum corrections. As simulation time evolves, H2, C2H6, and some other hydrocarbon species start to form. The variety of the hydrocarbon chains increases with shock speed (Figures 6a and 6b). We analyze trajectory frames that are 150 fs apart and apply the following bond distance criteria to determine the breaking and forming of bonds: C−H ≤ 1.57 Å, H−H ≤ 1.09 Å, and C−C ≤ 1.98 Å. We take the bond length cutoff to be given by the first minimum in the corresponding pair radial distribution function. Also, we performed shock simulations with a different γ value (1 THz), and no obvious difference in chemistry or final thermodynamic state was found. Figures 6a and 6b show different reaction onset points on the Hugoniot curves of QB-MSST and MSST. The higher temperatures achieved in the quantum case undoubtedly play a role in this difference. As shown in Figure 6a, even at the same temperature along the Hugoniot, methane shows higher reactivity in QB-MSST. Besides the temperature factor, does the nonuniform energy spectrum from quantum corrections also contribute? To explore this possibility, we consider three possible methane decomposition reactions: (1) CH4 → 1/2C2H6 + 1/2H2, where C2H6 represents a shorter-length hydrocarbon, (2) CH4 → 1/4C4H10 + 3/4H2, where C4H10 represents a longer-length hydrocarbon, and (3) CH4 → C (diamond) + 2H2. The three reaction zero-point energy changes (difference between ZPE of products and ZPE of reactant) are approximately −0.06, −0.11, and −0.40 eV/CH4 respectively. (All the zero-point energy data are taken from NIST Chemistry Webbook, except for the diamond zero-point energy datum (0.232 eV/atom) which was calculated at 200 10457



Article

Hugoniot is therefore primarily a result of the higher classical temperature in QB-MSST. Two separate laser-heated diamond anvil cell (DAC) experiments indicate the formation of diamond, molecular hydrogen, and hydrocarbon polymers at temperatures above 2500 K and pressures as low as 13 GPa40 and at temperatures of about 2000 and 3000 K and pressures between 10 and 50 GPa.35 Reactions forming saturated hydrocarbons containing 2−4 carbons, molecular hydrogen, and graphite are observed in another laser-heated DAC experiment at pressures above 2 GPa and temperature between 1000 and 1500 K.41 Gas-gun shock experiments on methane are interpreted to indicate that chemistry occurs above 23 GPa on the principal Hugoniot,33 but determination of the chemistry onset via the Hugoniot alone is challenging.42 Later gas-gun experiments combined with the electrical conductivity measurements of methane can be interpreted to suggest that the liquid is mostly unreacted at pressures around 30 GPa and below.42 It may be important to note that the gas-gun experiments achieve high pressures for shorter time periods (microseconds) than the DAC experiments (seconds). The different time scale between these two methods may play a role in the different reaction onset suggested by the two experimental methods. Several separate DFT-based MD studies predicted formation of hydrocarbon mixtures at 4000 K and 100 GPa34 and above 3420 K in pressure and 37 GPa in pressure.37 The hydrogen atom is predicted to start forming at 3243 K and 32 GPa.19 Also, a tight-binding MD study suggests a 43 GPa reaction threshold.43 The DFT-based calculations are limited to picosecond time scales. Picosecond time scale AIREBO empirical potential studies by Elert et al. showed carbon chemistry above 80 GPa and 6000 K.44 Longer time scale AIREBO potential studies by Bolesta et al. lasted 200 ps and found the methane decomposition threshold to be at a shock temperature of 4400 K.45 Bolesta et al. attributed the reason for the reaction onset difference between ref 44 and their work to the difference in time scale and compression methods. Neither of the AIREBO studies report evidence of diamond formation. In QB-MSST calculations with ReaxFF (red in Figures 6a and 6b), as the shock speed reaches 12.2 km/s, the highest value of m in the hydrocarbon CmHn increases with time to 11 at 240 ps behind the shock front. Similar to the AIREBO simulations, no diamond fragments are observed in the QBMSST calculations. Polymer hydrocarbon formation is suspected to be an intermediate step toward diamond formation, and a longer simulation time may lead to diamond formation eventually. The observed CH4 → H + CH3 threshold is approximately 2120 K in temperature and 20 GPa in pressure, and carbon chemistry is observed as low as 2540 K and 25 GPa. These thresholds appear to be in better agreement with experiments than other theories. The QB-MSST observed product types and relative abundances are similar to those shown in the AIREBO simulation mass spectra of ref 44.

Figure 6. (a) Numbers of remaining CH4 molecules and the primary products CH3, H, H2, and C2H6 per supercell. (b) Simulated mass spectra of product hydrocarbons CmHn at shock speed 9, 10, 11, and 12.2 km/s. These four shock speeds correspond to the highest T vs ρ and P vs ρ points in Figures 5a and 5b. All the chemical compositions are the averaged values for the last 3 ps at the end of each 200−300 ps simulation.

GPa from ref 39. We expect the zero-point energy for diamond at ambient pressure to be smaller than 0.232 eV/atom, making the zero-point energy change of reaction 3 even more negative.) For each reaction, the consideration of zero-point energy not only increases the reactivity of methane but also biases the products more toward diamond. This suggests that the nonuniform vibrational spectrum could play a role in the chemistry observed. To systematically study the roles of the quantum temperature shift and nonuniform spectrum, we run QM and classical constant-temperature simulations at the same density (ρ = 0.423 g/cm3) and the same classical temperature for 1 ns. In the QM simulation, the system is thermalized in a quantum bath, and TQM is set to 2300 K; the final classical temperature is 2556 K. The second simulation is classical MD at 2556 K using a Nosé−Hoover thermostat. For both runs, we find similar reaction products and similar kinetics. The quantum-corrected energy spectrum does not appear to play a dominant role under these thermodynamic conditions, but could under different conditions. The different onset of chemistry observed along the

■

CONCLUSIONS In conclusion, this work presents the new QB-MSST method for molecular dynamics simulations of shocked materials including quantum nuclear effects in a semiclassical fashion. This method has computational expense similar to classical MD yet provides quasi-harmonic quantum nuclear corrections in a self-consistent way. We have applied this method to shocked methane described by a ReaxFF energy model. The quantum corrections are found to shift the Hugoniot curves to higher10458



Article

(19) Goldman, N.; Reed, E. J.; Fried, L. E. J. Chem. Phys. 2009, 131, 204103/1−204103/8. (20) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (21) Wang, C. Z.; Chan, C. T.; Ho, K. M. Phys. Rev. B 1990, 42, 11276−11283. (22) Dammak, H.; Chalopin, Y.; Laroche, M.; Hayoun, M.; Greffet, J.-J. Phys. Rev. Lett. 2009, 103, 190601/1−190601/4. (23) Ceriotti, M.; Giovanni, B.; Parrinello, M. Phys. Rev. Lett. 2009, 103, 030603/1−030603/4. (24) Callen, H. B.; Welton, T. A. Phys. Rev. 1951, 83, 34−40. (25) Dammak, H.; Chalopin, Y.; Laroche, M.; Hayoun, M.; Greffet, J.-J. Phys. Rev. Lett. 2011, 107, 198902/1−198902/2. (26) Beren, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. J. Chem. Phys. 1983, 79, 2375−2389. (27) Barrat, J.-L.; Rodney, D. J. Stat. Phys. 2011, 144, 679−689. (28) Gygi, F.; Galli, G. Phys. Rev. B 2002, 65, 220102/1−220102/4. (29) Mattson, W. D.; Balu, R. Phys. Rev. B 2011, 83, 174105/1− 174105/7. (30) Reed, E. J.; Rodriguez, A. W.; Manaa, M. R.; Fried, L. E.; Tarver, C. M. Phys. Rev. Lett. 2012, 109, 038301/1−038301/5. (31) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (32) Gray, D. L.; Robiette, A. G. Mol. Phys. 1979, 37, 1901−1920. (33) Nellis, W. J.; Ree, F. H.; van Thiel, M.; Mitchell, A. C. J. Chem. Phys. 1981, 75, 3055−3063. (34) Ancilotto, F.; Chiarotti, G. L.; Scandolo, S.; Tosatti, E. Science 1997, 275, 1288−1290. (35) Benedetti, L. R.; Nguyen, J. H.; Caldwell, W. A.; Liu, H.; Kruger, M.; Jeanloz, R. Science 1999, 286, 100−102. (36) Gurvich, L. V.; Veyts, I. V.; Alcock, C. B. Thermodynamic Properties of Individual Substances, 4th ed.; Hemisphere Publishing Corp.: New York, 1989; Vols. 1 and 2. (37) Li, D.; Zhang, P.; Yan, J. Phys. Rev. B 2011, 84, 184204/1− 184204/7. (38) Radousky, H. B.; Mitchell, A. C.; Nellis, W. J. J. Chem. Phys. 1990, 93, 8235−8239. (39) Gao, G.; Oganov, A. R.; Ma, Y.; Wang, H.; Li, P.; Li, Y.; Iitaka, T.; Zou, G. J. Chem. Phys. 2010, 133, 144508/1−144508/5. (40) Zerr, A.; Serghiou, G.; Boehler, R.; Ross, M. High Pressure Res. 2006, 23−32. (41) Kolesnikov, A.; Kutcherov, V. G.; Goncharov, A. F. Nat. Geosci. 2009, 2, 566−570. (42) Nellis, W. J.; Hamilton, D. C.; Mitchell, A. C. J. Chem. Phys. 2001, 115, 1015−1019. (43) Kress, J. D.; Bickham, S. R.; Collins, L. A.; Holian, B. L. Phys. Rev. Lett. 1999, 83, 3896−3899. (44) Elert, M. L.; Zybin, S. V.; White, C. T. J. Chem. Phys. 2003, 118, 9795−9801. (45) Bolesta, A. V.; Zheng, L.; Thompson, D. L. Phys. Rev. B 2007, 76, 224108/1−244108/6.

temperature and higher-pressure phase regions, bringing the ReaxFF potential-based simulation results closer to the experimental measurements. We find that quantum effects are responsible for a significant shift in the onset of chemistry to lower pressures along the Hugoniot, and the additional classical temperature brought by the quantum heat capacity contributes the most to the lowered reaction onset predicted by QB-MSST.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work uses resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Some calculations were performed in part using the Stanford NNIN Computing Facility (SNCF), a member of the National Nanotechnology Infrastructure Network (NNIN), supported by the National Science Foundation (NSF). We also acknowledge the support from the NASA AIM laboratory.

■

REFERENCES

(1) Dlott, D. D. Annu. Rev. Phys. Chem. 2011, 62, 575−597. (2) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. Phys. Rev. Lett. 2003, 90, 235503/1−235503/4. (3) Reed, E. J.; Fried, L. E.; Henshaw, W. D.; Tarver, C. M. Phys. Rev. E 2006, 74, 056706/1−056706/9. (4) Reed, E. J.; Fried, L. E.; Manaa, M. R.; Joannopoulos, J. D. A Multi-Scale Approach to Molecular Dynamics Simulations of Shock Waves, in Chemistry at Extrmeme Conditions; Elsevier: New York, 2005; pp 297−326. (5) Reed, E. J.; Maiti, A.; Fried, L. E. Phys. Rev. E 2010, 81, 016607/ 1−016607/9. (6) Goldman, N.; Reed, E. J.; Kuo, I.-F.; Fried, L. E.; Mundy, C. J.; Curioni, A. J. Chem. Phys. 2009, 130, 124517/1−124517/6. (7) Reed, E. J.; Manaa, M. R.; Fried, L. E.; Glaesemann, K. R.; Joannopoulos, J. D. Nat. Phys. 2008, 4, 72−76. (8) Mundy, C. J.; Curioni, A.; Goldman, N.; Kuo, I. W.; Reed, E. J.; Fried, L. E.; Ianuzzi, M. J. Chem. Phys. 2008, 128, 184701/1−184701/ 6. (9) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. J. Am. Chem. Soc. 2009, 131, 5483−5487. (10) Goldman, N.; Reed, E. J.; Fried, L. E.; Kuo, I. F. W.; Maiti, A. Nat. Chem. 2010, 2, 949−954. (11) Reed, E. J. J. Phys. Chem. C 2012, 116, 2205−2211. (12) Zhakhovsky, V. V.; Nishihara, K.; Anisimov, S. I. JETP Lett. 1997, 66, 99−105. (13) Ravelo, R.; Holian, B. L.; Germann, T. C.; Lomdahl, P. S. Phys. Rev. B 2004, 70, 014103/1−014103/9. (14) Maillet, J.-B.; Mareschal, M.; Soulard, L.; Ravelo, R.; Lomdahl, P. S.; Germann, T. C.; Holian, B. L. Phys. Rev. E 2000, 63, 016121/1− 016121/8. (15) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. I. J. Phys. Chem. A 2008, 112, 1040−1053. (16) Goldman, N.; Fried, L. E. J. Phys. Chem. C 2012, 116, 2198− 2204. (17) Moriarty, J. A.; Benedict, L. X.; Glosli, J. N.; Hood, R. Q.; Orlikowski, D. A.; Patel, M. V.; Soderlind, P.; Streitz, F. H.; Tang, M. J.; Yang, L. H. J. Mater. Res. 2006, 21, 563−573. (18) Miller, W. H. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6660− 6664. 10459


Simulations of Shocked Methane Including Self-Consistent

Recommend Documents