Temperature and pressure effect on the HMX/Graphene mixture system thermal decomposition via ReaxFF molecular dynamic simulation

Study the detail mechanisms of energetic materials thermal decomposition process at high temperature, is helpful to understand the detail reaction information and mechanism, which are foundation for understanding reactivity of energetic materials, design of mixed explosives and safety. In this work, new model and method effect of temperature and pressure for graphene (Gr) based HMX crystal were determined by ReaxFF molecular dynamic simulation. The thermal decomposition process, HMX molecules absorption on Gr surface of perfect HMX crystal, different direction crystal planes and HMX/Gr mixed models were studied at high temperature and pressure. The different congurations of HMX molecules adsorbed on graphene surface in the mixed system were conrmed by theoretical calculation method. We observed 3, 5 and 3 HMX congurations absorption on the graphene in the (001)/Gr, (010)/Gr and (100)/Gr range from normal pressure to 31 GPa, respectively. The time-dependent curves of fragments evolution, intermediates and pyrolysis products were analyzed. The results show that the rate constants of HMX molecules during thermal decomposition process can be signicantly affect by adding graphene with HMX crystal. Gr is the most obvious inhibition effect on the thermal decomposition reaction of (010) surface of HMX/Gr, which indicates that graphene has a coupling between the temperature and thermal anisotropy effect, because of the NO 2 functional groups steric hindrance and the interaction between graphene and HMX molecules. Gr also affects the initial reaction pathways of homolytic cleavage of N-NO 2 bond forming nitro radical and HONO through the C=O, C-OH and C-OC bonds formation on graphene surface. in mixture explosive is the based approach for reveal the decreasing sensitivity and complex detonation reactions. ReaxFF-lg molecular dynamics simulations (RMDS) [14] , which has been proposed by Liu et al, has been studied the reactions of HMX [15] , ICM-102 [16] , TATB [17] , CL-20 [18] , and CL-20/TNT. The ne results are compared with quantum mechanics calculations and experiments. adsorbed on the graphene surface as shown in the gure 8b. This change shows that the interaction between HMX and graphene on different crystal planes is different. To further study the amount of HMX molecules adsorbed on the surface of graphene at different pressures, taking the supercell of HMX/Gr at normal pressure and room temperature as the starting point, the structure of supercell at 300K and different pressures was calculated. The number distribution of HMX molecules adsorbed on the surface of graphene was analyzed. The results show that the nitro functional groups on the A-axis and E-axis of single molecule HMX can be adsorbed on the surface of graphene. At the pressure range between 1GPa to 31GPa, HMX molecules with two different congurations are adsorbed on the graphene surface along the (001) surface. In one conguration shown in gure 9a, the one oxygen atom in nitro functional group is adsorbed on the graphene surface, which can be descripted as C 4 H 8 O 7 N 8 O(ads), and in the other conguration shown in gure 9b, one N atom and two O atoms in the nitro functional group are adsorbed on the surface with the chemical structural formula of C 4 H 8 O 6 N 7 N(ads)O(ads) 2 . In the same pressure range, along the (100) plane shown in gure 9c, the rst conguration of HMX molecule absorption on the Gr is that the N atom in the CN ring constituting of HMX ring bonding on the graphene surface with the chemical formula C 4 H 8 O 8 N 7 N(ads), the second conguration is that the one of O atom in nitro functional group in the HMX molecule adsorbed on the graphene surface shown in gure 9d, and two molecules of this conguration with the chemical formula of C 4 H 8 O 7 N 8 O(ads) are located above and below the graphene surface respectively.


