Theoretical Prediction of Structures and Properties of 2,4,6-Trinitro-1,3,5-Triazine (TNTA) Green Energetic Materials from DFT and ReaxFF Molecular Modeling

Nitryl cyanide, O2NCN, as a new high-energy molecule, has not yet been successfully synthesized. It has prompted us to conduct a theoretical study of its possible space structures and properties. The RESP charges and the most stable spatial structures demonstrate that crystal morphology is affected by both the main nonbonded interactions and the molecular arrangement. The crystal structure prediction indicated that there are seven structures, namely P1, P21, P212121, P21/c, Pna21, Pbca, and C2/c. The most stable space structure is likely to be Pna21 and the corresponding cell parameters are Z = 4, a = 8.69 Å, b = 9.07 Å, c = 9.65 Å, and α = β = γ = 90.0°. To further study the intermolecular interactions of TNTA, a series of theoretical analyses were employed, including Hirshfeld surface analysis and fingerprint plots. The pyrolysis mechanism and properties show that high temperatures can promote decomposition. The systematic search approach can be a new strategy to identify structures effectively and has the potential to provide systematic theoretical guidance for the synthesis of TNTA.


Introduction
The ab initio theory was used to assess the energy performance and other important properties of 2,4,6-trinitro-1,3,5-triazine (TNTA) in 1907 [1]. The explosive performance of TNTA is superior to that of cyclotetramethylene tetranitramine (HMX), and it has a higher density of 2.1 g·cm −3 [2]. Korkin, A. A. et al. [3] conducted a theoretical study of its possible derivatives with similar performance (exothermic decomposition) but with possibly increased stability and a higher density in its condensed state. The research found that the higher the stability of the six-membered ring structure, the higher the density and energy released when it decomposes into a stable species, determining that TNTA is a potential high-energy compound. Yang K. [4] studied potential synthetic routes using the MP2/6-311G(d,p)//B3LYP/6-31G(d) levels of theory. Values of the heat of formation in the solid phase were predicted from density functional theory calculations. Densities were estimated from a regression equation obtained by molecular surface electrostatic potentials for TNTA. This work also suggested that TNTA might be formed in a solution if the trimerization reaction is carried out in a concentrated solution of nitryl cyanide. However, TNTA was not successfully synthesized. This current study rapidly predicted the structure and pyrolysis mechanism of TNTA, which can provide theoretical guidance for the experimental preparation.
Some studies [14,15] have predicted the structures and properties of high-energy materials. Other studies [16][17][18][19] have studied the thermal decomposition mechanism of high-energy materials by ReaxFF reactive molecular dynamics simulations. In this study, firstly, partial RESP charges and unit cell parameters were used to predict the space group and packing of TNTA. Then, Hirshfeld surface analysis and fingerprint plots were used to study the molecular interactions of TNTA. Finally, the thermal decomposition of TNTA was researched by reactive molecular dynamics using ReaxFF-lg.

Computational Details
The RESP charges were generated as follows: All molecules were optimized at the B3LYP/6-311g(d,p) density functional theory (DFT) level using Gaussian 16 software [20]. Then, the RESP charges were fitted from the optimized geometry and wave function using Multiwfn software [21]. The RESP charges of the monomers were derived from the respective trimers optimized by DFT, which is similar to the reported treatment method [22].
The Hirshfeld surfaces [23] or isosurfaces of the electron density were mapped by Crys-talExplorer [24], which is a freely available software. Two-dimensional mapping [25,26], a simple color plot, was used to analyze the intermolecular interactions quantitatively and qualitatively. The distances to the nearest atoms outside, de, and inside, di [26][27][28], were defined as the points of the intermolecular contacts.
The ReaxFF-lg force field using the LAMMPS program package performed all the molecular simulations. The original cell structure was a CIF file downloaded from CCDC, and the single-crystal cell was expanded into 3 × 4 × 2 supercells as the research substrate for dynamic simulation. The molecular formula and atomic numbering of TNTA are shown in Figure 1. First, the canonical ensemble (NVT) and the Berendsen thermostat were applied to the molecular dynamics (MD) simulation with a total time of 10 ps at 300 K, which further relaxed the TNTA supercell. Then, ReaxFF-lg isothermal-isochoric MD (NVT-MD) simulations were performed for 300 ps at 1800, 2250, 2500, 3000, and 3500 K, respectively, controlled by the Berendsen thermostat based on the relax supercell. An analysis of the fragments was performed with a 0.3 bond order cutoff value for each atom pair to identify the chemical species. The information of the dynamic trajectory was recorded every 50 fs, which was used to analyze the evolution details of TNTA in the pyrolysis process.
theoretical guidance for the experimental preparation.
The strength of intermolecular and intramolecular nonbonded interactions can be measured by restrained electrostatic potential (RESP) [5][6][7]. Zhang et al. [8][9][10][11][12] provided a detailed summary of the intermolecular interactions that play a dominant role in the packaging structures. Predicting crystal structures remains a challenging problem [13]. Some studies [14,15] have predicted the structures and properties of high-energy materials. Other studies [16][17][18][19] have studied the thermal decomposition mechanism of highenergy materials by ReaxFF reactive molecular dynamics simulations. In this study, firstly, partial RESP charges and unit cell parameters were used to predict the space group and packing of TNTA. Then, Hirshfeld surface analysis and fingerprint plots were used to study the molecular interactions of TNTA. Finally, the thermal decomposition of TNTA was researched by reactive molecular dynamics using ReaxFF-lg.

