Effect of Cocrystal Behavior on Sensitivity and Thermal Decomposition Mechanisms of CL-20/HMX via Molecular Dynamics Simulations

Yizheng Fu,^*1,2 Penghui Zhao,^1,2 Luxia Yang,³ Ruizheng Miao,⁴ Congyun Zhang,^1,2 Zhanhu Guo⁵ and Yaqing Liu*^1,2

¹School of Materials Science and Engineering, North University of China, Taiyuan 030051, China

²Shanxi Province Key Laboratory of Functional Nanocomposites, North University of China, Taiyuan 030051, China

³Business College of Shanxi University, Taiyuan 030051, China

⁴Xi’ an North Huian Chemical Industries Co. Ltd. ,Xi’ an 710302,China

⁵Integrated Composites Laboratory (ICL), Department of Chemical & Biomolecular Engineering, University of Tennessee, Knoxville, Tennessee 37996, USA

Abstract 

Molecular simulations were performed to investigate the effect of cocrystal behavior on the properties of 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane/1,3,5,7-tetranitro-1,3,5,7 tetrazacy-clooctane (CL-20/HMX) explosive. The sensitivity, binding energy and the mechanical properties of the cocrystals were simulated by molecular dynamics (MD) with COMPASS force field, and the thermal decomposition mechanisms were studied by reactive molecular dynamics (RMD) simulations with ReaxFF/lg. The results showed the binding energy between CL-20 and HMX in cocrystal structure was larger than mixture structure, indicating the more excellent stability of the former. Compared with pure crystal of CL-20, the value of tensile (E), bulk (K), shear (G) moduli of CL-20/HMX cocrystals and the physical mixture all decreased. Thus, both mixture and cocrystal structure reduced the stiffness and increased flexibility of CL-20, and they showed better security. RMD simulation showed that the potential energy change of mixture was larger than that of the cocrystal at the earlier thermal decomposition stage. However, the potential energy drops faster in cocrystal than in mixture in the later process due to the interactions of the products from CL-20 decomposition. The main products were NO₂, N₂, NO, H₂O, HONO, HON and CO₂.

Table of Content

The sensitivity, binding energy,  mechanical properties and the thermal decomposition mechanisms of the CL-20/HMX cocrystals were investigated through molecular dynamics (MD).

 

 

Keywords:CL-20/HMX;Sensitivity;Mechanical properties;Thermal decomposition mechanism

---

