Decomposition and Reaction of Polyvinyl Nitrate under Shock and

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Decomposition and Reaction of Polyvinyl Nitrate under Shock and Thermal Loading: A ReaxFF Reactive Molecular Dynamics Study Md Mahbubul Islam, and Alejandro Strachan J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b06154 • Publication Date (Web): 06 Sep 2017 Downloaded from http://pubs.acs.org on September 9, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry C 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 43

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

Decomposition and Reaction of Polyvinyl Nitrate under Shock and Thermal Loading: A ReaxFF Reactive Molecular Dynamics Study 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


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 43

ABSTRACT We use molecular dynamics (MD) simulations with the reactive force field ReaxFF to investigate the response of Polyvinyl Nitrate (PVN), a high-energy polymer, to shock loading using the Hugoniostat technique. We compare predictions from three widely used ReaxFF versions and, in all cases, we observe shock-induced, volume-increasing exothermic reactions following a short induction time for strong enough insults. The three models predict NO2 dissociation to be the first chemical relatively similar final product populations; however, we find significant differences in intermediate populations indicating different reaction mechanisms due to discrepancies in the relative stability of various intermediate fragments. A time-resolved spectral analysis of the reactive MD trajectories enables the first direct comparison of shock-induced chemistry between atomistic simulations and experiments; specifically, ultra-fast spectroscopy on laser shocked samples. The results from one of the ReaxFF versions are in excellent agreement with the experiments both in terms of threshold shock strength required for the disappearance of NO2 peaks and in the timescale associated with the process. Such direct comparison between physical observables is an important step towards a definite determination of detailed chemistry for high-energy density materials.


Page 3 of 43

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 Shock-induced chemistry of high-energy (HE) materials and the shock to detonation transition involve complex and inter-related physical and chemical processes through which the energy in the shockwave is transferred to chemical bonds with characteristic lengths of Angstroms and vibrational periods of few femtoseconds. In many cases the microstructure of the material plays a key role in localizing the energy in the shock in space, and these hotspots1,2 initiate chemical reactions, deflagration and eventually detonation. At the same time, inter- and intra-molecular relaxation processes transfer energy from long wavelength, low-frequency modes of the shock to high-frequency bond vibrations, a process called up-pumping.3–5 To complicate matters, these processes occur under extreme conditions of temperature, pressure and strain rate and involve short timescales. It is, thus, not surprising that despite significant efforts in experiments and simulations a definite a molecular-level description of shockinduced chemistry is still lacking even in relatively simple materials like nitromethane, a liquid at room temperature, or polyvinyl nitrate (PVN), an amorphous polymer. Recent experimental efforts are yielding insight and quantitative information about the molecular level processes responsible for shock-induced chemistry.6–10 For example, ultrafast spectroscopy study of laser shocked PVN enabled the detection of chemical reactive events under shock loading7, broadband multiplex vibrational sum-frequency generation (SFG) technique allows real-time observation of shockfront interactions with molecular monolayers.11 The experiment on PVN revealed that 18 GPa shocks resulted in chemical reactions within approximately 150 picoseconds indicated by the disappearance of a peak associated with NO2. Despite such significant accomplishments, these experimental techniques are not without limitations. For example, experiments do not enable a direct characterization of chemical mechanisms nor the identification of all intermediates and products. Current spatio-temporal and analytical resolution preclude in-situ evaluation of the detailed behavior of molecules behind a detonation front.8 On the simulation side, ab initio and quantum-based molecular dynamics (MD) simulations have been used to identify initiation mechanisms, reaction barriers, and rates associated with the decomposition of various energy materials12–19. However, the computationally intensity of these simulations limits the spatial and temporal scales achievable. Reactive force fields, like ReaxFF,20–22 provide a less computationally intensive alternative and enable multi-million atom simulations and capturing the role of microstructure and defects with scales approaching those in experiments.23,24 Recent efforts provide an atomic picture of the shock to deflagration transition.23 While these simulations provide unparalleled resolution in space and time, they involve approximations whose effect on predictions remain poorly



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 4 of 43

understood. To different degrees, ab initio, tight-binding, and reactive potentials approximate atomic interactions and forces; in addition, the use of classical ionic dynamics is also an approximation. We note that shock-induced chemistry is an extremely challenging problem for an atomistic model as it needs to capture two interdependent and complex processes: i) The thermo-mechanics of shock loading: the amount of energy the shock deposits into the system and how this energy is distributed among different degrees of freedom via variety of relaxation processes: volumetric compression and heating as well as plasticity, fracture, phase transformations, void collapse and interfacial sliding that can lead to energy localization and the formation of hot spots. ii) The model should also capture the thermodynamics and kinetics of chemical decomposition at various conditions of pressure and temperature; including uni- and multi-molecular processes. Given the complexity of the problem, a definite understanding of shock-induced chemistry will likely require a synergistic combination of experimental and theoretical investigations25 where the experiments validate the simulations and simulations help interpret the experimental results. In this paper, we use MD simulations with three different versions of the ReaxFF14,26,27 force field to establish the initial decomposition of PVN following shock loading. We build on the successful use of ReaxFF in a family of high-energy molecular crystals RDX, HMX, PETN, and NM13,16,18,26,28,29 and apply it to an energetic polymer. This is motivated by the existence of spectroscopic data on shocked PVN which enables the first direct comparison of ReaxFF predictions of shock-induced chemistry against experiments. Importantly, the amorphous nature of PVN simplifies the comparison between experiments and simulations due to the lack of grain boundaries and other crystal defects that can significantly affect shock response. In addition, PVN is of interest as a binder in explosive and propellant formulations and understanding its physics and chemistry can contribute to their safe operation. A high energy density binder provides not just adhesion but also contributes to the total energy release.6 Despite the interest of the PVN polymer in explosive applications, only a limited number experiments have been performed to investigate the mechanistic details of its response to shock loading. Moreover, no atomistic simulation studies on PVN chemical dissociation mechanisms under shock and thermal loading conditions are available in literature. As such the detailed chemistry of PVN remains unknown. Therefore, we aim to provide a detailed characterization of the mechanical, chemical, and spectroscopic properties of the PVN polymer. We compute shock Hugoniots and investigate how shock and thermal loading dictates the chemical reaction pathways. We find that PVN chemical dissociation initiates via the formation of NO2 and the activation energy decreases with the increased system densities. For the first time, our simulations enabled a one-on-one comparison of the reaction initiation timescales with the ultrafast shock experiments via time-resolved spectra. In addition, each of the force fields predicted results are compared


Page 5 of 43

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