Introduction
Energy materials (EMs) like explosives, pyrotechnics, and propellants play an important role in many elds, etc. aeronautics, and defense industry [1] . Because of the high density, energy and tension, EMs are widely used applications ranging from engineering blasting, automobile air bags, oil exploitation to weapon design and rocket propellants preparation [2] . Unfortunately, the pure EMs have very high sensitivity, and should be prepared to mixed components by adding some binder, additive, desensitizer, which will be improving the thermal decomposition mechanism, tensile mechanical property and sensitivity. Considering numerous excellent properties of graphene (Gr) [3] , such as large speci c surface area, high thermal conductivity, and strong conductivity, the family of graphene-based materials (GBMs), such as Gr, graphene oxide (GO) and reduced graphene oxide (rGO) have become a focus in the application of EMs/GBMs [4] . Recently, GBMs are proposed for additive and desensitizer and are expected to replace the graphite and other compounds in the TATB based [5] and CL-20 based [6][7] plastic-bonded explosive (PBX).
These studies [6][7] have been sought to reveal the reduce mechanical sensitivity and thermal stability reasons when GBMs adding to CL-20 crystal. The results showed that high speci c surface areas of graphene and its derivatives make GBMs have large buffer space when the mixture of EMs/GBMs encounter mechanical shock and impact stimulation. The thermal decomposition process was also studied at low temperatures. Due to the good thermal conductivity of graphene, the heat ow will be transferred from graphene to the PBX surface quickly, resulting in the decomposition of PBX ahead of time and reducing its thermal stability; while the rich oxygen-containing functional groups on the surface of GO will rst decompose and release OH radicals under high temperature, and transfer part of the heat ow to PBX, which also leads to PBX decomposition ahead of time; the thermal conductivity of rGO is normal, and most of the oxygen-containing groups on the rGO surface have been reduced, the PBX thermal stability is not affected adding rGO powder.
The activation energy in the HMX/GO composites is signi cantly increased when the GO add to HMX crystal, and the thermal stability of HMX explosive is increased [8] because of the ame retardant effect, which is due to use many metal salts in the GO preparation process [9] . The desensitizing effect of GO sheets are better than fullerene and Carbon nanotubes (CNTs). Niu et al. [10] have been prepared the HMX/rGO/graphite via an in situ chemical reduction coating method. These results suggest that GO sheet and graphite can be used as co-desensitizer in HMX explosive. In a word, thermal reactions at high temperature and pressure depend heavily on preparation method under different conditions. In fact, the results of thermal stability of HMX and CL-20 are opposite when the GO adds to the HMX and CL-20. On the other hand, the detail information of thermal decomposition mechanisms of Gr and GO in EMs/GBMs are still elusive under extreme conditions because of lacking experiment data.
Theoretical methods, such as quantum mechanics and molecular dynamic simulation, were reported to investigate the mechanisms of EMs/GBMs and clarify the effects on the sensitivity of GBMs adding to the EMs. For researching the mechanisms of functionalized graphene sheets (FGSs) adding and dispersing into the liquid nitromethane (NM), ab initio molecular dynamics simulations were calculated by Liu et al [11] . The results have been shown that the NM molecules and thermal decomposition products are greatly accelerating with the FGSs, forming water, nitrogen, and carbon dioxide. The reactions pathways, such as exchange of protons or oxygens, have been happened in the system and the FGSs play a catalyst role in the thermal decomposition process. Zhang el al. [12] have been determined thermal stability of GO, ATRZ and GO-ATRZ, by using the ReaxFF-lg force eld molecular dynamics simulations. because of the GO strong space effect, the activation energy of GO-ATRZ can increase 16.1 kJ/mol comparison than the pure substance of ATRZ, indicating that the GO improve the thermal decomposition property of ATRZ. GO can be combined with ATRZ as a desensitizer candidate. The interaction interface of FOX-7-GO mixture composites [13] , which is an ideal prototype system, were studied by dispersion-corrected density functional approach. The results revealed that the interfacial charge transfer from FOX-7 to GO and hydrogen bonds are two important issues, leading to relatively strong interaction between FOX-7 and GO.
In the present work, molecular con gurations of single HMX molecule absorption on the Gr surface and HMX/Gr solid systems thermal decomposition at high temperature were studied using RMDS. This work focused on the thermal behaviors of perfect HMX crystal, HMX crystals with different crystal planes and HMX-Gr mixture systems. A deeply understanding of Gr effect in complex explosive chemical mechanisms and estimation of the different factors in uence on the degradation pathways are essential for the applications of preparation, storage, transportation and use of mixed explosives utilization.

