Reactive Molecular Dynamics Simulations to Investigate the Shock

Subscriber access provided by United Arab Emirates University | Libraries Deanship

C: Physical Processes in Nanomaterials and Nanostructures

Reactive Molecular Dynamics Simulations to Investigate the Shock Response of Liquid Nitromethane Md Mahbubul Islam, and Alejandro Strachan J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b11324 • Publication Date (Web): 10 Jan 2019 Downloaded from http://pubs.acs.org on January 12, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Reactive Molecular Dynamics Simulations to Investigate the Shock Response of Liquid Nitromethane

Md Mahbubul Islam, Alejandro Strachan*

School of Materials Engineering and Birck Nanotechnology Center, Purdue University, West Lafayette, Indiana 47907, United States

*Corresponding author: [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

Abstract We use molecular dynamics (MD) simulations with the ReaxFF reactive force field to investigate the thermo-mechanical, chemical, and spectroscopic response of nitromethane (NM) to shock loading. We simulate shocks using the Hugoniostat technique and use four different parameterizations of ReaxFF to assess the sensitivity of the results with respect to the force field. The predicted shock states, both for the unreacted and reacted materials, are in good agreement with experiments and two of the force fields capture the increase in shock velocity due to exothermic reactions in excellent agreement with experiments. The predicted detonation velocities with these two force fields are also in good agreement with experiments and differences in predicted values is linked to the differences in the reaction products. Across all force fields, NM decomposes predominantly via bimolecular reactions and the formation of nitrosomethane (CH3NO) is found as a dominant initiation pathway. We elucidate the mechanisms of secondary reactions leading to stable products, whose predicted populations with all four descriptions are in good agreement with experiments. We also calculated the timeresolved spectra from the trajectories during the shock processes to help correlate the underlying reaction mechanisms with spectral features and enable a one-to-one comparison with laserdriven shock experiments. This study demonstrates the potential of reactive molecular dynamics to describe the physics and chemistry of high-energy density materials under shock loading and complement experimental efforts to derive a definite, validated understanding.


Page 3 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I. INTRODUCTION A detailed understanding of the mechanical and chemical response of high energy (HE) materials to shock loading is crucial to evaluate their performance and safety. When HE materials are subjected to strong mechanical shocks, they undergo complex chemical and physical processes that can eventually lead to the detonation. Once a steady detonation is established, a leading shock compresses and heats up the material to the von Neumann spike where chemistry is initiated. The chemical reactions that sustain the detonation occur within the reaction zone that ends in the Chapman-Jouguet (CJ) point. Throughout the reaction zone, with length scales ranging from 10 μm to 1,000 μm and time scales of nanoseconds to microseconds,1 the material heats up and expands. The equations of state (EOS) of the unreacted explosive and CJ reaction products are key to understand detonation properties, and their determination either experimentally or theoretically remains challenging. The evaluation of the detailed chemical reactions that takes initial HE material to final products has proven even harder to characterize. Recent experimental breakthroughs are beginning to shed light into the shock-induced process at micron and picosecond timescales, tracking detonation chemistry using picosecond resolution spectroscopy2 or mechanical response of the materials using ultrafast dynamic ellipsometry.3 Other notable recent development includes utilization of Raman spectroscopy to probe vibrational energy redistribution,4 velocity interferometer system for any reflector (VISAR) to measure particle velocities directly,5 or coherent anti-Stokes Raman scattering technique to obtain vibration spectra.6,7 Similarly, recent advancements in reactive atomistic simulations are elucidating the details of molecular-level processes responsible for shock-induced chemistry both with electronic structure methods8–13 and reactive interatomic potentials. The reactive force field ReaxFF14,15 enabled large-scale shock simulations of realistic HE materials to investigate the unreacted and reacted states accessible via shocks (the Hugoniot EOS), chemical decomposition mechanisms, and the role of microstructures and defects16,17 on hotspot formation in HE materials under shock loading. Recently, multimillion atoms simulations were able to probe the shock to deflagration transition upon collapse of a nanoscale cylindrical pore,18,19 to elucidate the dynamic transition in the crystal structure during chemical reactions at shock front prior to detonation.20 Reactive 3 ACS Paragon Plus Environment


Page 4 of 49

molecular dynamics (MD) simulations have been used extensively to investigate the chemical decomposition and mechanical response of a range of HE materials such as RDX, HMX, PETN, TNT, and NM.16,17,21–28 However, despite remarkable progress in atomistic simulations, they are not without limitations. The parameterization of the force fields to describe the tremendously intricate coupling between the mechanics of shock loading and chemical decomposition of HE materials often entails tradeoffs in the description of different physical and chemical properties. Therefore, it is not surprising that various parameterizations of ReaxFF are available. Thus, it is imperative to make a systematic effort to evaluate differences in predictions and eventually quantify uncertainties that result from the use of various force field parameterizations to develop a comprehensive and accurate understanding of shock processes. Nitromethane (NM) is an insensitive explosive material. However, under certain experimental conditions of high shock pressure it undergoes rapid detonation. Often times the sensitivity of NM is increased by the addition of a weak base, such as diethylenetriamine (DETA)3 and ethylenediamine (EDA).29 In either condition—pure or sensitized—the reaction mechanism of shocked NM is subjected to debate. The limitations in the spatiotemporal resolution required for in-situ measurements prohibits the determination of the shock Hugoniots and reaction kinetics of shocked NM. Thus, only limited information about the shock properties and reaction mechanisms are available from experiments. Recently, ultrafast dynamic ellipsometry,3 UV-vis absorption spectra,30 time-resolved Raman,7 and optical absorption spectroscopy29 techniques have been utilized to understand molecular-level changes in the shocked NM. Brown et al.3 for the first time reported us-up data for overdriven NM using ultrafast dynamic ellipsometry. They found a volume-increasing chemical reaction for shocks above a piston velocity up = 2.4 km/s for neat NM. Recently, Bhowmick et al.31 also observed volume expanding reactions in shocked NM at up=~2.3 km/s in their tabletop, laser-driven flyer plate detonation experiments. Pangilian et al.7 reported frequency hardening and broadening of the CN, CH3, and NO2 stretch/CH3 bend at a shock pressure of 14 GPa believed to be precursors to shock-induced chemical reactions. Later, Winey et al.6 applied stronger shocks on NM sample in the range of 14-17 GPa. Using timeresolved Raman spectroscopy, after a shock pressure of 17 GPa they observed the disappearance of all peaks associated with CN, NO2, and CH3, which indicates chemical decomposition reactions. 4 ACS Paragon Plus Environment

Page 5 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Interestingly, despite its simple chemical structure, the detailed chemistry of shocked NM is not known with certainty. Atomistic simulation studies seeking to understand NM decomposition chemistry have been carried out using both electronic structure and reactive force field methods. Density functional theory (DFT) based MD simulations performed on dense NM (ρ=1.97 and 2.20 g/cm3) by Manaa et al.8 at high temperatures (3000 K and 4000 K) reported both inter and intramolecular proton transfer as initiation pathways. The intermolecular proton transfer leads to the formation of CH3NO2H+ and the aci-ion H2CNO2-, while the intramolecular proton transfer produces CH2NO2H, an isomer of the NM molecule. Pressure induced initiation of the NM via intermolecular proton transfer was also observed in Car-Parrinello MD (CPMD) simulations.32 On the other hand, at ambient pressure, C-N bond scission, along with the proton transfer, has been identified as initiation pathway in CPMD simulations33 on solid NM, that illustrates the possibility of unimolecular reactions at lower pressures. Despite this progress, the computationally intensity of quantum chemistry (QC) methods hinders simulation of shock process in HE materials. In this regard, computationally less intensive classical force fields offer a practical alternative to the limitations of electronic structure methods, enabling the simulation of larger systems and longer timescales. ReaxFF simulations34 on solid and liquid phase NM (ρ=1.97g/cm3) at high temperatures (2000-3000 K), revealed initiation reactions via both inter and intramolecular proton transfer, consistent with the DFT studies,8,32 as well as isomerization reaction to form CH3ONO. In contrast, Rom et al.21 showed that, at the densities near the CJ states, ReaxFF predicts the formation of CH3NO via N-O bond cleavage as the dominant initiation pathway, while at low densities the initiation route switches to the unimolecular C-N bond breaking in agreement with CPMD simulations.33 Overall, these studies indicate a significant pressure dependency of the initiation of NM under thermal loading. Complementing the thermal decomposition work, Periott et al.35 studied shock loading on NM using Hugoniostat simulations with density functional based tight binding (DFTB) and ReaxFF methods. The predicted pressure-volume and shock velocity vs. particle velocity (us-up) curves underestimate the experimental measurements, see Marsh et al.36 Furthermore, in contrast to the laser-driven shock studies from Brown et al.3 no shift in the us-up curves was observed following reactions. In summary, despite its simple molecular structure and 5 ACS Paragon Plus Environment


Page 6 of 49

lack of microstructure in its liquid form, several questions remain open regarding the shockinduced decomposition of NM. We believe that a detailed description of the chemical response will be obtained via a combination of ultra-fast spectroscopy and reactive MD simulations, here we present our first steps in this direction. Herein, we use the reactive force field ReaxFF19,21,37,38 to simulate the response of NM to shock loading. We report unreactive and reactive Hugoniots, Chapman-Jouguet (CJ) states, detonation velocities, chemical decomposition mechanisms, and spectral changes under shock loading. In order to evaluate uncertainties associated with the description of atomic interactions, we utilized four different versions of ReaxFF and performed systematic validation of each potential. We find close agreement among the four force fields in terms of unreacted and reacted us-up curves and good agreement with experiments. Importantly, two of the descriptions show a jump in the us-up curves corresponding to exothermic reactions in good agreement with experiments. We discuss differences in the detailed chemistry and the resulting spectral evolution between force fields that could be elucidated via ultra-fast spectroscopic experiments. The paper is organized as follows. Section 2 provides simulation details, and Section 3 focuses on the evolution of thermodynamic properties under shock loading, shock Hugoniots, and CJ state properties. Section 4 describes shock-induced chemical reactions, the influence of pressure on NM decomposition kinetics, and Section 5 delineates the spectral changes during shock loading. Finally, conclusions are drawn in section 6.

II. SIMULATION DETAILS II.A. Force Fields ReaxFF,14,39 an empirical force field method, is based on the formalism of partial bond order together with environment-dependent partial atomic charges. Differentiable and continuous functional forms describe all the bonded interactions and ensure a seamless bond breaking and forming process. van der Waals and Coulomb partial energy terms are calculated between all atoms without any connectivity dependent exceptions. These two features enabled ReaxFF to model a wide range of covalent and ionic materials.40–44 For a detailed description of the ReaxFF method, readers are referred to Refs. [15,39,45]. 6 ACS Paragon Plus Environment

Page 7 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Various versions of C/H/N/O ReaxFF force fields available in literature to simulate HE materials have their inherent strength and limitations. Therefore, it is imperative to evaluate their predictions under extreme loading conditions and validate them against shock experiments or high-level QC studies to establish their accuracy in describing shock physics and chemistry. We are interested in assessing how different descriptions of atomic interactions affect the predicted reaction mechanisms, kinetics of decomposition and reaction as well as thermodynamic properties like reacted and unreacted equations of state and detonation velocity. Therefore, this study uses four available ReaxFF C/H/N/O parameter sets, namely, ReaxFF-2014,37 ReaxFF2018,19 ReaxFF-IW (inner-wall),21 and ReaxFF-lg.38 ReaxFF-2014 force field has been developed especially to improve post-detonation combustion chemistry and spectroscopic properties. However, this parameterization resulted in a poor description of the lattice parameters of the condensed phase HEs. Thus, the ReaxFF-2018 force field adds the so-called low-gradient (lg) corrections and modifies other related van der Waals parameters resulting in improved density predictions. The ReaxFF-IW and ReaxFF-lg parameter sets also differ only by the van der Waals descriptions. A detailed description of the force fields can be found in Ref. [46,19]

II.B. Nitromethane Models NM crystal structure was obtained from the Cambridge Structural Database (CSD).47 The unit cell was replicated 6 x 6 x 6 times resulting in 6048 atoms. The structure was at first relaxed using the conjugate gradient minimization scheme. The minimized geometry then heated to 400 K under NVT (constant volume, temperature) and held at the temperature for 10 ps to equilibrate liquid structure. The liquid NM was cooled to room temperature, followed by an equilibration simulation for another 10 ps. We used a time step of 0.1 fs and periodic boundary conditions for all the simulations. The experimental density of 1.13 g/cm3 was considered in this study. A snapshot of the system is shown in Figure 1.



Page 8 of 49

Figure 1. Snapshot of the simulation cell with 864 Nitromethane molecules (color code: oxygen, red; carbon, gray; hydrogen, white; nitrogen, blue)

II.C. Shock Simulations Using the Hugoniostat Method All simulations are performed using LAMMPS48 and the Hugoniostat MD simulation technique to describe the shock loading. The Hugoniostat uniaxially compresses the system and controls total energy to satisfy the Hugoniot jump conditions at a desired applied shock pressure.49,50 The computationally intensity of this method is significantly less as compared to nonequilibrium MD simulations since herein only the final states of the shocked materials are simulated, circumventing the actual shock wave propagation through the materials. Thus, it enables longer simulations time scales. The Hugoniostat equations of motion are: p i  v p r i mi

r& i 

p& i  F i  (v p  vH  ) p i L&  v p L

&

vH [E  E H (t)]   H  B0V0

 

vp B0

(   P )   p )


