Formation and growth mechanisms of natural metastable twin boundary in crystalline β-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine: a computational study

Abstract

The twin boundary (TB), a typical planar defect, occurs naturally in molecular explosives and manipulates their sensitivities to external stimuli. We systemically studied the formation and growth mechanisms of the TBs in the β-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (β-HMX) crystal by using self-consistent charge density functional tight binding with dispersion corrections. The three TB species along the [010] zone axis based upon the experiments were considered. The (001)/[010] TB species are more favorable energetically than other TB species. The twinning in β-HMX most probably occurs on the (001) plane rather than on the (101) plane, which reproduces the naturally occurred twinned crystals in the experiments well. The TB-induced symmetry breaking alters not only the geometries but also on the electronic structures of the HMX molecules located at TB. The inner surface, outer surface, and intersection of the two surfaces are suggested to play vital roles in sensitizing the condensed phase β-HMX and to act as a trigger in initiating the chemical decomposition. The HMX molecule is most likely to be adsorbed in the concave site on (001) plane through either normal or twinning pathway in a competitive manner. After the grooves on (001) surface being filled, new grooves merge naturally for further adsorption.

---

1. Introduction

Although the molecular structural characteristics of high energy density materials are the intrinsic factors in determining their physical and chemical properties, the molecular microscopic environments including the size effects, physical phases, and imperfections (vacancies, internal voids, dislocations, and grain boundaries) often play a crucial role.^1–6 HMX^7–9 (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) is an important high-energy molecular explosive due to extensive applications in high explosives, polymer-bonded explosives, solid rocket propellants, etc. Experimentally, its β form, the most stable polymorph at ambient conditions among the four polymorphs of HMX, possesses a tendency to form twinned structures at suitable conditions than other forms due to its monoclinic crystal system. Both experimental and theoretical studies found that the existence of twin boundaries tends to significantly increase the shock sensitivity of β-HMX.^10–16 Shock initiation experiments¹⁵ indicated that both coarsely and finely twinned HMX crystals have higher shock sensitivities than the non-twinned HMX with the same apparent density. Molecular dynamic simulations with reactive force field (ReaxFF) and multiscale shock technique (MSST)¹⁶ reported that the decomposition and heat release of the twinned β-HMX crystal are faster than those of the non-twinned one. These studies indicate that the disordered molecular environments in the vicinity of the twin boundary (TB) in β-HMX can change decomposition pathways and decrease activation barriers, thus leading to the disparity of sensitivity to external stimuli loaded upon the HMX crystals. However, experimental identification of the microscopic picture and twinning tendency of the TBs is typically difficult and indirect. Therefore, a comprehensive understanding of the twin-induced microscopic instability to the external stimuli in β-HMX is still far from being clear.

In this work, we studied the formation and growth mechanisms of the TBs in the β-HMX crystal using the self-consistent-charge density functional tight binding methods with dispersion corrections (SCC-DFTB-D). Our main purpose here is to realize: (i) the most probable twinning orientation of TB in β-HMX; (ii) the most possible TB structure in β-HMX for explaining the experimental observations; (iii) the internal metastability of TB for sensitizing β-HMX. Our results may shed light on the TB-induced effects on the sensitivity of condensed phase molecular explosives.

2. Computational details

Our calculations were performed with the self-consistent-charge density functional tight binding method¹⁷ with Universal Force Field (UFF)-based Lennard–Jones dispersion corrections (SCC-DFTB-D) implemented in the DFTB+ code.^18–20 This method is based on a second-order expansion of the Kohn–Sham total energy in density functional theory (DFT) with respect to charge density fluctuations, which allows description of total energies, atomic forces, and charge transfer in a self-consistent manner. Since the SCC-DFTB-D method has considerable accuracy and efficiency for studying energetic materials like nitromethane, TATB, and HMX,^21–23 this method with both the pbc^24,25 and CHNO²⁶ Slater–Koster library parameter sets was considered in this study.