Computational Details
The RESP charges were generated as follows: All molecules were optimized at the B3LYP/6-311g(d,p) density functional theory (DFT) level using Gaussian 16 software [20]. Then, the RESP charges were fitted from the optimized geometry and wave function using Multiwfn software [21]. The RESP charges of the monomers were derived from the respective trimers optimized by DFT, which is similar to the reported treatment method [22].
The Hirshfeld surfaces [23] or isosurfaces of the electron density were mapped by CrystalExplorer [24], which is a freely available software. Two-dimensional mapping [25,26], a simple color plot, was used to analyze the intermolecular interactions quantitatively and qualitatively. The distances to the nearest atoms outside, de, and inside, di [26][27][28], were defined as the points of the intermolecular contacts.
The ReaxFF-lg force field using the LAMMPS program package performed all the molecular simulations. The original cell structure was a CIF file downloaded from CCDC, and the single-crystal cell was expanded into 3 × 4 × 2 supercells as the research substrate for dynamic simulation. The molecular formula and atomic numbering of TNTA are shown in Figure 1. First, the canonical ensemble (NVT) and the Berendsen thermostat were applied to the molecular dynamics (MD) simulation with a total time of 10 ps at 300 K, which further relaxed the TNTA supercell. Then, ReaxFF-lg isothermalisochoric MD (NVT-MD) simulations were performed for 300 ps at 1800, 2250, 2500, 3000, and 3500 K, respectively, controlled by the Berendsen thermostat based on the relax supercell. An analysis of the fragments was performed with a 0.3 bond order cutoff value for each atom pair to identify the chemical species. The information of the dynamic trajectory was recorded every 50 fs, which was used to analyze the evolution details of TNTA in the pyrolysis process. Simulated annealing of the global minimum energy, a thermodynamic problem, was performed by Monte Carlo (MC). The crystal was heated quickly to a high temperature and then cooled slowly to obtain the annealing structure. The packing groups were modified by the cell parameters by rotating and translating the rigid molecular units. For each space group, the low potential energy (E) packings were selected.

Crystal Structure Predict
The detailed COMPASS force field parameters and partial RESP charges are given in Table 1. The partial RESP charge of C1-C3 was 0.78. Meanwhile, the partial RESP charge of N1-N3 was −0.73. The greater the difference in the electronegativity between C1-C3 and N1-N3, the stronger the stabilizing bond-antibond interactions [5,6]. This demonstrates that the benzene ring-like structure is much more stable than other bonds. Moreover, the nonbonds of (C1-C3)/(N1-N3) have a significant effect on the intermolecular stacking [7]. Both the partial RESP charges of C1-C3 and N4-N6 were positive. The difference in the partial RESP charges between C1-C3 and N1-N3 was exactly the same as the value of the NO 2 fragment. This demonstrates that the bonds of C1-C3 and N4-N6 are relatively weak. The bonding of C1-C3 and N4-N6 is mainly due to the interaction between N4-N6 and O1-O6. This nonbonded interaction has an effect not only on the intramolecular but also on the intermolecular stacking. The nonbonded interactions of (C1-C3)/(N1-N3) and (N4-N6)/(O1-O6) determine the stacking of molecules. Table 1. Atom names, force field parameters, and RESP charges taken from COMPASS for prediction of crystal morphology.

