Next Article in Journal
Artemisia santolinifolia-Mediated Chemosensitization via Activation of Distinct Cell Death Modes and Suppression of STAT3/Survivin-Signaling Pathways in NSCLC
Next Article in Special Issue
Synthesis of Energetic 7-Nitro-3,5-dihydro-4H-pyrazolo[4,3-d][1,2,3]triazin-4-one Based on a Novel Hofmann-Type Rearrangement
Previous Article in Journal
Lentzeacins A-E, New Bacterial-Derived 2,5- and 2,6-Disubstituted Pyrazines from a BGC-Rich Soil Bacterium Lentzea sp. GA3-008
Previous Article in Special Issue
5-Nitrotetrazol and 1,2,4-Oxadiazole Methylene-Bridged Energetic Compounds: Synthesis, Crystal Structures and Performances
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Design of High Energy RDX-Based Derivatives: Property Prediction, Intermolecular Interactions, and Decomposition Mechanisms

Institute for Computation in Molecular and Materials Science, School of Chemistry and Chemical Engineering, Nanjing University of Science and Technology, Nanjing 210094, China
*
Author to whom correspondence should be addressed.
Molecules 2021, 26(23), 7199; https://doi.org/10.3390/molecules26237199
Submission received: 24 October 2021 / Revised: 23 November 2021 / Accepted: 25 November 2021 / Published: 27 November 2021
(This article belongs to the Special Issue Promising High-Energy-Density Materials)

Abstract

:
A series of new high-energy insensitive compounds were designed based on 1,3,5-trinitro-1,3,5-triazinane (RDX) skeleton through incorporating -N(NO2)-CH2-N(NO2)-, -N(NH2)-, -N(NO2)-, and -O- linkages. Then, their electronic structures, heats of formation, detonation properties, and impact sensitivities were analyzed and predicted using DFT. The types of intermolecular interactions between their bimolecular assemble were analyzed. The thermal decomposition of one compound with excellent performance was studied through ab initio molecular dynamics simulations. All the designed compounds exhibit excellent detonation properties superior to 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20), and lower impact sensitivity than CL-20. Thus, they may be viewed as promising candidates for high energy density compounds. Overall, our design strategy that the construction of bicyclic or cage compounds based on the RDX framework through incorporating the intermolecular linkages is very beneficial for developing novel energetic compounds with excellent detonation performance and low sensitivity.

1. Introduction

Energetic materials (EMs) are imperative parts in numerous areas due to their high energy. For example, they have applications in the military, mining, space, and other fields [1,2,3,4]. With the advancement of science and technology, there are continuing requirements for new EMs with fabulous performance and low sensitivity. Usually, a huge challenge is that the high energy and low sensitivity of EMs are mutually conditioned. Therefore, searching for new EMs with high-energy and insensitivity is in urgent need [5]. For this reason, the advancement of novel EMs has become one of the hot spots within the field of EMs.
As the most promising class of compounds as high energy density compounds (HEDCs), cyclic compounds have been pulled in significant consideration due to their favorable energy and low sensitivity [5,6,7,8,9,10,11,12,13]. For instance, the two famous explosives 1,3,5-trinitro-1,3,5-triazinane, called RDX, and 1,3,5,7-tetranitro-1,3,5,7-tetrazocane, called HMX (Figure 1). They have been broadly utilized for military and civilian applications for a long time. Specifically, RDX has good detonation properties and low impact sensitivity. Its impact sensitivity (h50) is 26 cm, detonation velocity is 8.70 km/s, and detonation pressure is 34.0 GPa [14]. The detonation performance is exceptionally near to the standard of HEDCs, whose detonation velocity is more than 9.0 km/s and detonation pressure is more than 40.0 GPa [15,16]. Accordingly, the skeleton of RDX could be regarded as a skeleton structure for the improvement of HEDCs.
It is found that joining nitro groups by means of C- or N-functionalization in the RDX skeleton may be a well-known and effective technique for achieving better energetic performances. At present, some new cyclic compounds have been successfully created as the potential candidates of HEDCs. Unfortunately, when these derivatives have more nitro groups, their energy performance has been greatly improved, but at the same time this also makes them very sensitive to external stimuli and very dangerous in experiments and industry. Furthermore, energetic compounds containing more nitro groups are difficult to be synthesized. Hence, further study to find new ways to design new HEDCs based on the RDX structure is needed.
In general, compared to monocyclic compounds, polycyclic and cage compounds such as 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20, Figure 1), etc. have huge ring or cage strain [4,17]. So, the polycyclic and cage compounds may have higher detonation properties. Pan et al. [18] reported a new method to build bicyclic and cage compounds on the skeleton of HMX and found that these novel candidates of HEDCs have both fabulous performance and low sensitivity. Therefore, building polycyclic and cage structures based on the system of RDX might lead to surprising increase in the detonation performance. Additionally, due to incomplete redox reaction, the low oxygen balance (OB) of RDX (−21.61%) leads to a lower amount of energy released. Of course, the introduction of nitro groups into newly designed RDX-based compounds to improve oxygen balance is also essential.
Within the design of HEDCs, their crystal structure must be considered [19]. However, the designed HEDCs have not been synthesized, so their crystal structure cannot be decided by tests. Therefore, it is very important to predict the crystal structure of designed HEDCs. Johnson et al. [20] demonstrated that the RDG strategy is a viable method to show the interactions, such as hydrogen bond, van der Waals force, and steric effect. Meanwhile, it is well-known that data on the decomposition process of the explosives are valuable to understand the behaviors. Recently, ab initio molecular dynamics (AIMD) strategy has been effectively employed to investigate the decomposition mechanisms [21,22]. Additionally, these studies give some valuable data to understand the decomposition process of the explosives at high temperature. For example, Xiang et al. [23] investigated the decomposition process of the 4,10-dinitro-2,6,8,12-tetraoxa-4,10-diazaisowurtzitane (TEX) crystal at high temperature of 2160 K by AIMD simulations. Wu et al. [24] performed AIMD simulations to examine the decomposition of 3,6-di(azido)-1,2,4,5-tetrazine (DiAT) at 3000 K. In spite of the fact that this research give important information to understand the decomposition reactions of the explosives, a comprehensive understanding of their decomposition process is still unclear.
In this work, seven novel high energy compounds (Figure 2) were designed using a new strategy. To begin, -N(NO2)-CH2-N(NO2)- and -N(NO2)- were selected as intermolecular linkages rather than substituents. At the same time, the -O- linkage can improve the oxygen balance effectively. In addition, three new bicyclic compounds R1, R2, and R3 were generated by connecting the C2 and C6 (or C4) atoms in the RDX structure using intermolecular linkages. After that, four new cage compounds were produced by adding -N(NH2)-, -N(NO2)-, and -O- linkages to connect the C atoms within the molecular structure of the bicyclic compounds. Then, their electronic structures, heats of formation (HOFs), energetic properties, OBs, and impact sensitivities of the designed compounds were studied carefully using density functional theory (DFT). The interaction energy, nitro charge change, electron density topology, reduction density gradient, and electron density difference analysis were considered to analyze the interactions between biomolecules and to discuss the stability of these high energy compounds. Their crystals were predicted by Polymorph and then the beginning chemical processes of the crystal of R3 at high temperature of 2500 K was simulated by AIMD. Our fundamental aim is to look for new high-energy insensitive energetic compounds.

