As an important energetic material, octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine (C₄H₈N₈O₈, HMX) has four known crystalline phases, namely, α, β, δ and γ. It has been observed that the δ form has 4 coplanar carbon atoms with a C₂ axis; α-HMX has a C₂ axis, and β-HMX has a center of inversion symmetry^[1]. The stabilities of the crystalline forms at room temperature are β>α>γ>δ, the same as the order of the densities^[2]. When heated to temperatures above 435 K, the monoclinic β-HMX converts to the hexagonal δ morphology^[3-4], the most reactive phase^[5-6]. For the transformation from β phase to δ phase, it involves major changes of both the crystal lattice and the molecular ring conformation^[3-6]. Accompanying the transition process is that a significant lattice expansion^[7] results in an appearance of cracks and large amount of hot spots in the material and the fast growth of reaction during shock compression^[8]. Those can be contributed to sensitivity increase of the δ phase^{[3-4, 8]}. The decomposition of the δ-HMX practically coincides with the β→δ phase transition and the possible decomposition of the β-phase.

Many experimental reports declare a high sensitivity of δ phase, atomistic theoretical investigations are scarce^[9-11]. Many factors, such as density, defects, electronic excitations, and particle sizes, suggest to contribute to the sensitivity of energetic materials. Among them, defects including voids, molecular vacancies, dislocations, surface, and interfaces etc. play a crucial role to initiate explosive decomposition reactions, referred to as sensitivity. At present, the investigations on defect effects mainly focus on the β-HMX^[12-14], and only a literature had studied the surface-accelerated decomposition of δ-HMX^[15]. These investigations have indicated that the surfaces, interfaces, voids, or vacancies can lower the activation barriers and accelerate kinetics of decomposition of HMX. Seeing that the special function of δ phase in the initial reaction of decomposition, dynamics characterization of δ-HMX containing defects are very important to further understand the reaction mechanism of HMX.

In our previous work^[14], we investigated the void defects contained in crystalline β-HMX. Void collapse observed with runtime was promoted by the change of the strain field and the increased system pressure due to applying canonical ensemble (NVT). In this paper, aimed at gaining a better understanding of molecular dynamics characterization of void defects, we further researched the void defects in crystalline δ-HMX applying NPT (normal pressure and temperature) simulations, including the activation barriers for an isolated molecule, void evolvement, and molecular conformations change along with simulation time. This investigation may provide some atomistic details for the observed sensitivity increase of the δ phase compared to the β phase and the degradation mechanism of HMX.

The structure of crystalline δ-HMX is hexagonal and belongs to the space group P6₁. The unit cell of the crystal is shown in Fig. 1. The initial structure was created using the crystallographic data from Ref. 4. In the δ-phase lattice, each C₄H₈N₈O₈ molecule has a boat conformation, having all four NO₂ groups on the same side of the C₄N₄ ring (Fig. 1). According to neutron diffraction experiments of β-HMX^[16], we constructed the unit cell model presented in Fig. 1. β-lattice belongs to monoclinic space group P2₁/c and contains 2 independent HMX molecules in unit cell. Each β-HMX molecule presents chair form and has a symmetry center. The β-HMX molecule with central symmetry has lower energy than δ-HMX molecule with C₂ symmetry. The super-cells used for simulation of void defects are composed of 5×5×2 unit cells for δ crystal and 5×5×6 for β one (300 molecules, 8400 atoms). After that, cohesive voids were introduced by removing 30 molecules from the 300-molecules super-cell, first one and its neighbors over an increasing distance; this corresponds to a 10% concentration of vacancies in the material.

图 1(Fig. 1)

Fig.1 Unit cell structures for β and δ phases and HMX molecular conformations in δ and β phases. Legend: carbon-gray, hydrogen-white, nitrogen-blue, oxygen-red

The active energy barriers were determined by the QST3 method, which searched for the transition structure using the Synchronous Transit-Guided Quasi-Newton (STQN) method. This method was developed by H. B. Schlegel and coworkers^[17-18], and it uses a linear synchronous transit or quadratic synchronous transit approach to get closer to the quadratic region around the transition state and then uses a quasi-Newton or eigenvector-following algorithm to complete the optimization. Density function theory (DFT) method combined with 6-31++G (d, p) basis set and B3LYP functional was applied. Convergence criterion was set to be "tight". Firstly, the geometry optimization was performed at the B3LYP/6-31++G (d, p) level. Secondly, the search for the transition state was carried out by QST3 method. Finally, the harmonic vibrational frequencies were calculated based on the optimized geometries. All quantum chemical calculations were completed through the Gaussian software^[19].