Page 9 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Where

    (

p i p i mi

 r i F ) / V is

the

internal

stress

tensor

and

1 EH (t )  E0  [ zz (t)  P0 ][V0  V(t)] , subscript 0 refers to the initial unshocked properties. σzz is 2

the normal component of the stress tensor in the shock direction. The Hugoniostat technique uses a thermostat and a barostat to apply the jump conditions at the imposed shock pressure. The room temperature structure described above was utilized as the initial conditions for the Hugoniostat shock simulations that were performed in the pressure range of 4 to 40 GPa.

II.D. Spectral Analysis The MD simulation trajectories were post-processed to obtain vibrational spectra via the Fourier transform of time-dependent atomic velocities. The power spectrum was obtained as:37 P   

 N

N 1

3N

m v  nt  e j

j 1

2  i 2 nt

j

n 0

where vj(t) are the atomic velocities at time t, τ=nt is the total sampling period, mj is the mass of atom j, Δt is the time increment, and N is the number of frames to be analyzed. We used fast Fourier transform to evaluate the power spectra. The atomic trajectories were sampled every 2 fs to resolve the high-frequency modes accurately. The time-resolved spectra calculations were carried out in 10 ps intervals, and each case was averaged over 5000 frames. We note that herein we are reporting power spectra including all modes, not just IR active ones. This is in contrast to FTIR spectra. IR signals can be obtained from atomistic trajectories from the time evolution of the dipole moment,51 however such analysis results in significant noise due to the short duration of the sampling period and small system sizes. While comparing total power spectra with FTIR results should be done with care, most of the modes of interest are IR active and this is particularly true as chemical decomposition occurs.