Force Field Types
Atom Type Description RESP Charges The impact sensitivity (IS) could largely be affected by the crystal packing structure [5][6][7]. Therefore, we investigated the variations in the crystal packing structures from a viewpoint of seven space groups, as shown in Figure 2. The space groups P2 1 2 1 2 1 , P2 1 /c, and C2/c are not affected by the intermolecular interactions on molecular stacking. So, they are not the most reasonable spatial arrangement. The P1 and P2 1 space groups are affected by the nonbonded interaction of (N4-N6)/(O1-O6). However, the strongest nonbonded interaction of (C1-C3)/(N1-N3) is ignored. The packing structure of the Pbca space group is the exact opposite to that of the P1 and P2 1 space groups. The random molecular arrangement is not conducive to impact sensitivity. The packing structure of Pna2 1 combines the main nonbonded interaction, as shown in Table 1, and the graphene-like structure of TATB [29]. The Pna2 1 space group is the most likely spatial arrangement.
We employed the COMPASS force field and the Polymorph module in Materials Studio (MS) 4.4 2008 to obtain the possible molecular packing in the crystal state. Statistical data [30][31][32][33] demonstrate that most crystals belong to seven space groups (P1, P2 1 , P2 1 2 1 2 1 , P2 1 /c, Pna2 1 , Pbca, and C2/c), and the global search was confined to these groups only. The space group with the lowest energy was selected as the most likely molecular packing. The input structure for the polymorph search came from the ground-state geometry calculation at the B3LYP/6-311G(d,p) level. Here, we applied the PBE and PW91 pseudopotentials to the GGA and the CA-PZ pseudopotential to the LDA functional to the input structures. Table 2 shows the packing energy and cell parameters of the seven space groups. It was found that the differences among the packing energies estimated by the GGA-PBE for the seven space groups were much smaller than those estimated by the GGA-PW91 and LDA-CA-PZ. The energies were in the range of −4.24 to −0.87 kJ·mol −1 and the Pna2 1 space group had the lowest energy. The lower the lattice energy, the more stable the crystal structure. Therefore, due to the lowest Gibbs free energy, the Pna2 1 space group is the most likely crystal structure for TNTA. The corresponding lattice parameters were Z = 4, a = 8.69 Å, b = 9.07 Å, c = 9.65 Å, and α = β = γ = 90.0 • .  We employed the COMPASS force field and the Polymorph module in Materials Studio (MS) 4.4 2008 to obtain the possible molecular packing in the crystal state. Statistical data [30][31][32][33] demonstrate that most crystals belong to seven space groups (P1, P21, P212121, P21/c, Pna21, Pbca, and C2/c), and the global search was confined to these groups only. The space group with the lowest energy was selected as the most likely molecular packing. The input structure for the polymorph search came from the ground-state geometry calculation at the B3LYP/6-311G(d,p) level. Here, we applied the PBE and PW91 pseudopotentials to the GGA and the CA-PZ pseudopotential to the LDA functional to the input structures. Table 2 shows the packing energy and cell parameters of the seven space groups. It was found that the differences among the packing energies estimated by the GGA-PBE for the seven space groups were much smaller than those estimated by the GGA-PW91 and LDA-CA-PZ. The energies were in the range of −4.24 to −0.87 kJ·mol −1 and the Pna21 space group had the lowest energy. The lower the lattice energy, the more stable the crystal structure. Therefore, due to the lowest Gibbs free energy, the Pna21 space group is the most likely crystal structure for TNTA. The corresponding lattice parameters were Z = 4, a = 8.69 Å, b = 9.07 Å, c = 9.65 Å, and α = β = γ = 90.0°.

