Decomposition Mechanisms of α-Octahydro-1,3,5 ... - ACS Publications

Mar 21, 2017 - For example, the nanoscale explosives possess low melting point(1, ... is lacking an atomistic picture on the thermal decomposition pro...
1 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCC

Decomposition Mechanisms of α‑Octahydro-1,3,5,7-tetranitro1,3,5,7-tetrazocine Nanoparticles at High Temperatures Zhichao Liu,† Weihua Zhu,*,† Guangfu Ji,‡ Kefeng Song,‡ and Heming Xiao† †

Institute for Computation in Molecular and Materials Science and Department of Chemistry, Nanjing University of Science and Technology, Nanjing 210094, China ‡ Laboratory for Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang 621900, China S Supporting Information *

ABSTRACT: Density functional tight-binding molecular dynamics simulations with dispersion corrections were performed to study the adiabatic initial decomposition processes of molecular explosive αoctahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (α-HMX) nanoparticles (NPs) with the diameters of 1.4−2.8 nm at high temperatures from 2400 to 3000 K. The results indicate that the global thermal decompositions of the HMX NPs present great dependence on the temperature and particle size. The initial decomposition process of the HMX NPs includes two sequential stages: (i) competition between rapid expansion and unimolecular decompositions at surfaces; (ii) subsequent complex uni- and bimolecular decompositions. The main decomposition pathway at low temperatures is the isomerization reaction of the HMX molecule, which is quite different from the N−NO2 homolysis with ring opening at high temperatures. A second-order rate model was established to rationalize the global decompositions of the HMX NPs at different temperatures as well as the elementary decomposition paths. Our findings suggest that the thermal decompositions of the HMX NPs are quite different from their solid phase decompositions in both decomposition mechanisms and global kinetics. These findings provide basic understandings of initial decompositions of nanosized explosives. friction sensitivity5,11 but to enhance the sensitivity to external heat.2,4−10 However, experimentally prepared explosive NP samples often have wide particle size distribution, which might still obscure the size effects on the thermal stability and performance. Hence a major obstacle to understanding the initial chemical events of the nanoscale explosives contributing to further ignition is lacking an atomistic picture on the thermal decomposition process. Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX), an outstanding representative of nitramine explosives, has been widely used in many fields of energetic materials. Over the past decades, the decomposition mechanisms of HMX in both gas and condensed phases have been extensively studied.12−18 It was thought that the decomposition of the nitramine explosives proceeded via unimolecular channels in both gas and solid phases and their global activation energies were similar in gas and solid phases.12,13,19−22 Recent density functional theory (DFT) studies on two phases of HMX proposed that unimolecular N−NO2 homolysis, HONO elimination, nitro− nitrite rearrangement, and concerted ring fission acted as the initial decomposition steps of HMX.12,13,19−22 Lewis et al.21

1. INTRODUCTION Explosive nanoparticles (NPs), as an intermediate state between a bulk material and its constituent molecule, possess unusual physical and chemical properties. For example, the nanoscale explosives possess low melting point1,2 and sublimation enthalpy3 due to large free surface energy and low crystal energy. Many experimental studies found that the nanosized explosives with increased specific surface area have lower critical ignition temperature,2,4,5 higher burn rate,1,6 and more efficient energy release and complete combustion7,8 than the microsized ones or bulk material. In addition, the particle size has remarkable influence on the sensitivity of the explosives to the external heat, impact, friction, and electrostatic discharge. In all, nanostructuring of the explosives has been becoming an attractive and fast developing field. To understand the evolution of various physical and chemical properties from the constituent molecule to bulk material, the size effects of the explosive NPs on their behaviors have attracted extensive attention. Recent experiments on various kinds of explosives (HMX, RDX, CL-20, TATB, NTO, etc.) with different particle sizes indicated that their NPs had relatively lower thermal decomposition temperature and activation energy than the microsized or macroscale particles.2−5,7−10 Additionally, decreasing the particle size tends to decrease the impact and © XXXX American Chemical Society

Received: February 4, 2017 Revised: March 18, 2017 Published: March 21, 2017 A

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Configurations of relaxed (a) α-HMX bulk and NPs with average diameters of (b) 1.4 nm, (c) 1.8 nm, (d) 2.2 nm, and (e) 2.8 nm. The C, H, O, and N are represented by gray, white, red, and blue balls, respectively.

2. COMPUTATIONAL METHODOLOGY 2.1. Models and DFTB Calculations. In all calculations, we used the third-order corrected self-consistent charge density functional tight binding method23−26 with the first-principles London dispersion correction D327 (SCC-DFTB-D3, abbreviated as DFTB in this work) as implemented in CP2K package.28 Here we conducted a series of DFTB calculations for the HMX NPs with the size range of 1.4−2.8 nm at initial temperatures of 2400, 2600, 2800, and 3000 K, respectively. The initial configurations of the NPs were cut spherically with specified average diameters from the infinite α-HMX crystal. The initial structure of crystalline α-HMX, stable from approximately 377 to 429 K, was taken from the experimental X-ray crystal structure29 with full relaxation. Recently the DFTB method has been applied to many organic materials and presents comparable results with DFT calculations.30−35 The DFTB method calculated the lattice parameters of α-HMX to be a = 15.16 Å, b = 23.87 Å, and c = 6.131 Å, which are in good agreement with the experimental values29 of a = 15.14 Å, b = 23.89 Å, and c = 5.913 Å. We optimized the geometries of the HMX NPs in a box of 50 × 50 × 50 Å3 with periodic boundary conditions to model the thermal decomposition confined within a relatively free space. The SCC tolerance is 10−5 au. The conjugate gradient method (CG) was used for the cell optimization. The limited-memory Broyden−Fletcher−Goldfarb−Shanno method (LBFGS)36 was used for the geometry optimization. The convergence criteria for the maximum geometry change was less than 3.0 × 10−3 bohr, the maximum force component less than 4.5 × 10−4 hartree·bohr−1, the root-mean-square geometry change less than 1.5 × 10−3 bohr, and the root-mean-square force less than 3.0 × 10−4 hartree·bohr−1. 2.2. DFTB Molecular Dynamics Simulations. The DFTB molecular dynamics simulations (DFTB-MD) were performed under constant energy and volume conditions (NVE ensemble)37 with a 0.5 fs time step. Initial temperature was set to be 2400, 2600, 2800, and 3000 K, respectively. All DFTBMD simulations were started from energy minimum structures with initial velocities randomly chosen at each initial temperature. The snapshots were sampled every 5 fs for postsimulation analysis. To identify the concentration of stable molecular species, the chemical species were defined based upon a set of bond-distance cutoff criteria30 (see Table S1). Two atoms were not considered to be bonded unless they stayed continuously within the bond-distance criteria during the minimum bond lifetime of 25 fs.34−38 The time evolution of the species identified by the bond-distance criteria has no effect on the qualitative conclusions in this work. Thus, the chemical