Page 10 of 49

III. SHOCK-INDUCED PHYSICAL RESPONSE III.A. Time Evolution During Shock Loading Figure 2 shows the time evolution of volume, pressure, and temperature for NM shocked at a variety of pressures and for the four force fields studied. We find that thermodynamic properties equilibrate to the unreacted shocked state within a timescale of approximately 1 ps. We note that for strong shocks this is a short-lived metastable state as chemical reactions follow the shock loading. The fast compression of the shocked material is followed by an induction time before the onset of rapid chemical reactions for strong enough shock pressures. We observe that the inter and intramolecular temperature of the NM molecules equilibrate within ~5ps, similar to the timescales observed in non-equilibrium simulations,52 see Figure S1 in the supplementary information (SI). Shock loading results in exothermic reactions within the 400 ps of simulation for pressures higher than or equal to 16 GPa, 12 GPa, 32 GPa, and 28 GPa for ReaxFF-2014, ReaxFF-2018, ReaxFF-IW, and ReaxFF-lg, respectively. The induction time for chemical reaction is strongly correlated to the temperature after the initial compression—due to the mechanical (PV) work done on the system. The threshold shock pressure for reaction predicted by ReaxFF-2014 (16 GPa), is in good agreement with the reported experimental values in neat NM of 17 GPa, 18 GPa, and 19 GPa by Winey,6 Bhowmick et al.,31 and Brown et al.3, respectively. ReaxFF-2018 result predicts a slightly lower threshold value and ReaxFF-IW and –lg force fields require higher pressures. The exothermic chemistry rapidly increases system temperature by ~3000 K in ReaxFF-2014 and 2018 and by ~2000K in ReaxFF-IW, -lg simulations. The higher sensitivity of the 2014 and 2018 force fields correlates with higher exothermicity. In comparing threshold shock strengths for reaction one should note that our simulations use classical (as opposed to quantum) ionic dynamics which results in an overestimation of the specific heat and underestimation of the shock temperature.53,54 Quantum statistics would result in an increase of unreacted shock temperature and a reduction of critical shock strength to induce chemistry.


Page 11 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Evolution of pressure, volume, and temperature during Hugoniot shock simulations of NM at different pressures using ReaxFF-2014, ReaxFF-2018, ReaxFF-IW, and ReaxFF-lg force fields. The pressure-density and temperature-density locus of all shock simulation, Figure 3, exhibit two distinct regimes – the nonreacted and reacted Hugoniots. Our predicted unreacted pressuredensity data is in good agreement with the LASL shock Hugoniot data.36 At the onset of chemistry, the nonreactive and reactive Hugoniots branch out for all force fields. The shock state shifts to the reactive Hugoniots due to the volume expanding exothermic reactions. Our predicted temperature-density data of reactive Hugoniots are comparable with the NM detonation modeling performed by Menikoff et al.55 For example, at 18 GPa shock pressure, they reported product temperature as ~4600K, while ReaxFF-2014 and -2018 simulations predict ~4300K. 11 ACS Paragon Plus Environment


Page 12 of 49

Figure 3. (a) Pressure vs. density and (b) Temperature vs. density for NM samples at various shock pressures using four different force fields and comparison with experimental data. Color code: red, ReaxFF-2014; cyan, ReaxFF-2018; green, ReaxFF-IW; blue, ReaxFF-lg; and black, experimental data.36

III.B. Shock Velocity vs. Particle Velocity Figure 4 shows the shock velocity us as a function of particle velocity up. The exothermic volume-increasing reactions lead to an upward shift in the us values. We find that the calculated unreacted us-up data is similar for all the force fields and in good agreement with experimental results.3,31,36 The sudden change in the us values due to the rapid chemistry is more apparent in ReaxFF-2014 and ReaxFF-2018 compared to other two force fields. Importantly, the ReaxFF-2014 and ReaxFF-2018 simulation predicted us-up data are in excellent agreement with the experimental data for both the reactive and nonreactive regimes. Ultrafast dynamic ellipsometry3 on laser-driven and photon Doppler velocimetry31 on laser-driven flyer plate shock studies predict a transition to the reacted Hugoniot for particle velocities between ~2.3 and 2.4 km/s. ReaxFF-2014 predicts the transition to reactive states at up~2.2 km/s, in good agreement with these experiments3,31 and mesoscale modeling using SESAME EOS55 data. ReaxFF-2018 predicts earlier chemistry and shows a shift in us values at up~1.7 km/s, and ReaxFF-lg and ReaxFFIG predict relatively higher values of ~3.3 and ~3.5 km/s, respectively. We also calculated sound speeds by extrapolating the unreacted us-up data to zero piston velocity. The ReaxFF-2014, ReaxFF-2018, ReaxFF-IW, and ReaxFF-lg predicted sound speeds are 12 ACS Paragon Plus Environment

Page 13 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.86±0.05 km/s, 1.67±0.03 km/s, and 1.58±0.01, and 1.35±0.04 km/s, respectively. These data agrees reasonably well with the experimentally reported sound speeds of 1.6531 and 1.3536 km/s.

Figure 4. us-up data for NM samples at various shock pressures using three different force fields and comparison with experimental data. Color code: red, ReaxFF-2014; cyan, ReaxFF-2018; green, ReaxFF-IW; blue, ReaxFF-lg, and black circle, Brown et al.3; black rectangle, Marsh et al.36; black cross, Bhowmick et al.31

