Theoretical studies of novel high energy density materials based on oxadiazoles

In this study, 32 energetic compounds were designed using oxadiazoles (1,2,5-oxadiazole, 1,3,4-oxadiazole) as the parent by inserting different groups as well as changing the bridge between the parent. These compounds had high density and excellent detonation properties. The electrostatic potentials of the designed compounds were analyzed using density functional theory (DFT). The structure, heat of formation (HOF), density, detonation performances (detonation pressure P, detonation velocity D, detonation heat Q), and thermal stability of each compound were systematically studied based on molecular dynamics. The results showed that the -N3 group has the greatest improvement in HOF. For the detonation performances, the directly linked -N=N- and -NH-NH- were beneficial when used as a bridge between 1,2,5-oxadiazole and 1,3,4-oxadiazole, and it can also be found that bridge changing had little effect on the trend of detonation performance, while energetic groups changing influenced differently. In general, the introduction of nitro groups contributes to the improvement of the detonation performance of the compounds. In this study, the compounds containing the highest amount of nitro groups were found to have better detonation performance than their counterparts and were not significantly different from RDX and HMX.


Introduction
In recent years, the design and synthesis of new high energy density materials is an important part of research in the field of energetic materials [1][2][3][4]. High energy density materials (HEDMs) are compounds of a given volume that can generate a considerable amount of energy simply through the creation and breaking of chemical bonds within the molecule. High energy density materials have a series of advantages such as high energy density, intensity, and stability. In addition, to complete the research on HEDMs, high nitrogen energetic materials as a new type of energetic material have also received more and more attention due to their considerable positive heat of formation, high density, excellent detonation properties (detonation velocity and detonation pressure), and admissible thermal stability [5].
1,2,5-oxadiazole (furazan), with the molecular formula C 2 H 2 ON 2 , is a five-membered ring. Zelinskii [6][7][8] Institute of Organic Chemistry, Russian Federation Academy of Science had synthesized a variety of furazan energetic compounds after more than 20 years of research. Chaoyang Zhang [9] reported nitro-furazan to increase in density of 0.06-0.08 g cm 3 and detonation velocity of about 300 m s −1 . 1,3,4-oxadiazole is an isomer of furazan with a low enthalpy of production; however, the synthesis of energetic compounds with high nitrogen content using 1,3,4-oxadiazole can reduce the sensitivity of energetic compounds and improve the oxygen balance. Shreeve et al. reported the use of 1,3,4-oxadiazole to stabilize the nitrofurazone ring [10]. Based on 1,3,4-oxadiazole and the nitrofurazone backbone, a series of high nitrogen energetic salts were synthesized with good detonation properties and satisfactory susceptibility and stability (detonation velocities >7493 m s −1 , detonation pressures >20. 4 GPa, IS = 15 J, FS = 120 N). Lu Ming et al. concluded that the compounds formed by incorporating azo 1,3,4oxadiazole with nitramino furazan and nitro-furazan have significant thermal stability (T d > 233°C, T d(RDX) = 204°C), and 1,3,4-oxadiazole is expected to contribute to the construction of thermally stable oxadiazole molecules [11]. Qian Wang combined isomers of oxadiazoles into bicyclic compounds using a simple synthetic strategy and obtained energetic materials with excellent detonation properties. (3,3′-dinitramino-4,4′-bifurazane features a detonation velocity of 9086 m s −1 ; 3′-dinitroxazafurazan, which includes two nitro-furazan fragments, has an excellent detonation performance, D = 9390 m s −1 ; P = 40.5 GPa.) In this work, we conducted systematic research on 1,3,4-oxadiazole and 1,2,5-oxadiazole as a backbone with different bridges and various energetic groups inserted (shown in Scheme 1). Density functional theory (DFT) and electrostatic potential (ESP) studies were used to analyze these compounds, along with their HOFs, densities, and detonation properties.