Method/ Functional
Space To obtain a better understanding of the crystal stacking, the intermolecular interactions [11,34] of single crystals were studied through Hirshfeld surface analysis using freeware [26,35]. In the Hirshfeld surface analysis, the red and blue areas represent the high and low probabilities of close contact with external molecules, respectively [36]. Two features can be seen on the surfaces in Figure 3. One is that the interaction takes place through both the external and internal atoms because the red dots are distributed on both the front and side faces. These irregular surfaces show that TNTA molecules are irregular, more uneven, and less planar. The spatial symmetry of the whole surface shape of Pna2 1 is relatively better than the other space groups. This shows that the crystal stacking of Pna2 1 is more spatially symmetrical. The other feature is that the red dots are concentrated around the oxygen atoms. O···O contacts become dominant in TNTA molecules. The more oxygen atoms exposed on the exterior of TNTA molecules, the more sensitive the molecules are. Figure 4 shows that the O···O contact percentage of the Pna2 1 space group is 39%, which is the lowest percentage of the seven space groups. These results demonstrate that the Pna2 1 space configuration is the most stable crystal stacking for TNTA molecules. The percentage of N···O contacts is the second highest among the seven space groups, at 43%, as shown in Figure 4. The C···O contacts are the other main intermolecular interactions. Looking at Figures 1 and 3, we can draw the conclusion that N···O and C···O contacts may contribute to planar conjugated molecular structures between the NO2 group and the heterocyclic group. This may be because the crystal stacking of Pna21 is much more stable than the others. Two-dimensional plots of these intermolecular contacts are shown in Figure 5. The narrow orange line denoting O···O contacts is much more obvious in the plot of the P1 space group, which suggests an increase in the O···O contacts. This demonstrates that the crystal stacking is much more sensitive in the P1 space group.  The percentage of N···O contacts is the second highest among the seven space groups, at 43%, as shown in Figure 4. The C···O contacts are the other main intermolecular interactions. Looking at Figures 1 and 3, we can draw the conclusion that N···O and C···O contacts may contribute to planar conjugated molecular structures between the NO2 group and the heterocyclic group. This may be because the crystal stacking of Pna21 is much more stable than the others. Two-dimensional plots of these intermolecular contacts are shown in Figure 5. The narrow orange line denoting O···O contacts is much more obvious in the plot of the P1 space group, which suggests an increase in the O···O contacts. This demonstrates that the crystal stacking is much more sensitive in the P1 space group. The percentage of N···O contacts is the second highest among the seven space groups, at 43%, as shown in Figure 4. The C···O contacts are the other main intermolecular interactions. Looking at Figures 1 and 3, we can draw the conclusion that N···O and C···O contacts may contribute to planar conjugated molecular structures between the NO 2 group and the heterocyclic group. This may be because the crystal stacking of Pna2 1 is much more stable than the others.
Two-dimensional plots of these intermolecular contacts are shown in Figure 5. The narrow orange line denoting O···O contacts is much more obvious in the plot of the P1 space group, which suggests an increase in the O···O contacts. This demonstrates that the crystal stacking is much more sensitive in the P1 space group.

Reactive Molecular Dynamics Using ReaxFF-lg
The RESP charges, space groups, and intermolecular interactions predicted that Pna21 has a much more stable crystal morphology. Therefore, the MD simulations were performed solely for the crystal structure of the Pna21 space group. The evolution of the potential energy (PE) of the system is shown in Figure 6. The curves, except the 1800 K curve, first decrease steadily, then reach an approximately horizontal line, indicating that TNTA breaks down rapidly and instantaneously at these extreme temperatures. The curve of 1800 K decreases steadily, implying that TNTA is not completely decomposed within 300 ps. The asymptotic value of PE increased with the temperature increase. The declining rate of PE demonstrates the accelerating release of heat as the temperature increased.

Reactive Molecular Dynamics Using ReaxFF-lg
The RESP charges, space groups, and intermolecular interactions predicted that Pna2 1 has a much more stable crystal morphology. Therefore, the MD simulations were performed solely for the crystal structure of the Pna21 space group. The evolution of the potential energy (PE) of the system is shown in Figure 6. The curves, except the 1800 K curve, first decrease steadily, then reach an approximately horizontal line, indicating that TNTA breaks down rapidly and instantaneously at these extreme temperatures. The curve of 1800 K decreases steadily, implying that TNTA is not completely decomposed within 300 ps. The asymptotic value of PE increased with the temperature increase. The declining rate of PE demonstrates the accelerating release of heat as the temperature increased.