III.C. Crussard Curves and Chapman-Jouguet (CJ) State In order to determine Chapman–Jouguet (CJ) states, we constructed Hugoniot curves for the products (known as Crussard curves) by performing Hugoniostat simulations on the reaction products of shocked NM. Previously, Guo et al.56 proposed a methodology (Rx2CJ) where they performed a set of NVT-MD simulations over a range of temperature as a function of density to obtain a family of Hugoniot function values (Hg). Next, from these curves, they identified Hg = 0 points to obtain the pressure and temperature for determining the Crussard curve. In our approach, we take the trajectories of the reaction products from shock simulations where products gases reach equilibrium. The pressure ranges considered for constructing Crussard curves using the four ReaxFF parameterizations are shown in Figure 5. The CJ point corresponds 13 ACS Paragon Plus Environment


Page 14 of 49

to the point determined by the tangent of the Crussard curve that passes through the initial state P0, V0 point; the line joining the CJ point and initial states is called Rayleigh line. Importantly, the detonation velocity is the square root of the slope of the Rayleigh line. To identify the CJ point from the Crussard curves, we fitted a power law function to the simulation data to express the evolution of pressure and temperature as a function of compression ratio (n=V/Vo) as presented in Figure 5 to obtain the tangent curves analytically. We describe the Crussard curve using a power law: P  n   ao n  a1 . At the CJ point, the slope of the Rayleigh line is equal to the derivative of the pressure with respect to the compression ratio.

P  Po dP(n) (CJ )  CJ dn nCJ  1

(1)

Assuming Po CH3NO+O reaction using B3LYP exchange and correlation function with a 6-311++G** basis set using ORCA.69 We considered both singlet and triplet states and the dissociation energies obtained from the minimum energy pathway is 94.2 kcal/mol. ReaxFF-2014, ReaxFF17 ACS Paragon Plus Environment


Page 18 of 49

2018, ReaxFF-IW, and ReaxFF-lg predict 50, 46, 92, and 92 kcal/mol, respectively [See Fig S3]. These data show that ReaxFF-2014 and -2018 can be expected to overpredict the scission of the N-O bond. We find that the qualitative distribution of intermediates produced across various force fields are comparable, while the magnitudes of their population are significantly lower in ReaxFF-IW and –lg simulations compared to the other two force fields. This can be attributed to the formation of molecular clusters followed by the onset of volume increasing reactions. In particular, the large molecular clusters are observed in the cases of ReaxFF-lg and –IW force fields that results in smaller quantities of the intermediates of interest. The details of the cluster analysis are described in Section V of the SI.

Figure 6. Molecular population distribution of NM at shock pressures of 18 and 24 GPa using ReaxFF-2014 force fields. Molecular fractions are normalized by the initial number of NM molecules.


Page 19 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. Molecular population distribution of NM at shock pressures of 12 and 16 GPa using ReaxFF-2018 force fields.



Page 20 of 49

Figure 8. Molecular population distribution of NM at shock pressures of 32 and 36 GPa using ReaxFF-IW force fields.

Figure 9. Molecular population distribution of NM at shock pressures of 28 and 32 GPa using ReaxFF-lg force fields.

IV.B. Final Products Following initiation, rapid chemistry causes the system temperature and volume to raise abruptly for sufficiently strong shocks. Stable chemical species such as H2O, CO2, N2, NH3, and H2 are produced at a pronounced rate (see Fig 6 to 9). The major secondary reaction pathways leading to the production of H2O, CO2, and N2 as predicted by ReaxFF-2014 and –lg force fields are discussed in Section V of the SI. Overall, we found H2O to be the dominant reaction product irrespective of simulation conditions and force fields used, consistent with NM detonation experiments.70 Figure 10 shows the populations of final products at the CJ point for all force fields as well as experimental data.70 The simulations nicely capture the overall trend in the relative quantity of the stable gases compared to the experiment. The absolute quantities of the ReaxFF-2014, 20 ACS Paragon Plus Environment

Page 21 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ReaxFF-IW, and ReaxFF-lg predicted H2O, CO2, N2, and NH3 and ReaxFF-2018 predicted CO2 and NH3 are in excellent agreement with the experimental data. This provides important validation of the force fields description of the thermodynamics of reaction products and, in a more indirect manner, the reaction mechanisms. We note discrepancies in the population of H2 predicted by ReaxFF-2014 and -2018, and H2O and N2 in ReaxFF-2018 simulations. Since we produce complete Crussard curves, we can assess how composition depends on pressure along the Hugoniot of the products. The mole fraction of stable gases produced as a function of pressure and experimental values from detonation calorimetric data (plotted at CJ pressure)70 are shown in Figure S6(a-e). One can see that the quantity of the product gases depends strongly on pressure along the products EOS. The reduction of the shock pressures results in an increase in the number of stable gases across all the force fields. The ReaxFF-IW and -lg simulations produce higher amount of H2O, N2, and H2 compared to the other two force fields. ReaxFF-2014 and -2018 results in only trace amounts of H2 molecules irrespective of the shock pressures. The lowest quantity of H2O and N2 are found in ReaxFF-2018 simulations. The CO2 production is found to be relatively insensitive to the force fields used.

Figure 10. Stable final products at the CJ pressure predicted by the respective force fields and comparison with the experiment.70 Finally, we investigated the presence of molecular clusters that are known to influence overall species distribution and detonation properties. While ReaxFF-2014 and ReaxFF-2018 predict no 21 ACS Paragon Plus Environment


Page 22 of 49

noticeable amount of cluster formation, we observe that about 15% of the system mass is retained in the clusters in the case of ReaxFF-lg and ReaxFF-IW simulations. These clusters are carbon-rich and impact the exothermicity of the reactions and specific volume of the products. We believe the propensity for cluster formation is linked with the lower detonation velocities in ReaxFF-lg and -IW as compared to the ReaxFF-2014 and -2018. Details of the cluster formation analysis are included in Section V of the SI.

IV.C. Shock Pressure Dependency on NM Decomposition Rates To investigate the effect of pre-reaction shock-induced pressure and temperature on NM decomposition, we extracted reaction rates assuming first-order reaction kinetics. The number of undissociated NM was fitted as a function of time using the following equation

ln

N  kt No

Where N is the number undissociated NM at time t, No is the initial number of NM, and k is the rate constant. The rate constant data were extracted via fitting unreacted NM profile as a function of time up to ~75% of remaining NM. The data are presented as a function of unreacted temperatures in Figure 11, for all the force fields. One can see that NM decomposition rate follows an exponential behavior as a function of both the pressure and temperature. The increase in shock strength resulted in higher temperatures of the shocked unreacted NM consequently lead to the increase in the reaction kinetics. It is also evident from the figure that the ReaxFF2018 and ReaxFF-2014 predicted reaction rates are significantly higher than the other two force fields.


Page 23 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11. NM decomposition rates as a function of the unreacted peak temperature due to the shock compression for all the ReaxFF parameter sets. k is the reaction rate constant in s-1. Color code: red, ReaxFF-2014; cyan, ReaxFF-2018, green, ReaxFF-IW; and blue, ReaxFF-lg.