investigated four possible decomposition pathways of gas phase α-HMX at the BLYP and B3LYP levels. They found that the N−NO2 homolysis (activation barriers in the range 40.5−41.8 kcal·mol−1) was energetically more favorable than other paths, i.e., the HONO elimination (42.9 kcal·mol−1), C−N bond scission (41.2−56.0 kcal·mol−1), and concerted ring opening (48.1−71.8 kcal·mol−1). Ab initio studies on the unimolecular decomposition of HMX at the B3LYP/6-31G(d) level including zero point energy corrections by Chakraborty et al.13 showed that the HONO elimination (44.6 kcal·mol−1) and N−NO2 homolysis (39.8 kcal·mol−1) were energetically more favorable than other pathways. Recent DFT calculations12,19,20 on the decomposition of the solid-state HMX by using the PBE functional and PAW method found that the ideal bulk HMX decomposed slower and required higher activation energies (43.67 kcal·mol−1 for the N−NO2 homolysis and 44.6 kcal· mol−1 for the HONO elimination) than the gas phase HMX. Sharia et al.19 found that the HMX molecules placed at surfaces and vacancies required lower activation energy for the N−NO2 homolysis (36.5−42.2 kcal·mol−1) and HONO elimination (38.1−43.1 kcal·mol−1) as compared to the bulk crystal. Also, Kuklja et al.20 proposed that the surface polarization in polar morphology of HMX significantly reduced the activation barriers and accelerated both the NO2 loss and HONO elimination. However, most molecules in the HMX NPs are exposed at the surfaces rather than in the interior of the crystal with bulk properties. Moreover, these molecules in the HMX NPs are different from those in the actual HMX surfaces. In the four polymorphs of HMX, the α-form is one of the polar polymorphs similar to the δ-form and was thought to have higher surface instability due to surface polarization. Therefore, the chemical reactivity and decomposition mechanisms of the α-HMX NPs with different sizes need to be clarified. In this work we report a quantum molecular dynamics simulation of the initiation chemistry of α-HMX NPs with the diameter range 1.4−2.8 nm at high temperatures (2400−3000 K) by maintaining constant energy and volume (NVE ensemble) to model adiabatic initiation process. A detailed physicochemical process during the initial decomposition is presented. We discuss the competing decomposition mechanism that depends on the system temperature and particle size from quantum molecular dynamics simulations. A second-order competing reaction model is used to rationalize the initiation decomposition of the α-HMX NPs. The results are carefully compared to previous theoretical and experimental results. B

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C reaction events and kinetic parameters were traced by the concentrations of the species. 2.3. DFT Calculations. To find the minimal energy path of each reaction and corresponding transition state, DFT calculations at different levels 22,30,39,40 (B3LYP/6-311g, B3LYP/6-311+g(d,p), and M06-2X/6-311+g(d,p)) as implemented in the Gaussian 09 package41 were performed to calculate accurate energies for the species involved in the decomposition as revealed by the DFTB-MD simulations. We performed the vibrational analysis of the structures at the same level to identify whether the optimized structures were local minima. The intrinsic reaction coordinates (IRC) using internal coordinates were used to guarantee that all transition states from DFT calculations were connected with the reactants and products. The activation barrier is defined by the difference between the total energies of the initial and transition states. The reaction energy is determined as the difference between the total energies of the initial and final configurations of the system. (See Figure 1.)

Figure 3. Snapshots of the 1.4 and 2.8 nm HMX NPs at starting temperature of 3000 K. The circles indicate where the initial decomposition takes place.

positions, which needs less than 5 ps for 2.8 nm@2400K_1 and 7.65 ps for 1.4 nm@2400K_1. Elevating the temperature will significantly accelerate this competing process, as shown in Figure 3. For the small NPs, the expansion seems to be easier to take place than the chemical decomposition. Also, note that the initial decompositions of the HMX NPs always start at the surfaces (see Figures 2b and 3b). This suggests that the surface molecules in the HMX NPs have higher reactivity than inner molecules. Figure 4 presents the initial decomposition time (t0) of the HMX NPs with different sizes at different temperatures. It is

3. RESULTS 3.1. Initial Decomposition Mechanisms. Here we studied the effects of different temperatures and sizes on the thermal decomposition of the α-HMX NPs using adiabatic DFTB-MD simulations with NVE ensemble. We performed 20 independent DFTB-MD simulations for each HMX NP (1.4 and 1.8 nm) at each temperature and one DFTB-MD simulation for each NP (2.2 and 2.8 nm) at each temperature, the total up to 124 independent simulations. The trajectory is denoted as 1.4 nm@2400K_n, where “1.4 nm” denotes the 1.4 nm NP, “2400K” denotes the starting temperature, and “n” denotes the independent simulation number (n = 1−20). Our simulations at different temperatures indicate that the initial decomposition process of the HMX NPs involves two sequential stages: (i) competition between rapid expansion and initial unimolecular decompositions at surfaces; (ii) subsequent complex uni- and bimolecular decompositions. 3.1.1. Competition between Expansion and Initial Decompositions. Figures 2 and 3 display the snapshots of the 1.4 and 2.8 nm HMX NPs at starting temperature of 2400 and 3000 K, respectively. It is seen in Figure 2 that the competing character between the expansion and chemical decompositions during the initial stage is remarkable for the large NPs at 2400 K. As the HMX NPs begin to sublimate, this competing stage is fast accompanied by chemical decom-

Figure 4. Initiation decomposition time (t0, ps) of different HMX NPs at various starting temperatures. t0 = 20 ps denotes that there are no reactions observed in the simulations. The circles denote the average values for the parallel simulations of the particular HMX NP.

found that the obtained t0 values display a clear dependence on both the stating temperature and particle size. The parallel tests for the 1.4 and 1.8 nm HMX NPs indicate that the initiation time presents certain statistical character in all the cases. The obtained t0 values for the small size HMX NPs display a wide and approximately uniform distribution, especially at low temperatures, where the fastest initiation is as short as the initiation at high temperature. But the slowest initiation is longer than the simulation time period (20 ps). Increasing the temperature tends to significantly shorten the initiation time, which leads to the narrower distribution of the t0 values. Moreover, it is seen in Figure 4 that the averaged t0 values from the parallel simulations show that the initial decomposition time for the large NPs is shorter than that for the small NP. Therefore, the particle size also affects the initial decomposition

Figure 2. Snapshots of the 1.4 and 2.8 nm HMX NPs at starting temperature of 2400 K. The circles indicate where the initial decomposition takes place. C

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C time. These findings further suggest the remarkable effects of the surface molecules on the decomposition chemistry of the HMX NPs. To understand the nature of expansion, we evaluated the surface energy (Esurf) and expansion enthalpy (ΔHsub) of the unreacted HMX NPs by the following equations:42 Esurf = A−1(Etotal − nHMX E HMX )