1. Introduction

     2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane, commonly referred to CL-20, with its rigid isowurtzitane cage and nitro groups being attached to the bridging nitrogen atoms, is the most powerful explosive material presently available, and has a promising prospect of application.^1-3 However, the highly energetic materials suffer from the relatively poor safety and high economic cost, which seriously hinder development and application of CL-20.^{4, 5}

       The cocrystallization of explosive, which is a multicomponent crystal of several neutral explosive molecules formed by non-covalent interactions including H-bond, electrostatic interaction in a defined ratio,^{6, 7} has drawn great attention because it can satisfy the demands for both power and safety somewhat. Recently a lot of cocrystal explosives have been synthesized and characterized, such as CL-20/1,3-dinitrobenzene (DNB)⁸, CL-20/2,4,6-trinitrotoluene (TNT),⁹ CL-20/Benzotrifuroxan (BTF),¹⁰ CL-20/1,3,5-triamino-2,4,6-trinitrobenzene (TATB)¹¹ and CL-20/caprolactam (CPL).¹²Bolton et al¹³cocrystallized CL-20 with the cheap and insensitive explosive HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazacyclooctane) to form a novel cocrystal explosive. The prepared cocrystal explosive with a lower sensitivity, makes up the deficiency of conventional energetic materials.

      Energetic materials may release a lot of energy in a short time through complex chemical reactions under extreme conditions, which would bring out some challenge in time and space. Thus, it is difficult to provide some microscopic information from the atomic or molecular scale in current experimental conditions. Molecular dynamics (MD) simulation is a newly analytical tool used to describe the detailed information at atomic level and femtosecond scale. Especially, the reaction molecular dynamics (RMD) simulation method based on the reactive force field, as the bridge to connect quantum chemistry and molecular dynamics methods, could provide a method to investigate property and chemical reaction of energetic materials in atomic level. They can not only facilitate the understanding of reaction mechanism of energetic materials in extreme conditions, but also help to obtain the energy release law.¹³

        Rom¹⁴studied the decomposition mechanism of hot liquid nitromethane at various compressions via ReaxFF molecular dynamics simulations. It has been demonstrated that there were differences of thermal decomposition between higher and lower compression. Han¹⁵showed that the reaction path was seriously influenced by the temperature through investigating thermal decomposition and subsequent reaction of the energetic material nitromethane (CH₃NO₂) using molecular dynamics with ReaxFF forcefield. Strachan¹⁶ reported the thermal induced chemistry in cyclotrimethylene trinitramine (RDX) at different temperatures and densities. With time evolution, the change of potential energy can be described well through a single exponential function that increases with declining density and shows an Arrhenius temperature dependence. Weismiller¹⁷ studied the chemical kinetics of ammonia borane oxidation. The results elucidated the pertinent chemical pathways, intermediate species and chemical kinetic mechanism. Strachan¹⁸demonstrated shock propagation of RDX through nonequilibrium MD simulations at different collision velocities, high impact velocities (>6 km/s). The RDX molecules decomposed and reacted to form all kinds of small molecules in very short time scales (<3 ps), and at lower velocities, only NO₂ was formed in agreement with experiment results. Zhang¹⁹ studied the initial thermal decomposition pathways of supercell structure and signal molecule of CL-20 explosive at various densities and temperatures by RMD. The results showed that the density and temperature affected the rate of pyrolysis, but exerted no effects on pyrolysis mechanism. Compared with the pure energetic materials, fewer references were reported in CL-20 eutectic of energetic materials by RMD simulation. Liu²⁰reported the initial thermal decomposition of condensed phase CL20/TNT cocrystal under high temperature. Products identification analysis showed that the main products were NO₂, NO, CO₂, N₂, H₂O, HON, HNO₃. Goddard²¹separately compared the thermal decomposition process of the TNT/CL-20 cocrystal with that of pure crystals of TNT, CL-20, a simple physical mixture of TNT and CL-20. They found that cocrystal exhibited a lower decomposition rate than pure CL-20, suggesting the decrease in sensitivity of the cocrystal. These all could provide us with a rational method to study the thermal decomposition mechanism of CL-20/HMX cocrystal.

     In this study, molecular dynamics simulation method was employed to study the pure CL-20, HMX, CL-20/HMX mixture and CL-20/HMX cocrystal with COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) forcefield²²，and the sensitivity and mechanical property of these structures were analyzed. Then, the thermal decomposition of these systems were studied by RMD with ReaxFF/lg forcefield.²³Therefore, the effect cocrystal behavior on the thermal decomposition mechanism of energetic materials was obtained.

 

2. Modelling and simulation details

2.1. Modelling construction

      CL-20 and HMX molecular structures are displayed in Figure 1. In this paper, Materials Studio 6.0 software was used to build all molecular structures, the initial unit cell structures of cocrystal CL-20/HMX,¹³CL-20²⁴and HMX²⁵ were taken from Cambridge Crystallographic Data Cenre. The CL-20/HMX cocrystalwith supercell 2×3×2 (Figure 2a), was constructed for avoiding the influence of temperature and heating rate fluctuation on reaction and obtaining more accurate statistical results. For comparison, pure CL-20 with supercell 3*2*2 (Figure 2b), pure HMX with supercell 2×3×2 (Figure 2c) and hybrid CL-20/HMX mixture (Figure 2d) with molar ratio 2:1 were constructed. The mixture was constructed by using Build Layers tool and put 3×4×2 CL-20 supercell and 4×3×1 HMX supercell into a supercell system.

Figure1. The molecular formulas for CL-20 (a) and HMX (b)