The X-ray crystal structure of β-HMX²⁷ at room temperature was used as the input structure to relax without any symmetry constraints. Then based upon this structure, we established a series of TB models along the [010] zone axis in HMX. A Monkhorst–Pack k-point mesh of 5 × 2 × 3 is for the HMX unit cell and 1 × 1 × 1 for the twinning models in all the calculations. The SCC tolerance is 10⁻⁸ a.u. During the geometry relaxation, the total energy of the system was converged less than 0.02 kcal mol⁻¹, the residual force less than 0.1 kcal mol⁻¹ Å⁻¹, and the displacement of atoms less than 0.001 Å.

3. Results and discussion

3.1 Models and validation

A reasonable structural model is prerequisite for investigating the effects of the twin defects in solid β-HMX. The relaxed β-HMX unit cell was sliced into different slab models along different crystal planes (see Fig. 1a) and then the slabs were combined with their mirror images to construct 1 × 1 × 1 models with a vacuum layer of 15 Å. These models were enlarged to build 2 × 2 × 1 super cells with six-molecule layers (see Fig. 1b). To construct the TBs, the upper half layer in the above models was shifted to parallel the ab plane with translation vector r (r = σa + λb) and kept apart from the lower half layer with normal interlayer spacing between respective crystal planes. Table 1 lists the σ and λ parameters varying in the range between 0 and 1. Thus, a series of TBs with the same [010] zone axis are obtained, named the (101̄)/[010] TB, (001)/[010] TBs, and (101)/[010] TB, respectively. The obtained models contain 72 HMX molecules (2016 atoms) for (101̄)/[010] and (001)/[010] TBs and 64 HMX molecules (1792 atoms) for (101)/[010] TB. For a comparison, we also constructed three surface models with the same size of corresponding TB models along the (101̄), (001), and (101) plane, respectively.

Fig. 1 (a) Side views of twinning orientations along (101̄), (001), and (101) crystal planes, respectively. (b) 2 × 2 × 1 super cells to construct the TB along (001) orientation as an example. (c) A HMX molecule. The C atoms are shown in gray, O in red, N in blue, and H in white.

Table 1 Translation vectors in modeling TBs along different crystal planes

Model	Twinning plane	Translation vector r (r = σa+ λb)
		σ	λ
(101̄)/[010] TB	(1 0 1̄)	1/6	1/4
(001)/[010] TB_1	(0 0 1)	0	1/2
(001)/[010] TB_2		1/4	1/2
(001)/[010] TB_3		1/2	1/2
(001)/[010] TB_4		3/4	1/2
(101)/[010] TB	(1 0 1)	1/4	1/4

---

We performed some test calculations to validate the SCC-DFTB-D method for studying organic molecular explosives with extensive intermolecular interactions. The relaxed lattice parameters and crystal densities of crystalline β-HMX, RDX, ε-CL-20, TATB, and PETN by using SCC-DFTB-D method with pbc and CHNO parameter sets together with the experimental data are shown in Table S1†.^27–31 Both the two parameter sets reproduce the experimental results well. For β-HMX, both the results are in excellent agreement with the experimental data with deviations of −1.68–1.05% for pbc parameter sets and of −1.38–1.59% for CHNO parameter sets. This confirms the reliability of the SCC-DFTB-D for studying this kind of nitramino explosives. Since the pbc-based calculations have about 3 times higher computational efficiency than the CHNO-based ones, the pbc parameter sets are selected for subsequent calculations.

3.2 Energetics of TBs along different twinning orientations

The tendency to form twins in the monoclinic β-HMX was found to be dependent on the solvents from which the TBs were grown.¹¹ Previous studies proposed that the twinning occurs on the (101) plane.^10,12 Thus, to get a better understanding of the most likely twinning orientation in β-HMX, we examined the core configurations in the (101̄)/[010], (001)/[010], and (101)/[010] TBs. Fig. 2 gives the side views of the molecular configurations of both unrelaxed and relaxed TBs. Since all the relaxed (001)/[010] TBs possess quite similar configurations at the TB core (see Fig. 2d), we used (001)/[010] TB_2 (see Fig. 2b) as a representative example to discuss their differences for brevity. Fig. 2d displays the side views of the relaxed (001)/[010] TBs with different translation vectors r. It is seen that all the three TB systems present strain relaxations at the interfaces of the two slabs within the thickness of one or two molecular layers. As seen in Fig. 2, for the (101̄)/[010] and (101)/[010] TBs, the sliced surfaces are very difficult to embed with each other due to the NO₂ groups stretching out to repel from each other. In addition, both the two TBs experience serious molecular deformations at the TB core. But for the (001)/[010] TB_2, only the flexible nitro groups were observed to rotate slightly under the influence of the internal stress and the whole configurations are kept unchangeably during the relaxations.