V. EVOLUTION OF VIBRATIONAL SPECTRAL FEATURES DURING SHOCK LOADING Establishing detailed chemistry of materials subjected to shock loading remains a grand challenge in condensed matter physics and chemistry. As mentioned above, ultra-fast spectroscopy is providing information about chemical processes with picosecond time resolution, which can contribute to addressing this challenge. Unfortunately, interpreting these experiments, specially deconvoluting spectral peaks, is challenging at high pressures and for fast exothermic processes. Reactive MD simulations can also predict the time evolution of spectra during shock loading enabling direct comparison with experiments and, given its atomistic nature, correlating spectral features to atomic-level processes. Thus, we expect these simulations to facilitate peak assignment of the experimental spectra and contribute to the identification of various chemical species and reaction pathways. The calculated spectra can effectively shed light into the explanation of reaction initiation mechanisms, where conventional experimental techniques for chemical species recognition cannot perform due to the extreme conditions. In order to assign vibrational peaks to vibrational modes, we performed normal mode analysis on 23 ACS Paragon Plus Environment


Page 24 of 49

an isolated NM molecule via diagonalization of the Hessian matrix. We also performed DFT calculations in ORCA using B3LYP functional and 6-311++G** basis set to obtain normal modes of an NM molecule. The normal modes calculated from ReaxFF-2014 and –lg force fields are presented in Fig S9, S10, and Table S2 compares the calculated data with experimental values.71 We find that the ReaxFF-2014 and ReaxFF-lg force fields quantitatively predict most of the normal mode frequencies, while a few deviate from the experiments and DFT data. We note that some of the deviations can be attributed to the mixed modes predicted by ReaxFF simulations. For example, the NO2 symmetric mode frequency is significantly blue-shifted, but, that can be due to the fact that this mode is combined with the CH2 stretching frequency. The peak assignment to our calculated time convoluted spectra is shown in Figure S11. As discussed in Section II, we computed the time evolution of the power spectra corresponding several different shock pressures—unreactive and reactive cases—using simulation trajectories from ReaxFF-2014 and ReaxFF-lg force fields. We find that ReaxFF-2014 and ReaxFF-2018 as well as ReaxFF-IW and ReaxFF-lg force fields exhibit qualitatively similar spectral behavior. The time evolution of the spectra calculated from the Hugoniostat simulations for various conditions are shown in Figure S8. However, comparing Hugoniostat simulations with experiments require some post-processing of the simulation results. In the Hugoniostat simulations the entire system is shocked simultaneously, this is in contrast with shock experiments where the shock wave gradually propagates through the system. Thus, the experimental results averaged over unshocked material and material that has been shocked different amounts of times. Therefore, in order to facilitate comparison with experimentally derived spectroscopic data, we performed a time convolution of the calculated time-dependent and unshocked materials’ spectra. The convolution is achieved via the weighted sum of the unshocked and the series of time-dependent spectra derived from Hugoniostat simulations at various proportions. We considered samples 2 μm in thickness, typical of laser shock experimentals,2,72, and the actual shock velocity corresponding to each case from the simulations. The time-convoluted spectra for samples shocked below the threshold for chemistry exhibit broadening and shifting of the signature peaks, but all initial peaks remain distinguishable and


Page 25 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


no new features appear, see left panels Figure 12. In contrast, chemical reactions induce striking changes in the spectra, see right panels in Figure 12.

Figure 12. Time convoluted full spectra as predicted using ReaxFF-2014 and ReaxFF-lg force fields: (a) 12GPa, (b) 18 GPa ReaxFF-2014; (c) 20 GPa, and (d) 28 GPa ReaxFF-lg force fields. The left panel is unreactive and the right panel is for the reactive cases. The negative time indicates unshocked materials. ν, δ, and ρ indicates stretching, bending, and rocking modes, respectively. The NO2 symmetric stretching mode (at ~1780 cm-1 for ReaxFF-2014 and 2300 cm-1 for ReaxFF-lg) is a good measure of the onset of rapid chemical reactions. The timescale of the disappearance of this peak is associated with the induction time for the shock-induced chemistry at the particular pressure. Thus, the disappearance of the NO2 peak is correlated to the decomposition of NM molecules. Our simulations indicate that C-N bond breaking is less favorable than the dissociation of N-O bonds; our analysis on the evolution of chemical bonds (see Figure S4(a)) corroborates the retainment of the C-N bonds while approximately 99% of the 25 ACS Paragon Plus Environment


Page 26 of 49

N-O bonds dissociate due to the exothermic chemistry. In Figure 12 (b), we observe the persistence of the C-N peak until the onset of rapid chemistry; then the peak is lost within a broad low-frequency feature. Interestingly, the disappearance of the NO2 peak is more evident in the spectra as the frequency region around with the peak is relatively less active with reaction intermediates and products. The formation of a large number of small molecules due to rapid exothermic reactions and the intermolecular collisions spawn a diffuse, broad peak at low and high frequencies that complicate precise molecular assignment; this exemplified the challenge in determining detailed chemistry from spectral analysis at extreme conditions.

Figure 13. Spectra of key intermediates and final products as extracted from the ReaxFF-2014 and ReaxFF-lg shock simulations at shock pressures of 18 GPa, and 28 GPa, respectively. We foresee that atomistic simulations, such as those reported in this paper, will help deconvolute the broad spectra typical of chemistry under shock conditions and, together with experiments, enable the identification of decomposition and reaction paths. To contribute to 26 ACS Paragon Plus Environment

Page 27 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


such efforts, we extracted the time evolution of key intermediates and final products form the overall shock simulation and performed spectral analysis on the individual species. Figure 13 shows the resulting spectra of key intermediates and products for both ReaxFF-2014 and ReaxFFlg force fields. We stress that these molecular spectra reflect the actual chemical environment, pressure and temperature during shock-induced reactions. The results indicate that the highfrequency peaks around 3000 cm-1 are due to the formation of gas species like H2O, N2, and NH3. We also observe the appearance of a carbon-oxygen peak near 1000 cm−1 (see peak assignment in Figure S12c) resulting from the evolution of CO2 species shown in Figs. 13(b) and (d). While we notice discrepancies between the two force fields, we find that N2 and NH3 result in sharper highfrequency peaks than the other product species. Further, we can derive insight into the reaction mechanisms from the spectra of the intermediates shown in Figs. 13(a) and (c). For example, the appearance of peaks in the vicinity of 1600 cm-1 is indicative to the production of HONO, HNO3 or HNO. In contrast, characteristic peaks at the range of 800-1000 cm-1 depict the intermediates with C-O bonds such as HOCN, HCHO, and HCO2. Such information will be valuable in discerning reaction paths observed experimentally. While these results include uncertainties associated with the force fields, the individual molecular contributions will be key to reconstruct the experimental spectra and identify possible decomposition pathways. We believe ab initio calculations of spectra at various times during the decomposition of HE materials will be important both to validate force field predictions (and tune force fields as needed) and to contribute to the interpretation of the experimental spectra. These computationally intensive ab initio simulations could build on force field simulations, for example, by using configurations under various degrees of reaction as starting point. Finally, in a recent laser-driven shock experiment, McGrane et al.73 reported preliminary spectroscopy data of shocked NM. Follow up work at Los Alamos under various shock strengths should provide data directly comparable to the results presented in this paper. In the experiments, changes in transmission was reported, and in order to enable a direct comparison with such data, we estimate our calculated spectra in terms of transmission and we compare them with the experimental data in the SI Section IX.