Model Building And Computational Approach
The perfect HMX crystal structure was obtained from experimental X-ray data [20] . Experiments showed that the shock sensitivity of HMX is different along different directions [21] . The compression results of β-HMX obtained from the diamond anvil cell have been shown that the lattice parameters of signal crystal are anisotropy [22] . In addition, the results of MD based on the COMPASS showed that the interactions between different HMX crystallographic plane and same binder are different [23] . The ReaxFF module of AMS software package with the 2019.304 version was carried in this work.
To research the crystal plane and Gr effect on the thermal decomposition of HMX/Gr, we built systemic models. Along x, y and z directions, 4×2×3 perfect HMX supercells and 8×8×1 graphite supercells were employed with expanding with unit cell of HMX and graphite. Then (001), (010) and (100) structures were obtained, by cleaving perfect HMX supercell along (001), (010) and (100) crystal planes, respectively. So, the total HMX molecules of different crystal planes were 48 (1344 atoms). The initial lattice parameters of perfect HMX, (001), (010), and (100) crystal planes with HMX and HMX/Gr, which can be replaced as P-HMX, P-001, P-010, P-100, M-001, M-010 and M-100 abbreviated forms, list in Table 1.  Energy minimization without the optimization lattices was applied to release the internal structures stress of seven supercells. Subsequently, the two step RMDS, which was 25 ps isothermal-isobaric (NPT) MD (0.25 fs timestep) and 25 ps isothermal-isochoric (NVT) MD (0.25 fs timestep), were calculated to relax the seven structures at room temperature and pressure range from 1atm to 31GPa with the Nose-Hoover chains (50 fs damping constant) and NPT NHCP anisotropic barostat (500 fs damping constant). Finally, based on the NVT ensemble, the total 100 ps RMDS (0.1 fs timestep) were calculated to research the thermal decomposition at 3000K with timestep of 0.1fs. By analyzing the evolution of RMDS, behaviors of HMX molecules and Gr molecule were studied to reveal the Gr effect on the mixture explosive.

Validation of ReaxFF-lg force eld
Con rming the applicability of the force eld to the calculation system, is the basic work of molecular dynamics simulations. Many works were studied the applicability of the ReaxFF force eld to the HMX and other EMs, reporting that the force eld can well describe the equation of state of HMX [14] and the reaction process under different loading conditions [15][16][17][18][19] , but the applicability of the ReaxFF force eld to the mixed system composed of HMX and graphene has never been studied.
To verify the ReaxFF-lg force eld tting to HMX and graphene, the seven HMX systems of 25 ps NPT-MD and 25 ps NVT-MD at 300 K and normal pressure were performed. Our results of the lattice parameters and densities of P-HMX unit cell are summarized in Table 2 and compared to the experiment and other theory methods. The density result is 1.808 g/cm 3 , which are agreed very well with earlier experiment [20][21][22][23][24][25] , DFT-D 2 [24] .
Expt. [20] 6.54 11.05 8.7 90 124.3 90 1.894 Expt. [21] 6.537 11 The lattice parameters of different HMX crystal planes and HMX/Gr mixture systems are summarized in Table 3, which were predicted by ReaxFF-lg MD simulations. The densities, which are the mainly factor effect on the thermal decomposition process, have been the same values scales about 1.8 g/cm 3 . As shown in Table 3, the values of a axis, b axis, α, β and γ become little changes, while the value of c axis change abruptly, indicating that the c axis is easily compression. Table 3 Lattice parameters of different crystal planes and HMX/Gr Noteworthy, the interaction between HMX/Gr mixture system can be expressed by binding energy, which has a positive correlation with the compatibility and stability. That is to say, the binding energy is bigger with the greater compatibility and stability of the mixture explosives. The binding energy can be calculated as: Where E bind is the binding value in the HMX/Gr system between HMX and graphene, and E interaction is the interaction energy between the HMX and graphene in HMX/Gr system. E total is the potential energy of HMX/Gr system, E HMX and E Gr are the potential energy of HMX and graphene in the HMX/Gr system with different crystal plane, which can be obtained by removing graphene and HMX. The largest and smallest binding energy systems prediction from RMDS studies are M-100 and M-010. He and coworkers [25] , investigated the interaction between different directions of HMX and Hydroxyl-terminated poly-butadiene (HTPB), which is a typical binder in plastic bonded explosives. They proposed that formation energy sequence between the HMX crystal and HTPB with different directions is the (010) (001) (100). The interaction sequences of binding energy between the HMX directions with graphene and HTPB mean that the strong interaction intensity, like hydrogen bonds, is forming between different crystal surfaces and graphene, especially (010).