Fig. 2 Molecular configurations of both unrelaxed (initial) and relaxed (final) (a) (101̄)/[010] TB, (b) (001)/[010] TB_2, and (c) (101)/[010] TB, (d) side views of the relaxed (001)/[010] TB species.

To evaluate the differences in forming the TBs along various orientations, we focused on the formation energetics of the TBs. First, we investigated the energetics of the surfaces before forming the twinned structures. The surface energy originated from the surface stress can affect the formation tendency of the twinning. The surface energies on different crystal planes, i.e. the (101̄), (001), and (101) plane, are calculated as:

γ = (E_Slab − nE_HMX)/2A (1)

where E_Slab is the total energy of the slab model of each surface with the same size of corresponding TB, A is the surface area, n is the number of HMX molecules in the system, and E_HMX is the normalized energy of the β-HMX unit cell per molecule. Table 2 lists the surface energies of different surfaces for β-HMX along with available theoretical³² and experimental^33,34 data. For comparison, periodic density functional theory with high precision was also used to calculate the surface energies different surfaces for β-HMX, as shown in Table 2. The SCC-DFTB-D results present the same variation trend in the surface energies as the GGA/PBE-D ones. For the (010) surface, both the dispersion-corrected (SCC-DFTB-D and GGA/PBE-D) values are quite larger than previous experimental value of 46 mJ m⁻², while the GGA/PBE value without dispersion corrections is comparable to the experiment. It is seen in Table 2 that among the five surfaces, the (001) surface possesses the lowest surface energy of 154.5 mJ m⁻². This suggests that the (001) surface is the most stable one and so has the largest tolerance to disorder (twin) among the five surfaces.

Table 2 Surface energies of different crystal planes for β-HMX compared with density functional theory with/without dispersion corrections. Units are in mJ m⁻²

Method	Crystal surface
	(101̄)	(001)	(101)	(100)	(010)
a The values were calculated at the GGA/PBE level with dispersion correction using Grimme scheme.b The values were calculated at the GGA/PBE level with no dispersion correction.c The experimental data were taken from ref. 33 and 34.
SCC-DFTB-D	225.7	154.5	231.2	243.6	179.1
GGA/PBE-D^a	170.4	118.7	177.6	199.8	138.7
GGA/PBE^b	63.1	49.7	73.5	83.9	57.5 (46.0)^c

---

Next we turn to study the twinning behaviors of these surfaces. The TB-induced strain can be evaluated by the specific interfacial energy ε, which is obtained as:

where E_TB is the total energy of the TB, E_Slab is the total energy of the surface with the same size of the TB model, and A is the surface area. Fig. 3 displays the specific interfacial energies of different TBs. In Fig. 3, the specific interfacial energies of different TBs decrease in the sequence: (101)/[010] TB > (101̄)/[010] TB > (001)/[010] TB_3 > (001)/[010] TB_1 > (001)/[010] TB_4 > (001)/[010] TB_2. This indicates that the (001)/[010] TB species with an average ε value of 71 mJ m⁻² are energetically more easy to form than other TBs with binding energies above 280 mJ m⁻². However, it is seen in Fig. 3 that the (101)/[010] TB is the most energetically unfavored than the (101̄)/[010] and (001)/[010] TBs. This result unexpectedly contradicts with previous experiments^10–12 that the (101) plane is considered as the twinning orientation.

Fig. 3 Specific interfacial energies of different TBs.

By the end, we compared the calculated TB models with the experimental observations. Due to the similarity of the twinned structures at the TB cores, the (001)/[010] TB_2 with the lowest ε value is chosen to compare with the twinned HMX crystal growing from cyclohexanone by optical microscopy (see Fig. 4).¹¹ The quasi-cross-like twinned HMX crystal is typical in many other experiments.^10–15 The calculated HMX molecules lie at a crossing angle of 84.2°, which is identical with the crossing angle of experimental twinned crystal. This confirms our conclusion that the twinned HMX crystal mostly occurs on the (001) surface rather than the (101) surface with the crossing angle of 149.5°.