Page 28 of 49

VI. CONCLUSIONS We investigated the detailed shock response behavior of liquid nitromethane using four different versions of the ReaxFF force field to understand the uncertainties in the predictions resulting in due to the choice of force fields. We found that ReaxFF-2014 predicted threshold pressure for shock initiation and the calculated us-up data are in excellent agreement with shock experiments. The ReaxFF-IW and –lg versions require stronger shocks for initiation, while ReaxFF2018 yields chemical reactions at a lower pressure. Our calculated detonation velocities using ReaxFF-2014 and -2018 are within a reasonable agreement with literature results, while ReaxFFIW, and –lg underestimates the data. We found that the formation of large clusters predicted by these two force fields can be linked to the lower detonation velocities and exothermicities. We find the Initial decomposition to be dominated by bimolecular reactions that produce nitrosomethane (CH3NO) as a primary intermediate species, consistent with the experimental observation of a decrease in unimolecular reactions at high pressures. We described the detailed decomposition mechanisms under shock loading for all the force fields and found relatively similar final products which are in excellent agreement with experiments. We find differences in the amounts of intermediates which can be attributed to the different reaction pathways predicted by the force fields. Especially, larger cluster formation in ReaxFF-IW and ReaxFF-lg force fields leads to the smaller quantity of the intermediates. A detailed analysis of the MD trajectories produces time-resolved spectra during decomposition that can be directly compared with ongoing laser-driven shock experiments towards the development of a definite understanding of reaction initiation mechanisms and timescales. From the simulations we extract the contribution of key intermediates and products to the spectra that could be used to interpret experimental data. As with any physics-based models, the MD simulations in this paper involve multiple uncertainties.74 Key among them is the description of atomic interactions (i.e. the force field used) and our use of four parametrizations aims at quantifying its effect on the various quantities of interest. Relating atomic interactions to quantities such as detonation velocity, kinetics or decomposition path in a quantitative manner is not possible, as these are emerging 28 ACS Paragon Plus Environment

Page 29 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


properties. Our results indicate correlations between properties that could be useful in future efforts to refine these potentials. For example, lower exothermicities correlates with a tendency to predict larger molecular clusters, lower detonation velocities and higher thresholds for initiation. The variability in the predictions of various properties provides a rough estimate of the model-form uncertainties; our team is currently assessing these uncertainties over a wide range of HE materials with the aim of quantifying the degree to which these simulations can predict relative properties in addition to absolute values.

AUTHOR INFORMATION Corresponding Author *AS ([email protected])

Notes The authors declare no competing financial interest. ASSOCIATED CONTENT Supporting Information The supporting Information is available free of charge on the ACS Publications website at DOI: Further details of the formation energies of key species, analysis of the evolution of chemical bonds, clusters formation, and resulted final products in the Hugoniostat simulations, reactions pathways leading to the key final products, as well as spectra and normal mode analysis. Acknowledgments This work was support by the US Office of Naval Research, Multidisciplinary University Research Initiatives (MURI) Program, Contract: N00014-16-1-2557. Program managers: Clifford Bedford and Kenny Lipkowitz.



Page 30 of 49

REFERENCES (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16)

(17)

Bdzil, J. B.; Stewart, D. S. The Dynamics of Detonation in Explosive Systems. Annu. Rev. Fluid Mech. 2007, 39, 263–292. McGrane, S. D.; Moore, D. S.; Funk, D. J. Shock Induced Reaction Observed via Ultrafast Infrared Absorption in Poly(Vinyl Nitrate) Films. J. Phys. Chem. A 2004, 108, 9342–9347. Brown, K. E.; McGrane, S. D.; Bolme, C. A.; Moore, D. S. Ultrafast Chemical Reactions in Shocked Nitromethane Probed with Dynamic Ellipsometry and Transient Absorption Spectroscopy. J. Phys. Chem. A 2014, 118, 2559–2567. Deàk, J. C.; Iwaki, L. K.; Dlott, D. D. Vibrational Energy Redistribution in Polyatomic Liquids: Ultrafast IR−Raman Spectroscopy of Nitromethane. J. Phys. Chem. A 1999, 103, 971–979. Hemsing, W. F. Velocity Sensing Interferometer (VISAR) Modification. Rev. Sci. Instrum. 1979, 50, 73–78. Winey, J. M.; Gupta, Y. M. Shock-Induced Chemical Changes in Neat Nitromethane: Use of Time-Resolved Raman Spectroscopy. J. Phys. Chem. B 1997, 101, 10733–10743. Pangilinan, G. I.; Gupta, Y. M. Time-Resolved Raman Measurements in Nitromethane Shocked to 140 Kbar. J. Phys. Chem. 1994, 98, 4522–4529. Riad Manaa, M.; Reed, E. J.; Fried, L. E.; Galli, G.; Gygi, F. Early Chemistry in Hot and Dense Nitromethane: Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 10146–10153. Manaa, M. R.; Fried, L. E.; Melius, C. F.; Elstner, M.; Frauenheim, T. Decomposition of HMX at Extreme Conditions: A Molecular Dynamics Simulation. J. Phys. Chem. A 2002, 106, 9024–9029. Schweigert, I. V. Ab Initio Molecular Dynamics of High-Temperature Unimolecular Dissociation of Gas-Phase RDX and Its Dissociation Products. J. Phys. Chem. A 2015, 119, 2747–2759. Schweigert, I. V. Quantum Mechanical Simulations of Condensed-Phase Decomposition Dynamics in Molten RDX. J. Phys. Conf. Ser. 2014, 500, 052039. Patidar, L.; Thynell, S. T. Quantum Mechanics Investigation of Initial Reaction Pathways and Early Ring-Opening Reactions in Thermal Decomposition of Liquid-Phase RDX. Combust. Flame 2017, 178, 7–20. Isayev, O.; Gorb, L.; Qasim, M.; Leszczynski, J. Ab Initio Molecular Dynamics Study on the Initial Chemical Events in Nitramines: Thermal Decomposition of CL-20. J. Phys. Chem. B 2008, 112, 11005–11013. Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409. Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; et al. The ReaxFF Reactive Force-Field: Development, Applications and Future Directions. Npj Comput. Mater. 2016, 2, 15011. Zhou, T.; Zybin, S. V.; Liu, Y.; Huang, F.; Iii, W. A. G. Anisotropic Shock Sensitivity for βOctahydro-1,3,5,7-Tetranitro-1,3,5,7-Tetrazocine Energetic Material under CompressiveShear Loading from ReaxFF-Lg Reactive Dynamics Simulations. J. Appl. Phys. 2012, 111, 124904. Zhou, T.-T.; Huang, F.-L. Effects of Defects on Thermal Decomposition of HMX via ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 278–287. 30 ACS Paragon Plus Environment