Evolution of HMX molecules
Decomposition reactions of solid HMX crystal at high temperature produced many small fragments and products by using RMDS [14] and experiments [10] . Fig.   3 compares the decay number of HMX molecules in the thermal decomposition reactions of HMX and graphene mixture systems at 3000K. HMX molecules have the same trends in the initial stage during the 0-0.3ps, then with the different rate constants, HMX molecules are gradually disappeared with several initial reactions pathways during 0.31 ps -1.8 ps.
By approaching the rst order phase transition nucleation theory and interatomic interaction model Monte Carlo simulation [26] , the result of graphene melting temperature is 4510 K, which is higher than the results of graphite. Especially, it is reported that the chemical bonds broken, and formation have been observed by using ab initio molecular dynamics calculations (AIMD), suggested that the graphene cannot melt when the system be heated up to 4500 K, leading a quasi-2D liquid state [27] . It means that in this work, the graphene is still no reactive during the thermal decomposition of HMX/Gr systems at 3000K. It can be con rmed that graphene and crystal face affect the decomposition process of HMX molecules in different calculation systems.
The rst stage of rate constant, which has been proposed to be critical component in the thermal decomposition, can be calculated by the number evolution of HMX molecule using the st order reaction rate equation [18] as where N(t) and N 0 are the number of HMX molecules at time t and t 0 , t and t 0 are the molecular dynamic simulations time. According to the equation, Fig. 4 shows the rate constant k 1 of HMX molecules, it can be suggested that the crystal plane and graphene have both been changed the rate constant compared with pure HMX. Parameters of HMX molecules reactions in P-001, P-100, M-001, M-010, M-100 are both smaller than in P-HMX.
The rate constants of HMX molecules in HMX/Gr is lower than with the same crystal plane of HMX. The crystal plane of (010) can accelerate the thermal decomposition reactions of HMX molecules in P-010, whereas the rate constant of k 1 is abruptly decreased when the graphene adds to the same crystal plane of HMX. The difference values between crystal plane and HMX/Gr of (001), (010) and (100) are 0.68 ps −1 ,5.17 ps −1 and 0.69 ps −1 , respectively. Thus, the graphene molecules have been shown great desensitized agent role in the HMX/Gr mixture systems, especially in the (010) crystal plane.
HMX have been exhibited the anisotropic responses under shock loadings by experiments and theory. Menikoff et al. [28] measured the wave pro les of HMX, which is an elastic-plastic material, and calculated nonlinear and transient wave behavior using a rate-dependent elastic-plastic model. It is reported that the effective yield strength values of (011) and (010) are 0.18 GPa and 0.31 GPa. Ge et al. [29] performed AIMD conjunction with multiscale shock technique of HMX under different lattice vectors. Compared with shock wave propagation along the lattice vectors of a, b and c, sliding rate of lattice vector b is the smallest, suggesting the lattice vectors b is main reason leading to anisotropic of HMX under shock loading. Interestingly, the results in our work suggest that graphene also has anisotropic reaction responses with the oriented of HMX crystal under thermal loading with different planes.

Important intermediates fragments
Decomposition pathways, intermediates and products of signal molecule and crystal structure under different loading have been reported by using DFT [30][31][32][33][34][35] , AIMD [29,36,37] and RMDS [38] . Initial pathways of HMX molecule are homolytic cleavage of N-NO 2 bonds producing the NO 2 radical, migration of H atom from the methylene forming HONO, breaking four C-N bonds in HMX molecule to form CH 2 N 2 O 2 , isomerization of NO 2 and then C-N bond fracture occurs to release N 2 O 2 fragments, C-H bond dissociation. The simulation of theory results [26][27] have also been suggested that the graphene molecule is not reactive lower than 4510 K. Thus, in this work of HMX-Gr systems, the origin fragments of intermediates are from HMX molecules.
The evolution of fragment number with simulation time is a commonly used method for studying the mechanisms of EMs under different conditions. We calculated the fragments in seven models, using the bond order with 0.3 at 3000K. The results mean that NO 2 , HONO, N 2 O 2 , HNO 3 and NO are the main intermediates in the thermal decomposition of pure HMX and HMX-Gr systems, which are completely consistent with other calculations [36][37][38] and experiments results [39][40] . The time evolution of NO 2 radical fragments is presented in Fig. 5  The rate constants of intermediates were obtained using equation 2 in the simulation time range from maximum formation populations to zero at 3000K, list in Tab. 5. As can been seen, the rate constants of HONO and N 2 O 2 are both increased in the HMX/Gr compared with P-HMX. While the NO rate constants are not changed in the thermal decomposition process in the seven system. Whereas NO 2 rate constants show complex behaviors in the different system.
Moreover, the rate constants of NO 2 in P-010 and P-100 changed abruptly compared with P-HMX. This trend has been changed when the graphene mixtures with HMX in the M-001 and M-010. Crystal planes and graphene molecules do not affect the initial decomposition pathways of HMX molecules in P-HMX and other systems, while the total amount populations of intermediates and the rate constants in HMX/Gr mixture systems are abruptly changed compared with the perfect HMX and different crystal planes supercells.