are selected based on the exposed crystal facets of all the studied HMX NPs. The surfaces of the α-HMX NPs and slab models are terminated by the NO2 and/or CH2 groups. Therefore, the exposed surfaces are polar. This results in the separation of positive and negative charges especially along the (001) facet. Among the (001), (011), (101), and (111) facets along different directions, the (001) facet possesses the highest surface tension. It is seen in Figure 5 that both Esurf and ΔHsub exhibit great dependence on the morphology of the HMX particle. The obtained surface energies of the HMX NPs in the range 160− 240 mJ·m−2 are close to those of the HMX crystal by our slab model calculations (see Figure 5a). This indicates that the HMX NPs have some characteristics of the HMX crystal surfaces. It is found in Figure 5b that the ΔHsub values of different sizes of HMX NPs in the range 1.4−4.0 nm are only 18.7−45.2% of that of the HMX bulk (24.8 kcal·mol−1). Unfortunately, our DFTB results significantly underestimate the experimental ΔHsub values of 44.16 and 41.80 kcal·mol−1.43 Similar underestimations on the sublimation energy using the DFTB method were also found in previous calculations.44,45 As the size of the HMX NP increases, its ΔHsub value gradually enhances. This is in agreement with the variation trend of the surface energy of the HMX NPs with the particle size. Our calculated results indicate that the HMX NPs have low sublimation enthalpies. Thus, they have lower melt and sublimation points than the bulk HMX. This will promote the diffusivity of the HMX molecules and decomposition products near the internal vacancies, voids, cracks, etc. even at low temperatures (lower than the temperature under the Chapman−Jouguet condition) as studied here. This suggestion is consistent with our DFTB-MD simulations on the HMX NPs. The high surface reactivity of the HMX NPs is further supported by the idea that microscopic gaps and voids play a key role in the formation of ejecta gas (the HMX molecules and

(1)

ΔHsub = −ΔE l − 2RT

(2)

where surface area of the HMX NP is A = πD in which D is the average diameter of the HMX NP, Etotal is the total energy, nHMX is the number of the HMX molecules in the system, EHMX is the normalized energy of the bulk HMX per molecule, and ΔEl is the lattice energy. For comparison we also calculated the surface energy of several crystal surfaces using the slab model, as listed in Figure 5a. The crystal facets summarized in Figure 5 2

Figure 5. DFTB-predicted surface energy and sublimation enthalpy of the HMX NPs with various particle sizes together with available theoretical and experimental results.

Figure 6. Time evolution of the main decomposition fragments of the 1.4 and 1.8 nm HMX NPs at different starting temperatures. D

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. A schematic diagram of the competing unimolecular decomposition channels for the HMX NPs. Two representative decomposition channels are listed in channel A, and one representative decomposition channel is listed in channel B.

HONO/HNO2, and HNO3 (see Figure 6) were observed at relatively later decomposition stages through subsequent radical reactions: quenching of two free radicals (NO2 + NO3 → N2O5) and trapping hydrogen atom from ring (NO2 + H−R → HONO/HNO2, and NO3 + H−R → HNO3). The DFTB-MD simulations show that the NO3 and NO2 possess high reactivity and can catalyze further the decomposition via a series of hydrogen transfer processes. Our predicted initial main products MN and N2O molecules and NO2 radical were supported by previous experimental results on gas phase HMX decomposition.13 The MN molecule was thought to undergo further decomposition reaction MN → N2O + CH2O as proposed in many studies.47,48 In all, the DFTB-MD simulations predicted that the decomposition mechanisms of the HMX NPs are governed by five unimolecular competition reaction channels: isomerization into linear isomer (channel A); N−NO2 homolysis with ring opening (channel B); concerted ring opening (channel C); HONO elimination (channel D); isomerization into tenmembered ring (channel E). A schematic diagram of the typical reaction pathways of the competing channels is illustrated in Figure 7. Other possible pathways with lower possibility are not presented for clarity since the first three channels make up a large majority (>95%) of the overall channels from the 124 independent trajectories (see below). For example, the HONO elimination at the studied temperatures was observed only 7 times and occurred only when the starting temperature was over 2800 K among all the 124 independent simulations. Our DFTB-MD predicted decom-

decomposition products) whose subsequent recompression magnifies the local temperature, which will enhance the sensitivity of the explosive HMX to initiation.46 3.1.2. Detailed Decomposition Mechanisms. After the competing stage, the global decompositions of the HMX NPs begin. Tracing all the decomposition channels from the 124 independent DFTB-MD simulations found that the decomposition of the HMX molecule is predominantly unimolecular. Their subsequent decompositions involve a complex set of uniand bimolecular reactions. Of course, this greatly depends on the system temperature. To demonstrate the temperature dependence, Figure 6 displays the time evolution of main chemical decomposition species for the 1.4 and 1.8 nm HMX NPs at various temperatures. For other size NPs, the time evolution of main chemical species is provided in Figure S1. At low temperature 2400 K, it is found that one of the dominant initial products is a linear isomer of the HMX molecule (Int1, see Figure 8), which to our knowledge has never been reported in previous studies. Our DFTB-MD simulations predict that the isomer Int1 possesses unexpected long lifetime (longest lifetime >17.85 ps), leading to a considerable concentration in early stage. Therefore, it may be suggested that the initial decomposition step of the HMX NPs at low temperature is the isomerization reaction of the HMX molecule by forming the less stable isomer (Int1), whose decomposition triggers the following degradation reactions. At high temperature, 3000 K, the most abundant products in the initial stage were found to be methylenenitramine (MN), NO2, and N2O. Less abundant intermediates like N2O5, E

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. Energy paths of (a) isomerization into linear isomer, (b) N−NO2 homolysis with ring opening, (c) concerted ring opening, (d) HONO elimination, and (e) isomerization into ten-membered ring observed in the DFTB-MD simulations by gaseous phase DFT calculations at various levels. The optimized configurations of transition states and intermediates were obtained at the UB3LYP/6-311+g(d,p) level. F

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 1. Activation Barriers (Ea, kcal·mol−1), Zero-Point Energy Corrected Activation Barriers (Ea,ZPE, kcal·mol−1), and Reaction Energies (ΔE, kcal·mol−1) of the Gas Phase Reactions observed in the Dominant Decomposition Channels of the αHMX NPs B3LYP/6-311G Ea

reaction channel A B

50.8 NO2 loss ring breaking

C D E e-NO2 loss a-NO2 loss a−e

39.2 (45.4,a 44.17,b 45.1d) 27.4 (23.03b) 70.0 (82.98e) 45.8 (47.29,b 41.5d) 43.7 (47.3d) 42.3 44.5

B3LYP/6-311+G(d,p) ΔE

Ea,ZPE 47.7

29.1

34.7 (39.8,a 39.72,b 40.5,c 40.3d) 25.1 (20.91,b 28.0c) 62.9 (71.8c) 40.8 (44.6,a 42.43,b 42.9,c 36.8d) 41.3 (44.8d) 37.4 40.1

