Molecular dynamics simulation of initial thermal decomposition mechanism of DNTF

In order to understand the thermal decomposition characteristics of 3,4-Bis(3-nitrofurazan-4-yl)furoxan (DNTF), the thermal decomposition reaction of DNTF at 300–4000 K temperature programmed and constant temperatures of 2000 K, 2500 K, 3000 K, 3500 K, and 4000 K was simulated by ab initio computational molecular dynamics method. The thermal decomposition mechanism of DNTF at different temperatures was analyzed from the aspects of product evolution, cluster, potential energy curve, and reaction path. The analysis of products shows that the initial small molecular products are NO, NO2, CO, CO2, and N2, and the final small molecular products are CO2 and N2. In the early stage, the ring-opening reaction of furoxan in DNTF structure is the main trigger reaction, and the C–C bond is broken at the initial stage of reaction. The carbon chain structure produced by decomposition forms various cluster structures in the form of C–N bond. In addition, it was found that temperature significantly affects the decomposition rate of DNTF, but does not change its initial decomposition path.


Introduction
The detonation process of energetic materials is a complex process involving both physical and chemical interactions and often occurs in a very short period of time, which make it hard to understand the reaction kinetics in the study of the detonation mechanism of energetic materials. The initial decomposition mechanism of energetic materials at high temperatures is closely related to their detonation phenomena and shock initiation characteristics. Through the study of the thermal decomposition mechanism of energetic materials in the process of detonation, it can help to reveal the energy release law of explosives and also help to provide theoretical guidance for the study of the thermal stability and sensitivity of explosives.
High energy density materials (HEDM) have become a hot topic of research in the field of energetic materials due to their excellent detonation performance. As a new type of potential HDEMs, furazan compounds have received a lot of attention from scholars at home and abroad because of their high nitrogen content [1][2][3]. 3,4-Bis(3-nitrofurazan-4-yl)furoxan (DNTF) is a hydrogen-free compound containing furazan and furoxan. Because of its high density, low melting point, and better detonation velocity and heat than HMX, it has attracted much attention in the field of energetic materials [4][5][6]. Kotomin et al. [7] proved that the detonation performance of DNTF is close to PETN. Ren [8] studied the fast thermolysis process of DNTF at different experimental temperatures and found that the ratio of NO and NO 2 in the gas phase products of DNTF pyrolysis gradually decreased with the increase of temperature. Sinditskii et al. [9] considered that the first stage of thermal decomposition of DNTF is the destruction of furoxan ring, and the second stage is the decomposition of nitrofurazan fragments at higher temperature, and its thermal stability is similar to HMX. Li et al. [10] studied the thermal decomposition characteristics of DNTF by differential scanning calorimetry (DSC) and found from DSC curve that the thermal decomposition process of DNTF is divided into two stages. In addition, Gao et al. [11] compared the thermal decomposition behavior of DNTF at different heating rates, and the results showed that the initial decomposition temperature of DNTF increased with the increase of heating rate. To sum up, the research mainly explores the thermal decomposition characteristics of DNTF and calculates the relevant kinetic parameters through the experimental means such as rapid temperature transition in situ pool, temperature jump Fourier transform infrared spectroscopy (T-Jump/FTIR) technology, and differential scanning calorimetry, but no relevant research has been done on the initiation reaction, intermediate products, and reaction path during the thermal decomposition of DNTF.
Lindsey et al. [12] used a machine learning approach to modify the standard density functional tight binding (DFTB) models. An approximate quantum mechanical model was obtained that can be used to describe the DNTF. The model is able to simulate the early decomposition of materials under extreme conditions on a longer time scale. They used the modified model to simulate the thermal decomposition of DNTF at high temperature and obtained the small molecule products at the early stage of DNTF thermal decomposition. In addition, they found that the large C x N y O z species produced by the decomposition of DNTF under shock could be the precursors of the experimentally observed carbon condensate. These findings have important guiding significance for understanding the decomposition mechanism of DNTF under extreme conditions. However, they did not analyze the formation mechanism of products and clusters of DNTF thermal decomposition and their evolution regulation with time. A complete study of the initial decomposition path of DNTF at high temperatures is also lacking. Therefore, it is necessary to investigate the thermal decomposition mechanism of DNTF such as initial decomposition path, products, and clusters at different temperatures. These studies not only help to understand the decomposition mechanism of DNTF under high temperature conditions but also enable to recognize the sensitive properties of DNTF at the molecular structure level.
Recently, Li et al. [13] used the efficient approximate density functional theory to study the decomposition mechanism of DNTF at the temperature programmed of 300-3000 K and 1800-2500 K constant temperatures. It is found that the thermal decomposition of DNTF is caused by the cleavage of furoxan ring. Important intermediate products, such as NO 2 , NO, and CO, were found during thermal decomposition. These studies are essential to understand the mechanism of DNTF thermal decomposition. However, the evolution mechanisms of final products and clusters during the thermal decomposition of DNTF have not been solved in this work. In addition, the thermal decomposition mechanism of DNTF at higher temperatures (close to detonation temperature) has not been explored. Therefore, it is necessary to further investigate the thermal decomposition mechanism of DNTF. Understanding the cluster structure of DNTF detonation process not only contributes to the development of thermodynamic calculation model but also can improve the change law of DNTF under the action of heat [14].
With the development of computational science, molecular dynamics and quantum chemistry have become effective methods to analyze the thermal decomposition of energetic materials [15][16][17]. The ReaxFF [18] and ReaxFF-lg [19] reaction force field molecular dynamics can simulate the complex chemical reaction changes in large systems over long time scales and can accurately describe the breaking and formation of chemical bonds, so it is widely used to study the thermal decomposition mechanism of energetic materials such as TNT [20], CL-20 [21,22], RDX [23,24], and TATB [25]. Although the ReaxFF reaction force field is optimized by hydrocarbon reaction force field training set, it has some restrictions for specific molecular structures. In the previous study, after optimizing the unit cell of DNTF by using this force field, it was found that the ring-opening of furoxan ring of DNTF molecule occurred, and the optimized unit cell parameters were quite different from the experimental values. Therefore, it is considered that this force field is not suitable for the simulation of DNTF molecules.
Compared with classical molecular dynamics, ab initio computational molecular dynamics, first proposed by Car and Parrinello in 1985, equates a polyatomic system to a multi-particle system composed of nuclei and electrons [26]. The density generalized function theory (DFT) approach based on ab initio algorithms does not require the provision of force field parameters; it describes the properties of the system in terms of the density distribution function of electrons and solves the Schrödinger equation by the Kohn-Shan equation to obtain the interatomic forces. Due to the higher computational accuracy of DFT compared to ab initio calculations based on Hartree-Fork self-consistent-field calculations, it has been successfully applied to study the thermal decomposition of explosives [27][28][29].
In this paper, VASP software package [30] is used for ab initio computational molecular dynamics simulation, which is a program for quantum mechanics-dynamics calculation using pseudopotential [31] and plane wave basis set [30]. The program adopts the first-principles calculation and iteratively solves Kohn-Shan equation based on density functional theory, which can make the properties of the simulated system closer to reality. Through the analysis of VASP output files, the initial decomposition path, intermediate products, potential energy evolution, the type and quantity of cluster molecules, and the evolution trend of species number under different temperature conditions were obtained. In addition, the influence mechanism of different temperature conditions on the thermal decomposition of DNTF was revealed.