Figure 2The molecular structures for CL-20/HMX cocrystal(a), Pure CL-20(b), Pure HMX(c) and CL-20-HMX mixture(d)

2. 2 Simulation details

    After that, each unit cell was subjected to the energy minimization through smart minimizing method to remove the local non-equilibrium structure. To further relax local hot-spots and achieve better equilibrium of the system, MD simulations were performed at 300 K for 500 ps in the NVT canonical ensemble, and at 300 K, 1 bar for 500 ps in the NPT ensemble, successively. Trajectories were saved every 1 ps and the final 50 ps configurations were used to analyze the bond length, binding energy and mechanical property. Temperature and pressure were controlled by the Andersen²⁶ and Berendsen’s²⁷ methods respectively. The Ewald summation²⁸ was used to calculate the Coulombic interactions. The atom-based summation²⁹ was applied to calculate van der Waals’s interactions with a cutoff distance of 0.95 nm, a spline width of 0.1 nm and a buffer width of 0.5 nm. The time step was set as 1 fs. During the whole simulation process, COMPASS forcefield was used for computing interatomic interactions. 

    Large-scale atomic/molecular massive parallel simulator (Lammps)software was chosen to investigatethe thermal decomposition of cocrystal, pure crystals, and simple mixtures. RMD simulations with ReaxFF/lg forcefield were performed at 300 K for 500 ps in NVT ensemble, and at 300 K, 1 bar for 10 ps in the NPT ensemble, successively. Finally, at 2000 K, RMD simulation was performed for 50 ps in NVT ensemble, and Berendsen method was used for the temperature and pressure controlling. Therein, periodic boundary condition was chosen, and coupling parameter andthe whole simulation time were set as 100 fs, respectively. Trajectories were saved every 50 fs and band orders were set as 0.3.

 

3. Results and discussion

3.1. The trigger bonds affected by cocrystalization and mixture

        Based on a comparison of the bond orders of trigger bonds, a criterion³⁰ can be used to judge the relative stability and security of energetic compounds. As the “principle of the smallest bond order” (PSBO) indicates, for the trigger bond of energetic compound, its bond order iscloselyrelated to the sensitivity of the compound. Namely, the smaller the order is, the bigger the sensitivity is. This is especially suitable for the compound with similar molecular structures. Generally, the decrease of the bond order will lead to the increase of the bond length, and vice versa. Since the bond order and length decide the strength of the energetic compound, the statistical distribution of bond length obtained through classical MD simulation can be employed to evaluate the compound stability. The longer the bond length is, the more likely the bond is to break. This will enhance the molecule’s activity and decomposition. The molecule with the maximum bond length has the greatest activity and can be most easily activated and decomposed.

       The bond length of  are calculated N–NO₂ trigger bond in CL-20 and N-NO₂ trigger bond in HMX molecules are calculated. Table 1 shows the average bond length (L_ave) and maximum bond length (L_max) of trigger bondsin different simulation systems. It is observed that CL-20/HMX with mixture and cocrystal modes exert a slight influence on the average bond length comparing to pure CL-20 or pure HMX, whereas have an obvious effect on the maximum bond length. The mixture and CL-20/HMX cocrystalmade the maximum bond length of trigger bond N-NO₂ in CL-20 become shorter and N-NO₂ in HMX elongated. The sensitivity is closely connected with the maximum bond length of trigger bond, so the mixture and CL-20/HMX cocrystalreduced the sensitivity of CL-20 and increased the sensitivity of HMX. Otherwise, effects of cocrystal on the trigger bonds are more obvious than that of mixture.

Table 1The L_max and L_ave of trigger bonds of N-NO₂ in CL-20 and N-NO₂ in HMX (Å)

System	Lave of N-NO2 of CL-20	Lave of N-NO2 ofHMX	Lmax of N-NO2 ofCL-20	Lmax of N-NO2 of HMX
Pure CL-20	1.395	----	1.551	----
Pure HMX	----	1.397	----	1.518
CL-20/HMX Mixture	1.395	1.399	1.548	1.534
CL-20/HMX Cocrystal	1.392	1.398	1.537	1.541

 