34.7 14.9 (10.04,b 14.1c) 64.2 (65.5e) −6.69 (4.0,a 10.24,b 10.8c) 11.0 (11.5,c 12.8d) 37.4 40.1

M06-2X/6-311+G(d,p)

Ea

Ea,ZPE

ΔE

Ea

Ea,ZPE

ΔE

50.6

47.6

23.0

59.0

56.8

36.6

36.6 25.6

32.3 23.4

32.3 23.4

46.8 32.4

43.1 29.9

43.1 16.1

63.6

57.4

45.8

74.6

69.2

81.7

44.2

39.4

−13.2

48.5

44.6

−8.97

43.5

41.3

48.5

46.9

8.65

39.8 40.6

34.8 36.3

53.0 52.6

48.0 48.6

9.78 34.8 36.3

48.0 48.6

Available theoretical data were taken from refs 13, 22, 21, 20, and 12, respectively.

stages of the HMX NPs at low temperatures as found in our DFTB-MD simulations. In addition, our DFTB-MD simulations indicate that the formation of activated transition complex TS1 with a hexatomic ring structure, during which the three bonds participate in the transition reaction, is a multistep process, as illustrated in Figure 9. The DFTB-MD simulations on 2.8 nm@3000K_1 indicate that the formation of Int1 requires 44 fs even at high temperatures. The decomposition of Int1 has two intramolecular isomerization pathways with the barriers of +20.8 kcal·mol−1 (via TS2) and +14.5 kcal·mol−1 (via TS3) by UB3LYP/6-311+g(d,p), respectively. Further decomposition of Int1 starts with a low barrier of +6.6 kcal· mol−1 (via TS4) to form four fragments: Int1−4, NO2, NO3, and N2O, with a reaction energy of only +0.4 kcal·mol−1. It is seen in Figure 8b that there are two concerted steps in channel B: the barrierless N−NO2 homolysis and the C−N bond dissociation in meta-position. Figure 9 displays the snapshots of the concerted process from exemplary simulation on 2.8 nm@3000K_1. The barrier of the N−NO2 homolysis was estimated to be +32.3 kcal·mol−1 at the UB3LYP/6311+g(d,p) level. This is in agreement with previous gas phase calculations.12,13,19−22 The following reaction is the ring opening with a barrier of +23.4 kcal·mol−1 to form a radical Int4 (shown in Figure 8b), which is unstable and continues to release MN molecules via a low barrier process (ca. + 20.4 kcal· mol−1). Considering the first step being barrierless, the consecutive two steps have a total barrier of +55.7 kcal·mol−1 by UB3LYP/6-311+g(d,p), which is 8.1 kcal·mol−1 higher than channel A. Our calculated results for channel B (N−NO2 homolysis, ring opening, and release of MN molecules) are in good agreement with previous studies by Chakraborty et al.13 using the B3LYP/6-31G(d) method as well as the results by Zhang et al.22 at the B3LYP/cc-pVDZ level. Chakraborty et al. found that the ring opening after the N−NO2 homolysis required the lowest energy (+28.0 kcal·mol−1) and subsequent elimination of a MN molecule needs a barrier of +21.8 kcal· mol−1. Our DFT results confirm that channel B has the higher activation barrier, which supports our DFTB-MD simulations that channel B was less favorable than channel A at low temperatures.

position mechanisms of the HMX NPs were different from previous reports12,13,19−22 that the N−NO2 homolysis and HONO elimination predominated in the initial decompositions of gaseous and solid phase HMX. Also, Kuklja et al.20 pointed out that the surface polarization had dramatic effects on the activation barriers and reaction rates of polar δ-HMX and nonpolar β-HMX from DFT calculations. In the following section we try to explain the differences by studying thermal decomposition kinetics of the HMX NPs based on the DFTBMD simulations. 3.2. Kinetics of Thermal Decompositions. 3.2.1. Activation Barriers and Reaction Energies from DFT Calculations. The initiation chemistry of the α-HMX NPs predicted by our DFTB-MD simulations was different from previous studies on gaseous or condensed phase HMX. To rationalize the predicted decomposition mechanisms, we calculated and compared the energy paths of the decomposition reactions observed from the DFTB-MD simulations. Figure 8 illustrates the energy paths of the typical decomposition channels of the HMX NPs using three density functional methods: UB3LYP/6-311g (cyan), UB3LYP/6311+g(d,p) (black), and M06-2X/6-311+g(d,p) (red). All the methods present similar variation trends of reaction energies for each pathway. The activation barriers and reaction energies of the dominant decomposition channels with available theoretical data are summarized in Table 1. The structures of intermediates and transition states are provided in Figure 8. It is seen in Figure 8a that the HMX molecule isomerizes to form its isomer (Int1) through the NO2 group transfer together with the ring opening. The energy barrier for this process at TS1 is +47.6 kcal·mol−1 at the UB3LYP/6-311+g(d,p) level (+56.8 kcal·mol−1 at the M06-2X/6-311+g(d,p) level). We calculated the activation barrier for axial and equatorial N−NO2 bond cleavages to be +34.8 and +36.3 kcal·mol−1 at the UB3LYP/6-311+g(d,p) level (+48.0 and +48.6 kcal·mol−1 at the M06-2X/6-311+g(d,p) level), respectively. This indicates that the isomerization of the HMX molecule has higher energy barrier than the N−NO2 homolysis (ca. 11.3−12.8 kcal·mol−1). Therefore, it seems to be concluded that the path a is more energetically unfavorable than the N−NO2 homolysis. However, channel A is still dominant in the early decomposition G

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 9. Snapshots of isomerization into linear isomer (a), N−NO2 homolysis with ring opening (b), concerted ring opening (c), HONO elimination (d), and isomerization into ten-membered ring (e) observed in the decompositions of the HMX NPs at starting temperature of 3000 K.

precludes it from being a primary decomposition pathway. Despite the high activation barrier, the concerted ring opening pathway still serves as a main source for producing the MN molecules (main products at high temperatures). We found that the UB3LYP method with different basis sets predicted well-defined transition structures. However, the UB3LYP method with the 6-311g basis set can predict the channel with a barrierless reaction, while the M06-2X method failed to predict the transition structure.

The concerted ring opening pathway to yield four MN molecules in sequence is found to be the third dominant decomposition channel observed in our DFTB-MD simulations (see snapshots in Figure 9). This reaction proceeds via the highest barrier of +57.4 kcal·mol−1 by UB3LYP/6-311+g(d,p) among the observed decomposition channels, but lower than 71.8 kcal·mol−1 reported by previous theoretical studies at the B3LYP level and 82.98 kcal·mol−1 using the CNEBM methods.13,21 This channel is highly endothermic with 45.8 kcal·mol−1 above the ground-state α-HMX molecule, which H

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 10. Time dependence of the temperature and potential energies for the 1.4 and 2.8 nm HMX NPs at starting temperatures of 2400 and 3000 K. T0 and Tave denote the initial temperature and average temperature of the systems during the chemical decomposition period, respectively. The potential energies are set to be zero for the initial configurations of the HMX NPs.

