Reactive molecular dynamics simulation of the high-temperature pyrolysis of 2,2′,2′′,4,4′,4′′,6,6′,6′′-nonanitro-1,1′:3′,1′′-terphenyl (NONA)

2,2′,2′′,4,4′,4′′,6,6′,6′′-Nonanitro-1,1′:3′,1′′-terphenyl (NONA) is currently recognized as an excellent heat-resistant explosive. To improve the atomistic understanding of the thermal decomposition paths of NONA, we performed a series of reactive force field (ReaxFF) molecular dynamics simulations under extreme conditions of temperature and pressure. The results show that two distinct initial decomposition mechanisms are the homolytic cleavage of the C–NO2 bond and nitro–nitrite (NO2 → ONO) isomerization followed by NO fission. Bimolecular and fused ring compounds are found in the subsequent decomposition of NONA. The product identification analysis under finite time steps showed that the gaseous products are CO2, N2, and H2O. The amount of CO2 is energetically more favorable for the system at high temperature or low density. The carbon-containing clusters are a favorable growth pathway at low temperatures, and this process was further demonstrated by the analysis of diffusion coefficients. The increase of the crystal density accelerates the decomposition of NONA judged by the analysis of reaction kinetic parameters and activation barriers. In the endothermic and exothermic stages, a 20% increase in NONA density increases the activation energies by 3.24 and 0.48 kcal mol−1, respectively. The values of activation energies (49.34–49.82 kcal mol−1) agree with the experimental data in the initial decomposition stage.


Introduction
Increases in the depth of oil and gas wells have led to corresponding changes in drilling and operating conditions. For the explosives employed, the increase in the temperature and pressure during production operations and trouble-shooting is of decisive signicance. In particular, super heat-resistant energetic materials maintain stability at temperatures of up to 250 C. 1-3 2,2 0 ,2 00 ,4,4 0 ,4 00 ,6,6 0 ,6 00 -Nonanitro-1,1 0 :3 0 ,1 00 -terphenyl (NONA, C 18 H 5 O 18 N 9 ) that is safe, reliable and stable at high temperature is currently recognized as a promising heatresistant explosive (decomposition temperature: 440-450 C). 4-6 NONA, as one of the organic polynitro compounds, releases large amounts of energy during thermal decomposition. Hence, it plays a crucial role in the exploration and exploitation of oil or space. Understanding the physical and chemical events is essential to obtain modied models for its combustion and detonation. However, the evolution details of the ignition and explosion of the condensed phase NONA remain unclear. The rapid reaction process of energetic materials under extreme conditions involves short time scales and complex chemical reaction kinetics, which poses signicant challenges to current experimental techniques. Reactive molecular dynamics (ReaxFF MD) simulations enable to directly probe the evolution processes of chemical reaction under the microscopic scale. The ReaxFF MD simulations have been widely used to explore the thermal decomposition mechanism of condensed phase explosives, such as nitromethane (NM), 7 TNT, 8,9 RDX, 10 HMX, 11 CL-20, 12 and 1,3,5-triamino-2,4,6trinitrobenzene (TATB). 13 Nitroaromatic explosives have higher activation energies in the condensed phase as can be seen in Table 1, and the major nal products are N 2 , CO 2 , and H 2 O.
Considering that NONA belongs to nitro-based organic explosives, we summarize the basic research and applied research on the thermal decomposition of a literature review. Rom et al. 7 used ReaxFF MD simulations to investigate the decomposition mechanism of hot liquid NM at various compressions. They observed two different initial thermal decomposition schemes: (1) unimolecular C-N bond cleavage followed by subsequent C-NO 2 steps forming eventually three NO 2 (g) molecules and $CH 3 radical at low density; (2) the formation of CH 3 NO via H-transfer and N-O bond rupture at high density. Rom et al. 8 used ReaxFF MD simulations to analyze the thermal decomposition of liquid phase TNT at high temperature. The dissociation of the C-NO 2 bond is putatively the initial mechanism of liquid phase TNT decomposition at low density, whereas dimer formation and decomposition at higher compressions. Furman et al. 9 elucidated the difference between activation energies in gas phase and condensed phase of TNT, and discussed the uni-and bimolecular reaction mechanism of thermal decomposition of condensed phase TNT. Zhang et al. 13 used the ReaxFF MD simulation to study the process of thermal decomposition of TATB to form carbon clusters. They found that TATB decomposes to form a large amount of carbon clusters (15-30% of the total mass of the system). Chen et al. 14 suggested two distinct mechanisms for initial decomposition of 2,2 0 ,4,4 0 ,6,6 0 -hexanitrostilbene (HNS): C-NO 2 bond dissociation and nitro-nitrite isomerization. The activation energy of thermolysis for NONA is a crucial parameter of performance. Zeman et al. 5 have shown the activation energy for thermal decomposition of NONA (56.48 kcal mol À1 ) by the Russian manometric method. Keshavarz et al. 15 predicted the activation energy of the thermolysis in the condensed state for polynitro arenes, and they obtained activation energy of NONA (51.14 kcal mol À1 ).
However, the detailed understanding of the transient evolution of condensed phase NONA as well as its interplay with the environmental conditions at initial stage of thermal decomposition is currently lacking. In this work, we performed the simulations to study the thermal decomposition of NONA. This paper is organized as follows: Computational details for both the simulation method and kinetic parameter analysis are described in Section 2. Results and discussions (decomposition pathways, exothermic phase, fragments analysis, and reaction kinetic parameters) are presented in Section 3. Finally, conclusions are summarized in Section 4.