3.2. The binding energy

      Binding energy (E_bind) is defined as the negative value of the intermolecular interaction energy (E_inter), E_bind =- E_inter, which can well reveal the interaction between the two components³¹. The intermolecular interaction energy was calculated based on the total energies of the whole system and the individual component energy in the system. E_bind between CL-20 and HMX was expressed as Eq. 1.

                                                   

Where E_total represents the total energy of the whole system, and E_CL-20 and E_HMX are the energies of CL-20 and HMX, respectively. Higher binding energy reveals stronger interaction between components.^{32, 33}

     Table 2 shows the binding energies (E_bind), the van der walls energies (E_vdW), and the electrostatic energies (E_elec) ofCL-20/HMX cocrystaland their mixtures. The binding between CL-20 and HMX of CL-20/HMX cocrystal is far stronger than CL-20/HMX mixture, which indicates that the CL-20/HMX cocrystal structure is more stable than that of their mixture. As shown in Table 2, with the COMPASS force field, binding energy E_bind is the summation of E_vdw and E_elec, and electrostatic portion is the major contributions to the binding energy. Therefore, the cocrystallization changed the interior constitution and structure of CL-20 crystal and resulted in the formation of strong intramolecular interaction of the two components,³⁴ which all increase the stability of cocrystal explosives, enhance the vibrationresistance to external forces, and improve the security of explosives.

Table 2E_bind, E_elec and E_vdW of different systems (kcal/mol)

System	Ebind	Eelec	EvdW
CL-20/HMX Cocrystal	861.57	810.28	13.08
CL-20/HMX Mixture	257.18	176.74	19.67

 

3.3. The mechanical property

     Material stress (σ) and strain tensors (ε) are applied to explosive materials to exhibit small deformations when subjected to external forces. Constant stress molecular dynamics simulations were then used to analyze the elastic moduli obtained from the stress-strain behavior of materials when subjected to an applied load. The dependence of the stress on the strain for elastic materials is shown in Eq. 2.

                                                                           （2）                                

C_ijkl (i, j, k, l=1, 2, 3) are the elements of the elastic symmetric matrix. To fully describe the stress-strain behavior of an arbitrary material, a maximum of up to 21 constants are required. The σ and ε are expressed by Eq. 3 and 4, where i, m_i, v_i and f_i are the number, mass, velocity of the particle and force acting on the particle, respectively. The parameter V₀ represents the volume of the system.

                                                                   

                                                                                         

       In Eq4, h₀ and h represent the matrix of the initial dimension of the structure and the dimension of the deformed structure, respectively. G is the metric tensor h^Th.

     The Reuss average is used to obtain the effective isotropic compliances of given materials averaged over all orientations. The effective bulk moduli (K) and shear moduli (G) are shown in Eq. 5 and 6, respectively, where the subscript R represents the Reuss average. The compliance matrix Sstands for the inverse matrix of the elastic symmetric matrix C, i.e., S=C^-1. The isotropic linear tensile moduli (E) and Poisson’s ratio (ν) can be obtained from Eq. 7.

                                                         

                               

       Tensile, bulk and shear moduli revealthe mechanical properties of pure CL-20, HMX, CL-20/HMX cocrystal and CL-20/HMX mixture, which were summarized in Table 3. E, K and G of CL-20/HMX cocrystal and mixture were lower than those of CL-20, which indicate that cocrystal and mixture decreased the stiffness and made CL-20 system more flexible and "softer". When the system were subjected to an external force, it can effectively buffer, disperse force and reduce the friction between the explosive particles, leading to the internal stress distribution more evenly. Thus, the cocrystal and mixture reduce the "hot spot" formation and sensitivity of CL-20 explosive, increasing the security of CL-20.³⁵

Table 3 The mechanical properties of different systems

Parameter	CL-20	HMX	Cocrystal	Mixture
Tensile modulus (E)	9.4	7.5	2.6	6.6
Poisson ratio (v)	0.3	0.3	0.4	0.3
Bulk modulus (K)	7.3	5.3	3.9	5.4
Shear modulus (G)	3.4	3.8	2.3	2.4
K/G	2.1	1.4	1.7	2.3

 