2. Results and Discussion

2.1. Property Prediction

2.1.1. Molecular Geometry and Electronic Structures

There are two important orbitals in molecular orbitals, the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO). Additionally, their energy gaps (ΔE) can reflect the activity. It was reported that B3LYP is an effective method used to calculate the energy gaps of the energetic compounds and can give reliable qualitative variation trend for the energy gaps [6,18], so the energy gaps for the designed compounds were calculated at this level. Figure 3 shows the 3D plots of HOMO, LUMO, and ΔE of the seven designed compounds, where the positive phase was appeared in red area and the negative was appeared in green area. As shown in Figure 3, all the designed compounds have energy gaps ranging from 0.200 to 0.227 eV. Additionally, R3 has the biggest energy gap, whereas R4 has the smallest one. This is due to large differences in the molecular structures between R3 and R4. R3 has smaller cage with fewer nitro groups, whereas R4 has bigger cage with more nitro groups.
Previous studies showed that the energetic compounds are more sensitive when they have greater positive potential in the molecules [25]. Figure 4 outlines the molecular electrostatic potentials (MEPs) of the 0.001 electron/bohr3 electron density isosurface of the title compounds. The red color represents negative potential and blue signifies positive potential with the range from −0.03 to +0.05 hartree. Based on the MEPs of the designed compounds, the positive potentials are at the center of the bicyclic or cage, whereas the negative potentials locate on the NO2. Among these compounds, R4, R5, R6, and R7 have more positive potential than the other three compounds, which may lead to higher sensitivity.

2.1.2. HOFs

To estimate the detonation velocity and pressure by theoretical calculations, the value of ΔHf needs to be found first. However, it is difficult to measure ΔHf through tests, so theoretical strategy is used to predict ΔHf by means of isodesmic reactions. Table S1 of Supporting Information lists the total energies (E0), zero-point energies (ZPE), thermal correction (HT), and HOFs of the reference compounds within the isodesmic reactions.
Table S2 of Supporting Information lists the total energies, ZPEs, thermal corrections, gas-phase HOFs, heat of sublimation, and solid-phase HOFs of the designed compounds. Clearly, all the designed compounds show highly positive gas and solid-phase HOFs. Among them, the gas and solid phase HOFs of R5 are the highest, which are 1556.46 kJ/mol and 1423.97 kJ/mol, respectively. It can be seen that the design of bicyclic or cage compounds by combining high-energy groups can significantly increase the HOFs. This is as the inherent strain existing in bicyclic or cage compounds can release more energy during combustion [4,26].

2.1.3. Detonation Properties

Table 1 lists the density (ρ), heat of detonation (Q), detonation velocity (D), detonation pressure (P) and oxygen balance (OB) of the designed compounds. All the designed compounds, save R1, have higher densities than RDX. Among them, R4, R5, R6, and R7 have higher densities than HMX. However, these compounds have lower densities than CL-20. Furthermore, R5 and R7 have the highest density (1.98 g/cm3) as they have the foremost nitro groups.
As appeared in Table 1, the OB values of all the compounds are in the range from –6.78% to +7.77%, which are higher than those of RDX, HMX, and CL-20, whose values are −21.62%, −21.62%, and −10.96%, respectively. In particular, the OB values of R2 and R6 rise to zero, making them possess excellent detonation properties as the energy is discharged completely within the explosion. Additionally, all the designed compounds have higher Q than their parent compound RDX. Among them, R2 has an extraordinarily high Q value of 2291.20 cal/g.
The detonation velocities of the designed compounds vary from 9.63 to 10.27 km/s, which are much higher than those of RDX, HMX, and CL-20, whose values are 8.70 km/s, 9.10 km/s, and 9.40 km/s, respectively. Additionally, R6 has the highest detonation velocity of 10.27 km/s among these compounds, which surpasses that of CL-20. The detonation pressures of the designed compounds are from 40.86 to 48.91 GPa, much higher than those of RDX, HMX, and CL-20, which are 34.00 GPa, 39.30 GPa, and 42.00 GPa, respectively. R6 has the highest detonation pressures of 48.91 GPa among the designed compounds, over that of CL-20. In addition, the -O- linkage within the compounds R3 and R7 decreases the detonation velocity and detonation pressure altogether. Along with the increasing numbers of the nitro groups, the detonation velocities and pressures of the designed compounds improve. Overall, all the designed compounds possess surprisingly high detonation performances over RDX, HMX, and CL-20 due to the bicycle and cage structure within the molecules. Therefore, the construction of the bicyclic or cage compounds based on the RDX system can successfully improve energy properties.

2.1.4. Impact Sensitivity

Impact sensitivity (h50) is an index to evaluate the sensitivity of energetic compounds [27]. In most cases, the higher h50 is, the more stable the explosive is.
As seen in Table 1, all designed compounds have higher h50 values than CL-20, lower than the parent compound RDX but for R1 and R3, and lower than HMX but for R1 and R3. It can be concluded that R1–R7 are more stable than CL-20. Among them, R1 and R3 are more insensitive to external impact than RDX and HMX. In particular, R3 has the highest h50 value of 40.73 cm. Subsequently, the joining of the oxygen linkage can successfully decrease the sensitivity of the bicyclic and cage compounds. This is consistent with the results of the previous molecular electrostatic potential analysis.
Overall, all the designed compounds are regarded as potential candidates of HEDCs. In particular, R1 and R3 display fabulous detonation properties superior to RDX, HMX, and comparable to CL-20, and lower impact sensitivity than RDX and HMX. Consequently, R1 and R3 may be exceptionally appealing potential HEDCs.

2.2. Molecular Packing

In order to determine the molecular packing of the designed compounds, 10 space groups which account for nearly 86.49% of organic molecular solids are mainly considered. They are P21/c of 36.59%, P 1 ¯ of 16.92%, P212121 of 11.00%, C2/c of 6.95%, P21 of 6.35%, Pbca of 4.24%, Pna21 of 1.63%, CC of 0.9%, C2 of 0.9%, and Pbcn of 1.01% [28,29]. These ten possible space groups are sufficient to find the most stable polymorphs for organic compounds. The energies and densities of the packing structures within the ten most common space groups are shown in Figure S1.
Seen from Figure S1, R1 and R6 are anticipated to have the CC space group for the possible packing, with the lowest energy and highest crystal density. R3, R4, and R5 have the P 1 ¯ space group for the most possible packing. R2 has the Pna21 space group, while R7 has the P212121 space group.

2.3. Intermolecular Interactions

In this section, we will study the intermolecular interactions in the bimolecular assembles for the 7 designed compounds R1–R7 to estimate their stability.

2.3.1. Geometrical Structures

The optimized structures of the bimolecular assembles by using B3LYP/6-311G (d, p) are displayed in Figure S2 of Supporting Information. It is seen in Figure S2 that the intermolecular interactions of the bimolecular assembles mainly exist among the H, N, and O atoms.
Table 2 lists the distances of the most intermolecular interactions for the bimolecular assembles. For the intermolecular interactions between H···O, the distances are in the range of 2.44~2.86 Å except for R6 and R7, belonging to the hydrogen bonding distance. The shorter the distance is, the stronger the hydrogen bonding between O and H atoms is. Thus, there are extensive hydrogen bonding interactions in the bimolecular assembles, but the interactions in R2 and R3 are stronger.