Simulation method
All the molecular simulations were performed with the ReaxFF-lg force eld using LAMMPS program package. The original unit cell obtained from the experiment was expanded into a 3 Â 4 Â 2 NONA periodic supercell, and the system contains 192 molecules (24 units, 9600 atoms). 6 The molecular formula, threedimensional structure and initial supercell of NONA are shown in Fig. 1. First, in order to minimize the energy, a conjugate gradient algorithm was used to relax the NONA supercell. The canonical ensemble (NVT) and the Berendsen thermostat were applied to the MD simulation with a total time of 10 ps at 300 K, which further relaxed the NONA supercell. The relaxed supercell is performed with isothermal-isobaric (NPT) MD simulation at 0 atm and 300 K for 20 ps, and the equilibrium structure was obtained. In order to verify the feasibility of the parameters, the Expt  equilibrium supercell structure was compared to the single crystal diffraction data. Then, the relaxed supercell (V 0 ¼ 4.646 Â 10 3 cm 3 , crystal density r ¼ 1.816 g cm À3 ) was used to generate the volumetrically expanded and compressed supercells and would be used in a high temperature pyrolysis reaction. By changing the lattice parameters of the xed fractional coordinates of the superelement, we obtained the corresponding 1.1V 0 (10% expansion, r ¼ 1.365 g cm À3 ) and 0.9V 0 (10% compression, r ¼ 2.492 g cm À3 ) supercells. To further reduce any articial pressure caused by altering unit cell parameters, the supercell is minimized by energy and further evolved by the microcanonical ensemble (NVE) for another 5 ps to obtain a fully relaxed supercell for thermal decomposition simulation. Then, the relaxed supercell was used as the initial structure for subsequent MD simulations. Then, ReaxFF-lg isothermal-isochoric MD (NVT-MD) simulations were performed for 200 ps at 1800, 2250, 2500, 3000 and 3500 K, respectively, controlled by the Berendsen thermostat. The damping constant was 100.0 fs. Periodic boundary conditions are applied to all simulations. Newton's equation of motion was calculated by using the velocity-Verlet algorithm with a time step of 0.1 fs. An analysis of the fragments was performed with a 0.3 bond order cutoff value for each atom pairs to identify chemical species. The information of dynamic trajectory was recorded every 50 fs, which was used to analyze the evolution details of NONA in the pyrolysis process.