Simulation method
The initial structure of DNTF used in this paper is derived from Cambridge Crystal Database (CCDC) [32]. The single cell of DNTF was first optimized by Gaussian 09 at the calculated level of M062X-D3/6-311 + G (d, p) to obtain the molecular structure of DNTF at the lowest energy. The bond length and Wiberg bond order of DNTF were analyzed. By comparing with the standard bond length and Wiberg bond order, we can judge the chemical bond which may be broken during the initial decomposition of DNTF.  The 2 × 1 × 1 DNTF unit (176 atoms in total) was optimized by VASP software package, and the canonical ensemble (NVT) was used to relax supercell 2 ps at 300 K, as shown in Fig. 1. The optimization process includes two nested processes: electronic iteration and ion relaxation. In the process of electronic iteration, the truncation energy of plane wave base is 240 eV, and the total energy and structural parameters have good convergence. Electronic selfconsistent convergence criterion is 1.0 × 10 −5 eV. Conjugate gradient method (CG) is used to calculate ion relaxation. The convergence standard of interatomic interaction force is − 1.0 × 10 −2 eV, the maximum number of steps of ion relaxation is 100 steps.
The lattice parameters under the minimum binding energy and equilibrium structure are obtained, and the optimized and relaxed atomic position parameters are used as the original parameters of POSCAR file. Molecular dynamics simulation is carried out by using periodic boundary conditions, in which the electron exchange correlation function is in the form of Perdew-Burke-Ernzerh (PBE) in generalized gradient approximation (GGA), and the energy in Brillouin zone is integrated with the point grid accuracy of 1 × 1 × 1. The purpose of choosing periodic boundary conditions is to make the forces on the boundary atoms in the supercell more comprehensive, thus, avoiding systematic errors caused by finite boundaries. The plane wave base truncation energy is 550 eV, and the electron self-consistent convergence criterion is 1.0 × 10 −5 eV, and the maximum number of electron self-consistent cycles is 60 steps The electronic iterative solution adopts residual minimization scheme (RMM) algorithm and uses Nose-Hoove thermostat to control the temperature. Dispersion correction was performed using the DFT-D3 method with Becke-Jonson damping. The molecular dynamics simulations of DNTF were carried out at 300-4000 K temperature programmed and 2000 K, 2500 K, 3000 K, 3500 K, and 4000 K constant temperatures for 10 ps using canonical ensemble (NVT). Then the thermal decomposition mechanism of DNTF under aforesaid temperature was analyzed. The choice of the temperature programmed conditions is to accelerate the collision between molecules and the thermal decomposition of DNTF, thus, shortening the simulation time. The constant temperature of 2000-4000 K is chosen because the explosive temperature of energetic materials is higher during detonation, and the high temperature environment is more in line with the thermal decomposition conditions of DNTF during detonation.