with the available literature to establish their performance in predicting the mechanistic details of the PVN polymer. II. SIMULATION METHODOLOGY A. ReaxFF and Force Field Description ReaxFF is a bond order (BO)30,31 empirical reactive force field which allows the description of chemical reactions during MD simulations. It uses the concept of partial BOs between pairs of atoms to describe covalent interactions including 2-, 3-, and 4-body bonded interactions. BOs are a continuous, many body, and functions of the atomic coordinates. The non-bonded van der Waals and Coulomb interactions are calculated between every pair of atoms, within the respective cutoffs, regardless of covalent interactions.32,33 In the non-bonded energy expressions, shielding parameters are introduced to circumvent excessive repulsion at short distances, and a seventh order taper function is used to eliminate any energy discontinuity.34,35 The electronegativity equalization method (EEM)36 is utilized to obtained environment dependent partial atomic charges that are updated at every step during the MD simulations. For a more detailed description of the ReaxFF method, see Refs. [20–22,37,38]. Over the last decade, a number of ReaxFF parameter sets have been developed to describe highenergy density materials. The accurate prediction of the complex chemistry of these materials at extreme conditions of pressure and temperature as well as the description of physical, thermodynamic, and spectroscopic properties are challenging. Therefore, several versions are currently available, each emphasizing different properties. Importantly, the uncertainties associated with different versions have not been characterized, and comparisons across different force fields are not common. Therefore, it is an evident necessity to systematically evaluate the performance of the available force fields in terms of their predictions of mechanical, chemical and spectroscopic properties. In this study, we used three widely used ReaxFF C/H/N/O force fields to formulate a detailed chemical and mechanical description of the PVN polymer. The original C/H/N/O ReaxFF force field was developed and applied for RDX19,39 and since then several updates have been proposed. The standard ReaxFF C/H/N/O force field resulted in a less than ideal description of lattice parameters and equations of states due to the inadequate description of the London dispersion. Rom et al.26 incorporated additional inner-core repulsions for the van der Waals interactions to improve the description of lattice parameters and updated the previous C/H/N/O force field.39 We will denote this force field as ReaxFF-IW (inner-wall). The additional inner core repulsion has the following form

(

Evdw ( IW ) = c1.exp c2 (1 − rij / c3 )

) 5 ACS Paragon Plus Environment


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 6 of 43

Where c1, c2, and c3 are force field parameters, and rij is the interatomic distance between atom i and j. To improve the bulk properties of the energetic materials while, at the same time, maintaining the description of the chemical reactions of the unmodified ReaxFF– in the ReaxFF-lg force field, a low gradient (lg) attractive term was incorporated to account for the long-range London dispersion using following expression27

Clg,ij

N

Elg = − ∑

6 ij ,i < j ij

r + dReij6

Where rij and Reij are the distance and equilibrium vdW radius between atom i and j, respectively, and Clg,ij is the dispersion energy correction parameter. It has been shown that the addition of the lg-parameters on the ReaxFF-IW force field substantially improves the equation of states of energetic materials.27 We note that the difference between ReaxFF-IW and ReaxFF-lg force fields are only the low-gradient parameters. More recently, Wood et al.14 developed an updated version of the C/H/N/O force field specifically to improve the spectroscopic properties and the post-detonation combustion chemistry (denoted here as ReaxFF-2014). B. Structure Preparation Amorphous PVN (molecular formula (-CH(-ONO2)CH2-)n) structures were prepared using the Polymer Modeler tool available at nanoHUB.40 In order to investigate the effect of chain molecular weight (MW) on the polymer density, four different structures were constructed. The systems contain 135 chains each 8 monomers long (MW=714 g/mol), 22 50-monomer chains (MW=4,452 g/mol), 11 100monomer chains (MW=8,902 g/mol), and 5 200-monomer chains (MW=17,800 g/mol). All systems are built using a continuous configurational-bias Monte Carlo method as implemented in Polymer Modeler40 at an initial density of 0.5g/cm3. The structures are relaxed via energy minimization using the DREIDING non-reactive force field, followed by thermalization at temperatures between 300 K and 800 K via isobaric, isothermal MD (NPT ensemble) using a time step of 1fs. The annealing simulation steps are comprised of i) heating the system from 300K to 800 K in 500 ps ii) holding at 800 K for 2ns and iii) cooling down to 300 K in 500ps. Next, the system was further equilibrated at 300 K and atmospheric pressure for additional 10 ns using NPT MD to ensure molecular relaxation and to obtain the equilibrium density at ambient condition. A nonreactive force field was used to generate well-relaxed amorphous polymer structures for two main reasons. First, large forces due to bad contacts are common in the structures resulting from Monte Carlo builds, polymer modeler uses a series of relaxations with scaled van der Waals parameters to mitigate the effect of bad contacts and the use of a reactive force field during


Page 7 of 43

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


these initial relaxation steps could result in unwanted chemical reactions. Second, non-reactive reactive force fields are computationally less intensive than reactive ones, enabling longer relaxation procedures. The 300 K DRIEDING equilibrated structures were used as starting configurations for ReaxFF simulations. We re-equilibrated the structures using all three ReaxFF force fields with NPT simulation at room temperature and pressure for 300 ps, enough for the density to achieve steady state. Finally, we carried out isochoric, isothermal MD (NVT ensemble) at 300 K for 100 ps to prepare the initial configurations for the shock and isothermal decomposition simulations. We used a time step of 0.2 fs in the ReaxFF relaxation simulations, the time step is reduced to 0.1 fs in all shock and thermal decomposition simulations. A snapshot of the simulation cell with 200 monomers and the densities predicted by the ReaxFF-2014, ReaxFF-IW, and ReaxFF-lg for all the polymer structures considered are shown in Fig. 1 together with the experimental value of 1.34 g/cm3 measured by gas pychnometry.41 One can see that the densities are converged with respect to molecular weight for values larger than 4,450 g/mol. The densities predicted by ReaxFF-2014 and ReaxFF-IW compare well with the experimental value41 of 1.34 gcm-3, while ReaxFF-lg overestimates the density. In this study, the bulk PVN geometry consisting of 200-monomer chains are used for all the shock and isothermal decomposition simulations.

Fig 1. (a) A snapshot of the simulation cell with 200 monomers; insert highlights a section of the polymer chain (Color code: Oxygen: red, Carbon: gray, Nitrogen: blue, and Hydrogen: white). (b) The calculated densities as a function of molecular weight of the PVN polymer chains. ReaxFF-2014, ReaxFF-IW, and ReaxFF-lg are red, green and blue, respectively.

C. Hugoniostat Simulations 7 ACS Paragon Plus Environment


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 8 of 43

Non-equilibrium MD simulations are often used for shock simulations; however, large system sizes (typically millions of atoms) are required for these simulations due to the fast propagation velocity of shock waves and the need to achieve steady shock front. Therefore, significant computational resources are spent in solving the equations of motion for the atoms which are not participating in the shock chemistry. Alternatively, equilibrium MD simulations have been proposed to simulate the final state of the system after the passage of shock, such as Hugoniostat42 and MSST.43 The Hugoniostat MD simulation technique, compresses (uniaxially) and heats up the system to a state corresponding to the desired shock by satisfying the Hugoniot jump conditions. These simulations are less intensive computationally as the actual shock wave propagation through the samples is circumvented. Thus, smaller system sizes can be simulated and longer time-scale simulations are possible with moderate computational resources, resulting in better statistics of the chemical decomposition reactions. The jump conditions imposed by the Hugoniostat are:

mass : ρ ous = ρ ( us − u p ) momentum : Pzz = P o + ρo us u p 1 energy : EH − Eo = ( Pzz + Po )(Vo − V ) 2 P −P us = zz o Vo Vo − V u p = ( Pzz − Po )(Vo − V ) Where subscript 0 refers to the unshocked initial properties, us and up are shock and particle velocities, respectively. Pzz is the normal component of the stress tensor in the shock direction and ρ=1/V. Given the desired shock pressure the Hugoniostat technique uses a thermostat and a barostat to satisfy the jump conditions. We performed Hugoniot shock simulations in the pressure range of 4 to 36 GPa. In order to investigate the influence of initial densities of the systems, we carried out two sets of simulations with the densities predicted by the respective force fields at ambient conditions (denoted ‘FF density’) and the experimental density (denoted ‘experimental density’). We note that for the case of ReaxFF-2014 force field both the densities are identical. D. Isothermal Decomposition Simulations The 300 K NVT equilibrated structures with MW=8,902 g/mol were used as starting configurations for studying the homogeneous isothermal decomposition of PVN at three different temperatures: 2400 K, 2800 K, and 3200 K. Three different densities were considered – experimental 8 ACS Paragon Plus Environment

Page 9 of 43

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


density (ρo), 1.2ρo, and 1.4ρo. The simulation cell with the experimental density was used to generate the volumetrically compressed systems by adjusting volumes to the target densities using deformation simulations. The compressed systems were further relaxed using 300 K NVT simulations for 10 ps. Next, the target temperatures were reached by rescaling the atomic velocities at a heating rate of 1000 K/ps using NVT ensemble. The Nose-Hoover thermostat was used with a coupling constant of 20 fs. The rapid heating rate is comparable to the timescales of shock loading. Once the target temperatures were reached, the systems were held at the temperature with an NVT ensemble in the range of 300-800 ps. The lower temperature cases required longer simulation time to attain steady-state. The isothermal simulations allowed us to extract Arrhenius parameters to study reaction kinetics of the PVN polymer.

The

characteristics time for the reaction rate calculations is considered as the time required to reach the half of the total exothermicity in each case. E. Analysis of MD Simulations The MD simulation trajectories were post-processed to obtain power spectra using the following formula:14

P (ω ) =

βτ N

2

N −1

3N

∑m ∑v ( n∆t ) e j

j =1

− i 2πω n∆t

j

n =0

where vj(t) are the atomic velocities at time t, β is defined as (kBT)−1, T is the temperature, τ is the sampling period, mj is the mass of atom j, ∆t is the sampling rate, 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 at every 2.5 fs to resolve the high-frequency modes accurately. The time-resolved spectra calculations were carried out at 10 ps interval, and each case was averaged over 4000 frames. The chemical bonds and molecular species were detected using a conservative bond order cutoff of 0.30 for all the simulations. The accuracy of on-the-fly species recognition was improved by sampling bond order values at every 0.50 fs and averaged over 10 samples. This approach instead of using instantaneous bond orders also abate identification of spurious bonds in high-density systems, typical for shocked compression. The molecular fractions represented in this study are normalized by the total number of monomers.

III. SHOCK STATES AND EXOTHERMIC REACTIONS A. Evolution of Pressure, Volume, and Temperature



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 10 of 43

The pressure, volume and temperature evolution of shocked PVN polymer starting at the experimental density are shown in Fig. 2. The evolution of all the thermodynamic properties of the ambient pressure simulations is presented in Fig. S1 of Supporting Information (SI). During the Hugoniostat shock simulations, system pressure, volume, and temperature equilibrate to their shocked state within 1 ps timescale for all three force fields. The high deformation rate is typical for shock conditions. This initial fast response is driven by the uniaxial compression in response to the loading. Interestingly, for strong enough shocks the systems exhibit a sudden increase in volume and temperature following the initial compression and an induction time. This is due to the exothermic, volume-increasing reactions induced by a combination of temperature and density increase due to the shock. During this phase, small molecules are produced (expansion stage). The simulations show a decrease in the induction time and reaction-induced volume expansion with increasing shock strength (and pressure). Interestingly, as can be seen in Fig. S1, ReaxFF-lg force field does not exhibit exothermic reactions when shocked from its ambient pressure density for the shock strengths up to 36 GPa, and timescales of 300 ps.


Page 11 of 43

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


Fig 2: Pressure, volume, and temperature evolution during shock loading of PVN samples starting at the experimental density systems using ReaxFF-2014, ReaxFF-IW, and ReaxFF-lg force fields. The shock pressure-density data extracted from the shock simulations for the three force fields and two types of initial conditions are presented in Fig. 3; Fig. 4 shows the corresponding temperature-density plots. The shock states achieved from the experimental densities, Figs. 3(a) and 4(a), exhibit two branches: the shock states follow the unreactive Hugoniot for relatively weak loading and transition to the reactive Hugoniot as reactions occur. In order to show the transition from non-reactive to reactive Hugoniot we included a sample case in Fig S2. We find that ReaxFF-2014 predicts rapid chemistry for shocks of 18 GPa and higher pressures. The reaction involves a reduction in density of approximately 10% and an increase in temperature of approximately 2000 K. The induction time required for exothermic reactions at 18 GPa is ~180 ps, see Figure 2. The critical shock strength and reaction timescales predicted by ReaxFF-2014 are in excellent agreement with experimental results of McGrane et al.41 Their ultrafast spectroscopy experiments show no reactions for shock pressures lower than 17 GPa and reactions within 150 ps for a shock strength of 17 GPa. ReaxFF-IW and ReaxFF-lg are slightly less reactive and require 24 GPa shocks to exhibit reaction within a timescale of 300 ps. It is evident from Figs. 3(a) and 4(a) that the nonreactive P-ρ Hugoniots and T- ρ obtained starting from the experimental density are very similar for all the force fields. Interestingly, the more sensitive ReaxFF-2014 exhibits a higher exothermicity than the other two force fields. Shock simulations starting from the ambient pressure conditions for each force field exhibit larger discrepancies among each other; see Figs 3(b) and 4(b). The ReaxFF-2014 and -IW force fields predicted volume-increasing reactions and pressures of 18GPa and 30 GPa respectively, while no significant reactions occurred in the ReaxFF-lg simulations for shocks up to a pressure of 36 GPa. The lower reactivity in these simulations can be attributed to the smaller amount of mechanical (PV) work done on the system during shock compression. The equations of state in Fig. 3 (b) show a significant difference in density for low pressures among the three force fields, but this difference reduces with compression. Thus, the overestimation of density in the ReaxFF-IW and -lg force fields resulted in a smaller amount of PV work and heating of the system as compared with the ReaxFF-2014 simulations. As such the first rapid chemistry observed in the ReaxFF-IW force field simulations is at a higher pressure of 30 GPa compared to the 18GPa in the case of ReaxFF-2014 force field.



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 12 of 43