Fig. 4 Side views of relaxed (001)/[010] TB_2 (σ = 1/4, λ = 1/2) compared with experimental twinned HMX crystal growing from cyclohexanone by optical microscopy.¹¹

3.3 Electronic structures

To understand how the TBs affect the electronic structure of the β-HMX solid, we compared the band structures of the representative example (001)/[010] TB_2 with the (001) surface with the same size of corresponding TB, as shown in Fig. 5. Both the valence band maxima (VBM, along the G–F and Q–Z line) and the conduction band minima (CBM, along the Z–G line) of the (001)/[010] TB_2 are smooth and slightly broaden to the lower energy region than those of the (001) surface. This slightly decreases the energy gap of TB to be 4.725 eV as compared with the gap of 4.804 eV for the surface. By and large, the (001)/[010] TB_2 keeps the band structure of the ideal bulk HMX. Thus, the interfacial interactions approach the interactions in the bulk system. This can also be supported by the band gaps as illustrated in Fig. 6.

Fig. 5 Band structures of (a) (001)/[010] TB_2 versus (b) (001) surface with the same size of corresponding TB.

Fig. 6 Band gaps of different TBs and surfaces with the same size of corresponding TBs.

Fig. 6 compares the band gaps of different TBs with corresponding surface models as well as the ideal bulk. The band gap of ideal bulk HMX was calculated to be 4.084 eV by the SCC-DFTB-D method. This is comparable with the experimental band gap (ca. 5.40 eV) of the ideal bulk HMX in acetonitrile.³⁵ The reduction in the band gaps of (101̄)/[010] and (101)/[010] TBs is found to be originated mainly from the surface relaxations and partly from the twinning. The disordered surfaces greatly influence both the VBM and CBM and so lead to the band gap reduction. But for the (001)/[010] TBs, the uncovering surface do not affect the frontier energy levels, but it is the molecular distortions at the TB core that contribute to the redistribution of the frontier energy levels. This suggests the TB core has higher activity than the (001) surface.

In Fig. 7, we focused on the intramolecular charge transfer in the HMX molecules under different molecular environments. Based on the calculated Mulliken charge analysis, the electrons were assigned to corresponding atoms. It is found that the HMX molecules in (001)/[010] TB_2 remain neutral, but the charge transfer exists inside these molecules at the TB core as well as on the outer surfaces accompanied with symmetry breaking and sequenced polarization. Fig. 7a lists the atomic charges of N and O atoms of the axial (aNO₂) and equatorial (eNO₂) nitro groups as a function of distance from the twinning plane in (001)/[010] TB_2. The horizontal lines give the comparative charges of the corresponding atoms/atom sets in the ideal bulk. The charge distribution along the c axis can be divided into three main regions: the TB core (middle), two surfaces (side), and two bulk regions (sandwiched). Anisotropy can be found to exist between the aNO₂ and eNO₂ groups of the HMX molecules in all the regions. In the bulk region, the HMX molecules possess good symmetry, while the atom pairs are ‘polarized’ at the TB core and the O atoms close to the TB are more positive and their counterparts close to the bulk region are more negative. This appears to be similar with the situation of the HMX molecules placed on the surface (see surface regions in Fig. 7). In addition, the more far the distance of the HMX molecule from the twinning plane is, the more positive the charge of the N1 atoms in the aNO₂ groups will be. But the charges of the N2 atom pairs in the eNO₂ groups present the opposite variation trend.

Fig. 7 Charge transfer in (a) nitro groups and (b) nitramino groups in (001)/[010] TB_2. Each segment denotes a pair of atoms/atom sets in symmetry in a HMX molecule. The horizontal lines give the charges of the corresponding atoms/atom sets in the ideal bulk. The vertical lines divide the TB model into three main regions: TB core (middle), two surfaces (side), and two bulk (sandwiched) regions.