Reactive Molecular Dynamics Using ReaxFF-lg
The RESP charges, space groups, and intermolecular interactions predicted that Pna21 has a much more stable crystal morphology. Therefore, the MD simulations were performed solely for the crystal structure of the Pna21 space group. The evolution of the potential energy (PE) of the system is shown in Figure 6. The curves, except the 1800 K curve, first decrease steadily, then reach an approximately horizontal line, indicating that TNTA breaks down rapidly and instantaneously at these extreme temperatures. The curve of 1800 K decreases steadily, implying that TNTA is not completely decomposed within 300 ps. The asymptotic value of PE increased with the temperature increase. The declining rate of PE demonstrates the accelerating release of heat as the temperature increased.  whole decomposition stage under the five extreme conditions. All the curves, except the 1800 K curve, initially rapidly increase and then decrease until they disappear over time, indicating that TNTA first breaks the C-NO2 bond to generate NO2, and then all of the NO2 fragments decompose into more stable products, such as N2, over time. As the temperature increased, the maximum value of the NO2 and the rates of increase and decrease in the amount of NO2 were higher, indicating that increases in temperature can accelerate the decomposition rate of TNTA. At 1800 K, the amount of NO2 first increased and then oscillated around the equilibrium value of 30. The results show that TNTA does not decompose completely in 300 ps. This conclusion is consistent with the results in Figures 6 and 7. The curves of the amount of NO initially rapidly increase and then decrease until reaching an equilibrium value, indicating that TNTA first converts nitro to nitroso, and then a part of the NO fragments decompose into more stable products over time. As the temperature increased, the maximum value of NO and the rate of increase in the amount of NO were higher, indicating that increases in the temperature can accelerate the decomposition rate of TNTA. However, the equilibrium value was lower with the temperature increase, demonstrating that higher temperatures can accelerate NO fragment decomposition into more stable products.
The curves of the amount of NO3 initially slightly increase and then decrease with time evolution. As the temperature increased, the maximum value of NO3 and the rate of increase in the amount of NO3 was lower, indicating that increases in the temperature can inhibit the production of NO3. The NO3 disappeared, except at 1800 K, it reached an equilibrium value. This demonstrates that higher temperatures can accelerate the complete decomposition of TNTA.  Figure 8 demonstrates the time evolution of the main fragments of NO 2 during the whole decomposition stage under the five extreme conditions. All the curves, except the 1800 K curve, initially rapidly increase and then decrease until they disappear over time, indicating that TNTA first breaks the C-NO 2 bond to generate NO 2 , and then all of the NO 2 fragments decompose into more stable products, such as N 2 , over time. As the temperature increased, the maximum value of the NO 2 and the rates of increase and decrease in the amount of NO 2 were higher, indicating that increases in temperature can accelerate the decomposition rate of TNTA. At 1800 K, the amount of NO 2 first increased and then oscillated around the equilibrium value of 30. The results show that TNTA does not decompose completely in 300 ps. This conclusion is consistent with the results in Figures 6 and 7.
The curves of the amount of NO initially rapidly increase and then decrease until reaching an equilibrium value, indicating that TNTA first converts nitro to nitroso, and then a part of the NO fragments decompose into more stable products over time. As the temperature increased, the maximum value of NO and the rate of increase in the amount of NO were higher, indicating that increases in the temperature can accelerate the decomposition rate of TNTA. However, the equilibrium value was lower with the temperature increase, demonstrating that higher temperatures can accelerate NO fragment decomposition into more stable products.
The curves of the amount of NO 3 initially slightly increase and then decrease with time evolution. As the temperature increased, the maximum value of NO 3 and the rate of increase in the amount of NO 3 was lower, indicating that increases in the temperature can inhibit the production of NO 3 . The NO 3 disappeared, except at 1800 K, it reached an equilibrium value. This demonstrates that higher temperatures can accelerate the complete decomposition of TNTA.
The decay rate of TNTA molecules increased with increasing temperature. The evolution of the amounts of the final products (N 2 , and CO 2 ) over time is shown in Figure 9. The final products were constantly produced, and their amounts were nearly stable at the end. The N 2 and CO 2 increased sharply, then reached equilibrium, except at 1800 K. However, the value of N 2 and CO 2 increased steadily over 300 ps. This demonstrates that TNTA decomposes totally except at 1800 K. The equilibrium values were larger as the temperature increased. It took less time to reach equilibrium as the temperature increased. This implies that higher temperatures can accelerate the decomposition of TNTA. An interesting phenomenon is that the variation trend of N 2 is consistent with that of CO 2 . At the same time, the amount of N 2 was slightly higher than the amount of CO 2 during the whole reaction time at all five high temperatures. The decay rate of TNTA molecules increased with increasing temperature. The evolution of the amounts of the final products (N2, and CO2) over time is shown in Figure 9. The final products were constantly produced, and their amounts were nearly stable at the end. The N2 and CO2 increased sharply, then reached equilibrium, except at 1800 K. However, the value of N2 and CO2 increased steadily over 300 ps. This demonstrates that TNTA decomposes totally except at 1800 K. The equilibrium values were larger as the temperature increased. It took less time to reach equilibrium as the temperature increased. This implies that higher temperatures can accelerate the decomposition of TNTA. An interesting phenomenon is that the variation trend of N2 is consistent with that of CO2. At the same time, the amount of N2 was slightly higher than the amount of CO2 during the whole reaction time at all five high temperatures.  The decay rate of TNTA molecules increased with increasing temperature. The evolution of the amounts of the final products (N2, and CO2) over time is shown in Figure 9. The final products were constantly produced, and their amounts were nearly stable at the end. The N2 and CO2 increased sharply, then reached equilibrium, except at 1800 K. However, the value of N2 and CO2 increased steadily over 300 ps. This demonstrates that TNTA decomposes totally except at 1800 K. The equilibrium values were larger as the temperature increased. It took less time to reach equilibrium as the temperature increased. This implies that higher temperatures can accelerate the decomposition of TNTA. An interesting phenomenon is that the variation trend of N2 is consistent with that of CO2. At the same time, the amount of N2 was slightly higher than the amount of CO2 during the whole reaction time at all five high temperatures. For each temperature, the system was heated to the target temperature within 0.3 ps and the C-NO 2 bond was broken to release NO 2 fragments. The life time was shorter as the temperature increased and the maximum number of NO 2 fragments was larger with the temperature increase, indicating that higher temperatures would promote the production of NO 2 and then promote its decomposition into more stable products. The appearance time of the fragments NO, NO 3 , CO 2 , and N 2 became shorter as the temperature increased. This demonstrates that a high temperature will accelerate the decomposition of TNTA. The maximum amounts of NO 2 , NO, CO 2 , and N 2 increased with the temperature increase. However, for NO 3 , the change rule was just the opposite. These results imply that higher temperatures can better promote the production of stable products. Table 3 demonstrates the relatively higher NF of the reactions at 2500 K and 3500 K. We selected the two temperatures (3500 K, the higher in the series) to illustrate the main trends because of the higher number of reactions.  4,6,12,16,18,21,26,29,30) that occurred due to the high concentration of CO 2 ; (b) the breakdown of fragments generating CO 2 (R1, 2,5,7,15,17,22,23,24,25,26,27) that occurred due to the stability of CO 2 . Conversely, the overall analysis of simulations shows that the other main gas product is N 2 . The chemical effect of N 2 can be expressed in terms of two types of reactions: R-C-N 2 → N 2 + R-C (R1, 9, 10, 20, 28) and R-N-N 2 → N 2 + R-N (R8, 13). Thus, in general, increases in CO 2 and N 2 production in the thermal decomposition of TNTA at high temperatures are expected.

Conclusions
We report predictions using quantum mechanics of the most stable crystal structures of TNTA, which is a promising green high-energy-density material. Firstly, the partial RESP charges were used with the COMPASS force field to find the main nonbonded interactions which may affect the space groups. Then, unit cell parameters, Hirshfeld surface analysis, and fingerprint plots of the seven space groups were used to predict the most promising space groups and packings. Finally, the thermal decomposition of TNTA was researched by reactive molecular dynamics using ReaxFF-lg. We found that the Pna2 1 space group with the corresponding cell parameters of Z = 4, a = 8.69 Å, b = 9.07 Å, c = 9.65 Å, and α = β = γ = 90.0 • was the most compatible with TNTA. The pyrolysis mechanism implies that increasing the temperature can accelerate the decomposition rate of TNTA. The two main products, CO 2 and N 2 , had the same variation trend, and the amount of N 2 was slightly higher than that of CO 2 during the whole reaction time.