Bond order analyze
The DNTF molecule was optimized by Guassian09 software at the calculated level of M062X-D3/6-311 + G(d, p) as shown in Fig. 2. The corresponding bond lengths and Wiberg bond orders are shown in Table 1. From Table 1 we can see that the bond orders between O2-N11, N16-C22, and N15-C17 are 0.943 Å, 0.903 Å, and 0.897 Å, respectively, all of which are smaller than the standard single bond order. In addition, the bond lengths of N16-C22 and N15-C17 are larger than other N-C bonds in the DNTF molecules. The bond length of O2-N11 is significantly larger than other O-N bonds in the DNTF molecule. Therefore, it is tentatively judged that these three bonds are more likely to break during the thermal decomposition of DNTF.

Feasibility analysis
The 2 × 1 × 1 DNTF cell (8 molecules, 176 atoms in total) was optimized by VASP software package, as shown in Fig. 1. The electronic exchange correlation function is in the form of PBE in generalized gradient approximation (GGA), and the energy of Brillouin zone is integrated with the point grid precision of 1 × 1 × 1 k points. The plane wave base truncation energy is 240 eV, and the electronic selfconsistent convergence criterion is 1.0 × 10 −5 eV. Then, the optimized supercell was equilibrated at 300 K for 2 ps to release the stress in the supercell. Table 2 summarizes and compares the bond length of optimized and relaxed DNTF structure with that of original cell molecules. It can be seen from Table 2 that the maximum error of bond length of DNTF molecule is only 5.50%, which indicates that the structure of DNTF simulated by this method is basically consistent with the experiment. Therefore, it is reasonable to use this pseudopotential to describe DNTF molecule. The optimized and relaxed structure was analyzed as the initial structure of molecular dynamics simulation. The plane wave base truncation energy is 550 eV, and the step size is controlled to 1 fs, and the same k-point distribution is adopted. The structures were subjected to molecular dynamics simulations for 10 ps at the temperature programmed of 300-4000 K and constant temperatures of 2000 K, 2500 K, 3000 K, 3500 K, and 4000 K, respectively, to compare the thermal decomposition mechanism of DNTF at different temperature conditions.