2.3.2. Interaction Energy Analysis

The intermolecular interaction energies E(int.) of seven bimolecular assembles were calculated at the B3LYP/6-311G (d, p) level. The results were corrected by basis set superposition error (BSSE) using counterpoise correction (CP). The intermolecular interaction energies E(int.) of seven bimolecular assembles are listed in Table 3.
The interaction energies of the seven bimolecular assembles increase in the order of R3 < R1 < R2 < R5 < R6 < R7 < R4. For the most part, the smaller the interaction energy is, the more stable the assemble is. Therefore, their stability position is R3 > R1 > R2 > R5 > R6 > R7 > R4. There is hydrogen bonding within the intermolecular interactions in the bimolecular assembles. This is the same as the conclusion drawn from the impact sensitivity that R3 is most stable in these assembles.

2.3.3. AIM Topological Analysis

AIM topological analysis at the basic point of the bond is used to examine the interactions between bimolecule assembles and to analyze the interaction [30]. The total electron density ρ(r) indicates hydrogen bond effect. The Laplace values of electron density ∇2ρ(r) and energy density H(r) can reflect the type of intermolecular interactions, where the energy density H(r) is equal to the whole of the potential energy density V(r) and Lagrangian kinetic energy density G(r) when the Laplace values of electron density and energy density are negative for the covalent bond interactions but are positive for the closed shell interactions (hydrogen bond, van der Waals interaction, etc.).
As shown in Table 2, ρ(r) is within the range of 0.00079~0.009488 a.u. ∇2ρ(r) and H(r) are both positive. This means that there are closed shell interactions in the bimolecule assembles. It is found based on the ρ(r) values of the seven assembles that the hydrogen bonding in R1, R2, and R3 is stronger. In addition, the interactions between O···O and O···N still exist in the bimolecule compounds, and their intensity is less than that of the hydrogen bonding. Therefore, the weak hydrogen bonding and van der Waals interactions exist in the bimolecule assembles.

2.3.4. Reduced Density Gradients

The reduced density gradients (RDG) of these seven bimolecular assembles were analyzed in Figure 5. The scatter plot and filled isosurface plots can visually analyze weak interactions. The denser the points in the scatter diagram are, the greater the electron density is, and the stronger the effect is. Within the isosurface filling map, different colored areas represent different kinds of interactions. The blue area represents a strong attraction, such as a hydrogen bond or strong halogen bond, the green area is van der Waals force, and the red area is strong repulsion, such as steric effect.
In Figure 5, the two peaks of sign(λ2)ρ in the range of −0.03~0 a.u. indicate two types of interactions in the molecules: weak hydrogen bonding (−0.03~−0.01 a.u.) and van der Waals interaction (−0.01~0 a.u.). In addition, hydrogen bonding is the main intermolecular force. The values of sign(λ2)ρ are over 0, so there is a steric effect and the red area is shown in the bimolecular assembles. For R1, R2, and R3, the dots on the graph are denser than those of other compounds, indicating that there are stronger weak hydrogen bonding and van der Waals interactions. This is the same as previous conclusions that R1, R2, and R3 are more stable than R4−R7. Within the RDG color-filled maps, there are green areas in the oxygen atoms of the nitro groups and the hydrogen atoms, showing weak hydrogen bonds and van der Waals interactions in the bimolecular assembles. In R2, R3, R6, and R7, the light green regions appear between the oxygen atoms of the nitro groups in the bimolecular assembles, indicating the existence of O···O interactions.

2.3.5. Electron Density Difference

The difference in electron density can intuitively describe the change in density between some atoms after the formation of the bimolecular assemble, showing the weakening and enhancement of corresponding bonds [31].
In Figure 6, green and blue correspond to the density difference in isosurfaces of 0.001 and −0.001, respectively, showing the main areas where the electron density increases or decreases. Due to the hydrogen bond and van der Waals interactions, the oxygen atom area shows green, and hydrogen atom presents blue, confirming that the electron density of the hydrogen atom is reduced when the hydrogen bond is formed. The green area near part of the nitro group indicates the transfer of electrons from the hydrogen atom to the initiation bond, thereby increasing the strength of the initiation bond. It can be found that R3 is the most stable due to the larger density difference, followed by R1, R2, and R4–R7.

2.4. Thermal Decomposition

Among the seven designed compounds, R3 has excellent detonation performance and lowest impact sensitivity. Therefore, we will stimulate its decomposition process at a high temperature of 2500 K.

2.4.1. Initial Decomposition Mechanisms

There are four different initiation reactions during the thermal decomposition of the R3 crystal. As shown in Figure 7, (a) the carbon-nitrogen ring was broken through C-N cleavage. The C-N bond length increases from 1.468 to 1.768 Å. This reaction occurs at 0.008 ps of the decomposition. (b) The oxygen ring was opened by breaking the C-O bond. The C-O bond length increases from 1.396 to 1.750 Å. (c) The N-NO2 bond broke and released NO2, and the bond length increases from 1.364 to 1.717 Å. (d) The breaking of the C-H bond released H radical, and the bond length changes from 1.092 to 1.314 Å.
In order to better understand the initial decomposition, the bond dissociation energy (BDE) of these initial reactions were further analyzed by DFT calculations. According to the DFT calculated results in Figure 7, the BDE of N-NO2 is 162.11 kJ/mol and is the smallest. This is the reason that the N-NO2 bond is first to decompose.

2.4.2. Subsequent Decomposition

Four initial reactions promoted consequent decomposition of the R3 crystal together. Particularly, the H radical in path (d) could directly attack other parts and lead to their decomposition. Xiang et al. [23] explained in detail that the H radicals are very active and may cause other reactions. This is consistent with our finding in the decomposition of the R3 crystal.
Different initial decomposition reactions result in different subsequent decompositions. More detailed decomposition processes happened after the four initiation reactions were displayed in Figure 8. It is found that the cage-opening mechanisms of the molecules are different: (1) The H radical attacked the C-N bond to open the cage structure. As more than one kind of C–N bond were broken, some isomers were produced. (2) The H radical attacked the C-O bond and then opens the cage structure. The cage breaking by the C-O bond cleavage could form the six-membered ring. After the cage structure was destroyed, the C-O and C-N bonds were broken to form long chains and release NO2. All the long chains are very unstable and rapidly decomposed into small radicals. (3) The H radical attacked the N atoms within the cage and caused the N-NO2 bond to crack. The catalytic behavior of the hydrogen radical is similar to that described in previous reports on the thermal decomposition of crystalline TEX [23] and furoxan [24] by AIMD simulations. Therefore, our results can provide information for a comprehensive understanding of the thermal decomposition process for the energetic materials.
For the subsequent decomposition process, there are two kinds of reaction paths. One is that long chains gradually decomposed into small radicals. Another is to form a new ring after the ring was opened, which was broken to a certain content. The long chain has undergone a rearrangement to break down into small molecules, such as NO2, NO, and CO. The snapshots of the decomposition mechanism (i) were displayed in Figure S3 of Supporting Information. Figure S4 of Supporting Information shows the snapshots of the decomposition mechanism (iii). After the six-membered ring was formed by the C-O bond cleavage, the new cage was produced by forming a new N-N bond. Through a series of bond cleavage and recombination reactions, a new cage was formed (Figure S4b). Afterwards, the new cage was broken, and the long chain decomposed into small fragments.

