Carbon-free energetic materials: computational study on nitro-substituted BN-cage molecules with high heat of detonation and stability

A new series of high-energy density materials (HEDMs) B6N6H6−n(NO2)n (n = 1–6) are studied at the M06-2X/6-311++G**, ωB97XD/6-311++G** and B3LYP/6-311++G** levels. Analysis of the structural changes caused by substituting the NO2 and the electronic structures, such as electron localization function (ELF), Wiberg bond index (WBI), charge transfer and bond dissociation energies (BDE), provide important insights into the essence of the chemical characteristics and stability. Moreover, the Born–Oppenheimer molecular dynamic (BOMD) simulation is performed to verify their stability, which suggests that only the BN-cage derivatives with one and two nitro groups bonding with boron atoms (NO2-1-1 and NO2-2-1) can remain stable under ambient conditions. To predict the detonation performance and sensitivity of these two stable BN-cage energetic molecules accurately, the density, gas phase enthalpy of formation, enthalpy of sublimation, detonation performance, impact sensitivity and BDE are calculated systematically. The calculation results show that both NO2-1-1 and NO2-2-1 have a higher heat of detonation, higher value of h50, and larger BDE of trigger bonds than CL-20.


Introduction
High-energy density materials (HEDMs) with both superior detonation performance and low sensitivity have always drawn the attention of research scientists. [1][2][3][4][5][6][7][8][9][10][11] The increased energy density oen comes at the expense of molecular stability. Seeking new HEDMS with a ne balance between high detonation performance and low sensitivity remains to be an interesting, but very challenging task.
Inspired by this, we attempted to replace the C atoms of CL-20 with B atoms to develop a new carbon-free BN-cage energetic system. The nitro (NO 2 ) groups are bonded to the BN-cage structure to increase the nitrogen and oxygen content, which can achieve a good oxygen balance. In this theoretical study, the BN-cage and its NO 2 derivates are systematically investigated to understand their stability, electronic structures, safety and detonation properties at the M06-2X/6-311++G**, uB97XD/6-311++G** and B3LYP/6-311++G** levels of density functional theory (DFT). Through the systemic study of the above compounds, it is predicted that two BN-cage derivatives, NO 2 -1-1 and NO 2 -2-1, will have a superior detonation performance and low sensitivity, and could therefore potentially be new HEDMs.

Theoretical methods
All quantum mechanical calculations in this paper were performed using the Gaussian-09 soware package 22 at the M06-2X, uB97XD and B3LYP level with the 6-311++G** basis set. The M06-2X meta-hybrid density functional developed by Zhao and Truhlar is an advanced method to calculate energies, which have shown a good performance for main group chemistry and kinetics studies. 23,24 For comparison, the theoretical methods of uB97XD and B3LYP were also employed with the basis set of 6-311++G**. The uB97XD method includes a 100% long-range exact exchange, a small fraction (about 22%) of the short-range exact exchange, a modied B97 exchange density functional for a short-range interaction, and the B97 correlation density functional and empirical dispersion correction. 25,26 uB97XD yields satisfactory accuracy for thermochemistry and kinetics. The DFT method of B3LYP combines Becke's threeparameter (B3) 27 functional with the Lee-Yang-Parr (LYP) 28 correlation functional, which has been considered to be capable of accurately predicting the structural parameters and frequencies of many nitro-substituted compounds. [29][30][31] Moreover, the 6-311++G** basis set can generally give satisfactory geometries. 32 All of the optimized molecular structures in this paper belong to local minima on their singlet spin state potential surface and have no imaginary frequency. Structures of the nitro-substituted BN-cage compounds are denoted as NO 2 -n-N in the present article, where n is the number of energetic groups, N orders the structures according to their relative energies using the M06-2X method.
Calculation of the electron localization function (ELF), deformation energies, net charge and charge transfer at the same level are used to analyze the electronic structure and stability of the NO 2 -n-N series. Moreover, Born-Oppenheimer molecular dynamic simulations (BOMD) are employed to exam the stability of the designed cage compounds to select the stable BN-cage derivates with a DFT in the level of M06-2X/6-311++G**. BOMD deals with the electronic and nuclear problems separately. In this method, the electronic structure in the ground state is calculated at each set of atomic positions, usually by optimization of the Kohn-Sham orbitals using an iterative method. 33 The BOMD method in general needs more CPU time than for other MD soware, but the method is more robust and stable.
Heat of detonation, detonation velocity and detonation pressure were calculated using EXPLO 5, 34 using the density and solid phase enthalpy of formation (D f H (s)). The gas phase enthalpy of formation (D f H (g)) and the enthalpy of sublimation ðDH sub Þ were used to calculate the solid phase enthalpy of formation (D f H (s)).
The following isodesmic reactions (1) combined with the computational Formula (2) were used to calculate the D f H (g) of the studied compounds.
In eqn (1) and (2), n is the number of nitro groups, DH 298 is the enthalpy change of the reaction at 298 K, DE is the change in total energy between the products and the reactants, DZPE is the change of zero-point energy between the products and the reactants, and DH T is the thermal correction from 0 to 298 K.
Electrostatic potential analysis (ESP) is employed to analyze the molecular surface at the M06-2X/6-311++G level of theory. In the eqn (3)-(7), the electrostatic potential on the surface is characterized by V s (r). V + s (r i ), V À s (r i ) and V s (r i ) represents the positive, negative and the overall value of V s (r) at any point r i on the surface. The V s þ , V s À and V s represent their averages, and s + 2 , s À 2 , and s tot 2 are the corresponding variances. n is the electrostatic balance parameter. m and n are the number of positions with positive and negative potentials on the molecular surface. The parameters mentioned above will be used to further calculate the enthalpy of sublimation, density and impact sensitivity in the corresponding equations.
The enthalpy of sublimation ðDH sub Þ can be evaluated using the developed method by Rice and Politzer et al. 35
The crystal density (r) of the studied compounds can be calcu- In order to understand the sensitivity and safety of the studied cage compounds, the impact sensitivity, 37 Wiberg bond index (WBI), 38,39 and bond dissociation energy (BDE) 40 were calculated. Impact sensitivity can be tested using the drop hammer test, which is measured by dropping a given mass (2.5 kg) upon the compound and recording the height (h 50 ) which has a 50% probability of producing an explosion. The value of h 50 is calculated using eqn (10). 37 In eqn (10), a 2 ¼ À0.0064, b 2 ¼ 241.42, and g 2 ¼ À3.43.
The BDE is calculated using eqn (11). In general, if the value of BDE is greater than 84 kJ mol À1 , the corresponding bond is considered to be stable.
3. Results and discussion