Figure 11. Time evolution of the decay of the HMX molecules (a), proportion of the dominant decomposition channels (b), and global reaction rate (c) for the 1.8 nm HMX NP at starting temperatures of 2400 and 3000 K.

compared to the former three decomposition channels. Figure 9 presents the snapshots of the formation of Int 7 from 1.4 nm@3000K_9. This channel is slightly endothermic (+11.0 kcal·mol−1) as revealed by all the methods, in good agreement with previous work with the reaction energies of +11.5 kcal· mol−1 at the B3LYP/level21 as well as +12.8 kcal·mol−1 with the PBE functional.10 Further ring opening of Int7 was also proposed to proceed via a low energy barrier (ca. + 20 kcal· mol−1) in previous calculations.20 3.2.2. Second-Order Kinetics of Thermal Decompositions from DFTB-MD Simulations. Figure 10 presents the temperature and potential energy of the HMX NPs during our adiabatic DFTB-MD simulations. In all the cases the system temperature displays a downward trend at the earliest stages. At low temperatures, the potential energy for each NP system increases gradually as the simulation proceeds. This indicates

The HONO elimination mechanism requires a relatively higher barrier (+39.4 kcal·mol−1) than the N−NO2 homolysis, but significantly lower than the former three decomposition channels A, B, and C (Figure 8). The results agree well with previous calculations at different theoretical levels.12,13,19−22 However, the HONO elimination was found to occur to a much lesser extent (7 times observed in all 124 independent simulations of the HMX NPs and only found at starting temperatures over 2800 K). All three methods predicted similar activation barriers and weak exothermicity for the HONO elimination pathway. In addition, a stable 10-memberd oxy-ring isomer (Int7) was observed to form via TS9 with a barrier of +41.3 kcal·mol−1 (see Figure 8b), which is over 1.9 kcal·mol−1 for the barrier of the HONO elimination with +39.4 kcal·mol−1. This process occurred 5 times in all the simulations and was still unfavorable I

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 2. Calculated Ratio of Decomposition Channel (ω, %), Time History of Decomposition Channel (N) at t = 20 ps, Average Temperature (Tave), Initial Decomposition Time (t0, ps), Rate Constants of the Dominant Decomposition Channels (k, [N]−1·ps−1), and Global Rate Constant (K, [N]−1·ps−1) for the Decompositions of the HMX NPs Using the Second-Order Rate Model channel A

channel C

NP

Tave (K)

t0

ωA

kA × 10

ωB

kB × 10

1.4 nm

1175 1278 1399 1408 1195 1399 1180 1420

2.695 0.9675 0.95 1.165 0.9925 0.2 1.1525 0.2

80.2 57.3 50.3 31.1 66.5 30.1 31.6 45.1

8.50 8.94 13.9 16.8 2.41 8.07 11.4 110.0

11.4 28.9 39.7 65.4 24.5 60.1 68.4 51.9

1.21 4.51 11.0 35.2 0.888 16.1 24.8 126.6

1.8 nm 2.8 nm a

channel B

The global rate constant K is derived from eq 2: K =

5

1 t − t0

(

1 NHMX



1 NHMX,0

that the HMX NPs present the endothermic behavior at the studied time range. At high temperatures, the potential energies of all the NP systems have a maximum at 4−5 ps. This is consistent with the analysis from our DFTB-MD trajectories that the first 5 ps is the competition stage of the expansion and initial decomposition for the HMX NPs. Both of them are endothermic. The following decomposition by this stage is a series of exothermic reactions of intermediates. Next we try to rationalize the dependence of the decomposition kinetics of the HMX NPs on the temperature and particle size by establishing a reaction rate model for the thermal decomposition paths observed from our DFTB-MD simulations. Based on the energy paths of the decomposition reactions observed from the DFTB-MD simulations, it is found that the rate-determining step is the initial decomposition step of the HMX NP. Thus, the global kinetics in all the cases is determined by the rate-determining steps of various decomposition channels. Here we take the decompositions of the 1.8 nm NP as an example. The decay of the HMX molecules and the assignment of different decomposition channels (N) are illustrated in Figure 11a, and corresponding results for other size NPs are illustrated in Figure S2. To determine the reaction orders of the rate-determining steps, we computed the proportion (ω) of the dominant decomposition channels of the HMX NPs. It is found that the variation of ω for the dominant channels tends to be stable as the decomposition proceeds (see Figure 11b for the 1.8 nm NP and Figure S3 for other size NPs). This indicates that the ratedetermining steps have the same reaction order. Thus, the relation between the rate constant and the proportion of the channels can be obtained by eq 3: kA :kB:k C = NA :NB:NC = ωA :ωB:ωC

ωC

kC × 105

K × 105a

8.39 11.5 4.02 1.31 4.40 5.22 0.00 1.13

0.889 1.79 1.11 0.706 0.160 1.40 0.00 2.76

10.6 15.6 27.6 53.9 3.63 26.8 36.2 243.9

5

). NHMX =

NA =

NHMX,0 NHMX,0K (t − t0) + 1

(4)

2 kA NHMX,0 K (t − t0) K NHMX,0K (t − t0) + 1

(5)

kA = ωA K ,

kB = ωBK ,

k C = ωCK

(6)

where NHMX,0 and NHMX are the number of total and reacted HMX molecules, respectively; K is the global rate constant (K ≈ kA + kB + kC). A simple relation between the elementary and global rate constants can be obtained by eq 6. The fitting parameters are presented in Table 2. The fitting curves are displayed in Figure 11a for the 1.8 nm HMX NP and in Figure S2 for other size NPs. The fitting quality was good in general when compared with the DFTB-MD results except for the relatively poor one at low temperature due to the limited decomposition data at the simulation time period. To validate the model, we computed the global reaction rate (K) for the 1.8 nm HMX NP at different temperatures by using the number of unreacted HMX molecules. The obtained K confirms that the decomposition kinetics of the HMX NP has second-order character at the earliest stages (see Figure 11c). In all, this conclusion is quite different from previous studies19,20,49 on the decompositions of condensed phase HMX that the unimolecular decomposition reactions have the characteristics of first-order kinetics. The kinetic parameters in the second-order rate model of the decompositions of the HMX NPs based on the DFTB-MD simulations can be obtained as follows. Table 3 lists the obtained reaction rates and activation energies of different channels for different sizes of HMX NPs under different

(3)

Table 3. Activation Barriers (Ea, kcal·mol−1) and PreExponential Factors (ln A, [N]−1·ps−1) of the Global Decomposition and Rate-Determining Steps of the Dominant Decomposition Channels in the Decompositions of the HMX NPs Predicted by DFTB-MD Simulations