2.4.3. Decomposition Products

To achieve full decomposition process intuitively, Figure 9 shows the changes in the number of main intermediates and products over time during the decomposition of the R3 crystal. As shown in Figure 9, the molecule in the R3 crystal has decomposed completely after 0.25 ps. It shows that the decomposition rate of the crystal is very fast at a temperature of 2500 K.
Seen from Figure 9, the whole decomposition process of the R3 crystal could be roughly divided into two stages. We chose 9 ps as a critical point for the decomposition of the R3 crystal. The decomposition before 9 ps is a main decomposition process, whereas that after 9 ps is a stable stage. Before 9 ps, NO2 and NO were the main decomposition products, following the mechanisms mentioned above. When the C-NO2 and N-NO2 bonds were broken, the amount of NO2 reached a maximum at about 1.5 ps, then decreased steadily, and finally disappeared after 9 ps. At this time, the amount of NO increased sharply. Most NO were produced by the N-O bond breaking in NO2. Similarly, the amount of NO gradually decreased after reaching the maximum and stabilized at around 2 molecules after 9 ps. NO existed longer than NO2. As we know, CO and N2 are the main products formed in the decomposition process. At the beginning, the number of the CO molecule increased and finally reached around 4 molecules after 6 ps. The number of N2 started to increase after 0.9 ps. At 4.4 ps, the amount of N2 reached a stable value. This stable amount of N2 and CO lasted in the remaining time. We also note that the change in the amount of N2O is related to the amount of N2. When the concentration of N2 reached a stable value, the amount of N2O began to decrease. This also proves that most of N2 came from N2O.
In the whole decomposition process, although some small molecules are not high in concentration, they cannot be ignored. NCO was formed after 4.2 ps and after that its concentration reached around to 0–4 molecules. OH was produced at 1.1 ps, then its concentration fluctuated around about 1–3 molecules, and finally disappeared after 3.1 ps. HONO was produced at 0.7 ps, then its concentration fluctuated around about 1–2 molecules, and finally disappeared after 2.5 ps. In addition, the change in the number of OH is related to that of HONO. This also indicates that most of OH was come from HONO by the O-N bond cleavage. HNCO was formed at 1.4 ps as its structure is more stable than HONO, so HNCO lasted for a long time.

3. Computational Methods

Previous studies have shown that the basis set DFT-B3LYP/6-311G (d, p) could accurately predict the structure and energy of organic compounds. Therefore, the DFT-B3LYP method with the 6-311G (d, p) basis set is used to optimize molecular structure and to predict properties [6,32,33,34].
The gas phase HOF can be estimated on the basis of the isodesmic reactions, in which the number of various bonds could be kept constant to reduce the calculation error [35,36,37]. The isodesmic reactions used to obtain the gas-phase HOFs of the designed compounds were shown in Figure S5 of Supporting Information.
For an isodesmic reaction, the heat of reaction ΔH298K at 298 K can be calculated from Equation (1).
Δ H 298 K = Δ H f , p Δ H f , R
where ΔHf,R and ΔHf,P are the HOFs of reactants and products at 298 K, respectively. At the same time, the following Formula (2) can be used to calculate ΔH298K.
Δ H 298 K = Δ E 298 K + Δ ( P V ) = Δ E 0 + Δ Z P E + Δ H T + Δ n R T
where ΔE0 is the total energy change between the products and the reactants at 0 K, ΔZPE is the difference between the zero-point energies (ZPE) of the products and reactants at 0 K, and ΔHT is thermal correction from 0 K to 298 K. The Δ(PV) value is PV work, for ideal gas reaction, it is equal to ΔnRT.
Generally, most energetic compounds are solid, so the calculation of explosive properties requires solid-phase HOF (ΔHf,solid) [38]. According to Hess’s law of constant heat summation [39], the gas-phase HOF (ΔHf,gas) and heat of sublimation (ΔHsub) can be used for evaluating solid-phase HOFs by the following expression (3).
Δ H f , solid = Δ H f , gas Δ H sub
ΔHsub can be calculated by the empirical expression suggested by Politzer et al. [24,40].
Δ H sub = a A 2 + b ( ν σ tot 2 ) 0.5 + c
where A is the surface area of the 0.001 electrons/bohr3 isosurface of the molecule electronic density, ν describes the balance between positive and negative potential on the isosurface, and σtot2 is a measure of the variability of the electrostatic potential on the molecular surface. The coefficients a, b, and c are 2.670 × 10−4 kcal/mol/A4, 1.650 kcal/mol, and 2.966 kcal/mol, respectively [41]. The descriptors A, ν, and σtot2 can be calculated using the computational procedures proposed by Bulat et al. [42]. This method has been proven to accurately estimate the heat of sublimation of many energetic compounds [43,44].
The crystal density can be obtained on the basis of the improved method proposed by Politzer et al. [45] in which the interaction index νσtot2 is introduced by expression (5).
ρ = α ( M V ( 0 . 001 ) ) + β ν ( σ tot 2 ) + γ  
where M is the molecular mass (g/mol) and V(0.001) is the volume of the 0.001 electrons/bohr3 electronic density of the molecule (cm3/molecule). The coefficients α, β, and γ are 0.9183, 0.0028, and 0.0443, respectively [46].
The detonation velocity and pressure are estimated by the Kamlet–Jacobs equations [47].
D = 1.01 ( N M ¯ 1 / 2 Q 1 / 2 ) 1 / 2 ( 1 + 1.30 ρ )
P = 1.558 ρ 2 N M ¯ 1 / 2 Q 1 / 2
where each term is defined as follows: D, detonation velocity (km/s); P, detonation pressure (GPa); N, the number of moles of detonation gases per gram explosive; M ¯ , the average molecular weight of these gases; and Q, heat of detonation (cal/g). According to the exothermic reaction principle, it can be evaluated by the differences of HOF between product and explosive; and ρ, the crystal density of explosives (g/cm3).
The impact sensitivity (h50, cm) can be predicted by the equation proposed by Pospíšil et al. [14,48].
  h 50 = m σ + 2 + n ν + q