Conguration
Fig . 1 shows an optimized structure of the designed BN-cage that was conrmed to be a global minimal without imaginary frequencies at the M06-2X level. The distances of B3-N7, B5-N6 and B16-N15 are 1.716Å, 1.716Å and 2.803Å, respectively. The corresponding WBI, determined using natural bond orbital (NBO) analysis were 0.5085, 0.5085 and 0.0049, respectively. The relatively long distance and small WBI between B16 and N15 suggests the absence of the corresponding bond. Additionally, the bond lengths (WBI) of B1-N8, B2-N9 and B4-N10 were found to be 1.370Å (1.2746), 1.403Å (1.0820) and 1.403Å (1.0820), respectively. The shorter bond lengths and larger WBI of these three bonds compared with the other B-N bonds are indicative that double bonds exist between B1-N8, B2-N9 and B4-N10. Also, the topological analysis of the ELF was performed to understand the bond character of the B-N bonds in the designed BN-cage. The corresponding color-lled maps and the curve maps of ELF for these bonds are shown in Fig. S1, † which further conrms the existence of the B3-N7, B5-N6, B1-N8, B2-N9, and B4-N10 covalent bonds, and the inexistence of the B16-N15 bond.
The H atoms in the designed BN-cage were replaced by nitro groups to construct a series of new energetic molecules. All of the optimized structures could retain the integrity of the BNcages. Moreover, the calculation results indicated that the energetic groups bond with boron atoms in the cage led to lower-energy structures than that with nitrogen atoms. Fig. 2 shows the optimized geometries of the NO 2 -n-1 (n ¼ 1-6).
The four mononitro-substituted BN-cage compounds shown in Fig. S1 and Table S2 † were found to lie less than 200 kJ mol À1 in energy above the lowest-energy structure NO 2 -1-1. The nitro group in the lowest-lying structures NO 2 -1-1 is bonded to the B5 atom by replacing H12 in the BN-cage. It should be pointed out that the distance between B3 and N7 in NO 2 -1-3 is 2.621Å, the long distance of the bond and the corresponding ELF analysis (shown in Fig. S3 †) suggests that the B3-N7 bond does not exist.
The ve dinitro-substituted BN-cage compounds shown in Fig. S2 and Table S3 † are found to lie less than 200 kJ mol À1 in energy above NO 2 -2-1. Similar to the mononitro-substituted compounds, the nitro groups bond with B3 and B5 in NO 2 -2-1 to obtain the lowest-energy structure with the symmetry of the C s . The NO 2 -2-2 lies at an energy of 27.2 kJ mol À1 above the NO 2 -2-1S at the M06-2X level, having the B5-NO 2 and B16-NO 2 bonds. However, the NO 2 -2-3, NO 2 -2-4 and NO 2 -2-5, with the N-NO 2 bonds have much higher total energies than the NO2-2-1 and NO2-2-2 due to the high energy of the N-N bonds.
As shown in Fig. S4 and Table S5, † three tetranitrosubstituted BN-cages are found to lie less than 200 kJ mol À1 in energy above NO 2 -4-1. Both of the rst two structures NO 2 -4-1 and NO 2 -4-2 have three B-NO 2 bonds and one N-NO 2 bonds, while the NO 2 groups in NO 2 -4-3 bond with two B atoms and two N atoms. It should be pointed out that the distances between the B3 and N7 atoms in NO 2 -4-2 and NO 2 -4-3 are 2.534 and 2.491Å, respectively. The relatively long distances and the corresponding ELF analysis shown in Fig. S7 and S8 † suggests the absence of corresponding bonds.
Only the two pentanitro-substituted BN-cages were found, which are shown in Fig. S5 and Table S6. † In NO 2 -5-1, the ve nitro groups bond with the three boron atoms and two nitrogen atoms. Additionally, the distance of B3/N7 and B16/N15 in NO 2 -5-1 are 2.608 and 2.629Å respectively, which are about 1Å longer than the other B-N bonds. This is indicative of the absence of these two bonds, combined with the color-lled map and the curve map of ELF (shown in Fig. S10 †).
Finally, only one hexanitro-substituted BN-cage NO 2 -6-1 is found, as shown in Fig. S6 and Table S7. † All of the H atoms in the designed BN-cage are replaced by nitro groups leading to the C s symmetry for NO 2 -6-1.

