It is particularly a challenge to obtain a clear picture on the decomposition mechanisms of condensed phase explosives under extreme conditions due to their complex behaviors and danger in experimental measurements The atomistic simulation is an effective way to model the physical and chemical properties of complex solids under extreme conditions at the atomic level as a complement to experimental work. Recently, the application of ab initio molecular dynamics (AIMD) method in studying the decomposition mechanisms and sequences of the explosives at realistic reaction temperatures and pressures has seen tremendous progress^[1-4].

1, 3, 5-Trinitro-1, 3, 5-triazocyclohexane (RDX) is one of the most important nitramine explosives widely used for both civilian and military applications. The decomposition of RDX has been investigated by a number of groups using both experimental and theoretical methods^[5-11]. These studies, which were conducted under a variety of ambient- and high pressure environments, discuss a number of decomposition mechanisms that RDX may undergo; these include the loss of NO₂ groups^[11], the cleavage of the C—N bond to create CH₂—N₂—O₂ species, as well as whether the decomposition is a unimolecular or bimolecular process^[12]. There is an obvious difference between the two researches, so it is necessary to study the initial and subsequence decomposition mechanisms of RDX under extreme conditions further. In this work, AIMD was used to systematically study the reactive events of the initial decomposition and following reactions of crystalline α-RDX under high temperature (3000 K, which is the flame temperature of α-RDX) coupled with high pressure (34.5 GPa, which is the detonation pressure of α-RDX).

All calculations were performed in the framework of DFT^[13-14] using norm-conserving pseudopotentials^[15] and a plane-wave expansion of the wave functions based on the CASTEP code^[16]. The exchange and correlation were treated with the generalized gradient approximation (GGA), using the functional form of Perdew-Burke-Ernzerhof^[17] (PBE) method. To correct DFT for neglecting van der Waals (vdW) interaction, Grimme (G06)^[18] and Tkatchenko and Scheffler (TS)^[19] correction to the PBE functionals (labeled as PBE-G06 and PBE-TS) were used. RDX belongs to the orthorhombic space group P2₁/n with experimental lattice constants (a=13.18 Å, b=11.57 Å, c=10.71 Å) ^[20]. One unit cell consists of four molecules (168 atoms). Both geometry and cell volume were optimized using periodic DFT with vdW correction (DFT-D). A plane wave cutoff of 500 eV for molecular dynamics and 750 eV for geometry and cell optimizations was utilized. Constant temperature AIMD simulations were performed using the Car-Parrinello scheme^[21] with Nosé-Hoover thermostats^[22] for nuclei. Time steps 1 fs and 2 fs were used for the integration of the equations of motion during production and equilibration runs, respectively. The results indicate that 1 fs is much more suitable. The time step of 1 fs is chosen to longer time simulation of 20 ps. The convergence criteria were set to a 1×10 ^-6 eV energy difference for geometry optimization and a 1×10 ^-5 eV energy difference for DFT-MD simulation. he structures for the α-RDX were optimized individually before AIMD simulations. The procedure for the DFT-MD simulations was as follows: First, the system was equilibrated at 298.15 K for 10 ps using NVT. Then, AIMD simulations were carried out using NPT at 3000 K coupled with 34.5 GPa based on this equilibrated system. The simulation time was 20 ps for the AIMD simulations.

Here, the PBE-G06 and PBE-TS functionals to fully relax α-RDX crystals at ambient pressure without any constraint were applied, respectively. Table 1 lists the experimental and relaxed cell parameters of α-RDX crystal. It is found that calculated lattice parameters by the PBE-TS functional are 2.35%, 0.86%, 1.68% and 4.93% over the experimental parameters. The PBE-G06 calculation errors are 0.99%, 0.26%, 0.56%, and 1.72% for lattice constants a, b, and c, respectively. This indicates that the lattice parameters estimated by the PBE-G06 here are much closer to the experimental value than those by the PBE-TS. Therefore, the PBE-G06 functional is much suitable for studying α-RDX.

表 1(Tab. 1)

Tab. 1 Comparison of relaxed lattice parameters of α-RDX with experimental data under ambient pressure conditions method a /Å b /Å c /Å cell volume/Å3 expt.[14] 13.18 11.57 10.71 1633.90 PBE-TS 13.49 (2.35) 11.67 (0.86) 10.89 (1.68) 1714.39 (4.93) PBE-G06 13.31 (0.99) 11.60 (0.26) 10.77 (0.56) 1662.06 (1.72) Note: Values in parentheses correspond to the percentage differences relative to the experimental data.

Tab.1 Comparison of relaxed lattice parameters of α-RDX with experimental data under ambient pressure conditions

At 3000 K coupled with 34.5 GPa, the initial decomposition step of α-RDX is the unimolecular C—H bond breaking to form hydrogen radical (Fig. 1). The H radicals then induced the cleavage of two C—N bonds in other α-RDX molecule. The bond lengths of the two C—N bonds increase to 1.857 Å and 1.762 Å at first 3 fs, respectively. This suggests that the H radicals play a catalytic role in promoting subsequent decompositions.

图 1(Fig. 1)

Fig.1 Initial decomposition mechanisms of α-RDX crystal at high temperature 3000 K coupled with detonation pressure 34.2 GPa

To identify the breaking of intermolecular C—H bond, the density of states (DOS) of C and H atoms in RDX was investigated. Fig. 2 displays the DOS of the C and H states before and after the decomposition of α-RDX. It is found that the variation trends for the DOS curves of the C and H states before the decomposition of α-RDX are similar. But they have a few changes after decomposition. Three main strong peaks occur almost at the same energy in the DOS of the C and H states before and after decomposition of α-RDX. Although their DOS peak height is almost the same before its decomposition, their peaks obviously reduce after its decomposition.