Page 31 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Wood, M. A.; Cherukara, M. J.; Kober, E. M.; Strachan, A. Ultrafast Chemistry under Nonequilibrium Conditions and the Shock to Deflagration Transition at the Nanoscale. J. Phys. Chem. C 2015, 119, 22008–22015. (19) Wood, M. A.; Kittell, D. E.; Yarrington, C. D.; Thompson, A. P. Multiscale Modeling of Shock Wave Localization in Porous Energetic Material. Phys. Rev. B 2018, 97, 014109. (20) Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C.; Goddard III, W. A. Dynamic Transition in the Structure of an Energetic Crystal during Chemical Reactions at Shock Front Prior to Detonation. Phys. Rev. Lett. 2007, 99, 148303. (21) Rom, N.; Zybin, S. V.; van Duin, A. C. T.; Goddard, W. A.; Zeiri, Y.; Katz, G.; Kosloff, R. DensityDependent Liquid Nitromethane Decomposition: Molecular Dynamics Simulations Based on ReaxFF. J. Phys. Chem. A 2011, 115, 10181–10202. (22) Furman, D.; Kosloff, R.; Dubnikova, F.; Zybin, S. V.; Goddard, W. A.; Rom, N.; Hirshberg, B.; Zeiri, Y. Decomposition of Condensed Phase Energetic Materials: Interplay between Uniand Bimolecular Mechanisms. J. Am. Chem. Soc. 2014, 136, 4192–4200. (23) Rom, N.; Hirshberg, B.; Zeiri, Y.; Furman, D.; Zybin, S. V.; Goddard, W. A.; Kosloff, R. FirstPrinciples-Based Reaction Kinetics for Decomposition of Hot, Dense Liquid TNT from ReaxFF Multiscale Reactive Dynamics Simulations. J. Phys. Chem. C 2013, 117, 21043– 21054. (24) Zhang, L.; Zybin, S. V.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A.; Kober, E. M. Carbon Cluster Formation during Thermal Decomposition of Octahydro-1,3,5,7-Tetranitro-1,3,5,7Tetrazocine and 1,3,5-Triamino-2,4,6-Trinitrobenzene High Explosives from ReaxFF Reactive Molecular Dynamics Simulations. J. Phys. Chem. A 2009, 113, 10619–10640. (25) Joshi, K. L.; Chaudhuri, S. Reactive Simulation of the Chemistry behind the CondensedPhase Ignition of RDX from Hot Spots. Phys. Chem. Chem. Phys. 2015, 17, 18790–18801. (26) Joshi, K.; Chaudhuri, S. Observation of Deflagration Wave in Energetic Materials Using Reactive Molecular Dynamics. Combust. Flame 2017, 184, 20–29. (27) Shan, T.-R.; Thompson, A. P. Shock-Induced Hotspot Formation and Chemical Reaction Initiation in PETN Containing a Spherical Void. In Journal of Physics: Conference Series; IOP Publishing, 2014; Vol. 500, p 172009. (28) Li, Y.; Vashishta, P. Multistage Reaction Pathways in Detonating High Explosives. Appl. Phys. Lett. 2014, 105, 204103. (29) Constantinou, C. P.; Winey, J. M.; Gupta, Y. M. UV/Visible Absorption Spectra of Shocked Nitromethane and Nitromethane-Amine Mixtures up to a Pressure of 14 GPa. J. Phys. Chem. 1994, 98, 7767–7776. (30) Winey, J. M.; Gupta, Y. M. UV−Visible Absorption Spectroscopy To Examine Shock-Induced Decomposition in Neat Nitromethane. J. Phys. Chem. A 1997, 101, 9333–9340. (31) Bhowmick, M.; Nissen, E. J.; Dlott, D. D. Detonation on a Tabletop: Nitromethane with High Time and Space Resolution. J. Appl. Phys. 2018, 124, 075901. (32) Citroni, M.; Bini, R.; Pagliai, M.; Cardini, G.; Schettino, V. Nitromethane Decomposition under High Static Pressure. J. Phys. Chem. B 2010, 114, 9420–9428. (33) Chang, J.; Lian, P.; Wei, D.-Q.; Chen, X.-R.; Zhang, Q.-M.; Gong, Z.-Z. Thermal Decomposition of the Solid Phase of Nitromethane: Ab Initio Molecular Dynamics Simulations. Phys. Rev. Lett. 2010, 105, 188302. 31 ACS Paragon Plus Environment


Page 32 of 49

(34) Han, S.; van Duin, A. C. T.; Goddard, W. A.; Strachan, A. Thermal Decomposition of Condensed-Phase Nitromethane from Molecular Dynamics from ReaxFF Reactive Dynamics. J. Phys. Chem. B 2011, 115, 6534–6540. (35) Perriot, R.; Negre, C. F. A.; McGrane, S. D.; Cawkwell, M. J. Density Functional Tight Binding Calculations for the Simulation of Shocked Nitromethane with LATTE-LAMMPS. AIP Conf. Proc. 2018, 1979, 050014. (36) Marsh, S. P. LASL Shock Hugoniot Data; Univ of California Press, 1980; Vol. 5. (37) Wood, M. A.; van Duin, A. C.; Strachan, A. Coupled Thermal and Electromagnetic Induced Decomposition in the Molecular Explosive ΑHMX; A Reactive Molecular Dynamics Study. J. Phys. Chem. A 2014, 118, 885–895. (38) Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. ReaxFF-Lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115, 11016–11022. (39) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040–1053. (40) Islam, M. M.; Ostadhossein, A.; Borodin, O.; Yeates, A. T.; Tipton, W. W.; Hennig, R. G.; Kumar, N.; Duin, A. C. T. van. ReaxFF Molecular Dynamics Simulations on Lithiated Sulfur Cathode Materials. Phys. Chem. Chem. Phys. 2015, 17, 3383–3393. (41) Neyts, E. C.; van Duin, A. C. T.; Bogaerts, A. Changing Chirality during Single-Walled Carbon Nanotube Growth: A Reactive Molecular Dynamics/Monte Carlo Study. J. Am. Chem. Soc. 2011, 133, 17225–17231. https://doi.org/10.1021/ja204023c. (42) van Duin, A. C.; Merinov, B. V.; Han, S. S.; Dorso, C. O.; Goddard Iii, W. A. ReaxFF Reactive Force Field for the Y-Doped BaZrO3 Proton Conductor with Applications to Diffusion Rates for Multigranular Systems. J. Phys. Chem. A 2008, 112, 11414–11422. (43) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91, 098301. (44) Onofrio, N.; Guzman, D.; Strachan, A. Atomic Origin of Ultrafast Resistance Switching in Nanoscale Electrometallization Cells. Nat. Mater. 2015, 14, 440–446. (45) Russo Jr., M. F.; van Duin, A. C. T. Atomistic-Scale Simulations of Chemical Reactions: Bridging from Quantum Chemistry to Engineering. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. At. 2011, 269, 1549–1554. (46) Islam, M. M.; Strachan, A. Decomposition and Reaction of Polyvinyl Nitrate under Shock and Thermal Loading: A ReaxFF Reactive Molecular Dynamics Study. J. Phys. Chem. C 2017, 121 (40), 22452–22464. (47) Groom, C. R.; Bruno, I. J.; Lightfoot, M. P.; Ward, S. C. The Cambridge Structural Database. Acta Crystallogr. Sect. B Struct. Sci. Cryst. Eng. Mater. 2016, 72, 171–179. (48) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (49) Maillet, J.-B. Uniaxial Hugoniostat: A Method for Atomistic Simulations of Shocked Materials. Phys. Rev. E 2000, 63. (50) Ravelo, R.; Holian, B. L.; Germann, T. C.; Lomdahl, P. S. Constant-Stress Hugoniostat Method for Following the Dynamical Evolution of Shocked Matter. Phys. Rev. B 2004, 70, 014103. 32 ACS Paragon Plus Environment