Fig. 7b illustrates the charge redistribution in the nitramino groups. Both the aNO₂ and eNO₂ groups possess similar trends: the electrons flow from the nitro groups to the amino N atoms in the twinning region as well as in the surface regions. From these results, it can be concluded that the TB-induced symmetry breaking alters not only the geometries but also the electronic structures, which directly destabilize the HMX molecules placed at the TB as well as on the external surfaces.

To give an intuitive illustration of the electronic structures in TB, we calculated the partial charge densities for the VBM and CBM of the (001)/[010] TB_2 versus the (001) surface (see Fig. 8a and b). Compared with the (001) surface where both the VBM and CBM are distributed on the surface regions, the CBM of (001)/[010] TB_2 were found to be localized on the upper surface, while the VBM is shifted towards the twinning plane. The VBM of (001)/[010] TB_2 is contributed by both the aNO₂ and eNO₂ groups, while the CBM is contributed dominantly by the aNO₂ groups. The distribution of the frontier energy levels reveals the activity of HMX molecules at the internal and external surface. This is in good agreements with previous studies that the surfaces in HMX (crystal surfaces,^4–6 internal vacancies,^3,36,37 and voids²) promote the chemical decomposition of energetic materials. Their studies indicate that the chemical decompositions are more likely to start on the inner and outer surfaces. For instance, Sharia et al.^4,5 found the N–NO₂ homolysis of an HMX molecule placed at a surface or at a vacancy requires lower activation barriers of 6–7 kcal mol⁻¹ than that placed in the bulk crystal. The HONO elimination also exhibits the reduced activation barrier due to the surface effect. Zhou et al.³ found that the molecular vacancies accelerate the decomposition of HMX by increasing the reaction rate constant and reducing activation barriers. They are in agreement with our calculations by the slab models. Our studies present that the partial charge densities are localized on the surfaces rather than homogeneously distributed in the whole materials. Meanwhile, the partial charge densities analysis also shows that the VBM of the (001)/[010] TB_2 model shift towards the twinning plane from the surface, which suggests that the TB core can be more active than the surface as well as the bulk crystal. This indicates that the twinning plane possesses higher reaction activity than the surface and bulk region.

Fig. 8 Isosurfaces and energy levels of (001)/[010] TB_2 compared with (001) slab model with the same size of corresponding TBs.

In Fig. 8a and b, it seems that the electron transition in TB is much more difficult between the VBM and CBM since there is a thick molecular layer of ca. 20 Å between them. But in consideration of the intersections of the inner surface and the outer surface in real crystals, the intersectional zone is expected to notably facilitate the electron transfer between VBM and CBM. The frontier energy levels of both unrelaxed and relaxed TB and (001) surface models are listed in Fig. 8c. It is found that the frontier energy levels of (001)/[010] TB_2 is quite similar with those of the (001) surface after relaxation and its energy gap is much more wider than other TB species (also seen in Fig. 6).

In all, a simple deduction from the finite systems in this work may be proposed. For large twinned HMX crystal, hardly rejected in experiments, the HMX molecules on the surface are more active than in the bulk and those in internal TB are more active than on the surface. Thus, the inner surface, outer surface, and intersection of the two surfaces are suggested to play vital roles in sensitizing the condensed phase β-HMX and to act as a trigger in initiating its chemical decomposition.

3.4 Twinning mechanism of (001)/[010] TB_2

Twinned crystals are often formed due to the interruption during crystal formation and growth or due to the deformation or sliding on the crystal that already exists. The formation of TB in the twinned crystal is the result of an interruption or change in the lattice during crystal formation or growth due to a possible deformation from large substitutes. The tendency for twinning in β-HMX may reveal the compatibility of crystal growth via the normal and twinning pathway. To this end, we investigated the competition details in forming the TB and two kinds of adsorption sites were considered on the (001) surface.