where h50 is measured with a hammer with a mass of 2.5 kg dropping upon a sample to determine the height from which there is a 50% probability of causing an explosion [49,50]. σ + 2 is the indicator the strengths and variability of the positive surface potentials, and ν is the balance of charges between positive potentials and negative potentials on the molecular surface. The coefficients m, n, and q are −0.0064, +241.42, and −3.43, respectively.
The Monte Carlo (MC) sampling was performed for each of the space groups using fixed molecular structures. The empirical Dreiding Force Field (FF) as implemented in Polymorph code was used to minimize the most stable structures with allowing the molecules and lattice to relax. The most possible molecular packing was predicted according to both the total energy and density.
The intermolecular interaction energy could be calculated by the energy differences between the bimolecular assemble and monomer. The BSSE was corrected by CP [51]. The intermolecular interaction energy of the bimolecular assemble can be calculated from the following Equation (9).
E ( int . ) = E ( A / A ) - E ( A ) - E ( A ) + E ( BSSE )
where E(int.) is the complex interaction energy, E(A/A) is the total energy of the bimolecular assemble, E(A) is the monomer energy, E(BSSE) is the CP correction energy, and the unit is kJ/mol.
In order to obtain more accurate information about hydrogen bonding interactions, a topological analysis on a variety of molecular clusters was performed using QTAIM [52]. The electron wave functions of the molecules in the clusters were calculated at the level of B3LYP/6-31G (d, p) in the Gaussian 09 package. The electron densities (ρ(r)), Laplacian of electron densities (∇2ρ(r)), Laplacian of kinetic energy (Gb(r)), potential energy density (Vb(r)), and energy density (Hb(r)) are derived from the wave functions of particular cluster structures and can be used to characterize the types and strengths of bonding. In addition, the RDG method has been successfully used to rationally explain the interactions in the explosives. RDG can display the interactions between molecules and help us analyze the interactions between molecules [20]. All the QTAIM and RDG analyzed in this work were generated by Multiwfn 3.0.
AIMD method in the DFT framework was used to simulate the thermal decomposition of the R3 crystal under the CASTEP package, using Troullier–Martins norm-conserving pseudopotentials and the plane-wave expansion of the wave functions [53]. The plane wave cutoff was set to 900 eV for structure optimization and 300 eV for AIMD simulations. Sampling in the Brillouin area adopts the Monkhorst-Pack scheme, and k-point grid is 1 × 1 × 1. The values of the kinetic energy cutoff and the k-point grid were determined to ensure the convergence of total energies [54].
The molecular, crystal, and supercell structures of R3 were shown in Figure S6 of Supporting Information. The structure was optimized by using the Broyden, Fletcher, Goldfarb, and Shannon (BFGS) method and Perdew-urke-Ernzerhof (PBE) functional in generalized gradient approximation (GGA) [55]. Considering long simulation time, a 2 × 1 × 1 supercell was selected as a model for AIMD simulations, as shown in Figure S6(c) of Supporting Information. An NVT ensemble was selected, and the system temperature was controlled by a Nosé thermostat [56]. The simulation procedure is as follows. First, the dynamic simulation was performed for the optimized crystal structure at 298.15 K for 10 ps. Then, based on the equilibrium system, another simulation was performed for 20 ps with a time step of 1.0 fs at 2500 K in the same way. We outputted one frame every 5 fs and sum accounts to 4000 frames. We analyzed one frame every 200 fs. The time of 200 fs is much longer than the barrier crossing time of about 100 fs in the transition-state theory. All of the major decomposition products appeared and were sufficient for us to describe the primary and secondary chemical reactions within the total simulation time. The bond-distance cutoff criteria used to identify the individual species during decomposition were provided in Table S3 [57]. A chemical bond was considered to be formed between two atoms if the two-atom spacing was less than the bond length distance criterion within a minimum bond lifetime of 200 fs. Therefore, we can derive the reaction rule of reaction products with the time as well as the decomposition mechanisms. The adoption of this criterion can present qualitative conclusions in this work.

4. Conclusions

In this work, three new bicyclic and four new cage compounds were designed on the basis of the RDX system. Then, their electronic structures, HOFs, energetic properties, impact sensitivity, and intermolecular interactions in the bimolecular assembles were studied using the DFT-B3LYP method. Finally, the thermal decomposition of R3 was stimulated by AIMD.
All the designed compounds have higher densities, better detonation properties, and better oxygen balances than the parent compound RDX. Especially, they all have very excellent detonation properties. In addition, R1 and R3 show excellent detonation properties superior to RDX and HMX, comparable with CL-20, and lower impact sensitivity than RDX. Thus, all of them are promising candidates for HEDCs.
An analysis of intermolecular interactions indicates that R3 is the most stable, followed by R1, R2, and R4–R7. The R3 crystal has four initial reaction mechanisms and two subsequent decomposition paths at high temperature. After a series of ring opening and recombination reactions, some long chains and complex carbon-rich heterocyclic rings were formed, and then gradually broke into small fragments: CO, N2, NO2, NO, HNCO, and OH.
Our design strategy that the construction of bicyclic or cage compounds based on the system of RDX by incorporating the intermolecular -N(NO2)-CH2-N(NO2)-, -N(NH2)-, -N(NO2)-, and -O- linkages has been verified to be a valuable approach to develop novel HEDCs with both excellent performance and low sensitivity.

Supplementary Materials

The following are available online. Table S1: Calculated total energy (E0), zero-point energy (ZPE), thermal correction (HT), and HOFs of the reference compounds. Table S2: Calculated total energy (E0), zero-point energy (ZPE), thermal correction (HT), gas-phase HOFs (ΔHf,gas), heat of sublimation (ΔHsub), and solid-phase HOFs (ΔHf,solid) for the designed compounds. Table S3. Bond type and bond cutoff radius in the fragment analysis of the AIMD trajectory. Figure S1: Relationship between crystal density and total energy of the packing for R1–R7. Figure S2: Optimized bimolecular structures and atom numbers. Figure S3: Snapshots of the decomposition mechanisms of (i). Figure S4: Snapshots of the decomposition mechanisms of (iii). Figure S5: The corresponding isodesmic reactions of seven designed compounds (R1−R7). Figure S6: Molecular structures of R3 (a) and a unit cell (b) and a 2 × 1 × 1 supercell (c) of R3 crystal. Gray, blue, red, and white spheres stand for carbon, nitrogen, oxygen, and hydrogen atoms, respectively.

Author Contributions

Conceptualization, L.T. and W.Z.; methodology, L.T. and W.Z.; validation, L.T. and W.Z.; formal analysis, L.T.; visualization, L.T.; writing—original draft preparation, L.T.; writing—review and editing, W.Z.; supervision, W.Z.; project administration, W.Z.; funding acquisition, W.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 21773119).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 21773119).

Conflicts of Interest

The authors declare no conflict of interest.

Sample Availability

Samples of the compounds are not available from the authors.