Fig. 3: Pressure vs. density plots using three different ReaxFF force fields for PVN samples shocked starting from (a) experimental density, (b) force field densities. The non-reactive densities are not shown for the cases where volume-increasing reactions have occurred. Color code: Red: ReaxFF-2014, Green: ReaxFF-IW, and Blue: ReaxFF-lg.

Fig. 4: Temperature vs. density plots using three different ReaxFF force fields: (a) experimental density, (b) force field densities. The non-reactive temperatures are not shown for the cases where volume increasing reactions have occurred. Color code: Red: ReaxFF-2014, Green: ReaxFF-IW, and Blue: ReaxFF-lg.

B. Shock Velocity-Particle Velocity Hugoniots The shock and particle velocity was calculated using Hugoniot relations (Section IIC). The shockvs. particle velocity data for the experimental density simulations are shown in Fig. 5(a) together with the 12 ACS Paragon Plus Environment

Page 13 of 43

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


non-reactive experimental curve from Ref. [10]. The exothermic reactions lead to a second branch of the us-up plots and a sudden increase in the shock velocities. However, the slope of the reactive us-up curve is lower, and upon increase in shock pressure, the reactive and extrapolation of the unreactive us-up curves approach each other. This behavior was observed in nitromethane shock experiments.7 The predicted us-up curves slightly underestimate the experimental results, but overall they exhibit reasonable agreement. The shift in the us-up plots is indicative to the chemical reactions. While the nonreactive us-up data are analogous for all the force fields, the differences in the us-up plots at the reactive zone is due to the fact that ReaxFF-2014 predict chemical reactions at 18 GPa in contrast to 24 GPa for the other two force fields. However, at the pressure of 24 GPa or higher, us-up values from all the force fields are in close agreement. The ReaxFF-2014 predicts the transition to fast chemical decomposition at up=2.2 km/s and us=6.5 km/s, while in ReaxFF-IW and –lg simulations it occurs at approximately up=2.7 km/s and us=6.8 km/s. The sound speed (value of us at up=0) predicted by the ReaxFF-2014, ReaxFF-IW, and ReaxFF-lg force fields starting from the experimental densities are 1.75±0.05, 1.55±0.05, and 1.05±0.15 km/s, respectively. These values are smaller compared with the experimental data10 of 3.2±0.3 km/s. However, the FF density simulations predict a relatively higher sound speed values, that is, 2.16±0.02 and 3.14±0.04 km/s for ReaxFF-IW and ReaxFF-lg, respectively.

Fig. 5. The us-up plots for the PVN at various shock pressures using three different ReaxFF force fields and comparison with the experimental data: (a) Experimental density, (b) FF densities. Color code: Red: ReaxFF-2014, Green: ReaxFF-IW, and Blue: ReaxFF-lg, Black: experimental data.

The effect of initial systems densities is evident in the results presented in Fig. 5(b). In the FF density simulations, different force fields exhibit a wider variation. The computed shock-particle velocities follow the trend of densities predicted by each of the force fields. The higher density resulted in a stronger shock velocities. Therefore, the highest shock velocities are predicted by the ReaxFF-lg force 13 ACS Paragon Plus Environment


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 14 of 43

field, consistent with the highest density predicted by this force field. The two regimes in the us-up plots have been observed for ReaxFF-2014 and ReaxFF-IW force fields. The ReaxFF-IW force field simulations, the volume-increasing reactions occurred at up=3.0 km/s and us=7.5 km/s. Since there is no experimental or computational data available for the reactive regimes of the shocked PVN polymer, we cannot make any comparison. C. Shock-Induced Chemical Decomposition The processes of shock-induced decomposition of PVN can be categorized into three steps: i) initial non-reactive compression and heating when energy transfers from the shock wave to the material, ii) a shock-induced onset of chemical reactions iii) followed by rapid exothermic reactions. Therefore, investigation of chemical decomposition mechanisms is required to elucidate the coupling between mechanical loading and chemistry. The decomposition pathways are contingent upon the nature of the applied load, such as temperature and pressure conditions. With increasing shock strength, a larger amount of energy will be transferred to vibrational degrees of freedom of the molecules. If this results in sufficiently high temperature, chemical reactions will be induced within short timescales.44 In this study, we are primarily interested in probing the initial chemical reactions and subsequent pathways towards stable species formation.


Page 15 of 43

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


Fig 6. Molecular population distribution of PVN polymer at shock pressures of 18 and 22 GPa using ReaxFF-2014 force field. Molecular fractions are normalized by the total number of monomers. Possible reaction initiation mechanisms of PVN are NO2 fission, HONO elimination, C-C bond breaking, and hydrogen migration. These paths have been observed as the initiation of reactions in related compounds such as RDX,23,24 HMX,29 CL-20,45 and PETN.28 During the initial compression stage of 12ps, all systems were unreacted. As reaction initiates, the homolytic cleavage of the O-NO2 bond was observed as the first reactive event for all three force fields. The time evolution of the major intermediates and stable species are shown in Figures 6, 7, and 8 as determined using ReaxFF-2014, ReaxFF-IW, and ReaxFF-lg force fields, respectively. For clarity, only species with populations greater than 2% of the total number of monomers are shown. For each force field, we show results for the weakest shock at which we observe chemistry and a stronger shock to highlight how pressure affects chemical decomposition. The rapid exothermic chemistry leads to the generation of final stable gaseous products, most notably H2O, CO2, N2, and NH3 for all the force fields. As expected, ReaxFF-IW and ReaxFF-lg that only differ in their description of van der Waals interactions result in very similar chemical reactions. The final population of products is similar in all three force fields; the main difference being the higher population of CO2 in ReaxFF-2014. In addition, we observe a significant amount of H2 formation in the ReaxFF-IW and -lg simulations. H2O was found as the most abundant species irrespective of simulation conditions and force field used. Unfortunately, we are not aware of any experimental data on the distribution of final products of shocked PVN available for the comparison.



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 16 of 43

Fig. 7. Molecular population distribution of PVN polymer at shock pressures of 24, and 32 GPa using ReaxFF-IW force field (experimental density). The evolution of intermediate populations can shed light on the reaction mechanisms and also explain differences in the relative amount of final products predicted by each of the force fields. We observe that the primary intermediate molecules reach a maximum and then react away to form final, exothermic products. The reaction initiation yields a significant amount of NO2 and creates radical sites, followed by the generation of other intermediates. In ReaxFF-2014, at 18 GPa shock pressure, the NO population increases steadily up to ~180 ps, then suddenly drops at the onset of rapid exothermic chemistry. The population of the most dominant intermediate, HCO2, increases to 25% of the total number of monomers. CH2O is generated at a relatively later stage. However, the decline of the both HCO2 and CH2O occurs concurrently at the transition-to-fast exothermic reactions. ONH radicals are observed at timescales similar to CH2O, but in a lower quantity. The formation of the radicals proceeded via breaking of the backbone C-C bond of the PVN polymer. The rapid chemistry leads to the generation of CNOH, and N2O2 intermediates and stable gases are formed at a pronounced rate. At the onset point of rapid chemistry, ~80% of the total C-C bonds of the PVN chains are dissociated. In contrast to ReaxFF2014, a large quantity of HNO3, and HONO are produced in the ReaxFF-IW and -lg simulations. HONO elimination has been observed in the previous experiment,46 QM,15,16 and ReaxFF17,29,39 studies on energetic materials. 16 ACS Paragon Plus Environment