Computational methods
In this work, all quantum mechanical calculations were done using Gaussian 16 software [12]. Structural optimization of the designed compounds was at B3LYP/ 6-311G (d,p) level [13][14][15], with all structures reaching "no imaginary frequency" to achieve stable structures. The frontier molecular orbitals (FMOs) and the electrostatic potential were both calculated at the identical theoretical level. After that, the total energy of the molecules was determined at the B3LYP/def2-TZVPP level. The reactions and associated schemes designed to anticipate the gas-phase HOFs of the designed compounds were as follows (shown in Scheme 2): The advantage of designing isodesmic reaction calculations was that the types and numbers of various bonds of the reactants and products in the isodesmic reaction are correspondingly the same, which can effectively reduce errors in the calculation of HOFs. HOFs at a certain temperature (ΔH 298K ) was calculated by the following equation: For the ΔH 298K , HOFs in Eqs. (1) and (2) were calculated. ΔH f,P and ΔH f,R are the HOFs of the reactants and products, respectively. ΔE 0 is the single-point energy, changes in energy from product to reactant; ΔZPE is the difference between the zero-point energy (ZPE) of the target product and reactant; ΔH T is thermal correction starting with 0 will 298 K; n is energy group in numbers; Δ(PV) equals ΔnRT.
Most of the energetic compounds are not in the gas state but the solid state. In the calculation of detonation performance (heat of detonation, detonation velocity, detonation pressure), the calculation of the solid-phase HOFs (ΔH f,solid ) was the first step. According to Hess's law, the relationship between ΔH f,solid and ΔH f,gas can be expressed by the following formula [16]: The relational equation proposed by Politzer et al. can be used to calculate the ΔH sub [17,18]: where a, b, and c are constants of 0.000267 kcal mol −1 Å −4 , 1.650 kcal mol −1 , and 2.966 kcal mol −1 , respectively. A denotes the surface area of a molecule with an electron density of 0.001e Bohr −3 equivalents. ν is a measure of the balance between positive and negative regions on the surface, σ 2 tot represents a measurement of the electrostatic potential variability of the molecular surface, which can be obtained by Multiwfn [19]. The computational formula proposed by Politzer et al. [20]. can be used to calculate the density of the detonation velocity and pressure.
α, β, and γ are constants with values of 0.9183, 0.0028, and 0.0443, correspondingly. M represents the molar mass of the molecule (g mol −1 ) and V represents the volume of the molecule (m 3 mol −1 ). ν σ 2 tot indicate the same as above. Detonation performance (detonation velocity and detonation pressure) is calculated according to the Kamlet-Jacobs equation [21].
Here, D stands for blast velocity (km s −1 ), and the P indicates the burst pressure (GPa). N, M , and Q are the number of moles of initiating gases per gram of explosive (mol g −1 ), the average combined molecular weight of these gases (g mol −1 ), and the heat of detonation (cal g −1 ), respectively.

Electronic structures
The frontier molecular orbital, i.e., the highest occupied molecular orbitals (HOMO) and the lowest unoccupied molecular orbitals (LUMO), can gather suitable information on optical polarizability, dynamic stability, and reactivity [22][23][24][25]. The HOMO, LUMO, and their energy gaps (ΔE LUMO-HOMO ) of the designed structures are shown in Table 1.
From Table 1, it is not difficult to find that for HOMO, when inserting different energetic groups, the four series of A, B, C, and D have similar trends in energy level changes. Particularly for each series, HOMO energy increases with the substitution of -NH 2 and -NHNH 2, while the HOMO energy decreased when -CH(NO 2 ) 2 , -C(NO 2 ) 3 , and -NO 2 were inserted. The ranking of the contribution of each group to HOMO is as follows: -NH 2 ≈ -NHNH 2 > -N 3 > -ONO 2 > -NHNO 2 > -CH(NO 2 ) 2 >-NO 2 >-C(NO 2 ) 3 . The same is true for the LUMO energy level. Overall, the bridge has little effect on the trend of HOMO and LUMO energy changes for each designed compound. However, comparing to other bridging links, D series (-NH-NH-) has the highest LUMO energy level except for D6, which meant the backbone played a proactive role in LUMO.
The trend of ΔE LUMO-HOMO for each series can also be seen in Table 1. It is easy to find that three series A, B, and C have similar trends in ΔE LUMO-HOMO , in which the gap is relatively small when the energetic groups are -N 3 , -NHNH 2 , and -ONO 2 , the structure is relatively less stable, and the intensity of interatomic interactions is higher. It indicates a shift towards lower frequencies in their electronic absorption spectra. The difference is that in the D series, the ΔE LUMO-HOMO changes drastically, yet the D7 (4.26 eV) and B7 (4.29 eV) have similar gap values. Figure 1 shows the varied trends of the ΔE LUMO-HOMO of A series (directly linked) in the designed compounds. Clearly, all compounds of A series have comparatively larger energy gaps from 4.55 to 5.38 eV, meaning that these molecules possess good chemical stabilities. The findings indicate that the parent structure is the main influence of the ΔE LUMO-HOMO variation in series A, B, and C, the order of these series in average ΔE LUMO-HOMO can be written as follows: directly linked >-C=C->-N=N-; the energetic group is the main influence of the D series. According to previous studies, the smaller the ΔE LUMO-HOMO in a similar structure, the more unstable the structure is, which implies that the polarization rate is larger when the energy difference between occupied and non-occupied orbitals is smaller [25]. In the case where the bridges are both conjugated, the introduction of the -N=N-makes ΔE LUMO-HOMO smaller than when -C=C-is introduced. From all designed compounds, D2 (5.66 eV) has the highest ΔE LUMO-HOMO , while C2 (3.35 eV) has the lowest. In other words, the volatility of compound C2 is more evident than other compounds.