Small molecular products
Analysis of the trajectory files shows that the intermediate small molecule products of the thermal decomposition of DNTF include NO, NO 2 , N 2 O, CO, CO 2 , and N 2 . Figure 3 demonstrates the evolution curves of DNTF products at different temperatures. Figure 3 demonstrates that the higher the system temperature is, the greater the generation rate of NO and NO 2 is, and the earlier the peak time of evolution curve of NO and NO 2 arrives. In addition, the maximum content of NO at different temperatures is much higher than that of other products. This is consistent with the simulation results of Li [13] at 2500 K. It is worth mentioning that intermediate products such as C 2 N 2 O 2 and C 2 N 2 O were observed by Li [13], and the amount of both decreased with increasing temperature. However, these products were not found in our research. This may be due to the difference in temperature as well as the number of molecules in the system. As the reaction proceeds, NO 2   significantly compared with other temperature conditions. This is because the increase of temperature will promote the thermal decomposition reaction of DNTF and accelerate the fracture of C-C bond in the DNTF molecule. From the evolution curve of CO at a constant temperature of 4000 K, it can be inferred that CO is also involved in the reaction as an intermediate product.
The initial generation of NO 2 fragments is due to the denitration of DNTF molecules, while the initial decomposition of DNTF has different reaction paths besides denitration. The shed NO 2 fragments will also react with other small molecules, which limit the amount of NO 2 in the decomposition process under 10 and gradually disappears with the progress of the reaction. However, CO 2 and N 2 still show an upward trend at 10 ps, which suggests that CO 2 and N 2 are the final small molecular products of DNTF decomposition. The evolution of decomposition products of DNTF at constant temperature of 3500 K and 4000 K proves this conjecture. Lindsey [12] found that the small molecules of DNTF in the early stage of thermal decomposition at high temperature were mainly N x O y , but no "oxygen trapping" ion NO 3 was observed. The small molecule products in the late stage of thermal decomposition are mainly CO x and N 2 . This is basically consistent with our research results. At lower temperature and in the initial stage of reaction, the small molecular products of DNTF thermal decomposition are mainly NO and NO 2 . At higher temperature, the main intermediate products of thermal decomposition of DNTF are NO, NO 2 , and CO. This is probably because the high temperature accelerates the decomposition rate of DNTF, which causes a large number of C-C bond breaks in the DNTF molecule. The detached carbon-containing small molecule fragments undergo a series of reactions to produce CO. Comparing and analyzing the evolution trend of small molecular products of DNTF decomposition under temperature programmed and constant temperature, it can be preliminarily inferred that temperature will affect the thermal decomposition rate of DNTF, but has no significant effect Fig. 6 Structure of the largest cluster molecule at different temperatures on the evolution trend of its initial thermal decomposition small molecular products. Figure 4 illustrates the variation curve of cluster number during the thermal decomposition of DNTF. The clusters in this paper are structures with C atoms greater than 6. Since almost no clusters were produced within 10 ps at 2000 K, subsequent analyses of clusters were ignored. It can be seen from Fig. 4 that during the initial decomposition process of DNTF at constant temperatures of 2500 K, 3000 K, and 3500 K, the maximum number of clusters is greater than the maximum number of clusters generated at 4000 K. Except for 4000 K, the higher the temperature, the earlier the maximum number of clusters appears, which is because the increase of temperature will accelerate the destruction of DNTF molecular structure, which cause more species and a faster cluster formation. However, under the high temperature of 4000 K, DNTF is easier to decompose into small molecular products, resulting in a significant reduction in the maximum number of clusters.