The dynamics characterization of void defectwas investigated by a MD simulation. A large number of literatures have indicated that the COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) force field^[20] can be effectively applied to the theoretical simulations for the crystalline HMX^{[14, 21-24]}. Therefore, COMPASS force field was used in our simulations. Experimental, δ-HMX is stable in the temperature range of 433-553 K^[25] and β-HMX under room temperature at an atmospheric pressure^{[2, 25-26]}. Besides, the melting point of β-HMX is in the range of 540-550 K due to different determined methods. So, the simulation temperatures were set to be 200, 300 K and 500 K, containing the stable and unstable temperatures for two forms. The Coulombic and van der Waals interactions were calculated by the standard Ewald method with quality "fine"^[27]. The charge was calculated using the QEd method^[28]. The basis of this method is the equilibration of atomic electrostatic potentials with respect to a local charge distribution. The NPT ensemble was applied with a time step of 1 fs. Temperature was controlled using the Nose thermostat and pressure by Berendsen method. The system was simulated over a total time of 1ns and the trajectory was output every 5000 time steps. The simulation was run with the commercial molecular modeling software package Materials Studio 3.0^[29].

Four unique low-energy conformers were found for HMX molecule in vacuum, two whose conformational geometries correspond closely to those found in HMX polymorphs, β and α, and two additional, lower energy conformers that are not seen in the crystalline phases, boat-chair (BC) and boat-boat (BB)^[30]. The relative energies of all conformers of an isolated molecule can be found in literature [30], namely, α>β>BB>BC. These conformations are known to exist within a narrow interval of energies.

In this work, the energies of transition states for the conformational transitions were obtained by the QST3 method at the level of B3LYP/6-31++G(d, p). The calculated active barriers were listed in Table 1 with the correction of quantum zero-point vibrational effects. The frequency calculations indicate that there is one and only one imaginary frequency for all the transition structures, and the four conformers, namely, α, β, BB, and BC, have no any imaginary frequencies. Therefore, we can determine the structures obtained by the QST3 method corresponding to the transition states.

表 1(Tab. 1)

Tab. 1 Active barriers of the conformation transitions (in kJ·mol-1) transitions β⇆α β⇆BB β⇆BC α⇆BB α⇆BC BB⇆BC forward 32.72 33.22 10.33 21.46 -3.64 4.85 backward 21.46 25.56  9.58 25.56  5.44 9.58

Tab.1 Active barriers of the conformation transitions (in kJ·mol-1)

From the active barriers listed in Table 1, we can find that the transitions from other conformations to BC have lower active barriers. Specially, the barrier from α to BC is negative, which means this conversion is very favorable in energy. In all possible transitions, the value of the highest energy barrier for the β→BB conversion is 33.22 kJ·mol^-1. Secondly, the barrier from β transits to α is 32.72 kJ·mol^-1. The two barriers are approximately 7 times higher than RT at 500 K. Other barriers are much lower than the two values. Based on the Arrhenius equation, we predict that all conversions may easily take place thermodynamically.

For the MD simulations at different temperatures, we have observed the different evolvement trends of voids inserted into two different forms, δ and β (see Fig. 2). For δ polymorph, the void has been completely collapsed at three temperatures. At 200 K and 300 K, the molecules collapsed into the void present the characteristics of partial order, but when the temperature increases up to 500 K, not only the molecules near the void but also the whole lattice melt into the liquid state. This can be further verified by the mean square displacement (MSD) analysis later (see Fig. 3). Under low temperatures of 200 K and 300 K, except for the molecules located at or near the surface of a void, the whole lattice has also large change, losing strictly periodic structure seen from Fig. 2. However, for β phase, at the temperatures of 200 K and 300 K, the void is still hold on, only the shrink of void and some shift of molecules toward the center of void occur due to the change of the strain field. When the temperature is up to 500 K, the case is the same as the δ phase, namely, the whole lattice completely collapsed into liquid state.

图 2(Fig. 2)

Fig.2 Final geometries for δ and β phases containing a void defect after 1 ns NPT simulations at 200, 300 K and 500 K, respectively

图 3(Fig. 3)

Fig.3 Calculated MSD for δ-HMX system containing a void defect at different temperatures