References

  1. Thottempudi, V.; Gao, H.X.; Shreeve, J.M. Trinitromethylsubstituted 5-Nitro- or 3-Azo-1,2,4-triazoles: Synthesis, Characterization, and Energetic Properties. J. Am. Chem. Soc. 2011, 133, 6464–6471. [Google Scholar] [CrossRef] [PubMed]
  2. Wu, Q.; Tan, L.H.; Hang, Z.S.; Wang, J.Y.; Zhang, Z.W.; Zhu, W.H. A New Design Strategy on Cage Insensitive High Explosives: Symmetrically Replacing Carbon Atoms by Nitrogen Atoms Followed by the Introduction of N-oxides. RSC Adv. 2015, 5, 93607–93614. [Google Scholar] [CrossRef]
  3. Wu, Q.; Zhu, W.H.; Xiao, H.M. A New Design Strategy for High-energy Low-sensitivity Explosives: Combining Oxygen Balance Equal to Zero, A Combination of Nitro and Amino Groups, and Noxide in One Molecule of 1-amino-5-nitrotetrazole-3N-oxide. J. Mater. Chem. A 2014, 2, 13006–13015. [Google Scholar] [CrossRef]
  4. Sikder, A.K.; Sikder, N. A Review of Advanced High Performance, Insensitive and Thermally Stable Energetic Materials Emerging for Military and Space Applications. J. Hazard. Mater. 2004, 112, 1–15. [Google Scholar] [CrossRef] [PubMed]
  5. Guo, Y.Y.; Chi, W.J.; Li, Z.S.; Li, Q.S. Molecular Design of N-NO2 Substituted Cycloalkanes Derivatives Cm(N-NO2)m for Energetic Materials with High Detonation Performance and Low Impact Sensitivity. RSC Adv. 2015, 5, 38048–38055. [Google Scholar] [CrossRef]
  6. Wang, F.; Wang, G.X.; Du, H.C.; Zhang, J.Y.; Gong, X.D. Theoretical Studies on the Heats of Formation, Detonation Properties, and Pyrolysis Mechanisms of Energetic Cyclic Nitramines. J. Phys. Chem. A 2011, 115, 13858–13864. [Google Scholar] [CrossRef] [PubMed]
  7. Qiu, L.; Gong, X.D.; Ju, X.H.; Xiao, H.M. Substituent Effect on the Molecular Stability, Group Interaction, Detonation Performance, and Thermolysis Mechanism of Nitroamino-substituted Cyclopentanes and Cyclohexanes. Sci. China. Ser. B-Chem. 2008, 51, 1231–1245. [Google Scholar] [CrossRef]
  8. Qiu, L.; Xiao, H.M.; Gong, X.D.; Ju, X.H. Theoretical Studies on the Thermodynamic Properties and Detonation Performances of Bicyclic Nitramines: TNAD Isomers. Chin. J. Struc. Chem. 2006, 25, 1309–1320. [Google Scholar]
  9. Pan, Y.; Zhu, W.H.; Xiao, H.M. DFT Studies on Trinitromethylor Dinitromethyl-modified Derivatives of RDX and β-HMX. Comput. Theor. Chem. 2013, 1019, 116–124. [Google Scholar] [CrossRef]
  10. Du, M.R.; Wang, X.G.; Guo, Z.R. Theoretical Design of Bicyclo[2.2.1]heptane Derivatives for High-energy Density Compounds with Low Impact Sensitivity. Comput. Theor. Chem. 2016, 1095, 54–64. [Google Scholar] [CrossRef]
  11. Jin, X.H.; Hu, B.C.; Lu, W.; Gao, S.J.; Liu, Z.L.; Lv, C.X. Theoretical Study on a Novel High-energy Density Material 4,6,10,12-Tetranitro-5,11-bis(nitroimino)-2,8-dioxa-4,6,10,12-tetraaza-tricyclo-[7,3,0,03,7]dodecane. RSC Adv. 2014, 4, 6471–6477. [Google Scholar] [CrossRef]
  12. Qiu, L.; Xiao, H.M.; Gong, X.D.; Ju, X.H.; Zhu, W.H. Theoretical Studies on the Structures, Thermodynamic Properties, Detonation Properties, and Pyrolysis Mechanisms of Spiro Nitramines. J. Phys. Chem. A 2006, 110, 3797–3807. [Google Scholar] [CrossRef]
  13. Shen, C.; Wang, P.C.; Lu, M. Molecular Design and Property Prediction for a Series of Novel Dicyclic Cyclotrimethylene Trinitramines (RDX) Derivatized as High Energy Density Materials. J. Phys. Chem. A 2015, 119, 8250–8255. [Google Scholar] [CrossRef] [PubMed]
  14. Pospíšil, M.; Vávra, P.; Concha, M.C.; Murray, J.S.; Politer, P. A possible Crystal Volume Factor in the Impact Sensitivities of Some Energetic Compounds. J. Mol. Model. 2010, 16, 895–901. [Google Scholar] [CrossRef]
  15. Talawar, M.B.; Sivabalan, R.; Mukundan, T.; Muthurajan, H.; Sikder, A.K.; Gandhe, B.R.; Rao, A.S. Environmentally Compatible Next Generation Green Energetic Materials (GEMs). J. Hazard. Mater. 2009, 161, 589–607. [Google Scholar] [CrossRef] [PubMed]
  16. Xu, X.J.; Xiao, H.M.; Ju, X.H.; Gong, X.D.; Zhu, W.H. Computational Studies on Polynitrohexaazaadmantanes as Potential High Energy Density Materials. J. Phys. Chem. A 2006, 110, 5929–5933. [Google Scholar] [CrossRef] [PubMed]
  17. Ghule, V.D.; Jadhav, P.M.; Patil, R.S.; Radhakrishnan, S.; Soman, T. Quantum-Chemical Studies on Hexaazaisowurtzitanes. J. Phys. Chem. A 2010, 114, 498–503. [Google Scholar] [CrossRef] [PubMed]
  18. Pan, Y.; Zhu, W.H. Theoretical Design on a Series of Novel Bicyclic and Cage Nitramines as High Energy Density Compounds. J. Phys. Chem. A 2017, 121, 9163–9171. [Google Scholar] [CrossRef]
  19. Thottempudi, V.; Shreeve, J.M. Synthesis and Promising Properties of a New Family of High-Density Energetic Salts of 5-Nitro-3-trinitromethyl-1H-1,2,4-triazole and 5,5′-Bis(trinitromethyl)-3,3′-azo-1H-1,2,4-triazole. J. Am. Chem. Soc. 2011, 133, 19982–19992. [Google Scholar] [CrossRef]
  20. Johnson, E.R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A.J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498–6506. [Google Scholar] [CrossRef] [Green Version]
  21. Badgujar, D.M.; Talawar, M.B.; Asthana, S.N.; Mahulikar, P.P. Advances in Science and Technology of Modern Energetic Materials: An Overview. J. Hazard. Mater. 2008, 151, 289–305. [Google Scholar] [CrossRef]
  22. Mattsson, A.E.; Schultz, P.A.; Desjarlais, M.P.; Mattsson, T.R.; Leung, K. Designing Meaningful Density Functional Theory Calculations in Materials Science-A Primer. Modell. Simul. Mater. Sci. Eng. 2005, 13, R1. [Google Scholar] [CrossRef]
  23. Xiang, D.; Zhu, W.H. Thermal Decomposition of Isolated and Crystal 4,10-Dinitro-2,6,8,12-tetraoxa-4,10-diazaisowurtzitane According to Ab Initio Molecular Dynamics Simulations. RSC Adv. 2017, 7, 8347–8356. [Google Scholar] [CrossRef] [Green Version]
  24. Wu, Q.; Zhu, W.H.; Xiao, H.M. Ab Initio Molecular Dynamics Study of Thermal Decomposition of 3,6-Di(azido)-1,2,4,5-tetrazine. Phys. Chem. Chem. Phys. 2014, 16, 21620–21628. [Google Scholar] [CrossRef]
  25. Rice, B.M.; Hare, J.J. A Quantum Mechanical Investigation of the Relation between Impact Sensitivity and the Charge Distribution in Energetic Molecules. J. Phys. Chem. A 2002, 106, 1770–1783. [Google Scholar] [CrossRef]
  26. Zhang, J.Y.; Du, H.C.; Wang, F.; Gong, X.D.; Huang, Y.S. DFT Studies on a High Energy Density Cage Compound 4-Trinitroethyl-2,6,8,10,12-pentanitrohezaazaisowurtzitane. J. Phys. Chem. A 2011, 115, 6617–6621. [Google Scholar] [CrossRef]
  27. Li, J.S. Relationships for the Impact Sensitivities of Energetic C-Nitro Compounds Based on Bond Dissociation Energy. J. Phys. Chem. B 2010, 114, 2198–2202. [Google Scholar] [CrossRef]
  28. Mendoza-Cortes, J.L.; An, Q.; Goddard, W.A.; Ye, C.C.; Zybin, S. Prediction of the Crystal Packing of Di-tetrazine-tetraoxide (DTTO) Energetic Material. J. Comput. Chem. 2016, 37, 163–167. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Liu, Z.C.; Liu, Y.; Li, J.S.; Zhu, W.H. Prediction of Supramolecular Synthons and Crystal Packings of Supermolecular HMX/solvent Assemblies. RSC Adv. 2017, 7, 55482–55488. [Google Scholar] [CrossRef] [Green Version]
  30. Espinosa, E.; Molins, E.; Lecomte, C. Hydrogen Bond Strengths Revealed by Topological Analyses of Experimentally Observed Electron Densities. Chem. Phys. Lett. 1998, 285, 170–173. [Google Scholar]
  31. Scheiner, S.; Kar, T. Red-versus Blue-Shifting Hydrogen Bonds: Are There Fundamental Distinctions. J. Phys. Chem. A 2002, 106, 1784–1789. [Google Scholar] [CrossRef]
  32. Pan, Y.; Li, J.S.; Cheng, B.B.; Zhu, W.H.; Xiao, H.M. Computational Studies on the Heats of Formation, Energetic Properties, and Thermal Stability of Energetic Nitrogen-rich Furazano-[3,4-b]pyrazine-based Derivatives. Comput. Theor. Chem. 2012, 992, 110–119. [Google Scholar] [CrossRef]
  33. Fan, X.W.; Ju, X.H. Theoretical Studies on Four-membered Ring Compounds with NF2, ONO2, N3, and NO2 Groups. J. Comput. Chem. 2008, 29, 505–513. [Google Scholar] [CrossRef] [PubMed]
  34. Wei, T.; Zhu, W.H.; Zhang, X.W.; Li, Y.F.; Xiao, H.M. Molecular Design of 1,2,4,5-Tetrazine-Based High-Energy Density Materials. J. Phys. Chem. A 2009, 113, 9404–9412. [Google Scholar] [CrossRef] [PubMed]
  35. Wei, T.; Zhu, W.H.; Zhang, J.J.; Xiao, H.M. DFT Study on Energetic Tetrazolo-[1,5-b]-1,2,4,5-tetrazine and 1,2,4-Triazolo-[4,3-b]-1,2,4,5-tetrazine Derivatives. J. Hazard. Mater. 2010, 179, 581–590. [Google Scholar] [CrossRef] [PubMed]
  36. Pan, Y.; Zhu, W.H.; Xiao, H.M. Design and Selection of Nitrogen-rich Bridged Di-1,3,5-triazine Derivatives with High Energy and Reduced Sensitivity. J. Mol. Model. 2012, 18, 3125–3138. [Google Scholar] [CrossRef] [PubMed]
  37. Pan, Y.; Zhu, W.H.; Xiao, H.M. Theoretical Studies on the Structures, Heats of Formation, Energetic Properties and Pyrolysis Mechanisms of Nitrogen-rich Difurazano[3,4-b:3′,4′-e]piperazine Derivatives and Their Analogues. Struct. Chem. 2013, 24, 1071–1087. [Google Scholar] [CrossRef]
  38. Politzer, P.; Murray, J.S. Some Perspectives on Estimating Detonation Properties of C, H, N, O Compounds. Cent. Eur. J. Energy Mater. 2011, 8, 209–220. [Google Scholar]
  39. Atkins, P.W. Physical Chemistry; Oxford University Press: Oxford, UK, 1982. [Google Scholar]
  40. Politzer, P.; Lane, P.; Murray, J.S. Computational Characterization of a Potential Energetic Compound: 1,3,5,7-Tetranitro-2,4,6,8-tetraazacubane. Cent. Eur. J. Energy Mater. 2011, 8, 39–52. [Google Scholar]
  41. Byrd, E.F.C.; Rice, B.M. Improved Prediction of Heats of Formation of Energetic Materials Using Quantum Mechanical Calculations. J. Phys. Chem. A 2006, 110, 1005–1013. [Google Scholar] [CrossRef]
  42. Bulat, F.A.; Toro-Labbé, A.; Brinck, T.; Murray, J.S.; Politzer, P. Quantitative Analysis of Molecular Surfaces: Areas, Volumes, Electrostatic Potentials and Average Local Ionization Energies. J. Mol. Model. 2010, 16, 1679–1691. [Google Scholar] [CrossRef] [PubMed]
  43. Jaidann, M.; Roy, S.; Abou-Rachid, H.; Lussier, L.S. A DFT Theoretical Study of Heats of Formation and Detonation Properties of Nitrogen-rich Explosives. J. Hazard. Mater. 2010, 176, 165–173. [Google Scholar] [CrossRef]
  44. Wang, F.; Du, H.C.; Zhang, J.Y.; Gong, X.D. Comparative Theoretical Studies of Energetic Azo s-Triazines. J. Phys. Chem. A 2011, 115, 11852–11860. [Google Scholar] [CrossRef] [PubMed]
  45. Politzer, P.; Martinez, J.; Murray, J.S.; Concha, M.C.; Toro-Labbé, A. An Electrostatic Interaction Correction for Improved Crystal Density Prediction. Mol. Phys. 2009, 107, 2095–2101. [Google Scholar] [CrossRef]
  46. Rice, B.M.; Byrd, E.F.C. Evaluation of Electrostatic Descriptors for Predicting Crystalline Density. J. Comput. Chem. 2013, 34, 2146–2151. [Google Scholar] [CrossRef] [PubMed]
  47. Kamlet, M.J.; Jacobs, S.J. Chemistry of detonation. I. A Simple Method for Calculation Detonation Properties of C–H–N–O Explosives. J. Chem. Phys. 1968, 48, 23–35. [Google Scholar] [CrossRef]
  48. Pospíšil, M.; Vávra, P.; Concha, M.C.; Murray, J.S.; Politzer, P. Sensitivity and the Available Free Space per Molecule in the Unit Cell. J. Mol. Model. 2011, 17, 2569–2574. [Google Scholar] [CrossRef]
  49. Keshavarz, M.H. Simple Relationship for Predicting Impact Sensitivity of Nitroaromatics, Nitramines, and Nitroaliphatics. Propellants Explos. Pyrotech. 2010, 35, 175–181. [Google Scholar] [CrossRef]
  50. Li, J.S. A Multivariate Relationship for the Impact Sensitivities of Energetic N-nitrocompounds Based on Bond Dissociation Energy. J. Hazard. Mater. 2010, 174, 728–733. [Google Scholar] [CrossRef]
  51. Boys, S.F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced errors. Mol. Phys. 2006, 19, 553–566. [Google Scholar] [CrossRef]
  52. Bader, R.F.W. A Quantum Theory of Molecular Structure and Its Applications. Chem. Rev. 1991, 91, 893–928. [Google Scholar] [CrossRef]
  53. Segall, M.D.; Lindan, P.J.D.; Probert, M.J.; Pickard, C.J.; Hasnip, P.J.; Clark, S.J.; Payne, M.C. First-principles Simulation: Ideas, Illustrations and the CASTEP code. J. Phys. Condens. Matter 2002, 14, 2717–2744. [Google Scholar] [CrossRef]
  54. Wu, Q.; Zhu, W.H.; Xiao, H.M. Structural Transformations and Absorption Properties of Crystalline 7-Amino-6-nitrobenzodifuroxan under High Pressures. J. Phys. Chem. C 2013, 117, 16830–16839. [Google Scholar] [CrossRef]
  55. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. [Google Scholar] [CrossRef] [Green Version]
  57. Manaa, M.R.; Fried, L.E.; Melius, C.F.; Elstner, M.; Frauenheim, T. Decomposition of HMX at extreme conditions: A molecular dynamics simulation. J. Phys. Chem 2002, 106, 9024–9029. [Google Scholar] [CrossRef]