Fig. 9 presents the energy diagram of the HMX molecules adsorbed on the convex site and on the concave site of the (001) surface via normal and twinning paths, respectively. The concave site is more energetically favorable for adsorption than the convex site. The adsorption of the HMX molecules on the convex site prefers the twinned pathway than the normal one. The adsorption energy of the first molecule requires −0.66 eV by the twinning pathway, less than −0.25 eV by the normal pathway. Further adsorption on same site becomes easier and only needs the energy of ca. −0.97 eV. However, on the concave site, the situation is just the opposite. The twinned molecules adsorbed on the concave site need higher energies (ca. 0.27 eV) than the normal one. This means that there are very small energetic differences in both the two pathways. Therefore, it may be concluded that the HMX molecule is most likely to be adsorbed on the concave site through either normal or twinning pathway in a competitive manner. After the grooves on (001) surface being filled, new grooves merge naturally for further adsorption.

Fig. 9 Energy diagram of HMX molecules adsorbed in the convex site and concave site on the (001) surface via normal path and twinning path. The value in parenthesis represents the band gap of the relaxed system. All units are in eV.

4. Conclusions

In summary, the dispersion-corrected SCC-DFTB was used to systemically study the formation and growth mechanisms of the TBs in the β-HMX crystal. The three TB species along the [010] zone axis based upon the experiments were considered. Among the six surfaces here, the (001) surface possesses the lowest energy cost to uncover this surface and the largest tolerance to disorder (twinning) on this plane. The (001)/[010] TB species with an average specific interfacial energies value of 71 mJ m⁻² are more favorable energetically than other TB species with binding energies above 280 mJ m⁻². But the (101)/[010] TB is more unfavorable energetically than the (101̄)/[010] and (001)/[010] TB. The relaxed (001)/[010] TB species match the naturally occurred twinned crystals in experiments well. Thus, it may be concluded that the twinning in β-HMX most probably occurs on the (001) plane rather than on the (101) plane.

The representative (001)/[010] TB_2 keeps the band structure of the ideal HMX bulk except for a small reduction of the band gap. The TB-induced symmetry breaking alters not only the geometries but also on the electronic structures of the HMX molecules located at TB. The HMX molecules on both the inner and outer surfaces are more active than those in the bulk. The inner surface, outer surface, and intersection of the two surfaces are suggested to play vital roles in sensitizing the condensed phase β-HMX and to act as a trigger in initiating the chemical decomposition.

The normal and twinning pathways in the crystal growth along the (001) surface are energetically competitive. This suggests that the HMX molecule is most likely to be adsorbed in the concave site on (001) plane through either normal or twinning pathway in a competitive manner. After the grooves on (001) surface being filled, new grooves merge naturally for further adsorption. Further in-depth research on TB defects in energetic materials such as twinning-induced plasticity to the external stimuli is highly desirable for developing a reliable model for chemical decompositions of HEDMs.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 21273115), the project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, and the Innovation Project for Postgraduates in Universities of Jiangsu Province (Grant No. CXZZ13_0199) for partial financial support.

References