图 2(Fig. 2)

Fig.2 DOS curves of the C and H states in the C—H bond before and after the initial decomposition of α-RDX

This may be due to the C—H bond cleavage in the decomposition process. Fig. 3 displays the DOS of the C and N states before and after the decomposition of α-RDX. It is seen that the variation tendency for the DOS of the C and N states before the decomposition of α-RDX is similar. All of the peaks occur at the same energy in the DOS of the C and N atoms states before its decomposition. This indicates that the two atoms are strongly bonded. However, after the decomposition of α-RDX, some main peaks do not occur at the same energy. This clearly validates that the C—N bond is broken.

图 3(Fig. 3)

Fig.3 DOS curves of the C and N states in the C—N bond before and after the initial decomposition of α-RDX

After the initial decomposition, three subsequent decomposition processes take place, as displayed in Fig. 4. Three main interesting decomposition paths include: (1) the C—N bond hemolysis triggers the other C—N bonds of this ring to break; (2) the dissociation of N—NO₂ bond releases NO₂ gas; (3) the H radical attacks the O atom to release O—H radical by forming O—H bond. However, the decomposition processes are different from reported four different primary decomposition pathways^[23-24]. Although our primary decomposition pathways do not include successive HONO elimination^[25], there exist the HONO fragments in decomposition products of α-RDX. In our simulation, it is not observed that an oxygen migrates from NO₂ to its neighboring C tom, which was reported by Botcher et al^[24]. There are a large number of subsequent generated intermediates in the unimolecular decomposition of α-RDX. Among them, Int 1, Int 2, and Int 3 are main intermediates. Compared with Int 2 and Int 3, the incidences of Int 1 are very low. This indicates that the N—NO₂ fragments are very active. There are several different heterocyclic compounds observed in the decomposition. These heterocyclic products can transform to form chain fragments. The heterocyclic C₃N₃ is most likely to occur. This demonstrates that the heterocyclic C₃N₃ is much more stable than the N—NO₂ fragments.

图 4(Fig. 4)

Fig.4 Three subsequent decomposition paths after the initiation of α-RDX crystal

Fig. 5 shows the time evolution of the population of key decomposition products during the cleavage of α-RDX. The RDX molecules delayed quickly at the beginning of the decomposition. All the nitro groups released fast at about 0.05 ps and the number of the nitro groups reaches a maximum of 10. Then the NO₂ radicals decayed quickly within about 0.1 ps. At the same time, the number of HONO radicals increased to reach a maximum of 5 at about 0.05 ps. As the HONO is a very active intermediate, it decomposed rapidly to yield HO and NO fragments. The concentration of HONO was consumed practically during the time of about 0.1 ps. After the NO₂ and HONO radicals disappeared, nitric oxide accumulated quickly to reach a maximum concentration of 5 molecules within about 0.1 ps. Afterwards they slowly reacted to form N₂, N₂O, and NOH. Fig. 5 provides a proof for the NO₂ and HONO intermediates decomposing to release NO. NO₂ can react with H to form NO. At the decomposition end of the NO₂ radicals, the concentration of N₂O₄ reaches the maximum. But they strongly fluctuated due to high oxidization activity. N₂O₄, a prototypical species of weakly bonded systems, has been the subject of intense theoretical interest^[26-28]. NO₂ radical is known to exist in the crystal as a dimer N₂O₄ and forms via addition of two NO₂ radicals^[29].

图 5(Fig. 5)

Fig.5 Time evolution of the number of RDX molecule, H2O, N2, O—H, NO2, CO2, H, NO, HONO, N2O2 and N2O during the whole decomposition stage of RDX crystal

At about 2 ps, the concentration of nitrous oxide reaches the maximum, but they strongly fluctuates due to high oxidization activity. The main intermediates during the decomposition of RDX contributed to the production of nitrogen gas. The concentration of nitrogen gas gradually increases, reaching 15 molecules at 14 ps after the initiation of RDX, then decreases 12 molecules during the whole simulation time. Among the five products, N₂O is the most unstable one and its mean lifetime survived in the system is about 11.8 ps. N₂O is induced by the hydrogen radical to form N₂. The first CO₂ was released at about 1.9 ps. The variation tendency for the concentration of CO₂ is stable at the whole decomposition stage. CO₂ is found to be in significantly smaller amounts than other species described above. Although many C—O bonds are formed, no more CO₂ is generated in the subsequent simulations. At the beginning of the simulation, the O—H fragment is formed. Then the concentration of O—H increases to 17 molecules. O—H is induced by the hydrogen radical to form H₂O. The concentration of H increases at the maximum of 4 molecules. H₂O forms fast at the beginning and the number of H₂O reaches the maximum of 10 molecules at 6 ps. Then the concentration H₂O quickly decreases and reaches stable after 12 ps.

(1) The decomposition of α-RDX at 3000 K coupled with 34.5 GPa is triggered by the homolysis of the C—H bond.

(2) The HONO elimination reaction during the decomposition is observed.

(3) The primary reactions of producing NO₂, NO, N₂O, N₂, HONO, N₂O₄, H, O—H, H₂O, and CO₂ occur at very early stages

(4) After the initiation of RDX, there are three different decomposition pathways, their decomposition processes can be described by Fig. 4 shown in the text.