3.4. Thermal decomposition

      To study the effect of cocrystal on the thermaldecomposition  mechanism of CL-20, we employ ReaxFF/lg in RMD simulations to study the thermaldecomposition  mechanism in 2000 K. Fig. 3(a) exhibits the evolution of the potential energy of different systems with the change of simulation time. Clearly the potential energy of HMX system shows no significance descending tendency with the comparison of the others at initial simulation time (10 ps). For CL-20 and its mixture and cocrystal, with the increase of time, the total potential energy decreased sharply in the whole of simulation times because of the higher thermal sensitivity and thermal decomposition of CL-20. The potential energy curve of pure CL-20 crystal system shows a steeper tendency than those of the CL-20/HMX mixture and cocrystal. The reason is that some of CL-20 molecules in the mixture and cocrystal are isolated by relatively stable HMX molecules, which make the sensitive CL-20 molecules less likely to interact with other CL-20 molecules for the following chemical reactions and restrict the rapid decomposition. In contrast to CL-20/HMX mixture, CL-20/HMX cocrystal need a fusion process of cocrystal at the initial stage of the induction (endothermic) decomposition due to endothermic NO₂ dissociation, which result in the higher potential energy. After 20 ps, the CL-20/HMX cocrystal decomposed rapidly, and released plenty of energy.

Figure3. The potential energy (a) and total species (b) changes along with the time evolution. (The initial potential energy is set to be zero as a reference)

       Fig. 3 (b) shows the total number of species in the reaction process along with the simulate time evolution. HMX has a lower thermal sensitivity, thusthere is not any thermal decomposition reaction within the first 10 ps, resulting in the few number of species. While, CL-20 crystal has generated a large number of products within 5 ps, then fluctuated around 20 ps. The pure CL-20 and HMX system have less atom number (1152 atoms for CL-20and 512 atoms for HMX) in the initial state, so the mixture andCL-20/HMX cocrystalproduce more species than CL-20 and HMX system. Strong intermolecular interactions will promote the thermal decomposition reaction of the components in the system, whereas the specific mechanical chemical coupling mechanism needs to be further studied

3.5. Species analysis

     Fig. 4 demonstrates the detailed species analysis of the CL-20/HMX cocrystal, CL-20/HMX mixture, CL-20 and HMX at the temperature of 2000 K along with the time evolution. The formation of several key products and intermediates of CL-20 and HMX decomposition are observed, such as NO₂, N₂, NO, H₂O, HONO, HON and CO₂, which are in accordance with the experimental results.³⁶

 

Figure4. The distribution of main products in (a) cocrystal, (b) mixture of CL-20 and HMX, (c) pure CL-20 and (d) HMX along with the time evolution

       The primary decomposition products of NO₂, were formed faster in larger quantities. NO₂ was derived from the dissociation of N-NO₂ bond owing to the thermal decomposition of CL-20, in agreement with the characteristic initiation thermolysis reaction of cyclic nitramines.³²N-NO₂ is also the initiation bond of HMX thermolysis, which is a route to produce NO₂. Subsequently, much NO₂ would be consumed away by the secondary reaction, so the number of NO₂ decreased rapidly and NO, HONO, HNO₃ etc. would be generated.

      NO is also the main product of thermolysis, and the formation paths are shown as follows. N₂O₄ was created by two NO₂ molecular fragments, then it was decomposed into NO₃ and NO. N-N bond in CL-20 and HMX was dissociated to generate C-ONO group, which arranged and acquired a proton to form HONO, leading to the generation of NO. It agrees well with the results of Isayev.³⁷With similar to NO₂, NO fission shaped in the early stages of thermal decomposition, then generated N₂ due to the secondary reaction.

      The generation paths of H₂O are described as follows:at initial stage, NOOH forms through nitro-group; then, it captures a proton to form -NOOH-H; and 1 mol of H₂O is pulled off. Second, NO₂ is arranged in forms of ONO-. Then it seizes a proton to form HONO. HONO is thermally decomposed to HO- and promotes the generation of H₂O. Third, HNO₃ is thermally decomposed to H₂O.

      Few HONO is produced at the early stage in the CL-20 decomposition, which are successively dominated by NO₂fission and ring-breaking. No HONO is formed during the later stages of the reaction, probably due to the abstraction of a hydrogen atom of NO₂ from other free radical and intermediates. In the condition, the HONO molecule is not stable and rapidly dissociated into NO and OH radicals.

      Although the underpinning mechanism of decomposition reaction remains unclear, based on the analysis of the decomposition products and previously reported work,^37,38we can infer that the possible decomposition reactions may undergo three-steps, (1) unimolecular endothermic decomposition of the CL-20 or HMX molecular by breaking N–NO₂ trigger bond in CL-20 and N-NO₂ trigger bond in HMX molecules. (2)unimolecular decomposition of primary fragments into such intermediates as N₂O, HONO, HNO₂, etc. (3) gas-phase reaction to produce final stable products (N₂, NO, H₂O, etc.).