Heat of formation
The heat of formation (HOFs) is generally employed as an indicator with regard to the "energy content" of a HEDM [26]. Additionally, HOF is a key parameter that predicts the detonation properties (especially the heat of detonation) of energetic materials. In order to predict accurate HOFs, calculations are usually performed by atomization reactions (mainly for small molecules) or isodesmic reactions (mainly for complex compounds). Table 2 presents calculated total energies (E 0 ), zero-point energies (ZPE), and thermal corrections (H T ) for the reference compounds in the isodesmic reactions. Table 3 lists the total energies, ZPEs, thermal corrections, ΔH f,gas , A, v, σ 2 tot , ΔH sub , and ΔH f, solid of designed materials. Except for A8 (−97.87 kJ mol −1 ) and B8 (−67.23 kJ mol −1 ), it is shown that almost all compounds have positive ΔH f,gas range from A5 (18.95 kJ mol −1 ) to C1 (1248.67 kJ mol −1 ). About . It is noted that for each designed series, -N 3 energetic group has the most important effect in improving HOFs; however, -CH(NO 2 ) 2 and -ONO 2 both play a negative role to HOFs. Overall, after comparing with solid-phase HOFs of common HEDMs (HMX 272.6 kJ mol −1 ), 12 compounds are higher than HMX. These high HOFs can give a great contribution to the detonation properties (heat of detonation, detonation pressure, detonation velocities). Figure 2 shows the gas-phase heat of formation for the designed structures. It can be observed that the HOFs of each series change due to the variation of the high-energy groups. The four series of compounds, A, B, C, and D, have similar trends in HOF changes, and it can be

Detonation property
Oxygen balance (OB), density (ρ), the heat of detonation (Q), detonation velocity (D), and detonation (P) are the detonation properties and are shown in Table 4. Clearly, after calculating detonation performances of common energetic materials (TNT, RDX, HMX), the results are in agreement with the experimental data, proving that the calculation method is feasible. OB represents the degree to which all carbon atoms in the molecule are oxidized to carbon dioxide; hydrogen atoms are oxidized to water, and the increased production of CO 2 The densities of all designed compounds range from 1.5 to 2.0 g cm −3 . When -CH(NO 2 ) 2 , -C(NO 2 ) 3 , -NO 2 , and -ONO 2 are inserted, the density of the compounds is comparable to that of TNT, RDX, etc. However, the compounds with the substituent -C(NO 2 ) 3 have the highest densities.
Detonation velocity (D) and detonation pressure (P) are two main parameters to determine detonation performance of high-energy materials [34]. Table 4 lists detonation performance of the designed structures and representative HEDMs. Four series of designed structures have good performance in detonation velocity (between 6.10 and 9.41 km s −1 ) After a comparison, the velocity values of RDX and HMX are found to have a close value than 5 compounds (A8, B6, D6, C6, and A6), and it can be regarded that -C(NO 2 ) 3 gives the greatest contribution to the detonation velocity. This is also consistent with the results of the OB analysis. In the four series of designed objects, it can be seen that the -C=C-bond has a small detonation velocity when joined, and the -NH-NH-and -N=N-join have a small difference in detonation velocity (except D5).
Compounds with a substitution base of -NH 2 perform poorly in detonation velocity and are also relatively low in detonation pressure. The detonation pressure of the designed structure is basically in the range of 21.32-41.86 GPa, removing B4 (14.82 GPa), B2 (14.86 GPa), D2 (18.91 GPa), C2 (19.10 GPa), and A2 (19.28 GPa). Similarly, the detonation pressure is greatest at the insertion of the -C(NO 2 ) 3 energetic group. It can also be seen from Fig. 3d that the order of effect of substituents in the pressure is as follows: -C(NO 2 ) 3 >-ONO 2 ≈ -CH(NO 2 ) 2 ≈ -NHNO 2 ≈ -NO 2 >-N 3 >-NHNH 2 >-NH 2 . For detonation pressure, different from the similar level of three series A, C and D, the level of B is smaller particularly.
When there are more bi-atomic molecular gases CO and H 2 in the product, it indicates that less CO 2 and H 2 O will be produced, which will lead to a decrease in M ave and a decrease in Q [22]. It is clear from Fig. 3b that the structure with the -C(NO 2 ) 3 or -N 3 substituents in the same series has the largest Q and the smallest -NH 2 . When a single bond (A series) is added, the heat of detonation is relatively larger than the other bridging.