Electronic structure and thermal stability
To study the change caused by substituting the H of the designed BN-cage with NO 2 , the deformation energies of the skeleton (DE cage ) and the relative (DBDE) are calculated and summarized in Table 1. DBDE can be calculated using eqn (12), which is the relative BDE of cage-NO 2 bonds to the corresponding cage-H bonds in the designed BN-cage without substitution.
DE cage corresponds to the change of the BN-cage skeleton. The DE cage shown in Table 1 indicates that the substitution of the nitro groups leads to the deformation of the cage skeleton.
Adding nitro groups to the nitrogen atoms of the BN-cage brings more changes to the skeleton compared with the boron atoms.
It should be pointed out that NO 2 -1-3 has the most skeleton changes among the four isomers, caused by the breaking of the B3-N7 bond. Negative values of DBDE suggest that the bond strengths of all the cage-NO 2 bonds are lower than that of the corresponding cage-H bonds.
To further study the stability of the mononitro-substituted BN-cage isomers, the deformation energies of the BN-cage skeletons ðDE * cage Þ; the deformation energies of the nitro groups ðDE * NO2 Þ; and the relative BDE of the cage-NO 2 bonds (DBDE*) to NO 2 -1-1 are calculated and the results are shown in Table 2.
As shown in Table 2, comparing with NO 2 -1-1, the decrease of the BDE of the B-NO 2 bond in the NO 2 -1-2 structure leads to an increase of its energy, and the deformation of its skeleton could be ignored due to its small DE * cage value. For the Nmononitrosubstituted isomers NO 2 -1-3 and NO 2 -1-4, the difference of the B-NO 2 and N-NO 2 bond is the primary reason leading to their energy increase. Compared to the global  value and the change of the NO 2 could be ignored, as shown in Table 2.
The charge transfer calculated using the NBO analysis is also employed to analyze the stability differences among the different substitution positions. The valence-bond structure of the nitro group, shown in Fig. 3, indicates that the nitro is an electrophilic group, in which the nitrogen atom tends to combine with the boron atom in the cage to provide a stable conguration. Fig. 4 and Table 3 describe the molecular structures, net charge, charge transfer and BDE of the cage-NO 2 bonds for the mononitro or dinitro-substituted cages. As listed in Table 3, the results show that all the nitro groups in the B(N)mononitro or dinitro-substituted cages are electron accepters, and the charge transfers of the B-mononitro or dinitrosubstituted cage is always larger than those of the Nmononitro or dinitro-substituted cage, suggesting the higher dissociation energies of the B-NO 2 bonds compared with the N-NO 2 bonds.
The above calculations suggest that the B-nitrosubstituted cages have a relatively better stability than the Nnitrosubstituted cages. As such, the NO 2 -1-1, NO 2 -2-1 and NO 2 -3-1 structures were selected for further study, as discussed in the subsequent sections.