4. Conclusions

     Molecular dynamics simulation method was applied to study the structure and property of pure CL-20, HMX, CL-20/HMX mixture and cocrystal.The sensitivity, mechanical property and thermal decomposition process were analyzed. The mixture andCL-20/HMX cocrystalmade the maximum bond length of trigger bond N-NO₂ in CL-20 become shorter and N-NO₂ in HMX elongated. The mixture and CL-20/HMX cocrystaldecreased the sensitivity of CL-20.The binding between CL-20 and HMX of CL-20/HMX cocrystal was far stronger than that of CL-20/HMX mixture, and the CL-20/HMX cocrystal structure was more stable than their mixture. Cocrystal and mixture decreased the stiffness and made CL-20 system more flexible and "softer". The cocrystal and mixture reduce the "hot spot" formation and sensitivity of CL-20 explosive, increasing the security of CL-20. The several key products and intermediates of CL-20/HMX mixture and cocrystal decomposition are NO₂, N₂, NO, H₂O, HONO, HON and CO₂.

Conflict of interest

There are no conflicts to declare.

Acknowledgements

Project is supported by the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No.20131420120004). Science Foundation for Young Scientists of Shanxi Province（Grant No.201701D221103）.

References

1. 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.
2. Y. A. Bogdanova, S. A. Gubin, B. L. Korsunskii and V. I. Pepekin, Combust. Explo. Shock, 2009, 45, 738-743.
3. P. A. U. R. L. Simpson, D. L. Ornellas, G. L. Moody, K. J. Scribner, and D. M. Hoffman, Propellants Explos. Pyrotech., 1997, 249-255.
4. J. P. Agrawal, Energy Cornbust., 1998, 24, 1-30.
5. A. K. Sikder and N. Sikder, J. Hazard. Mater., 2004, 112, 1-15.
6. J. M. Lehn, Angew. Chem. Int. Ed. Engl., 1988, 27, 89-112.
7. F. Lara-Ochoa and G. Espinosa-Petrez, Supramol. Chem., 2007, 19, 553-557.
8. Y. P. Wang, Z. W. Yang, H. Z. Li, X. Q. Zhou, Q. Zhang, J. H. Wang and Y. C. Liu, Propellants Explos. Pyrotech., 2014, 39, 590-596.
9. O. Bolton and A. J. Matzger, Angew. Chem. Int. Ed. Engl., 2011, 50, 8960-8963.
10. Z. W. Yang, H. Z. Li, X. Q. Zhou, C. Y. Zhang, H. Huang, J. S. Li and F. D. Nie, Cryst. Growth Des., 2012, 12, 5155-5158.
11. H. F. Xu, X. H. Duan, H. Z. Li and C. H. Pei, Rsc Adv., 2015, 5, 95764-95770.
12. C. Y. Guo, H. B. Zhang, X. C. Wang, J. J. Xu, Y. Liu, X. F. Liu, H. Huang and J. Sun, J. Mol. Struct., 2013, 1048, 267-273.
13. O. Bolton, L. R. Simke, P. F. Pagoria and A. J. Matzger, Cryst. Growth Des., 2012, 12, 4311-4314.
14. N. Rom, S. V. Zybin, A. C. van Duin, W. A. Goddard, 3rd, Y. Zeiri, G. Katz and R. Kosloff, J. Phys. Chem. A, 2011, 115, 10181-10202.
15. S. P. Han, A. C. van Duin, W. A. Goddard, 3rd and A. Strachan, J. Phys. Chem. B, 2011, 115, 6534-6540.
16. A. Strachan, E. M. Kober, A. C. van Duin, J. Oxgaard and W. A. Goddard, J. Chem. Phys., 2005, 122, 54502.
17. M. R. Weismiller, M. F. Russo, A. C. T. van Duin and R. A. Yetter, Proc. Combust. Inst., 2013, 34, 3489-3497.
18. A. Strachan, A. C. van Duin, D. Chakraborty, S. Dasgupta and W. A. Goddard, 3rd, Phys. Rev. Lett., 2003, 91, 098301.
19. Z. Li, C. Lang, W. Chen and W. Jun-ying, Chin. J. Explos. Propellants, 2012, 4, 002.
20. H. Liu, Q. K. Li and Y. H. He, Acta Phys Sin.-Ch Ed, 2013, 62.
21. D. Z. Guo, Q. An, S. V. Zybin, W. A. Goddard, F. L. Huang and B. Tang, J. Mater. Chem. A, 2015, 3, 5409-5419.
22. H. Sun, J. Phys. Chem. B, 1998, 102, 7338-7364.
23. L. Liu, Y. Liu, S. V. Zybin, H. Sun and W. A. Goddard, 3rd, J. Phys. Chem. A, 2011, 115, 11016-11022.
24. D. I. A. Millar, H. E. Maynard-Casely, A. K. Kleppe, W. G. Marshall, C. R. Pulham and A. S. Cumming, CrystEngComm, 2010, 12, 2524-2527.
25. H. H. Cady, A. C. Larson and D. T. Cromer, Acta Crystallogr., 1963, 16, 617-623.
26. H. C. Andersen, J. Chem. Phys., 1980, 72, 2384-2393.
27. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684-3690.
28. P. P. Ewald, Ann. Phys., 1921, 369, 253-287.
29. M. P. Tosi, Solid State Phys., 1964, 16, 1-120.
30. W. Zhu, J. J. Xiao, J. Zheng, X. B. Zhao, Z. E. Chen and H. M. Xiao, Acta Chim. Sinica, 2008, 66, 2592-2596.
31. T. Su, R. Peng, Z. D. Hood, M. Naguib, I. N. Ivanov, J. K. Keum, Z. Qin, Z. Guo and Z. Wu, Chemsuschem, 2018, 11, 688-699.
32. X. Xu, J. Xiao, H. Huang, J. Li and H. Xiao, Sci. China Ser. B-Chem., 2007, 50, 737-745.
33. Z. Qin, T. Su, H. Ji, Y. Jiang, R. Liu and J. Chen, Aiche J., 2015, 61, 1613-1627.
34. T. Su, Q. Shao, Z. Qin, Z. Guo and Z. Wu, Acs Catal., 2018, 8, 2253-2276.
35. T. Sun, Q. Liu, J. J. Xiao, F. Zhao and H. M. Xiao, Acta Chim. Sinica, 2014, 72, 1036-1042.
36. J. C. Oxley, A. Kooh, R. Szekeres and W. Zheng, J. Phys. Chem., 1994, 98, 7004-7008.
37. O. Isayev, L. Gorb, M. Qasim and J. Leszczynski, J. Phys. Chem. B, 2008, 112, 11005-11013.
38. L. Z. Zhang, S. V. Zybin, A. C. T. van Duin, S. Dasgupta, W. A. Goddard and E. M. Kober, J. Phys. Chem. A, 2009, 113, 10619-10640.