Rate constants analysis
The reaction rate of NONA in the decomposition process can be divided into three stages (endothermic stage, exothermic stage, and generation of gaseous products). This method has already been used to calculate the rate constants of TNT, 8 CL-20, 12 and HNS 14 decomposition. The reaction rate of NONA at each stage can be described by the Arrhenius equation, shown below where t, T, R and a are time, temperature, gas constant and reaction progress, respectively. A is the pre-exponential factor, and E a is the activation energy. The A and E a can be obtained by linearly tting. f(a) is the reaction model, and k(T) is the temperature-dependent rate constant. The linear form of Arrhenius law relates the rate constant to activation energy and temperature, calculated by the following equation The three stages of the decomposition are: Fig. 2 Comparison of the RDF from ReaxFF-lg (300 K, 0 atm) and from experiment for NONA.   (1) The initial decomposition reaction uses the rst-order decay exponent N(t) to t the reduction in the number of NONA molecules from t 0 to t max . t 0 and t max are the times of the initial decomposition and the maximum potential energy, respectively.
where N 0 is the initial number of NONA molecules; t 0 is the time at which NONA begins to decompose; k 1 is the primary reaction rate.
(2) The intermediate decomposition reaction is an exothermic stage. The potential energy (PE) of the exothermic stage changes with time using the decay expression U(t) where DU exo ¼ U(t max ) À U N , U N is the equilibrium value of potential energy when t approaches innity; U(t max ) is the maximum potential energy, and k 2 is the secondary reaction rate.
(3) The nal product of the thermal decomposition of NONA supercell is tted with time to obtain the nal rate constant of product formation.
where C N is the asymptotic value when the amount of the product approaches equilibrium; t p is the time at which the product begins to form; k 3 is the nal product formation rate.

Validation of ReaxFF for NONA
To verify the feasibility of the ReaxFF-lg force eld, we compared the lattice parameters, density and bond length of the relaxed NONA supercell at 300 K and 0 atm with the corresponding parameters obtained by X-ray single crystal diffraction. 6 The calculation results are shown in Tables 2 and 3

Evolution of the potential energy
Evolution of PE is shown in Fig. 3. The variation of potential energy for the system is clearly divided into endothermic (seen in insets of Fig. 3) and exothermic stages. At the initial pyrolysis, the system quickly absorbs heat to increase the potential energy of the system. Subsequent reactions produce a series of intermediate and nal products and release a large amount of heat, thus a signicant decrease in system potential energy. The change of PE of each system is closely related to the chemical bond cleavage and formation in thermal decomposition reaction. The higher the temperature, the shorter the time required for the system to reach the dynamic equilibrium of the chemical reaction.