Dynamic simulation
Based on the above calculation results for the stability of these designed BN-cage derivatives, molecular dynamic calculations of the B-substituted NO 2 -n-1 (n ¼ 1, 2, 3) structures were carried out to exam their dynamic stability. BOMD were employed to perform the molecular dynamic simulation. The static congurations calculated at the M06-2X/6-311++G** level were run as seeds for the BOMD trajectories using the Gaussian-09 program package. Equilibration steps were performed over 1 ps and production runs are 3 ps long with the time step of 1 fs. Fig. 3 The valence-bond structures of the nitro group. The simulation results of the NO 2 -n-1 (n ¼ 1, 2, 3) structures show that the NO 2 -1-1 and NO 2 -2-1 remain stable at 298 K during the whole simulation, while the nitro groups in NO 2 -3-1 cannot remain stable at 298 K, as shown in Fig. 5. As can be seen from Fig. 5, the nitro groups in NO 2 -3-1 decompose from the cage with the cage-NO 2 bond breaking. Oxygen atoms in the nitro groups in the NO 2 -3-1 bond with the cage directly at 298 K, with B-O bond lengths of 1.452 and 1.316Å, while the NO groups decompose from the BN-cage compounds, with the distances between O and NO groups being $1.55Å. Moreover, deformation of the cage skeleton occurs with the distance between B5 and N6 being 2.821Å. Thus, the NO 2 -3-1 structure cannot remain stable at room temperature. Above all, the NO 2 -1-1 and NO 2 -2-1 structures have exhibited an attractive potential for application in the areas of energetic materials and the details are discussed in the following study on detonation performance.

Detonation performance
In order to quantitatively evaluate the detonation performance of the NO 2 -1-1 and NO 2 -2-1 structures, the predicted density (r), enthalpy of sublimation ðDH sub Þ; solid phase enthalpy of formation (D f H (s)), oxygen balance and detonation parameters (Q, D, P) were calculated systematically using EXPLO 5. The results, together with the related information for CL-20, are listed in Table 4 for comparison.
The results in Table 4 show that the r of NO 2 -1-1 and NO 2 -2-1 are 1.808 g cm À3 (1.789 g cm À3 , 1.732 g cm À3 ) and 1.931 g cm À3 (1.871 g cm À3 , 2.059 g cm À3 ) respectively at the M06-2X (uB97XD, B3LYP)/6-311++G** level. In the meantime, the calculated densities of CL-20 are 1.954 g cm À3 at the M06-2X level, 1.907 g cm À3 at the uB97XD level and 2.022 g cm À3 at the B3LYP level, which indicates the relatively high densities of these two newly proposed designed molecules. The Q of NO 2 -2-1 is 31.25% (M06-2X), 31.26% (uB97XD) or 24.79% (B3LYP) higher than that of CL-20 (6314 kJ kg À1 ), with a D of over 6800 m s À1 . The superior detonation heat of both the NO 2 -1-1 and NO 2 -2-1 structures indicate the great potential of these two energetic compounds to be HEDMs.

Sensitivity and safety
To test the sensitivity and safety of the above NO 2 -1-1 and NO 2 -2-1 structures, two methods have been employed in the present work: (1) impact sensitivity; and (2) BDE of the trigger bonds.
3.5.1 Impact sensitivity. The impact sensitivity is a crucial parameter responsible for effectiveness. It reects whether the explosive is stable enough to be handled or stored under typical conditions. In this paper, the theoretical impact sensitivity of NO 2 -1-1, NO 2 -2-1 and CL-20 were calculated at the M06-2X, uB97XD and B3LYP level to compare their safety and practicability (see Table 5). The calculation results show that the h 50 of the NO 2 -1-1 and NO 2 -2-1 structures were 49.498 cm (47.745 cm, 50.555 cm) and 34.459 cm (34.089 cm, 36.115 cm) at the M06-2X (uB97XD, B3LYP)/6-311++G** level, respectively, which are higher than that of CL-20. Such a result illustrates their high value of application.
3.5.2 Bond dissociation energy. The BDE of the trigger bond is oen a key factor in investigating the thermal stability and pyrolysis mechanism of energetic compounds. Generally, the smaller the energy needed for breaking a bond, the weaker the bond.

Conflicts of interest
There are no conicts of interest to declare. Table 4 Theoretical density (r), enthalpy of sublimation ðDH sub Þ; solid phase enthalpy of formation (D f H (s)), detonation velocity (D), detonation pressure (P), heat of detonation (Q) and oxygen balance of the NO 2 -1-1 and NO 2 -2-1 structures. Experimental parameters of CL-20 are labeled in parentheses (the theoretical densities of CL-20 are 1.954 g cm À3 at M06-2X, 1.907 g cm À3 at uB97XD and 2.022 g cm À3 at B3LYP)