where kA, kB, and kC are the rate constants for channels A, B, and C, respectively; NA, NB, and NC are the time histories of the rate-determining steps for channels A, B, and C, respectively; ωA, ωB, and ωC are the ratios of channels A, B, and C to the total channels, respectively. Several rate competing models were fitted to the systems, and their fitting parameters are provided in Table S2. It is found that the second-order rate model possesses higher fitting quality than others. Thus, the DFTB-MD predicted decomposition kinetics of the HMX NPs can be fitted by a second-order rate competing model:

1.4 nm

J

1.8 nm

2.8 nm

channel

Ea

ln A

Ea

ln A

Ea

ln A

global A B

22.9 9.57 47.6

12.2 6.24 20.6

32.6 19.7 47.2

15.0 9.17 19.8

26.5 31.7 22.6

14.9 15.8 12.9

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

ps). The DFTB-predicted expansion enthalpies for the HMX NPs with different particle sizes in the range 1.4−4.0 nm are less than 45.2% of that for the ideal bulk. The smaller the size of the HMX NP is, the easier it tends to sublimate. Thus, the first stage is dominated by the competition between the rapid expansion and the chemical initiation of the HMX molecules. At this stage, we find that the expansion of the smaller NPs was prior to the initiation decomposition. For the larger NPs, the competition is more remarkable, where the chemical initiation seems to be faster than the expansion. It is interesting to note that the initial decompositions occur preferentially at the surfaces of the HMX NPs (see 1.15 ps from 2.8 nm@2400K_1 in Figure 2 and 0.24 ps from 2.8 nm@3000K_1 in Figure 3). This is in good agreement with previous DFT studies that the surface HMX molecules possess lower activation barriers and higher reaction rates. The second stage is onset of the uni- and bimolecular competing decompositions of the HMX molecules. Our DFTBMD simulation presents that the dominant initial decomposition step of the HMX NP at lower temperature is an isomerization reaction by forming linear isomer (Int1 in channel A). However, our DFT calculations indicate that the activation barrier for this isomerization reaction is +47.6 kcal· mol−1, higher than that for the N−NO2 homolysis (+34.8−36.3 kcal·mol−1) or HONO elimination (+39.4 kcal·mol−1), as shown in Figure 8. This indicates that our DFTB-MD results are different from our DFT calculations and previous studies on the decompositions of gaseous and condensed phase HMX.12,13,19−22 Next we try to explain their differences. As observed from the DFTB-MD simulations, the formation of the activated complex TS1 is a multistep process rather than a concerted reaction as obtained from our DFT calculations. It is found in Figure 9 that the multiple processes to form Int1 last as long as 44 fs even at high temperatures (ca. 1500 K). The snapshots from 2.8 nm@3000K_1 in Figure 9 present that there hardly exists a well-defined transition structure as TS1 from DFT calculations. This indicates that the isomerization is more complex dynamically. In addition, the DFTB-MD simulations display that the isomer Int1 is capable of producing large amounts of smaller gas products (e.g., NO2, NO3, N2O, MN, and etc.) (see Figure 8), which leads to more than 4 times of volumetric expansion. Therefore, the local molecular packing and external conditions for the decompositions of the HMX NPs at low temperatures studied by our DFTB-MD simulations are completely different from our DFT calculations on gaseous HMX molecules and previous reports on both gas phase and condensed phase HMX. This may be main cause for leading to a difference in the decomposition mechanism between the HMX NPs and other (gas and solid) phase HMX. At higher starting temperatures 2800−3000 K, our DFTBMD simulations reveal that a consecutive process of the N− NO2 homolysis and subsequent ring opening prevails over channel A. The activation barriers for the N−NO2 homolysis and the ring opening are +32.3 and +23.4 kcal·mol−1 at the UB3LYP/6-311+g(d,p) level, respectively. However, the HONO elimination pathway is found to take place to a much less extent, which is observed only 7 times at higher starting temperatures over 2800 K. This suggests that both the N−NO2 homolysis and HONO elimination are more favorable at higher temperatures but seem to be unfavorable at the lower temperature range 2400−2800 K. The initial gas products are MN, N2O, and NO2, which are in agreement with previous

temperatures. Figure 12 presents the logarithm of the reaction rates of global decomposition and several decomposition

Figure 12. Logarithm of the reaction rates of global decompositions and several decomposition channels (inset) for different HMX NPs as a function of 1000/T.

channels (inset) for different HMX NPs as a function of 1000/T. It is seen in Table 3 that the global reaction rates for different sizes of HMX NPs are both size- and temperaturedependent.50 For a particular NP, as the starting temperature increases from 2400 to 3000 K, the global reaction rates increase by 5.1−6.7 times. For different sizes of HMX NPs, the size effect is also significant. The 2.8 nm NP has faster reaction rate than smaller NPs at the same starting temperature (3.4− 4.5 times faster than that of the 1.4 nm NP and 9.1−10.0 times faster than that of the 1.8 nm NP). It is found in Figure 12 that the decomposition rate of channel A for the HMX NPs increases as the temperature or the particle size increases. The same is true of other different decomposition channels. The global activation energies in all the cases predicted by the DFTB-MD simulations lie in the range 22.9−32.6 kcal·mol−1, which are lower than DFT predicted results (39.4−57.4 kcal·mol−1) and relatively lower than the gas phase experiments of 32.1−52.9 kcal·mol−1.15 The DFTB-MD predicted activation energies for channels A and B in Table 3 show that channel B has relatively higher activation energy than channel A with an exception for the 2.8 nm HMX NP. This is mainly due to the limited data for the lowtemperature decomposition of the HMX NP. Overall, the average activation energy for channel A is +20.3 kcal·mol−1 and for channel B is +39.1 kcal·mol−1, which is in qualitative agreement with our gas phase calculations.