Evolution of the chemical species
To obtain the decomposition process of NONA, a series of C++ scripts were used to obtain the fragment statistics during the decomposition process of NONA. Fig. 4 shows the change of molecular fragments with time at 2250 and 2500 K. NONA molecules have been decomposed at 11.64 ps and 7.13 ps, respectively. At 2250 K, the initial reactions generate NO 2 and NO occur at 1.96 ps and 3.54 ps, respectively, whereas at 2500 K, the initial generation times of these two species are reduced to 0.87 ps and 3.03 ps. In this respect, the reaction is accelerated by increasing temperature. NO appears slowly aer NO 2 , indicating that NO 2 generation is more likely to occur. The simulations reveal that two possible initial decomposition processes: (i) N-N homolysis and (ii) nitro-nitrite (NO 2 / ONO) isomerization followed by the NO ssion. The number of the abundance species NO 2 and NO at ve temperatures were summarized in Fig. 5 as time goes on the initial stage. More NO 2 and NO are produced as temperature increases, while the difference in the maximum values of NO 2 and NO is decreasing. Especially at 3000 K, the maximum value of NO exceeds that of NO 2 , indicating that the rearrangement of C-NO 2 moiety to C-ONO followed by O-NO homolysis is a thermodynamically more favorable approach. In addition, the concerted HONO elimination is deemed to preferentially occur in the condensed phase energetic materials. 9,16-19 However, it should be noted that the HONO eliminations is less favorable for the initial decomposition process of NONA. The decomposed products of NO 2 , NO and HONO reach their peak values at 9.15, 13.54 and 14.91 ps, respectively, corresponding to the species numbers of 334, 301   This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 5507-5515 | 5511 and 34. The amount of HONO generated is much smaller than those of NO 2 and NO (as shown in Fig. 4). Furman et al. 9 also proposed that the rst pathway of intermolecular H-transfer is possible only if H atoms are available on the ring substituents (i.e., -CH 3 , -NH 2 ), since the cleavage of a benzene ring bound hydrogen requires much more energy. Hence, for NONA with only nitro substituents, the intermolecular H-transfer is not favorable to participate in the HONO eliminations since the numbers of HONO and HNO are so small.

Subsequent decomposition pathway
Aer the homolysis of NO 2 and NO, the simulation results show that there are macromolecular fragment products containing C 18 and C 36 . Stability of terphenyl radicals might be responsible for their detection in the condensed phase. The cyclization reaction of terphenyl derivative occurs through the formation of N-O-C and O-C between the O and NO on the benzene ring and C radical in adjacent benzene ring (Fig. 6). The veand sixmembered rings are the derivatives of oxole and benzoxazine derivatives. In addition, the bimolecular mechanisms, especially for condensed phase, are a crucial decomposition process. However, most of the current researches are focusing on H transfer between bimolecules. 20,21 Fig. 7 shows docking modes for the bimolecular of NONA. Bimolecular is linked together by a C atom in one terphenyl derivatives removed from NO 2 group and an O atom of another terphenyl derivative. The mechanisms of ve initial reactions are illustrated in Fig. 8. Each path is independent in the analysis of reaction mechanism. (1) Path 1: the initiation of thermal decomposition of NONA occurs through C-NO 2 bond cleavage in intermediate benzene ring (t ¼ 4.71 ps), leading to the formation of NO 2 .
Then the another C-NO 2 bond in side benzene ring of NONA breaks immediately aer rst NO 2 dissociation (t ¼ 5.03 ps), leading to the formation of the carbon radical on each benzene ring. These benzene rings fuse with its neighbor nitro groups aer the twist of C-C single bond (t ¼ 5.10 ps). The planes of the two benzene rings approach parallel, and then the O atoms in the nitro group and the C of the adjacent benzene ring form a C-O bond (5.24 ps). A similar reaction process occurs at 7.61 ps, and nally forms a fused ring compound (INT3). (2) Path 2: the N-O in benzoxazine structure breaks leading to the formation of C-O and C-NO (9.34 ps). The CO bonding with external NO 2 radical to form nitrate ester C-ONO 2 (11.67 ps).

Evolution of the gaseous and cluster products
The simulation results show that the nal products of NONA decomposition are mainly N 2 , CO 2 and H 2 O. Here, Fig. 9 shows the evolution of the nal products with time at 2250 and 2500 K. The nal products are produced continuously, and their amounts are nearly equilibrium at 180 ps. The order of the amounts of the nal products are N 2 > CO 2 z H 2 O at 180 ps. In addition to forming gaseous products, some clusters come into being in thermal decomposition of NONA. Fig. 10 shows a snapshot of the NONA supercell at 200 ps for 2250 K, in which the formation of both small and large products, such as N 2 , CO 2 , H 2 O and large clusters, can be observed. As the -NO 2 and -NO bonds cleavage, the C-containing groups continue to polymerize to form large molecular clusters, which lead to a rapid increase in the molar mass of the carbon-containing clusters, followed by secondary chemical reactions. Large clusters are in a dynamic equilibrium by the cleavage of various small moieties and the polymerization of small clusters. The appearance of clusters is an important aspect to understand the properties of energetic materials. Experimental and theoretical researches have shown that carbon-rich explosives can form larger carbon-containing clusters in detonation process. 22-28