Compared to initial model, the density and volume of system have changed after 1ns MD simulations. Under low temperatures of 200 K and 300 K, the lattice shrink of δ-HMX results in the decrease of volume and the increase of density. For example, the density increases from 1.67 g·cm^-3 to 1.70 g·cm^-3 at 300 K and 1.73 g·cm^-3 at 200 K. The volume change accompanying void collapse is 1.53% and 3.20%, respectively. However, when the simulated temperature goes to 500 K, the lattice expands approximately 10.5%. Accordingly, the density reduces down to 1.51 g·cm^-3. Different from δ-HMX, the volume of β-HMX is on the increase with the simulated temperature. When the temperature changes from 200 K to 500 K, the volume grows from 2.1% to 16.4%. At 500 K, despite of β or δ phase, the system completely melts into liquid state and has the same density -1.15 g·cm^-3.

The different evolvement trends result from the different lattice arrangements and intermolecular interactions. The experimental value for heat of sublimation ΔH_sub is 184.76 kJ·mol^-1 and 175.90 kJ·mol^-1 for the β and δ phases, respectively^[31-32]. Structural stability and sublimation enthalpy are very useful for understanding the safety of HMX in different crystal phases. More sublimation enthalpy means that the system needs more energy to overcome the lattice bounding and inter-and intra-interactions when the solid energetic materials convert into small gaseous molecules. Besides, the formation energy of a void defect with size n, E_f(n), is also illustrated the different evolvement trends of δ and β phases. E_f(n) was calculated according to the method in Ref [33]. The value for E_f (30) in crystalline δ-HMX is calculated to be 6833.35 kJ·mol^-1, and the formation energy per molecule removed, E_f(n)/n, is 227.78 kJ·mol^-1. Applying the same method, this value of E_f (30) for β-HMX is 7009.66 kJ·mol^-1, E_f (30)/30 equals to 233.66 kJ·mol^-1. These quantitative values indicate the difference in the crystal lattice, the molecular conformation, and the microstructure of two crystalline phases. At the same time, it has also illustrated that β-HMX is more stable than δ phase. Our MD simulations show that the β-HMX containing a void defect is still more stable than δ-HMX with the same size void. Generally, the so-called "hot spots" in energetic materials are mainly the local regions of a crystal (crystal defects or deformations). Voids, molecular vacancies, dislocations, pore, impurities, and other types of defects play a crucial role in initial reactions. According to hot spots theory, the initiation of combustion in energetic materials is associated with "hot spots", in which collapse of voids may be one of the most important mechanisms the leads to ignition in pressed explosives. So, the different evolvement trends also give some atomistic details for the observed sensitivity increase of the δ phase compared to the β phase and the initial reaction mechanism of HMX.

The behavior of MSD as a function of time is presented in Fig. 3 and Fig. 4. The data in the range of 600-1000 ps were used to plot the MSD curve. The behavior of MSD as a function of time can be used to discriminate between a solid and a liquid phase. The MSD of the solid system oscillates about a mean value. For a fluid, without underlying regular structure, the MSD gradually increases with time. So, from Fig. 3, we can see that the δ-HMX crystal with a void defect presents part liquid-state feature at 200 K and 300 K, which is agreement with the geometries shown in Fig. 2. When simulated temperature rises to 500 K, the system exhibits normal linear diffusive behavior, and this means that the lattice has completely collapsed into liquid state. Diffusion coefficient (D) is calculated by the well-known Einstein's relation, which obtains diffusion coefficients as a tangent of the linear region of the displacement function. The Einstein's relation is expressed as:

图 4(Fig. 4)

Fig.4 Calculated MSD for β-HMX system containing a void defect at different temperatures

where r_i(t) denotes the position vector at time t, r_i(t₀) the position vector at time origin, the angular bracket denotes averaging over all choices of time origin within a dynamics trajectory, and N is the total number of the diffusion particles. The diffusion coefficient at 500 K was calculated to be 1.2×10^-10 m²·s^-1, significantly above the diffusion coefficient of baseline molecular of the order of 10^-12 m²·s^-1. This value is observed for a molecule in the defect-free lattice and restricted to thermal fluctuations around its equilibrium site.

For β phase, the MSD curves are shown in Fig. 4. At the temperature of 200 K, MSD oscillates about a mean value, showing an obvious solid-phase characteristic. At 300 K, we can still observe the oscillation phenomenon, but the MSD value increases slowly with runtime, indicating the lattice binding slightly reduces. When temperature rises up to 500 K, the case is very similar with the δ-HMX. Namely, the system finally changes into liquid state after 1 ns MD simulation. The difference is that normal linear region appears only after approximately 800 ps. This illustrates the rate of lattice collapsing is slower than that of δ phase. Besides, the diffusion coefficient is determined using the data of linear region 800-1000 ps, and the value is 1.5×10^-10 m²·s^-1, almost the same as that of δ phase.