Page 17 of 43

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


Fig. 8. Molecular population distribution of PVN polymer at shock pressures of 24, and 32 GPa using ReaxFF-lg force field (experimental density). To better understand the population of products and intermediate evolution for the three descriptions, we computed their formation energy using the three force fields and compared them to experiments in Figure 9. ReaxFF-2014 over stabilizes HCO2, CH2O, and CNOH radicals, compared to the other two force fields and experimental data. This stabilization leads to the production of these species in larger quantities, while they are almost absent in other two versions. Similarly, the high stability of HNO3 predicted by ReaxFF-IW and -lg force fields causes it to be one of the most dominant intermediate species (see Fig. 7 and 8). The formation of HONO is endothermic in ReaxFF-2014, as such this pathway has rarely been observed in ReaxFF-2014 simulations. Thus, our simulations show that the stability of intermediates and reaction paths have a direct effect on the final products, at least during the initial hundreds of picoseconds during decomposition. The major difference in the population of final products across various force fields is the amount of CO2 and NH3. ReaxFF-2014 correctly stabilizes CO2 with respect to the other force fields. The higher preponderance of CO2 production in ReaxFF-2014 simulations can also be attributed to the dissociation of the HCO2 radical pathway. The distribution of the



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 18 of 43

final products shows that ~20% of C atoms in ReaxFF-2014 are contained in small product molecules; this number if only 5% in the ReaxFF-lg and ReaxFF-IW force fields. This indicates that a significant amount of carbon remains as intermediate species or other final products. We further investigated the other carbon containing species, and they can be presented under the general formulas of HCxOy (x=1-4, y=2-5), H2CxOy (x=1-4, y=2-6), and H3CxOyN (x=1-2, y=2-3). A list of these representative species are included in the Table S1. We note that in our simulations we have not observed any carbon cluster formation.

Fig. 9. Formation energies of key intermediate and stable species as calculated using three different ReaxFF force fields and comparison with the experimental data [47]. In these calculations, reference states are considered as follows: (i) atomic energies are used for the formation energies of N2, O2, and H2, and (ii) N2, H2, O2, and graphite are used as references for all other cases. One can also see that the distribution of both intermediates and stable species indicates that the reaction initiation mechanisms are relatively insensitive to the magnitude of the applied shock. Overall, the kinetics of the exothermic reactions are relatively slower, and the exothermic heat release is lower in ReaxFF-IW and -lg simulations compared to ReaxFF-2014. D. Spectral Changes During Shock Simulations Although detailed reaction mechanisms are inaccessible in shock experiments, ultrafast spectroscopy provides important insight into some of the reactions under shock loading. Importantly, 18 ACS Paragon Plus Environment

Page 19 of 43

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


spectroscopy information can be extracted directly from MD simulations as discussed in Section II and this can provide a direct comparison between simulation results and experimental measurements. We aim to compare reactive MD simulations with published and ongoing spectral analysis for a direct validation to our simulations and, at the same time, to help interpret spectral data to detect various reaction pathways and associated timescales. Thus, we performed spectral analysis using atomistic velocities at various shock conditions and for the three force fields. The predicted per element spectra of a single polymer chain using ReaxFF-2014 is shown in the supplementary material, Fig. S3 (a). This enables the assignment of peaks to vibrational modes. In order to compare with the laser shock experiments by McGrane et al.41, we considered changes in the N-O stretching modes as a measure of chemical reactions in the system. The ReaxFF-2014 predicts NO2 symmetric and asymmetric stretching modes are at 800 and 1900 cm-1, respectively. These values are shifted from the experimental values of 1270 and 1624 cm-1.41 However, the discrepancy between experimental and our calculated vibrational modes does not significantly influence our comparison with experimental observations regarding the timescale or character of reaction initiation. The C-H stretching mode is in good agreement with literature data.

Fig. 10. (a) ReaxFF-2014 predicted time-resolved full spectra at 18 GPa shock simulation, (b) Expanded view of the disappearance of NO2 vibrational mode at a timescale of ~150ps. (c) time-dependent spectra as computed from a ReaxFF-lg 24 GPa experimental density shock simulation.



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 20 of 43

During shock simulations, time-resolved spectral data provides evidence of dissociation of various chemical groups, hence timescales of the initiation of related reactions pathways. In our 18 GPa shock simulation, using ReaxFF-2014, the νas(NO2) peak disappeared at a timescale of ~150 ps, the point at which rapid chemistry initiates, Figs. 10a-b. The calculated timescale of the NO2 peak disappearance is in excellent agreement with the experimental data McGrane et al.41. The calculated spectra also elucidate other reaction pathways. For example, the distribution of the intermediate species at 18 GPa in Fig. (6) depicts the production of a significant amount of NO starting at ~50 ps, which reaches a maximum and completely disappeared at ~200 ps. The peak appeared in our calculated spectra at ~1600 cm-1 nicely captures the entire NO evolution history. The timescale of the presence of the peak matches that of NO production. Thus our simulation predicted spectra enabled elucidation of various reaction pathways, at the same time, the timescale of the reactions. The rapid chemistry causes the disappearance of the major peaks, as the final gaseous products are produced. As a result, low-frequency modes denoting diffusive processes start to build up. The broad nature of the peaks with frequencies near 3000 cm-1 stems from the combination of O-H, N-H, and C-H stretching modes due to the formation of stable gaseous molecules. The spectra of a single PVN chain also calculated using ReaxFF-lg force field as shown in the Fig S3(b). It can be seen that the ReaxFF-lg predicted NO2 stretching mode is significantly blue shifted to ~2400 cm1

. During a 24 GPa shock simulation at the experimental density using ReaxFF-lg, the time evolution of