Diffusion analysis
To quantitatively describe the interactions between the NONA molecules, the three-dimensional diffusion coefficients for all atoms can be calculated as the derivative of its mean square displacement (MSD) at four temperatures (Fig. 11). At 1800 K, the C, H, O, and N atoms are estimated from the slope of the MSD to be 3.38 Â 10 À9 m 2 s À1 , 3.13 Â 10 À9 m 2 s À1 , 2.31 Â 10 À9 m 2 s À1 , and 3.53 Â 10 À9 m 2 s À1 , respectively. The order of the MSD of other four atoms were N > H > O > C, which is not consistent with the order of atom mass O > N > C > H. This indicates that there are other factors that affect the MSD. On the one hand, N belongs to NO 2 and participates in more reactions, whereas the amount of HONO and HNO are signicantly less in the initial pyrolysis reaction. It is found that the intermediate of NONA retains the H atom in Fig. 6 and 7. On the other hand, a large amount of C atoms participate in the formation of clusters and make its diffusion more difficult. As the temperature increases, the diffusion rate of all atoms increases, and the diffusion rate of H atoms gradually exceeds that of N atoms. This indicates that the temperature has a signicant acceleration on the diffusion coefficient of H atoms. Fig. 12 shows the evolution of PE for the system of 0.9V 0 and 1.1V 0 . When the volume of NONA supercells is reduced, and the corresponding atomic spacing and molecular spacing are also reduced. The equilibrium distance between molecules is reduced, and the repulsive force between molecules is dominant. On the contrary, the increase in the volume of NONA supercell makes the molecule and the atom appear attractive. The curve of PE for the different density of NONA shows an initial sharp increase followed by a slow decrease with time. The system of 0.9V 0 , V 0 and 1.1V 0 have a t max of 4.16, 4.87 and 5.59 ps at 3500 K, respectively, indicating that the increases in the volume of NONA supercells cause the endothermic phase to be prolonged. We subtract t max from the reach equilibrium time of t eq . The time required for the 0.9V 0 , V 0 , and 1.1V 0 systems to reach equilibrium during the exothermic phase were 12.15, 18.04, and 23.73 ps, respectively. This indicates that the increases in the volume of NONA supercells also prolong the exothermal time.

Thermal decomposition of evolution under two densities
To evaluate the degree of NONA decomposition for volumetrically expanded and compressed supercells of 1.1V 0 and 0.9V 0 at 1800 and 2250 K, the variations of remaining HMX and major decomposition products with time are showed in Fig. 13. These ReaxFF-lg MD simulations elaborate that the N-N homolysis and nitro-nitrite (NO 2 / ONO) isomerization followed by the NO ssion dominate the initial stage of NONA dissociation in the system of 1.1V 0 and 0.9V 0 . The decomposition schemes for NONA in 1.1V 0 and 0.9V 0 are similar to that for V 0 system except the amount and formation time of NO 2 and NO.
At 2250 K, NO 2 begins to appear at 5.14, 12.83, and 14.11 ps for the system of 0.9V 0 , V 0 , and 1.1V 0 , respectively, and the corresponding amount of NO 2 were 201, 339, and 454. NO 2 increases more quickly with the increase in density. Similarly, NO occurs at about 10.02, 19.24, and 30.63 ps for the system of 0.9V 0 , V 0 , and 1.1V 0 , respectively, and the corresponding amount of NO were 137, 270 and 354. Density has the same effect on the amount of NO.
The onset of various secondary exothermic reactions leads later to the formation of many nal products, such as N 2 , CO 2 , and H 2 O (Fig. 14). The equilibrium time required for N 2 and H 2 O is shortened to a low-density system of NONA, while the amount of N 2 and H 2 O fragments remains stable. However, the formation of CO 2 is favored at a low-density system of NONA. The lowdensity system of NONA is not conducive to the formation of C clusters, and more C detaches from clusters and forms CO 2 .