Thermal stability
There is an essential problem in the design and synthesis of new high-energy materials: whether HEDMs have sufficient dynamic stability to be of practical interest. The bond dissociation energy (BDE) can be regarded as an indispensable indicator in understanding thermal stability and decomposition process of high-energy materials [35]. Moreover, previous studies [36,37] on the BDE of nitro compounds especially nitroaromatic and nitramine molecules showed that the BDE of energetic material is related to its sensitivity and stability directly, and the smaller the BDE is, the more sensitive the compound is. Generally speaking, the smaller the bond dissociation energy required for bond dissociating, the weaker the bond is, which is also one of the ways to find out the trigger bond. Bond order is one of the quantitative descriptions of chemical bonds, which can be used for understanding the molecular electronic structure and predicting the molecular reactivity and stability Therefore, certain relative bonds are selected as the dissociating bond to calculate BDE according to the bond overlap population Laplacian bond order (LBO) [38]. Table 5 shows the BDE of almost designed compounds.
It is interesting that as listed in Table 5, most of the trigger bonds in the B series (-C=C-) are chemical bonds between the parent and the energetic group, and the value of BDE is relatively high (range from 249.89 to 468.45 kJ mol −1 ); thus, series B is more stable. The analysis also shows that most of the trigger bonds in series D (-NH-NH-) are bridges. For series A (directly linked) and series C (-N=N-), the chemical bonds in substitution groups become destroying bonds. Overall, it can indicate Fig. 2 The variation trends of ΔH f,gas of the designed compounds that the different bridges have a great influence on the determination of trigger bonds. Chung et al. [39] proposed that the BDE of a stable HEDM compound should be over the barrier 83.60 kJ/mol which means that almost all of the designed compounds possess suitable thermal stability (except for A3, A8, C3, C8, D2, D3, D4, D6).

Electrostatic potential
According to previous studies [23,40,41], there is a connection between the sensitivities of energetic compounds and the electrostatic potential (ESP) that is on the surface of the molecule. For ordinary organic  [31] molecules, the positive electrostatic potential will cover a larger area than the negative potential, but the latter tends to be a little stronger. However, in the case of the introduction of energetic groups (e.g., -C(NO 2 ) 3 , -NO 2 ) into the molecule, the positive potential will be strengthened and the negative potential will be weakened, so that the dominant position may become a positive potential. Figure 4 shows the ESP and surface area of representative designed compounds (series D), where the ESP is calculated based on an electron density of 0.001 a.u. (election/Bohr 3 ) and a lattice spacing of 0.25 Bohr. The local maximum and minimum values of the significant surface in the figure are shown in red and blue, respectively. It was observed that the red region (positive potential) was mainly concentrated on the parent, while the blue region (negative potential) was mainly distributed at the edges of the molecule, represented by the carbon atom on the parent oxadiazole and the oxygen atom on the nitro. The positive potential areas of the D1-D8 compounds were calculated as 149 Å2 (ratio 61%), 114 Å2 (ratio 54%), 145 Å2 (ratio 55%), 136 Å2 (ratio 58%), 153 Å2 (ratio 48%), 249 Å2 (ratio 68%), 117 Å2 (ratio 50%), and 165 Å2 (ratio 62%). Previous studies [42] reported that relatively sensitive molecules have regions of high electron indeed on covalent bonds throughout the molecular internal backbone, which are shown on the graph as a larger positive electrostatic potential on the surface area histogram. Easily, D1, D6, and D8 were the most sensitive. However, after comparing the electrostatic potential of the four series, it can be found when the substituent groups of the four series are the same, the positive potential areas of directly connected (A series) and the C and D series are not much different, while B series are slightly lower. And the comparison results also correspond to the results of the previous HOF and detonation performance. Detonation performances of the A, B, and C series are superior, the sensitivities are larger, and the situation is the opposite for the D series. Overall, -C(NO 2 ) 3 is the most effective strategy to improve the detonation performance of the design compounds (A6, D = 9.41 km s −1 , P = 41.86 GPa, Q = 1572.25 cal g −1 ). Except for B4, B2, C2, D2, A2, and B1, the detonation velocities and detonation pressures of all the design compounds were higher than those of the conventional energetic material TNT.
In summary, the designed compounds are forecast to have outstanding detonation properties and affordable sensitivity, information that would probably prove valuable to researchers interested in these compounds. The findings of this research should contribute to the development of other new high explosives designs and discoveries.