Figure 1. Molecular structures of RDX, HMX, and CL-20.
Figure 1. Molecular structures of RDX, HMX, and CL-20.
Molecules 26 07199 g001
Figure 2. Molecular frameworks of seven designed compounds (R1−R7).
Figure 2. Molecular frameworks of seven designed compounds (R1−R7).
Molecules 26 07199 g002
Figure 3. HOMO and LUMO energy levels and energy gaps of seven designed compounds.
Figure 3. HOMO and LUMO energy levels and energy gaps of seven designed compounds.
Molecules 26 07199 g003
Figure 4. Electrostatic potentials mapped of the designed compounds. Color coding for MEPs range from red (negative) to blue (positive).
Figure 4. Electrostatic potentials mapped of the designed compounds. Color coding for MEPs range from red (negative) to blue (positive).
Molecules 26 07199 g004
Figure 5. Scatter diagram and isosurface graph of RDG for configurations of seven bimolecular assembles.
Figure 5. Scatter diagram and isosurface graph of RDG for configurations of seven bimolecular assembles.
Molecules 26 07199 g005
Figure 6. Electron density differences of seven bimolecular assembles.
Figure 6. Electron density differences of seven bimolecular assembles.
Molecules 26 07199 g006
Figure 7. Four different initiation reactions in the thermal decomposition of the R3 crystal.
Figure 7. Four different initiation reactions in the thermal decomposition of the R3 crystal.
Molecules 26 07199 g007
Figure 8. Subsequent decomposition process of the R3 crystal.
Figure 8. Subsequent decomposition process of the R3 crystal.
Molecules 26 07199 g008
Figure 9. Time evolution of the number of main intermediates and products during the decomposition of the R3 crystal.
Figure 9. Time evolution of the number of main intermediates and products during the decomposition of the R3 crystal.
Molecules 26 07199 g009
Table 1. Calculated heats of detonation (Q), densities (ρ), oxygen balances (OB), detonation velocities (D), detonation pressures (P), and impact sensitivities (h50) for the designed compounds.
Table 1. Calculated heats of detonation (Q), densities (ρ), oxygen balances (OB), detonation velocities (D), detonation pressures (P), and impact sensitivities (h50) for the designed compounds.
Compoundρ (g/cm3)Q (cal/g)D (km/s)P (GPa)OB 1h50 (cm)
R11.782152.919.6340.86−4.5233.16
R21.882291.2010.2047.330.0025.27
R31.822089.759.6541.67−6.7840.73
R41.942138.0510.2348.46−4.1917.32
R51.982019.7310.1948.677.7713.87
R61.942223.8410.2748.910.0013.90
R71.982006.8710.1448.134.3519.68
RDX1.77
(1.82 2)
1559.99
(1501.00 3)
8.78
(8.70 2)
33.80
(34.00 2)
−21.6226.00
HMX1.84
(1.91 2)
1532.78
(1498.00 3)
9.01
(9.10 2)
36.58
(39.30 2)
−21.6230.00
CL-201.94
(2.04 2)
1671.87
(1738.20 3)
9.44
(9.40 2)
41.33
(42.00 2)
−10.9611.94
1 Oxygen balance (%) for CaHbOcNd: 1600(c-2a-b/2)/MW; MW = molecular weight of the designed compounds. 2 The experimental values were taken from Refs. [15,24]. 3 The calculated values were taken from Ref. [6].
Table 2. Topological parameters of seven bimolecular assembles at the (3, −1) critical points: calculated hydrogen bond length (d), electron densities (ρ(r)), Laplacian of electron densities (∇2ρ(r)), Laplacian of kinetic energy (Gb(r)), potential energy density (Vb(r)) and energy density (Hb(r)).
Table 2. Topological parameters of seven bimolecular assembles at the (3, −1) critical points: calculated hydrogen bond length (d), electron densities (ρ(r)), Laplacian of electron densities (∇2ρ(r)), Laplacian of kinetic energy (Gb(r)), potential energy density (Vb(r)) and energy density (Hb(r)).
AssembleInteractiondρ(r) × 1022ρ(r) × 102Gb(r) × 102Vb(r) × 102Hb(r) × 103
R1-R1H(34)···O(17)2.700.50232.00300.4067−0.31260.9407
H(30)···O(47)2.710.51772.12240.4342−0.33780.9640
H(29)···O(47)2.600.76532.73230.5946−0.50610.8850
R2-R2H(42)···O(9)2.440.94883.23140.7340−0.66010.7390
H(18)···O(32)2.500.84222.90140.6479−0.57040.7749
R3-R3H(40)···O(15)2.500.88632.92180.6633−0.59620.6713
H(20)···O(36)2.500.88642.92220.6634−0.59630.6713
R4-R4H(37)···O(23)2.780.40441.76580.3425−0.24360.9893
H(5)···O(55)2.780.40431.76540.3424−0.24350.9892
R5-R5O(44)···O(22)3.250.35141.68550.3526−0.28390.6875
O(52)···O(14)3.250.35131.68510.3526−0.28380.6872
R6-R6H(33)···O(21)2.440.78502.87780.6371−0.55470.8237
H(28)···O(52)2.860.27661.30430.2419−0.15770.8421
R7-R7H(5)···O(47)2.860.39551.78190.3393−0.23321.0616
Table 3. Intermolecular interaction energies E(int.) of seven bimolecular assembles.
Table 3. Intermolecular interaction energies E(int.) of seven bimolecular assembles.
AssemblesR1R2R3R4R5R6R7
E(int.) (kcal/mol)−8.19−7.02−9.46−4.94−6.18−6.16−5.42
E(int.) (kJ/mol)−34.27−29.37−39.58−20.67−25.86−25.77−22.68
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tang, L.; Zhu, W. Computational Design of High Energy RDX-Based Derivatives: Property Prediction, Intermolecular Interactions, and Decomposition Mechanisms. Molecules 2021, 26, 7199. https://doi.org/10.3390/molecules26237199

AMA Style

Tang L, Zhu W. Computational Design of High Energy RDX-Based Derivatives: Property Prediction, Intermolecular Interactions, and Decomposition Mechanisms. Molecules. 2021; 26(23):7199. https://doi.org/10.3390/molecules26237199

Chicago/Turabian Style

Tang, Li, and Weihua Zhu. 2021. "Computational Design of High Energy RDX-Based Derivatives: Property Prediction, Intermolecular Interactions, and Decomposition Mechanisms" Molecules 26, no. 23: 7199. https://doi.org/10.3390/molecules26237199

APA Style

Tang, L., & Zhu, W. (2021). Computational Design of High Energy RDX-Based Derivatives: Property Prediction, Intermolecular Interactions, and Decomposition Mechanisms. Molecules, 26(23), 7199. https://doi.org/10.3390/molecules26237199

Article Metrics

Back to TopTop