Cluster structure analysis
The N/C and O/C values of the initial structure of DNTF are both 1.33. Figure 5 shows the distribution curve of the ratio of the number of N, O, and C atoms in the cluster molecules produced within 10 ps of the thermal decomposition of DNTF. And from Fig. 5, the N/C and O/C values are less than 1.33 for the four temperature conditions, indicating that a large number of N and O atoms break away from the DNTF molecule at the initial stage of the reaction. This is due to the large amount of small gas molecules such as NO and NO 2 produced at the beginning of the reaction. In addition, the O/C values are smaller than N/C for all four temperature conditions, indicating that the number of O atoms detached from the DNTF structure is greater than the number of N atoms. Figure 6 demonstrates the maximum cluster structure of DNTF within 10 ps under different temperature conditions. It can be seen that there is almost no Fig. 7 Evolution of potential energy with time at different temperatures Fig. 8 Evolution curve of the total number of species in the system with time within 10 ps long chain structure with six carbon atoms in the cluster structure, which indicates the C-C bond cleavage in the early stage of thermal decomposition of DNTF. After that, the carbon chain structure was recombined with other intermediate products to form a variety of cluster structures in the form of C-N bond.

Evolution of species number
Potential energy (PE) is the result of the interaction of two molecules, which is related to the strength of the interaction between two molecules, and the interaction between two molecules is related to the relative position of two molecules. Therefore, the evolution of PE can reflect whether the chemical reaction has reached the equilibrium state. It can be seen from Fig. 7 that under the constant temperature conditions of 2500 K, 3000 K, 3500 K, and 4000 K, the PE curve suddenly rises and then falls, which indicates that DNTF experienced endothermic process first and then exothermic decomposition. The higher the temperature, the faster the decay rate of PE curve; that is, the greater the initial decomposition rate, and the sooner the system reaches equilibrium. Under the temperature programmed condition of 300-4000 K, with the increase of system temperature, In order to study the degree of influence of different temperatures on the thermal decomposition of DNTF, the total number of species of DNTF within 10 ps under different temperatures with time was counted, as shown in Fig. 8. The total number of species at this point is the total number of free atoms and molecules in the system at a given time point. From Fig. 8, it is found that DNTF hardly decomposes at 10 ps at 2000 K, which is consistent with the trend of PE evolution curve. Under the temperature conditions of 2500 K, 3000 K, 3500 K, and 4000 K, there is a sudden rise in species number within 2 ps, and the higher the temperature, the greater the slope of the curve and the faster the growth rate of species number. The final number of species is about 42 at 2500 K, about 50 at 3000 K, and about 55 at 4000 K. Therefore, the higher the temperature, the more complex the decomposition pattern of DNTF and the higher the number of species generated. During the temperature programmed process of 300-4000 K, the species number increased rapidly after 5.6 ps, and the final species number was stable at about 45.