Fragments formation on graphene surface
The absorption properties of HMX molecules on Gr in HMX/Gr system at room temperature and pressure range from 1atm to 31GPa were obtained. We found that there are 3, 5 and 3 different HMX signal molecule con gurations bonding on the Gr. At normal temperature and 1atm, the HMX molecules bonding on the graphene can be described as two types as shown in Fig Figure 11 shows the number of HMX molecules adsorbed on the graphene surface in the HMX/Gr mixed system along different crystal planes at atmospheric pressure to 31 GPa. It can be seen that, except for 19 GPa, 21 GPa and 23 GPa, in the remaining pressure range, the number of molecules adsorbed on the graphene surface by HMX molecules along (010) plane is less than that adsorbed on the graphene surface along (001) plane and (100) plane.
The chemical process in HMX/Gr in condensed phase at 3000K is more complex than in pure HMX supercell, the initial reaction pathways of HMX molecules experienced difference. Moreover, the fragments in HMX/Gr are obvious changed with crystal plane. Table. 6  the action of high temperature, the fragments generated by HMX explosive react with the C atoms on the surface of graphene to form different kinds of products, which further reduce the reaction rate constants of HMX molecules, prolong the decomposition time of HMX molecules, and improve the sensitivity of HMX explosive after adding graphene.  Figure 12 shows the value changes of O (ads) in HMX/Gr with simulation time. As can been seen, three system of O (ads) fragments have the same trends. The results of fragments evolution on the graphene surface in HMX/Gr with different crystal plane, show that the initial decomposition pathways of HMX/Gr are the C atoms in graphene react with the intermediates, forming C = O bonds, C-OH bonds and C-CO bonds between C atoms and nitro radicals, hydroxyl radical and CO at 3000K. Comparing the products distribution of three different crystal planes, it shows that the type and quantity of graphene products are the largest along the (010) crystal plane, and the decomposition reaction of HMX/Gr at high temperature is directional and selective, which provides theoretical guidance for the formulation design and safety improvement of graphene based explosives.

Final products of small fragments
Final products, like N 2 , H 2 O, CO 2 , and CO, in the thermal decomposition process were analyzed using bond order cutoff of 0.3. The products evolutions during the simulation time are presented in Fig. 14. It is suggested that small products are appeared quickly at the beginning of the simulation between 0 to 20 ps, then kept a plateau status between 20.1ps to 100 ps. This implies that graphene and crystal do not affect the kinetic parameters of the nal products.

Conclusions
In this work, ReaxFF molecular dynamic simulations were performed to study the initial thermal decomposition reactions of HMX and graphene mixture systems with different crystal planes. By analyzing the lattice parameters, binding energy, rate constants, evolution of intermediates and products fragments, the thermal stability of HMX is strongly dependent on the crystal plane and graphene; sensitive for thermal stability of HMX and HMX/Gr with the (010), but less sensitive for thermal stability of HMX/Gr with the (001) and (100)  Con icts of interest/Competing interests: I declare that I do not have any commercial or associative interest that represents a con ict of interest in connection with the work submitted.
Availability of data and material: The data used to support the ndings of this study are available from the corresponding author upon request.
Code availability: The ReaxFF module of AMS(Amsterdam Modeling Suite) software package with the 2019.304 version was carried in this work. Figure 1 The   Reaction rate constants of HMX molecules at 3000K Figure 5 Evolutions of NO 2 radical fragments at 3000 K Evolution of intermediates at 3000 K Number of HMX molecules bonding on the Gr in mixture systems from normal pressure to 31GPa Evolutions of HMX/Gr HO (ads) and OC(ads) at 3000 K Figure 14 Evolutions of HMX/Gr nal products at 3000 K

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. abstractgraphic.docx