For collapsed liquid-state molecules, tracing the single molecule along with runtime, it can be observed the frequent transitions of molecule between different conformations. To illustrate this, we cut out a section of snapshots from MD simulation with a time period of 20 ps (see Fig. 5). During 20 ps runtime, we have observed the obvious conformation transitions, namely from BB to α. Besides, we captured the transition state or intermediated structures undergone by this conversion, presenting the reversion of two downward NO₂ groups and adjustment of ring structure. When the system hasn′t completely melted into liquid state, for example the molecules collapsed into the void in δ crystal lattice under 200 K and 300 K, we can still observe the conformational transition, but the conversion rate is much slower due to the effects of the crystal field.

图 5(Fig. 5)

Fig.5 Change of molecular conformations during 20 ps from the MD simulation for liquid-state system at 500 K

The frequent transitions can be explained by the active barriers and the formation energy of the void defect. In section 3.2, we had already calculated the formation energies of voids inserted into β and δ crystal lattices, and they are 7009.66 and 6833.35 kJ·mol^-1, respectively. Combined with the active barriers listed in Table 1, we can predict that the formation energy is enough to induce conformational changes, rotations, and center-of-mass translations of the molecules in system, and even the melting of crystal lattice.

At the simulation temperature of 500 K, no matter δ or β phase, the system containing a void has collapsed into liquid state going through 1 ns simulation, and has the same density of 1.51 g·cm^-3. So, we have thought the conformation distribution of final liquid system to be roughly similar for two crystal phases. And that is true through our conformational analysis. We find that some molecules are in stable conformations(about 40%), but a large number of molecules with transition state or intermediate structures(approximately 60%). In the different conformations, the ratio of BC is the largest(approximately 18%), next BB(approximately 12%), and the proportion of α and β is relatively low(approximately 4% and 6%, respectively). The reason may be the combination of conversion barrier and the conformation stability. Such as, all the active barriers conversing to BC are relatively lower (see Table 1), and the energy of BC is the global minimum, which result in BC with highest ratio. Besides, a larger number of transition state and intermediate structures also make the collapsed system higher energy and more active. It can provide some evidence for the view that the collapse of voids seems to be one of the most important mechanisms that can lead to ignition in pressed explosives.

Under lower temperatures of 200 K and 300 K, for the molecules collapsed into the void in δ phase, we observe the conformation of most molecules is not change, still keeping the α conformation, owing to the effect of crystal field. Secondly, the transition state and intermediate are dominant in the rest of collapsed molecules. Besides, minority molecules are found with BC or BB structure.

For β phase, we have observed that the molecules on the surface of a void shift to the center of void at 200 K and 300 K from the Fig. 2. The anisotropic interactions with others in their vicinities lead to variations in numbers and strengths of interactions as well as possible conformational changes. We find the conformational transition is mainly from β to BC for the shifted molecules. It agrees with the lowest active barrier of β→BC among all transitions from β to other conformations (10.33 kJ·mol^-1, Table 1). Of course, the transition state and intermediates structures are also observed on the surface of void. These molecules have higher energy than that of β-HMX molecule, and will show more activity to initiate the degradation reaction of HMX. It had been demonstrated by DFT calculations coupled with transition state theory ^[12].

MD simulation for crystalline δ-HMX with a void defect was performed to study the dynamics characteristic of void evolvement and molecular conformation change, compared with β phase. The calculations for active energies of conformational transitions for the isolated molecules indicate that the higher barriers are 33.22 and 32.72 kJ·mol^-1, and others are much lower than the two values. This means that all conversions can easily take place in energy.

MD simulation shows that the crystalline system containing 10% vacancy concentration completely collapses into liquid phase after 1 ns NPT simulation at 500 K, no matter δ or β phase. Compared to initial models, the lattice expands result in the system density down to 1.51 g·cm^-3. For the liquid system, four conformations, namely, α, β, BB, and BC, have been observed. Among them, the proportion of BC is the largest. Besides, a large number of transition state and intermediate molecules are observed. Tracing the single molecule along with the runtime, frequent transitions between conformational states are found due to low energy barrier.

At 200 K and 300 K, different evolvement trends are observed for δ and β phases. For δ phase, the void has completely collapsed, and the whole crystal lattice has also lost the rigid periodicity. However, void and crystal lattice of β phase have been hold on, and there only occur the shrink of void and some shift of molecules toward the center of void. Contrary to the volume shrink of δ phase, lattice expansion is observed for β phase. Besides, for the molecules collapsed into void or on the surface of void, except for the stable conformations, transition state or intermediate molecules have been found. Different evolvement behaviors of the voids contained in δ and β phases may provide some important information to explain the sensitivity difference of two crystalline phases and the initial reaction mechanism of HMX.