Page 33 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(51) Boulard, B.; Kieffer, J.; Phifer, C. C.; Angell, C. A. Vibrational Spectra in Fluoride Crystals and Glasses at Normal and High Pressures by Computer Simulation. J. Non-Cryst. Solids 1992, 140, 350–358. (52) Jaramillo, E.; Sewell, T. D.; Strachan, A. Atomic-Level View of Inelastic Deformation in a Shock Loaded Molecular Crystal. Phys. Rev. B 2007, 76. (53) Lin, K.-H.; Holian, B. L.; Germann, T. C.; Strachan, A. Mesodynamics with Implicit Degrees of Freedom. J. Chem. Phys. 2014, 141, 064107. (54) Qi, T.; Reed, E. J. Simulations of Shocked Methane Including Self-Consistent Semiclassical Quantum Nuclear Effects. J. Phys. Chem. A 2012, 116, 10451–10459. (55) Menikoff, R.; Shaw, M. S. Modeling Detonation Waves in Nitromethane. Combust. Flame 2011, 158, 2549–2558. (56) Guo, D.; V. Zybin, S.; An, Q.; Iii, W. A. G.; Huang, F. Prediction of the Chapman–Jouguet Chemical Equilibrium State in a Detonation Wave from First Principles Based Reactive Molecular Dynamics. Phys. Chem. Chem. Phys. 2016, 18, 2015–2022. (57) Marquardt, D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. J. Soc. Ind. Appl. Math. 1963, 11, 431–441. (58) Levenberg, K. A Method for the Solution of Certain Non-Linear Problems in Least Squares. Q. Appl. Math. 1944, 2, 164–168. (59) Mochalova, V. M.; Torunov, S. I.; Utkin, A. V.; Garanin, V. A. 14th International Detonation Symposium, Idaho, April 11-16, 2010. (60) Tarver, C. M.; Shaw, R.; Cowperthwaite, M. Detonation Failure Diameter Studies of Four Liquid Nitroalkanes. J. Chem. Phys. 1976, 64, 2665–2673. (61) Yoo, C. S.; Holmes, N. C.; Souers, P. C. Time-Resolved Temperatures of Shocked and Detonating Energetic Materials. AIP Conf. Proc. 1996, 370, 913–916. (62) Bouyer, V.; Darbord, I.; Hervé, P.; Baudin, G.; Le Gallic, C.; Clément, F.; Chavent, G. Shockto-Detonation Transition of Nitromethane: Time-Resolved Emission Spectroscopy Measurements. Combust. Flame 2006, 144, 139–150. (63) Engelke, R. Effect of a Physical Inhomogeneity on Steady-State Detonation Velocity. Phys. Fluids 1979, 22, 1623–1630. (64) Bardo, R. D. Theoretical Calculations of Rate-Determining Steps for Ignition of Shocked, Condensed Nitromethane. Int. J. Quantum Chem. 1986, 30, 455–469. (65) Piermarini, G. J.; Block, S.; Miller, P. J. Effects of Pressure on the Thermal Decomposition Kinetics and Chemical Reactivity of Nitromethane. J. Phys. Chem. 1989, 93, 457–462. (66) Mueller, K. H. The Thermal Decomposition of Nitromethane at High Pressures. J. Am. Chem. Soc. 1955, 77, 3459–3462. (67) Taylor, H. A.; Vesselovsky, V. V. The Thermal Decomposition of Nitromethane. J. Phys. Chem. 1935, 39, 1095–1102. (68) Cottrell, T. L.; Graham, T. E.; Reid, T. J. The Thermal Decomposition of Nitromethane. Trans. Faraday Soc. 1951, 47, 584–590. (69) Neese, F. The ORCA Program System. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 2, 73−78. (70) Ornellas, D. L. Heat and Products of Detonation of Cyclotetramethylenetetranitramine, 2,4,6-Trinitrotoluene, Nitromethane, and Bis[2,2-Dinitro-2-Fluoroethyl]Formal. J. Phys. Chem. 1968, 72, 2390–2394. 33 ACS Paragon Plus Environment


Page 34 of 49

(71) Smith, D. C.; Pan, C.; Nielsen, J. R. Vibrational Spectra of the Four Lowest Nitroparaffins. J. Chem. Phys. 1950, 18, 706–712. (72) Moore, D. S.; McGrane, S. D. Comparative Infrared and Raman Spectroscopy of Energetic Polymers. J. Mol. Struct. 2003, 661–662, 561–566. (73) McGrane, S.; Bowlan, P.; Powell, M.; Brown, K.; Bolme, C. Broadband Mid-Infrared Measurements for Shock-Induced Chemistry. AIP Conf. Proc. 2018, 1979, 130004. (74) Alzate-Vargas, L.; Fortunato, M. E.; Haley, B.; Li, C.; Colina, C. M.; Strachan, A. Uncertainties in the Predictions of Thermo-Physical Properties of Thermoplastic Polymers via Molecular Dynamics. Model. Simul. Mater. Sci. Eng. 2018, 26, 065007.


Page 35 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

ACS Paragon Plus Environment

Page 36 of 49

Page 37 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 38 of 49

Page 39 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 40 of 49

Page 41 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 42 of 49

Page 43 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 44 of 49

Page 45 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 46 of 49

Page 47 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47





Page 48 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47


vas ( NO2 )

vs (CN)


DOS (a.u.)

Page 49 of 49

Reactive Molecular Dynamics Simulations to Investigate the Shock

Recommend Documents