the spectra (as shown in Fig. 10(c) indicates the disappearance of the NO2 vibrational mode at a timescale of ~100 ps. However, the evolution of the key peaks is not as pronounced as in the case of ReaxFF-2014 simulations.

IV. THERMAL-INDUCED CHEMICAL DECOMPOSITION As mentioned above, shocked-induced chemical reactions depend both on the thermo-mechanics involved in shock loading and the decomposition at the shock temperature and pressure. To decouple these two effects and compare the description of the chemistry of the three force fields, we performed additional simulations where PVN reaction and decomposition is studied under isothermal, isochoric conditions. A. Decomposition Reaction Rate Constants The isothermal decomposition simulations are performed following the methodology as described in Section IIC. The time evolution of potential energies at two different densities (ρo and 1.4ρo) are shown in Fig. S4 for all the three force fields. The first step in the analysis of the reaction rate constants is finding the characteristic time scales of reactions. We define it as the time required for the potential 20 ACS Paragon Plus Environment

Page 21 of 43

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


energy to reach half of the total exothermicity at a particular temperature and density. The inverse of the characteristic times, the reaction rates, are shown in Fig. 11 (a) as a function of temperature in an Arrhenius plot for all the densities and force fields. We observe that both temperature and compression of the structure plays a major role in the reaction rates. The increase in both temperature and density increases the reaction rates. For example, the reaction rate increases by factors of 60, 81, and 79 for the experimental density when system temperature is increased from 2400 to 3200 K, using ReaxFF-2014, ReaxFF-IW and ReaxFF-lg, respectively. The effect of compression of the bulk structure was also very pronounced. For example, at 2400 K, compressing to 1.4ρο leads to 32, 49, and 51 fold increase in the reaction rates with ReaxFF-2014, ReaxFF-IW and ReaxFF-lg, respectively.

Fig. 11. (a) Arrhenius plot for PVN dissociation at isothermal conditions (solid line: experimental density (ρο) dash; 1.2ρο and dotted; 1.4 ρο). (b) Overall activation energies associated with decomposition at different densities and for the three force fields. The slopes of the exponential regression analysis of Fig. 11(a) can be associated with an overall activation energy for decomposition at isochoric conditions; see Fig. 11(b). Similar approach has been used in previous studies to calculate activation energies.12,26,39 Consistent with the shock results, we find ReaxFF-2014 to exhibit lower activation energy than the other two force fields. The increased reaction rates with the volumetric compression are also evident from the decrease in activation energies as shown in Fig. 11(b). The ReaxFF-2014 force field calculated global activation energies are 16.0 15.9, and 15.5 kcal/mol for experimental (ρο), 1.2ρο, and 1.4ρο, respectively. The ReaxFF-2014 predicted activation energies of the PVN are significantly lower than the Wood et al.48 reported data for NM, HMX, and PETN. This indicates a higher reactivity of the PVN compared to these high-energy materials. In contrast to ReaxFF-2014, ReaxFF-IW force field predicts around 40% higher activation energies for all the cases. 21 ACS Paragon Plus Environment


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 22 of 43

We note that ReaxFF-lg calculated activation energies are slightly lower than the ReaxFF-IW at all the densities. B. Detailed Mechanism of Dissociation Chemistry Thermal decomposition of the PVN polymer at the early stages results in the formation of a large variety of radical intermediates. Subsequently, the radicals undergo chain-like secondary reactions to generate stable final products. The evolution of the species of experimental density samples at 3200 K is presented in Fig. 12. The results from the 2400 K experimental density and 2400K and 3200 K, 1.4ρο density simulations are included in the Fig S5-S7. The results presented in Fig. 12 reveals that a major initiation pathway is the production of NO2 via the dissociation of O-NO2 bond linked to the (-C-C-) backbone of the polymer; this is similar to the shock-induced decomposition. The NO2 peak reaches a maximum value at the very early stage of decomposition and then rapidly dwindles as the secondary reaction pathways consume this species. NO2 was found as a predominant initial reaction product in the thermal decomposition of other nitramines such as RDX,49 HMX,29 and CL-20.45 The NO2 evolution mechanisms are similar for all the force fields. However, the production of other intermediates varies significantly with the force fields. The ReaxFF-2014 generates a larger quantity of CH2O and HCO2 along with other radical intermediates such as NO, ONH, CNOH, and a trace of HONO. In contrast, in ReaxFFIW and –lg simulations, a significant amount of HNO3 and HONO has been observed. The secondary dissociation reactions of the intermediate radicals lead to the formation of CO2, H2O, N2, NH3, H2, and CO. The major secondary reactions that produce CO2 and H2O are given in the Table S2 and S3. The CO2 formation pathways tracked from the ReaxFF-2014 simulations exhibit that about 50% of the CO2 are evolved from the HCO2. The dissociation of the larger population of the HCO2, causes CO2 to be the most abundant species. However, we have not observed this CO2 reaction pathway in the ReaxFF-lg and ReaxFF-IW force field simulations. Across the temperature and pressure range considered in this study, H2O is the most dominant species in the case of ReaxFF-IW and –lg force fields. We observe differences in the proportion of reaction initiation pathways in between shock and isothermal simulations. The concentration of the NO2 produced (Fig. 12) in isothermal simulations is much higher than the shock simulations. To further illustrate the variation in the NO2 production with the simulation conditions, we calculated the O-NO2 bond breaking barrier at two different densities of the PVN. ReaxFF-2014 force field predicts the barriers as 43 and 60 kcal/mol for the densities of 1.40 and 2.40 gcm-3, respectively. We used a bond-restraint method50–52 to calculate the minimum energy pathway to dissociate O-N bond. Thus, with increasing density, the N-O bond breaking barrier increases that result in a smaller amount of NO2 production in shock simulations. The higher density caused by the shock loading reduces the distances among the neighboring monomer units, as such NO2 scission become less 22 ACS Paragon Plus Environment

Page 23 of 43

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


favorable. This observation is consistent with the reports on the decomposition of various nitroaminos explosives under shock and thermal loading.12,26,45,53 The reduction of NO2 scission due to the higher shock velocity or pressure has been observed in CL-20,45 RDX,19 and HMX.53 Moreover, a FTIR study on the thermal decomposition of RDX and HMX also corroborate our observation on the effect of pressure. It was found that an increase in system pressure reduces NO2 scission.49 To further understand the effect of shock vs. thermal loading, we performed isothermal simulations at equivalent thermodynamics conditions such as temperature, and density of the final state of the shocked material. The comparison of the species evolution in the cases of equivalent shock and thermal loading conditions are shown in Fig. S8-S10 and a summary of the distribution of the final gaseous products are presented in Fig. S11. The details of the thermodynamic conditions are given in the captions of Fig. S8-S10. From these figures, it is evident that the nature of the loading has a significant influence on the evolution of the intermediate species, as such reaction pathways. We observe a pronounced effect of mechano-chemistry due to the shock loading on the initial reaction pathways. The shock loading resulted in a decrease in the amount of the intermediates—those are evolved in a reasonable quantity—compared to the thermal loading for all the force fields. This implies that the PVN dissociates into many smaller fragments, the individual quantities of the fragments were small thus we ignored them in the population analysis. However, despite the variation in the amount of the intermediates, the distribution of final gaseous products is almost identical for both the simulation cases (see Fig. S11).



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 24 of 43

Fig. 12. Species evolution during isothermal decomposition simulations at 3200 K using three different force fields at experimental density. The molecular fractions are normalized by the total monomer units. We observe that higher compression and temperature alters the reaction pathways and final products distribution in the isothermal simulations. The increase in density reduces the amount of NO2 production for all the cases. In ReaxFF-2014 simulations, the amount of CHNO radicals produced increases with the system pressure. At 3200 K simulations, ReaxFF-2014 resulted in a reduction of CH2O and HCO2 and a substantial increase in NO and H2O population. At 2400 K, irrespective of the system densities, the amount of CO2 produced is notably higher than the other stable products; however, with the increasing temperature and density such as at 3200 K and 1.4ρo (see Fig. S6), this dominance dwindles, and CO2 and H2O are generated in a comparable quantity. This is due to the reduction of CO2 and an increase of H2O at higher temperatures. In contrary, higher temperature decreases the H2O population in 24 ACS Paragon Plus Environment

Page 25 of 43

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


both ReaxFF-IW and -lg simulations; however, the effect is less discernible at higher densities. At experimental density ReaxFF-IW and -lg simulations, CO2, H2O, and H2 reach at similar asymptotic quantities at 3200 K as a result of the reduced amount of H2O population. The quantity of the H2 production exhibits an increasing trend with the increasing system temperature. We note that H2 formation pathways are absent in ReaxFF-2014 force field simulations. One can also observe that the amount of N2 produced is insensitive to the change in temperature and densities. V. CONCLUSIONS In conclusion, the first atomistic study of the chemical decomposition of PVN polymer under shock and thermal loading conditions provides a detailed description of mechanical, chemical, and spectroscopic properties. The predicted shock vs. particle velocity Hugoniots using ReaxFF are in reasonable agreement with available experimental data. More importantly, this study enabled the first direct comparison between reactive MD simulations and experiments for the decomposition of a material under shock loading. We found that ReaxFF-2014 predicts the onset of chemical reactions and the timescale associated with the formation of NO2 in excellent agreement with shock experiments. The ReaxFF-IW and –lg versions required stronger shocks for initiation. The combination of reactive MD simulations and ultrafast spectroscopy opens the possibility of both the validation of ReaxFF at extreme conditions and can contribute to the interpretation of the experimental data relating changes in spectral features to atomic processes. Unfortunately, the narrow spectral region of existing experiments only enables a partial validation of the predictions. Broad spectral analysis of laser-shocks, currently under development, will enable a more complete validation of atomistic simulations and a detailed understanding of chemical decomposition paths. Our investigation of the reaction initiation pathways reveals a significant dependency on the loading type. NO2 formation is predominantly found as a reaction initiation step, however, the relative amount of the species is contingent on the nature of the applied load. The shock and volumetric compression results in the reduction of NO2 production, consistent with observations in other high energy explosives. The final product populations predicted by the three force fields are similar, the main difference being a higher population of CO2 in ReaxFF-2014. This can be attributed to the, correct, relative stability of this molecule, see Figure 9, but also by the overstabilization of the HCO2 intermediate by ReaxFF-2014. While the comparison experimental spectroscopic studies with the MD simulations indicates that ReaxFF-2014 provides a more accurate description of threshold for chemical initiation and its timescales, the lack of bradband data does not enable a comprehensive validation of the three force fields and the detailed chemistry they predict. We expect such comparisons to be possible in the near future. 25 ACS Paragon Plus Environment


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 26 of 43

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 temporal evolution of thermodynamics quantities during FF density shock simulations, per element spectra of a single polymer chain, evolution of potential energies and species in isothermal decompositions. 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. REFERENCES (1) (2)

(3) (4)

(5) (6) (7)

(8)

Tarver, C. M.; Chidester, S. K.; Nichols, A. L. Critical Conditions for Impact- and Shock-Induced Hot Spots in Solid Explosives. J. Phys. Chem. 1996, 100, 5794–5799. McGuire, R. R.; Tarver, C. M. Chemical-Decomposition Models for the Thermal Explosion of Confined Hmx, Tatb, Rdx, and Tnt Explosives; UCRL-84986; CONF-810602-15; Lawrence Livermore National Lab., CA (USA), 1981. Dlott, D. D. Shocked Molecular Solids: Vibrational up Pumping, Defect Hot Spot Formation, and the Onset of Chemistry. J. Chem. Phys. 1990, 92, 3798–3812. Wen, X.; Tolbert, W. A.; Dlott, D. D. Multiphonon up-Pumping and Molecular Hot Spots in Superheated Polymers Studied by Ultrafast Optical Calorimetry. Chem. Phys. Lett. 1992, 192, 315–320. Dlott, D. D. Theory of Ultrahot Molecular Solids: Vibrational Cooling and Shock‐induced Multiphonon up Pumping in Crystalline Naphthalene. J. Chem. Phys. 1990, 93, 1695–1709. Moore, D. S.; McGrane, S. D. Comparative Infrared and Raman Spectroscopy of Energetic Polymers. J. Mol. Struct. 2003, 661–662, 561–566. 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. Dlott, D. D. New Developments in the Physical Chemistry of Shock Compression. Annu. Rev. Phys. Chem. 2011, 62, 575–597. 26 ACS Paragon Plus Environment

Page 27 of 43

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


(9)

(10) (11) (12)

(13)

(14)

(15) (16) (17) (18) (19)

(20) (21) (22)

(23)

(24)

(25) (26)

(27)

(28)

Bassette, W. P.; Dlott, D. D. High Dynamic Range Emission Measurements of Shocked Energetic Materials: Octahydro-1,3,5,7-Tetranitro-1,3,5,7-Tetrazocine (HMX). J. Appl. Phys. 2016, 119, 225103. Moore, D. S.; McGrane, S. D.; Funk, D. J. Ultrafast Spectroscopic Investigation of Shock Compressed Energetic Polymer Films. AIP Conf. Proc. 2004, 706, 1285–1288. Patterson, J. E.; Lagutchev, A.; Huang, W.; Dlott, D. D. Ultrafast Dynamics of Shock Compression of Molecular Monolayers. Phys. Rev. Lett. 2005, 94, 015501. 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. Shan, T.-R.; Wixom, R. R.; Mattsson, A. E.; Thompson, A. P. Atomistic Simulation of Orientation Dependence in Shock-Induced Initiation of Pentaerythritol Tetranitrate. J. Phys. Chem. B 2013, 117, 928–936. 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. Schweigert, I. V. Quantum Mechanical Simulations of Condensed-Phase Decomposition Dynamics in Molten RDX. J. Phys. Conf. Ser. 2014, 500, 052039. 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. Joshi, K. L.; Chaudhuri, S. Reactive Simulation of the Chemistry behind the Condensed-Phase Ignition of RDX from Hot Spots. Phys. Chem. Chem. Phys. 2015, 17, 18790–18801. Manna, M. R.; Reed, E. J. Early Chemistry in Hot and Dense Nitromethane: Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 10146–10153. 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. 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. 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. 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. 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. Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C. T.; Goddard, 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. Li, Y.; Vashishta, P. Multistage Reaction Pathways in Detonating High Explosives. Appl. Phys. Lett. 2014, 105, 204103. Rom, N.; Zybin, S. V.; van Duin, A. C. T.; Goddard, W. A.; Zeiri, Y.; Katz, G.; Kosloff, R. Density-Dependent Liquid Nitromethane Decomposition: Molecular Dynamics Simulations Based on ReaxFF. J. Phys. Chem. A 2011, 115, 10181–10202. 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. Budzien, J.; Thompson, A. P.; Zybin, S. V. Reactive Molecular Dynamics Simulations of Shock Through a Single Crystal of Pentaerythritol Tetranitrate. J. Phys. Chem. B 2009, 113, 13142– 13151. 27 ACS Paragon Plus Environment


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

(29)

(30) (31) (32) (33)

(34)

(35)

(36) (37)

(38)

(39) (40) (41) (42) (43) (44) (45) (46)

(47) (48) (49)

Page 28 of 43

Zhou, T.; Song, H.; Liu, Y.; Huang, F. Shock Initiated Thermal and Chemical Responses of HMX Crystal from ReaxFF Molecular Dynamics Simulation. Phys. Chem. Chem. Phys. 2014, 16, 13914–13931. Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett. 1988, 61, 2879–2882. Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458–9471. Ashraf, C.; van Duin, A. C. T. Extension of the ReaxFF Combustion Force Field toward Syngas Combustion and Initial Oxidation Kinetics. J. Phys. Chem. A 2017, 121, 1051–1068. Islam, M. M.; Zou, C.; van Duin, A. C. T.; Raman, S. Interactions of Hydrogen with the Iron and Iron Carbide Interfaces: A ReaxFF Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2016, 18, 761–771. Van Duin, A. C.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. ReaxFFSiO Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803– 3811. 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. Mortier, W. J.; Ghosh, S. K.; Shankar, S. Electronegativity-Equalization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315–4320. 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. Islam, M. M.; Kolesov, G.; Verstraelen, T.; Kaxiras, E.; van Duin, A. C. T. eReaxFF: A Pseudoclassical Treatment of Explicit Electrons within Reactive Force Field Simulations. J. Chem. Theory Comput. 2016, 12, 3463–3472. Strachan, A.; Kober, E. M.; van Duin, A. C.; Oxgaard, J.; Goddard III, W. A. Thermal Decomposition of RDX from Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122, 054502. Haley, B. P.; Li, C.; Wilson, N.; Jaramillo, E.; Strachan, A. Atomistic Simulations of Amorphous Polymers in the Cloud with PolymerModeler. ArXiv Prepr. ArXiv150303894 2015. 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. 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. Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. A Method for Tractable Dynamical Studies of Single and Double Shock Compression. Phys. Rev. Lett. 2003, 90, 235503. Tokmakoff, A.; Fayer, M. D.; Dlott, D. D. Chemical Reaction Initiation and Hot-Spot Formation in Shocked Energetic Molecular Materials. J. Phys. Chem. 1993, 97, 1901–1913. Xue, X.; Wen, Y.; Zhang, C. Early Decay Mechanism of Shocked ε-CL-20: A Molecular Dynamics Simulation Study. J. Phys. Chem. C 2016, 120, 21169–21177. Maharrey, S.; Behrens, R. Thermal Decomposition of Energetic Materials. 5. Reaction Processes of 1,3,5-Trinitrohexahydro-S-Triazine below Its Melting Point. J. Phys. Chem. A 2005, 109, 11236–11249. NIST WebBook http://webbook.nist.gov/ (accessed Jun 7, 2017). Wood, M. A.; Strachan, A. Nonequilibrium Reaction Kinetics in Molecular Solids. J. Phys. Chem. C 2016, 120, 542–552. Oyumi, Y.; Brill, T. B. Thermal Decomposition of Energetic Materials 3. A High-Rate, in Situ, FTIR Study of the Thermolysis of RDX and HMX with Pressure and Heating Rate as Variables. Combust. Flame 1985, 62, 213–224.


Page 29 of 43

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


(50)

(51) (52)

(53)

Rahaman, O.; van Duin, A. C.; Goddard III, W. A.; Doren, D. J. Development of a ReaxFF Reactive Force Field for Glycine and Application to Solvent Effect and Tautomerization. J. Phys. Chem. B 2010, 115, 249–261. van Duin, A. C. T.; Damsté, J. S. S. Computational Chemical Investigation into Isorenieratene Cyclisation. Org. Geochem. 2003, 34, 515–526. Ostadhossein, A.; Cubuk, E. D.; Tritsaris, G. A.; Kaxiras, E.; Zhang, S.; Duin, A. C. T. van. Stress Effects on the Initial Lithiation of Crystalline Silicon Nanowires: Reactive Molecular Dynamics Simulations Using ReaxFF. Phys. Chem. Chem. Phys. 2015, 17, 3832–3840. Ge, N.-N.; Wei, Y.-K.; Ji, G.-F.; Chen, X.-R.; Zhao, F.; Wei, D.-Q. Initial Decomposition of the Condensed-Phase β-HMX under Shock Waves: Molecular Dynamics Simulations. J. Phys. Chem. B 2012, 116, 13696–13704.



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 30 of 43

TOC Graphic


Page 31 of 43

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


D

E

ACS Paragon Plus Environment


9ROXPH

5HD[)) ,:

5HD[))

3UHVVXUH

5HD[)) OJ

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


Page 32 of 43

7HPSHUDWXUH

Page 33 of 43

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


D

([SHULPHQWDO GHQVLW\

E


)) GHQVLW\


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

D ([SHULPHQWDO GHQVLW\

E )) GHQVLW\


Page 34 of 43

Page 35 of 43

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


D

E

([SHULPHQWDO GHQVLW\


)) GHQVLW\


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

6KRFN 5HD[)) ,QWHUPHGLDWHV

Page 36 of 43

6WDEOH 6SHFLHV

*3D

+2

+&2 &2 +&+2 12

1

12 &12+ +212

*3D

21+


&2

1+

Page 37 of 43

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




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


Page 38 of 43

Page 39 of 43

)RUPDWLRQ (QHUJ\ NFDO PRO

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


&+ 2 +&2

+212 &12+ +12

+2

&2

1+

1

+21

5HD[))

5HD[)) ,:


5HD[)) OJ

([SW

+

2

F 5HD[)) OJ


E 5HD[))

'HQVLW\ RI VWDWHV D X

D 5HD[))

Page 40 of 43


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



Page 41 of 43

30

(a)

25 ΔΕa (kcal/mol)

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


(b)

20

ρo

1.2ρo 1.4ρo

15 10 5 0 ReaxFF-2014


ReaxFF-IW

ReaxFF-lg


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

, KK,


Page 42 of 43

Page 43 of 43

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



Decomposition and Reaction of Polyvinyl Nitrate under Shock and

Recommend Documents