Reaction kinetic parameter analysis
To understand the effects of temperature and pressure for the decomposition process of NONA, the reaction kinetic parameters were employed and evaluated for the three different reaction stages.
(1) Endothermic stage: the system costs about 1.4 ps to heat to the target temperature. Hence, for endothermic stage, we should study the reduction process of NONA from 1.4 ps to t max . Table 4 lists the values of t max for the system of V 0 , 0.9V 0 and 1.1V 0 . Using eqn (3), we obtain the initial decomposition rate constants (k 1 ) at 1800, 2250, 2500, 3000, and 3500 K of the initial decomposition of NONA with different crystal density. The k 1 becomes larger with increasing crystal density and/or temperature, indicating that density and temperature can accelerate the decomposition of NONA. The pre-exponential factors and activation energies can be obtained by establishing a linearly tted eqn (2) (Fig. 15). These values are summarized in Table 5. The tted activation energy in the simulation is the apparent activation energy for overall reactions. The E a values (38-41 kcal mol À1 ) are lower than the reported experimental values (51.14 kcal mol À1 ) because some highly reactive radicals and intermediates take part in the reactions in addition to NONA decomposition when the temperatures are in the range of 1800 K to 3500 K. 15 (2) Exothermic stage: the second decomposition rate constants (k 2 ) are obtained by tting the attenuation of PE curve with eqn (4). The results of linear tting, as seen in Fig. 16, are extracted with eqn (2). The data of pre-exponential factors and activation energies for various systems is presented in Table 5. A higher temperature leads to a higher secondary decomposition rate, whereas the greater density results in smaller activation energy. The calculated activation energies agree well with the experimental data from Russian manometric (isothermal) method. 5 Fig. 15 Logarithm of NONA initial decomposition rate constant ln(k 1 ) against (1/T) in range of 1800-3500 K. The solid point and solid line represent calculated and linearly fitted values for three density systems. Orange, magenta, and dark cyan represent the systems of 0.9V 0 , V 0 , and 1.1V 0 , respectively.   (3) Generation of gaseous products: the gaseous products are mainly N 2 , CO 2 , and H 2 O in NONA decomposition process. The rate constants (k 3 ) of the gaseous products are obtained by tting these data of gaseous products into eqn (5). As shown in Fig. 17, the k 3 is increased from 1800 to 3500 K, and CO 2 with faster formation rate is observed at the NONA decomposition stage as the temperature increases.

Conclusions
We simulated and analyzed the thermal decomposition reaction of NONA crystals at different temperatures by using ReaxFF-lg reactive force eld. The primary initial reaction pathways are obtained by analyzing the evolution of the chemical species. We propose that two distinct initial decomposition mechanisms are the homolytic cleavage of C-NO 2 bond and the rearrangement-homolysis pathway of C-NO 2 / C-ONO followed by the NO ssion. The fused ring clusters and bimolecules were found in the subsequent decomposition of NONA. The latter is a thermodynamically more favorable pathway than the C-NO 2 homolysis at high temperature. The identication analysis of nal products showed that the gaseous products were CO 2 , N 2 , and H 2 O. The amount of CO 2 will be more favorable energetically for the system at high temperature or low density, and the clusters are a favorable growth pathway in low temperatures, and this process was further demonstrated by the analysis of diffusion coefficients. Crystal density has a signicant acceleration on the decomposition of NONA through increasing reaction kinetic parameter (k 1,2 ) and reducing activation barriers.

Conflicts of interest
There are no conicts to declare.