1. R. W. Armstrong, H. L. Ammon, W. L. Elban and D. H. Tsai, Thermochim. Acta, 2002, 384, 303–313.
2. S. Boyd, J. S. Murray and P. Politzer, J. Chem. Phys., 2009, 131, 204903.
3. T. T. Zhou and F. L. Huang, J. Phys. Chem. B, 2011, 115, 278–287.
4. O. Sharia and M. M. Kuklja, J. Am. Chem. Soc., 2012, 134, 11815–11820.
5. O. Sharia, R. Tsyshevsky and M. M. Kuklja, J. Phys. Chem. Lett., 2013, 4, 730–734.
6. M. M. Kuklja, R. V. Tsyshevsky and O. Sharia, J. Am. Chem. Soc., 2014, 136, 13289–13302.
7. Q. P. Rahul, G. Wang, G. Liu and S. de, Phys. Chem. Chem. Phys., 2014, 16, 19972–19983.
8. Q. Wu, G. Xiong, W. Zhu and H. Xiao, Phys. Chem. Chem. Phys., 2015, 17, 22823–22831.
9. G. Y. Wang, Q. Peng, G. R. Liu and S. de, RSC Adv., 2015, 5, 55892–55900.
10. S. J. P. Plamer and J. E. Field, Proc. R. Soc. London, Ser. A, 1982, 383, 399–407.
11. A. E. D. M. van der Heijden and R. H. B. Bouma, Cryst. Growth Des., 2004, 4, 999–1007.
12. H. H. Cady, MRS Proceedings, Cambridge University Press, 1992, vol. 296, p. 243.
13. R. W. Armstrong, H. L. Ammon, Z. Y. Du, W. L. Elban and X. J. Zhang, MRS Proceedings, Cambridge University Press, 1992, vol. 296, p. 227.
14. D. H. Liebenberg, R. W. Armstrong and J. J. Gilman, Materials Research Society, Pittsburgh, PA, 1993, vol. 296, p. 227.
15. H. Li, R. Xu, B. Kang, J. Li, X. Zhou, C. Zhang and F. Nie, J. Appl. Phys., 2013, 113, 203519.
16. Y. Wen, X. Xue, X. Zhou, F. Guo, X. Long, Y. Zhou, H. Li and C. Zhang, J. Phys. Chem. C, 2013, 117, 24368–24374.
17. G. Seifert, D. Porezag and T. Frauenheim, Int. J. Quantum Chem., 1996, 58, 185–192.
18. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7260.
19. G. Seifert, J. Phys. Chem. A, 2007, 111, 5609–5613.
20. B. Aradi, B. Hourahine and T. Frauenheim, J. Phys. Chem. A, 2007, 111, 5678–5684.
21. M. R. Manaa, E. J. Reed, L. E. Fried and N. Goldman, J. Am. Chem. Soc., 2009, 131, 5483–5487.
22. D. Margetis, E. Kaxiras, M. Elstner, Th. Frauenheim and M. R. Manaa, J. Chem. Phys., 2002, 117, 788–799.
23. N. N. Ge, Y. K. Wei, G. F. Ji, X. R. Chen, F. Zhao and D. Q. Wei, J. Phys. Chem. B, 2012, 116, 13696–13704.
24. E. Rauls, J. Elsner, R. Gutierrez and T. Frauenheim, Solid State Commun., 1999, 111, 459.
25. C. Koehler, Z. Hajnal, P. Deak, T. Frauenheim and S. Suhai, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 085333.
26. M. J. Dewar, E. G. Zoebisch, E. F. Healy and J. J. Stewart, J. Am. Chem. Soc., 1985, 107, 3902–3909.
27. C. S. Choi and H. P. Boutin, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1970, 26, 1235–1240.
28. P. Hakey, W. Ouellette, J. Zubieta and T. Korter, Acta Crystallogr., Sect. E: Struct. Rep. Online, 2008, 64, 1428.
29. A. T. Nielsen, A. P. Chafin, S. L. Christian, D. W. Moore, M. P. Nadler, R. A. Nissan, D. J. Vanderah, R. D. Gilardi, C. F. George and J. L. Flippen-Anderson, Tetrahedron, 1998, 54, 11793–11812.
30. H. H. Cady and A. C. Larson, Acta Crystallogr., 1965, 18, 485–496.
31. E. A. Zhurova, A. I. Stash, V. G. Tsirelson, V. V. Zhurov, E. V. Bartashevich, V. A. Potemkin and A. A. Pinkerton, J. Am. Chem. Soc., 2006, 128, 14728–14734.
32. The first-principle calculations were performed using the DFT method with Vanderbilt-type ultrasoft pseudopotentials and a plane-wave expansion of the wave functions as implemented in the CASTEP code. The generalized gradient approximation with the Perdew–Burke–Ernzerhof functional (and dispersion correction using the Grimme scheme), namely GGA/PBE(-D), was employed. The three-molecule-layer slab models along different planes were used in all calculations. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868; S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
33. W. L. Elban, J. Mater. Sci., 1979, 14, 1008–1011.
34. R. Y. Yee, A. Adicoff and E. J. Dibble, Surface Properties of HMX Crystal. The 17th JANNAF Combust. Meeting. NASA, Langley Res. Center, 1980, 2, 461–468.
35. J. K. Cooper, C. D. Grant and J. Z. Zhang, J. Phys. Chem. A, 2013, 117, 6043–6051.
36. M. M. Kuklja and A. B. Kunz, J. Phys. Chem. B, 1999, 103, 8427–8431.
37. Z. Liu, Q. Wu, W. Zhu and H. Xiao, Phys. Chem. Chem. Phys., 2015, 17, 10568–10578.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra15324g

---