Initial decomposition path analysis
From Table 3, the constant temperature conditions of 2000 K, 2500 K, 3000 K, 3500 K, and 4000 K, NO 2 fragments generated by C-NO 2 bond cleavage in DNTF structure almost occur within 1.5 ps. Especially at 2000 K, 2500 K, and 3000 K, the denitration frequency of DNTF reached the maximum. At 3500 K and 4000 K, high temperature accelerates the cleavage of C-C bond in DNTF, thus, reducing the frequency of C 6 O 8 N 8 → C 6 O 6 N 7 + NO 2 reaction. It can be seen from the decomposition reaction equation of DNTF under the temperature programmed condition of 300-4000 K that large amounts of NO 2 and NO are generated in the initial stage of thermal decomposition of DNTF, which significantly increases the amount of NO 2 and NO in the system. This phenomenon is consistent with the evolution trend of NO 2 and NO in Fig. 3, and so is that with the thermal decomposition path of DNTF under constant temperature, indicating that the initial thermal decomposition path of DNTF under the temperature programmed is consistent with that under constant temperature.
In addition to direct NO 2 removal, DNTF molecules are accompanied by NO removal and carbon chain cleavage, which leads to the increase of reaction types and the decrease of frequency. This is one of the reasons why the frequency of all the reactions in Table 3 is not high. At 2000 K, a relative lower temperature, the initial decomposition of DNTF is mainly denitration, and the small molecular product produced is mainly NO 2 , as shown in Fig. 3. Under the constant temperature conditions of 2500 K, 3000 K, 3500 K, and 4000 K, the production amount and production rate of small molecular product NO will increase significantly, which indicates that the increase of temperature will significantly promote the reaction of NO production.
From the previous analysis, it is known that the C-NO 2 bond in DNTF molecule is easy to cleavage, and the O2-N11 bond in its furoxan ring is also a weak chemical bond. However, the reaction equations in Table 3 cannot explain this phenomenon. Therefore, it is necessary to analyze the trajectory file of DNTF reaction deeply.
The 0-0.5 ps trajectory file was split into 100 frames and analyze the frequency of ring-opening between furoxan and furazan and nitro detachment in turn, as shown in Fig. 10. Figure 9 shows the case of DNTF furazan ring Fig. 9 Schematic diagram of DNTF furazan ring opening and furoxan ring opening opening and furoxan ring opening. It can be seen from Fig. 10 that the ring-opening reaction frequency of furoxan is extremely significant under other temperature conditions except 4000 K, and all of them are before DNTF denitration reaction. For example, at 3000 K, the first ringopening reaction of furoxan appeared in the third frame (0.015 ps) of DNTF. Therefore, the ring-opening reaction of furoxan is the main trigger reaction of its thermal decomposition. Of course, we can also find that within 0.5 ps, these reactions are to some extent reversible, and the lower the temperature, the smaller the decomposition degree of DNTF molecules, which is obvious under the constant temperature condition of 2000 K and the temperature programmed process of 300-4000 K. With the increase of temperature, for example, at 4000 K, the C-C bond of DNTF molecule cleavages a lot, which leads to the decrease of the ring-opening reaction frequency of DNTF, so only the reaction frequency within 9 frames is counted. In addition, the ring opening of DNTF molecule furoxan precedes the ring opening of furazan and the denitration during the temperature programmed reaction, which is consistent with the decomposition mechanism of DNTF at constant temperature. To sum up, the thermal decomposition of DNTF firstly occurs the ring-opening reaction of furoxan, then the ring-opening reaction of furazan and denitration reaction, and finally the thermal decomposition reaction of molecular fragments. This is consistent with the findings of Sinditskii [9] and Li [13].

Conclusion
The reaction process of DNTF thermal decomposition at 300-4000 K temperature programmed and constant temperature 2000 K, 2500 K, 3000 K, 3500 K, and 4000 K was analyzed, and the initial reaction path of DNTF thermal decomposition was obtained. It is found that the ring-opening reaction of furoxan occurs first in the thermal decomposition of DNTF. Then, furazan ring opening and the C-NO 2 bond cleavage of DNTF molecule and furoxan ring-opening were the main trigger reaction at the initial stage of the reaction. By comparing the evolution trend of species number of DNTF under six temperature conditions, it can be seen that the higher the temperature, the faster the decomposition rate of DNTF and the more species generated.
The initial stage of the thermal decomposition reaction of DNTF involves the C-C bond cleavage, and the carbon chain structure resulting from the cleavage combines with other intermediate products in the form of C-N bonds to form a variety of cluster structures. At lower temperature and in the initial stage of thermal decomposition reaction, the small molecular products are mainly NO and NO 2 , while the thermal decomposition intermediate products at high temperature are mainly NO, NO 2 , and CO. With the progress of the reaction, these small molecules will react with some short chain structures to form the final products N 2 and CO 2 .
By comparing and analyzing the evolution trend of small molecule products, reaction path, and trigger reaction of DNTF thermal decomposition at 300-4000 K temperature programmed and different constant temperatures, it is found that temperature affects the initial thermal decomposition rate of DNTF, but does not change the initial thermal decomposition mechanism of DNTF.