4. DISCUSSION To understand the differences between the thermal decomposition of the explosive NP and its bulk crystal, the early thermal decomposition stages of the nanoscale α-HMX particles were studied using the DFTB-MD simulations. The thermal decompositions of the HMX NPs include (i) competition between rapid expansion and unimolecular initial decomposition at surfaces; (ii) subsequent complex uni- and bimolecular decomposition reactions. Both the DFTB-predicted sublimation enthalpies and DFTBMD simulations for the HMX NPs indicate that the quantum size effect causes a significant reduction in the sublimation enthalpies of the NPs, which results in the rapid sublimation for the NPs during the earliest thermal decomposition stages (0−5 K

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

particle size accelerates the reaction rate. The global activation energies in the three main channels were calculated to be in the range 22.9−32.6 kcal·mol−1. The average activation energies for channels A and B are +20.3 and +39.1 kcal·mol−1, respectively. Although the results are lower than corresponding DFT results, their qualitative trend is consistent. Our findings suggest that the thermal decompositions of the HMX NPs are quite different from their solid phase decompositions in both decomposition mechanisms and global kinetics. These findings provide basic understandings of initial decomposition of nanosized energetic materials. Also, they provide further insights into the comprehensive understanding of the thermal decomposition of energetic materials.

experiments that the N2O, CH2O, and NO2 were found to be main gas products for the decompositions of HMX.13 The decomposition kinetics of the HMX NPs indicates that all the rate-determining steps of the three main channels have the characteristics of second-order kinetics. The results are quite different from previous studies on the decompositions of gaseous and solid phase HMX that all the N−NO2 homolysis, HONO elimination, and NONO rearrangement had welldefined transition states and were first-order reactions. However, our DFTB-MD simulations reveal that all the ratedetermining steps in predominated decomposition mechanisms are multistep processes. Thus, it may be inferred that the decomposition kinetics of the HMX NPs is different from the solid phase decomposition in the studied temperature range. Appling the second-order rate model to the DFTB-MD simulations found that the global reaction rates for different HMX NPs are both size- and temperature-dependent. Increasing either the temperature or particle size accelerates the decomposition reaction rate. The global activation energies in the three channels A, B, and C are predicted by the secondorder rate model to be in the range of 22.9−32.6 kcal·mol−1, smaller than the DFT-predicted results (39.4−57.4 kcal·mol−1 for the main decomposition channels) and relatively lower than the gas phase experiments of 32.1−52.9 kcal·mol−1.15 The average activation energies for channels A and B are +20.3 and +39.1 kcal·mol−1, respectively. Although the results are lower than the DFT results of +47.6 kcal·mol−1 for channel A and of +55.7 kcal·mol−1 for channel B, their qualitative trend is consistent. The overall mechanism, to the best of our knowledge, is quite different from previously proposed decomposition mechanisms for solid HMX or similar nitroamine explosives.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b01136. Bond-distance cutoff criteria, regression coefficient of linearity of different fitting models, and time evolution studies (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Weihua Zhu: 0000-0002-0295-8744 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the NSAF Foundation of National Natural Science Foundation of China and China Academy of Engineering Physics (Grant No. U1530104) and the Science Challenging Program.

5. SUMMARY AND CONCLUSIONS In summary, comprehensive DFTB-MD simulations in combination with DFT calculations have been performed to study the thermal decompositions of different α-HMX NPs at high temperatures. The global thermal decompositions of the HMX NPs present great dependence on the temperature and particle size. In the early decomposition stage, the expansion of the HMX NPs is remarkable in competition with the initial chemical decomposition. The smaller NPs have the lower sublimation enthalpies. The following thermal decompositions are greatly dependent on the system temperatures. At low temperatures, the predominating decomposition is the isomerization of the HMX molecule, while at high temperatures, the N−NO2 homolysis accompanied by ring breaking dominates. The concerted ring opening and HONO elimination occurs to a much less extent in the studied temperatures. Our simulations indicate that there are a large number of activated complexes during the low temperature decomposition channel, which leads to more than 4 times of volumetric expansion. This might be unfavorable in condensed phase HMX but can be available in real densely packed nanoscale HMX due to their differences in local molecular packing and external conditions. The subsequent expansion and recompression can magnify the local temperature and further accelerate their decompositions. The rate-limiting steps of the dominating decomposition channels at the earliest decomposition stages exhibit a clear second-order reaction character. It is found that the global reaction rates for different HMX NPs are both size- and temperature-dependent. Increasing either the temperature or



REFERENCES

(1) Spitzer, D.; Comet, M.; Baras, C.; Pichot, V.; Piazzon, N. Energetic Nano-materials: Opportunities for Enhanced Performances. J. Phys. Chem. Solids 2010, 71, 100−108. (2) Fathollahi, M.; Mohammadi, B.; Mohammadi, J. Kinetic Investigation on Thermal Decomposition of Hexahydro-1,3,5trinitro-1,3,5-triazine (RDX) Nanoparticles. Fuel 2013, 104, 95−100. (3) King, W. P.; Saxena, S.; Nelson, B. A.; Weeks, B. L.; Pitchimani, R. Nanoscale Thermal Analysis of an Energetic Materials. Nano Lett. 2006, 6, 2145−2149. (4) Akkbarzade, H.; Parsafar, G. A.; Bayat, Y. Structural Stability of Nano-sized Crystals of HMX: A Molecular Dynamics Simulation Study. Appl. Surf. Sci. 2012, 258, 2226−2230. (5) Fathollahi, M.; Pourmortazavi, S. M.; Hosseini, S. G. Particle Size Effects on Thermal Decomposition of Energetic Material. J. Energ. Mater. 2007, 26, 52−69. (6) Reid, D. L.; Kreitz, K. R.; Stephens, M. A.; King, J. E. S.; Nachimuthu, P.; Petersen, E. L.; Seal, S. Development of Highly Active Titania-Based Nanoparticles for Energetic Materials. J. Phys. Chem. C 2011, 115, 10412−10418. (7) Reid, D. L.; Russo, A. E.; Carro, R. V.; Stephens, M. A.; LePage, A. R.; Spalding, T. C.; Petersen, E. L.; Seal, S. Nanoscale Additives Tailor Energetic Materials. Nano Lett. 2007, 7, 2157−2161. (8) Armstrong, R. W.; Baschung, B.; Booth, D. W.; Samirant, M. Enhanced Propellant Combustion with Nanoparticles. Nano Lett. 2003, 3, 253−255. (9) Liu, R.; Yu, W.; Zhang, T.; Yang, L.; Zhou, Z. Nanoscale Effect on Thermal Decomposition Kinetics of Organic Particles: Dynamic

L

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

β-HMX under Shock Loadings via Molecular Dynamics Simulations in Conjunction with Multiscale Shock Technique. J. Phys. Chem. B 2014, 118, 8691−8699. (32) Ge, N.; Wei, Y.; Ji, G.; Chen, X.; Zhao, F.; Wei, D. Initial Decomposition of the Condensed-Phase β-HMX under Shock Waves: Molecular Dynamics Simulations. J. Phys. Chem. B 2012, 116, 13696− 13704. (33) Qi, T.; Bauschlicher, C. W., Jr.; Lawson, J. W.; Desai, T. G.; Reed, E. J. Comparison of ReaxFF, DFTB, and DFT for Phenolic Pyrolysis. 1. Molecular Dynamics Simulations. J. Phys. Chem. A 2013, 117, 11115−11125. (34) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. NitrogenRich Heterocycles as Reactivity Retardants in Shocked Insensitive Explosives. J. Am. Chem. Soc. 2009, 131, 5483−5487. (35) Manaa, M. R.; Fried, L. E.; Melius, C. F.; Elstner, M.; Frauenheim, Th. Decomposition of HMX at Extreme Conditions: A Molecular Dynamics Simulation. J. Phys. Chem. A 2002, 106, 9024− 9029. (36) Byrd, R. H.; Lu, P.; Nocedal, J.; Zhu, C. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM J. Sci. Comput. 1995, 16, 1190−1208. (37) Wood, M. A.; van Duin, A. C. T.; 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) Wu, C. J.; Fried, L. E.; Yang, L. H.; Goldman, N.; Bastea, S. Catalytic Behavior of Dense Hot Water. Nat. Chem. 2009, 1, 57−62. (39) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Oother Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (40) Kiselev, V. G.; Gritsan, N. P. Unexpected Primary Reactions for Thermolysis of 1,1-Diamino-2,2-dinitroethylene (FOX-7) Revealed by ab Initio Calculations. J. Phys. Chem. A 2014, 118, 8002−8008. (41) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al.; Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (42) Bisker-Leib, V.; Doherty, M. F. Modeling the Crystal Shape of Polar Organic Materials: Prediction of Urea Crystals Grown from Polar and Nonpolar Solvents. Cryst. Growth Des. 2001, 1, 455−461. (43) Lyman, J. L.; Liau, Y. Thermochemical Functions for Gas-Phase, 1,3,5,7-Tetranitro-1,3,5,7-tetraazacyclooctane (HMX), Its Condensed Phases, and Its Larger Reaction Products. Combust. Flame 2002, 130, 185−203. (44) Yuan, T.; Larsson, K. Theoretical Study of Size Effects on Surface Chemical Properties for Nanoscale Diamond Particles. J. Phys. Chem. C 2014, 118, 26061−26069. (45) Qu, X.; Latino, D.; Aires-de-Sousa, J. A Big Data Approach to the Ultra-Fast Prediction of DFT-Calculated Bond Energies. J. Cheminf. 2013, 5, 34−47. (46) Holian, B. L.; Germann, T. C.; Maillet, J. B.; White, C. T. Atomistic Mechanism for Hot Spot Initiation. Phys. Rev. Lett. 2002, 89, 285501. (47) Brill, T. B. Multiphase Chemistry Considerations at the Surface of Burning Nitramine Monopropellants. J. Propul. Power 1995, 11, 740−751. (48) Morgan, C. U.; Beyer, R. A. Electron-Spin-Resonance Studies of HMX Pyrolysis Products. Combust. Flame 1979, 36, 99−101. (49) Zhou, T.; Huang, F. Effects of Defects on Thermal Decomposition of HMX via ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 278−287. (50) 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.

Vacuum Stability Test of 1,3,5-Triamino-2,4,6-trinitrobenzene. Phys. Chem. Chem. Phys. 2013, 15, 7889−7895. (10) Kumar, R.; Siril, P. F.; Soni, P. Preparation of Nano-RDX by Evaporation Assisted Solvent-Antisolvent Interaction. Propellants, Explos., Pyrotech. 2014, 39, 383−389. (11) Risse, B.; Schnell, F.; Spitzer, D. Synthesis and Desensitization of Nano-β-HMX. Propellants, Explos., Pyrotech. 2014, 39, 397−401. (12) Sharia, O.; Kuklja, M. M. Modeling Thermal Decomposition Mechanisms in Gaseous and Crystalline Molecular Materials: Application to β-HMX. J. Phys. Chem. B 2011, 115, 12677−12686. (13) Chakraborty, D.; Muller, R. P.; Dasgupta, S.; Goddard, W. A., III Mechanism for Unimolecular Decomposition of HMX (1,3,5,7Tetranitro-1,3,5,7-tetrazocine), an ab Initio Study. J. Phys. Chem. A 2001, 105, 1302−1314. (14) Tarver, C. M.; Tran, T. D. Thermal Decomposition Models for HMX-Based Plastic Bonded Explosives. Combust. Flame 2004, 137, 50−62. (15) Brill, T. B.; Gongwer, P. E.; Williams, G. K. Thermal Decomposition of Energetic Materials. 66. Kinetic Compensation Effects in HMX, RDX, and NTO. J. Phys. Chem. 1994, 98, 12242− 12247. (16) Wu, C. J.; Fried, L. E. Ab Initio Study of RDX Decomposition Mechanisms. J. Phys. Chem. A 1997, 101, 8675−8679. (17) Long, G. T.; Vyazovkin, S.; Brems, B. A.; Wight, C. A. Competitive Vaporization and Decomposition of Liquid RDX. J. Phys. Chem. B 2000, 104, 2570−2574. (18) Chakraborty, D.; Muller, R. P.; Dasgupta, S.; Goddard, W. A. The Mechanism for Unimolecular Decomposition of RDX (1,3,5Trinitro-1,3,5-triazine), an ab Initio Study. J. Phys. Chem. A 2000, 104, 2261−2272. (19) Sharia, O.; Kuklja, M. M. Rapid Materials Degradation Induced by Surfaces and Voids: Ab Initio Modeling of β-Octatetramethylene Tetranitramine. J. Am. Chem. Soc. 2012, 134, 11815−11820. (20) Kuklja, M. M.; Tsyshevsky, R. V.; Sharia, O. Effect of Polar Surfaces on Decomposition of Molecular Materials. J. Am. Chem. Soc. 2014, 136, 13289−13302. (21) Lewis, J. P.; Glaesemann, K. R.; VanOpdorp, K.; Voth, G. A. Ab Initio Calculations of Reactive Pathways for α-Octahydro-1,3,5,7tetranitro-1,3,5,7-tetrazocine (α-HMX). J. Phys. Chem. A 2000, 104, 11384−11389. (22) Zhang, S.; Nguyen, H. N.; Truong, T. N. Theoretical Study of Mechanisms, Thermodynamics, and Kinetics of the Decomposition of Gas-Phase α-HMX (Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine). J. Phys. Chem. A 2003, 107, 2981−2989. (23) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (24) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678−5684. (25) Elstner, M. SCC-DFTB: What Is the Proper Degree of SelfConsistency? J. Phys. Chem. A 2007, 111, 5614−5621. (26) Elstner, M. The SCC-DFTB Method and Its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (27) Brandenburg, J. G.; Grimme, S. Accurate Modeling of Organic Molecular Crystals by Dispersion-Corrected Density Functional Tight Binding (DFTB). J. Phys. Chem. Lett. 2014, 5, 1785−1789. (28) CP2K Developers Group. CP2K Code; http://cpk2.berlios.de. (29) Cady, H. H.; Larson, A. C.; Cromer, D. T. The Crystal Structure of α-HMX and a Refinement of the Structure of β-HMX. Acta Crystallogr. 1963, 16, 617−623. (30) An, Q.; Liu, W.; Goddard, W. A.; Cheng, T.; Zybin, S. V.; Xiao, H. Initial Steps of Thermal Decomposition of Dihydroxylammonium 5,5′-bistetrazole-1,1′-diolate Crystals from Quantum Mechanics. J. Phys. Chem. C 2014, 118, 27175−27181. (31) Ge, N.; Wei, Y.; Song, Z.; Chen, X.; Ji, G.; Zhao, F.; Wei, D. Anisotropic Responses and Initial Decomposition of Condensed-Phase M

DOI: 10.1021/acs.jpcc.7b01136 J. Phys. Chem. C XXXX, XXX, XXX−XXX