*Article* **A Study of the Shock Sensitivity of Energetic Single Crystals by Large-Scale Ab Initio Molecular Dynamics Simulations**

**Lei Zhang 1,2,\*, Yi Yu <sup>1</sup> and Meizhen Xiang <sup>2</sup>**


Received: 6 August 2019; Accepted: 27 August 2019; Published: 3 September 2019

**Abstract:** Understanding the reaction initiation of energetic single crystals under external stimuli is a long-term challenge in the field of high energy density materials. Herewith, we developed an *ab initio* molecular dynamics method based on the multiscale shock technique (MSST) and reported the reaction initiation mechanism by performing large-scale simulations for the sensitive explosive benzotrifuroxan (BTF), insensitive explosive triaminotrinitrobenzene (TATB), four polymorphs of hexanitrohexaazaisowurtzitane (CL-20) pristine crystals and five novel CL-20 cocrystals. A theoretical indicator, *tinitiation*, the delay of decomposition reaction under shock, was proposed to characterize the shock sensitivity of energetic single crystal, which was proved to be reliable and satisfactorily consistent with experiments. We found that it was the coupling of heat and pressure that drove the shock reaction, wherein the vibrational spectra, the specific heat capacity, as well as the strength of the trigger bonds were the determinants of the shock sensitivity. The intermolecular hydrogen bonds were found to effectively buffer the system from heating, thereby delaying the decomposition reaction and reducing the shock sensitivity of the energetic single crystal. Theoretical rules for synthesizing novel energetic materials with low shock sensitivity were given. Our work is expected to provide a useful reference for the understanding, certifying and adjusting of the shock sensitivity of novel energetic materials.

**Keywords:** BTF; TATB; CL-20; cocrystal; energetic materials; shock sensitivity; large-scale ab initio molecular dynamics simulations

#### **1. Introduction**

Energetic materials (EMs) such as explosives, oxidizers and propellants are of significant importance in aerospace, oil-well drilling and other military and civilian applications. In this field, understanding the sensitivity of single crystals under shock or impact has long been a challenge. In engineering, shock sensitivity can be identified by the shock initiation threshold pressure, *P*90, which is obtained from the gap test and produces a detonation 50% of the time. *P*<sup>90</sup> results are generally reproducible and reliable. However, because of the high complexity of the test, the measurements have been performed for only limited types of energetic single crystals [1,2] and the test is not easily applicable to newly synthesized EMs. On the other hand, the drop-weight impact test, which characterizes the impact sensitivity of EMs by height *h*50%, is easy to implement. Therefore, there are abundant values of *h*50% in the literature compared to *P*<sup>90</sup> values. However, the *h*50% value is generally not reproducible as the results significantly vary depending on the conditions under which the tests are performed. For example, for benzotrifuroxan (BTF), the reported *h*50% values vary from 21 [3] to 50 [4] cm; for hexanitrohexaazaisowurtzitane/trinitrotoluene (CL-20/TNT) cocrystal, the values vary

from 30 [4] to 99 [5] cm. Therefore, the *h*50% values derived from the same experiments are comparable, while those from different equipment can only be used for a quantitative comparison of the mechanical sensitivity among various energetic single crystals.

With the rapid increase of computational capability and the development of material modeling methods, large-scale atomistic simulation becomes a powerful tool for understanding the physical processes of materials under extreme conditions [6,7]. However, in recent decades, there has been a strong tendency in the literature to elucidate the sensitivity of energetic single crystals by the electron density properties in a separate molecule of EM, such as electrostatic potential, molecular electronegativities, partial atomic charges, molecular weights, vibrational states, oxygen balance of the molecules, detonation gas concentrations and heats of detonation [1]. These quantities ignore the deformation of the chemical bonds, the motion of the molecules, the inner/inter molecular chemical reactions and the symmetrical structure of the crystals, and thereby cannot comprehensively characterize the reaction initiation of an energetic single crystal under shock.

To this end, we developed an *ab initio* molecular dynamics method [8] and extended its computational capability by improving the code's parallel calculation efficiency. We performed shock wave simulation tests on 11 types of energetic single crystals, with each simulated model composed of ~1000 atoms. On the basis of the calculations, we proposed a theoretical indicator to characterize the shock wave sensitivity of energetic single crystals, which is expected to be useful for the evaluation and adjustment of the shock sensitivity of novel EMs. We also revealed the shock reaction initiation mechanism, found factors that can inhibit the shock sensitivity of EMs and provided theoretical rules for synthesizing novel EMs with low shock sensitivity.

#### **2. Methodology**

#### *2.1. Multiscale Simulation Method of Shock Wave Tests*

Figure 1a schematically shows the simulated dynamical shock process in a single crystal. At the beginning when the shock wave reaches the single crystal, the simulation region starts to undergo lattice and molecular deformation. With the increase of simulation time, the simulation region gradually leaves the wave front relative to the unshocked material. The molecules in the simulated region are then compressed to react, eventually reaching the Chapman–Jouguet (CJ) state when the shock wave propagates steadily. We implemented the multiscale shock technique (MSST) [9,10] in the High Accuracy atomistic Simulation for Energetic Materials (HASEM) package [8], so as to capture the atomic motions and the chemical reactions of inner/inter molecules on the density functional theory (DFT) accuracy level.

#### 2.1.1. Continuum Theory Description of the Shock Wave Structure

The shock wave propagation was modelled using the one-dimensional Euler equations for compressible flow [9,10], which were represented by the conservation of mass, momentum, and energy, respectively, in the regions before and after the shock wave interface.

$$
\mu = v\_{\text{shock}} (1 - \frac{\rho\_0}{\rho}),
\tag{1}
$$

$$p - p\_0 = \upsilon\_{\text{shock}} 2\rho\_0 (1 - \frac{\rho\_0}{\rho}),\tag{2}$$

$$
\varepsilon - \varepsilon\_0 = p\_0(\frac{1}{\rho\_0} - \frac{1}{\rho}) + \frac{\upsilon\_{\text{shock}}^2}{2} (1 - \frac{\rho\_0}{\rho})^2,\tag{3}
$$

where *vshock* is the shock wave speed, *u* is the particle velocity, ρ is the material density, *e* is the energy per unit mass and *p* is the negative component of the stress tensor along the shock propagation direction. The quantities with 0 subscript describe the states of the unreacted material.

**Figure 1.** Dynamics simulation method of shock wave tests. (**a**) Schematics of the simulation method. (**b**) Validation of the code by comparing the Hugoniot curve of a single crystal model of octogen (HMX) single crystal to the reported experimental and calculational results in the literature. (**c**) Accuracy verification of the code by comparing the lattice lengths of the studied eleven energetic single crystals to the experimental values obtained by X-ray crystallography.

#### 2.1.2. Molecular Dynamics Description of the Atomic Motions

In the simulated region, the atomic motions were simulated using the molecular dynamics (MD) method. The Lagrangian per unit mass is

$$L = T\_c \left( \left\langle \stackrel{\cdot}{r}\_1 \right\rangle \right) - V\_c \left( \left\langle \stackrel{\cdot}{r}\_1 \right\rangle \right) + \frac{1}{2} Q \dot{\nu}^2 + \frac{1}{2} \frac{v\_{\text{shock}}}{\nu\_0^2} (\nu\_0 - \nu)^2 + p\_0 (\nu\_0 - \nu), \tag{4}$$

where υ = 1/ρ is the specific volume; *Te* and *Ve* are kinetic and potential energies, respectively, and their sum equals to *e*; *Q* is a parameter related to the mass of the simulated cell. When the volume of the simulated cell was fixed, that is, when . υ = 0, the Lagrangian expression is equivalent to the continuum Hugoniot relation of Equation (3).

The atomic position and velocity at each time step were obtained from the control equation of the simulated cell volume:

$$Q\ddot{\upsilon} = \frac{\partial T}{\partial \upsilon} - \frac{\partial V}{\partial \upsilon} - p\_0 - \frac{\upsilon\_{\text{shock}}}{\upsilon\_0^2} (\upsilon\_0 - \upsilon)\_{\text{\prime}} \tag{5}$$

which degenerates to the Rayleigh line of Equation (2) when the cell volume changes uniformly [9,10]. Thus, the system is restrained to fit the shock Hugoniot and the Rayleigh line of the material by changing the volume of the simulated cell.

#### 2.1.3. Density Functional Theory Description of the Electronic Structure

The atomic force of each time step was updated according to the DFT calculations of the electronic structure using HASEM software. The generalized gradient approximation was used for the exchange-correlation functional in the Perdew–Burke–Ernzerhof form. Norm-conserving pseudopotentials specialized for EM crystals were used to replace the core electrons. The valence electrons, described by linear combinations of numerical pseudoatomic orbitals, were calculated on a three-dimensional real-space grid. The reliability of the DFT calculations to describe the structures, energetics, dynamics, mechanical properties, detonation performance and sensitivity of EM crystals has been extensively confirmed in previous work [8,11–15].

In order to improve the computational capability of the dynamics simulation method, we reconstructed the HASEM software based on the J parallel adaptive structured mesh applications infrastructure (JASMIN), which has successfully accelerated many parallel programs for large scale simulations of complex applications on parallel computers [16]. Through this, the calculation efficiency of HASEM software was improved by one order of magnitude. Simulations of large-scale systems containing ~1000 atoms can thereby be achieved by using extended central processing units (CPUs) on supercomputers.

#### 2.1.4. Verification and Validation of the Dynamics Simulation Method

We constructed a single crystal model of octogen (HMX) and performed shock wave tests using the newly developed dynamics simulation method. A series of shock waves, with a speed smaller than 5 km/s, were applied on the HMX model. As shown in Figure 1b, the obtained Hugoniot curve satisfactorily agreed with the experiments [17–19] and other calculations [20], thereby confirming the reliability of the current method.

#### *2.2. Simulation Models of Eleven EM Single Crystals*

CL-20 has been proven to show excellent performance since it was first synthesized by the Naval Air Warfare Center China Lake 30 years ago [21]; however, it has not been widely used until now because of the sensitivity problems of its ε, β, γ and ζ polymorphs [21]. Cocrystallization, which mixes several components on a molecular level, has been considered a promising technique to obtain advanced EMs with good detonation performance and low sensitivity to accidental initiation. Under that circumstance, five novel cocrystals—CL-20/H2O, CL-20/TNT, CL-20/1,3-dinitrobenzene (CL-20/DNB), CL-20/N-methyl-2-pyrrolidone/H2O (CL-20/NMP/H2O) and CL-20/HMX, have been recently synthesized.

Herewith, we studied the four CL-20 polymorphs and the five novel CL-20 cocrystals, as well as the sensitive explosive BTF and the insensitive explosive triaminotrinitrobenzene (TATB). For the 11 EMs studied, we optimized the crystal geometries using the conjugate gradient method on a DFT level, with the initial inputs taken from the lattice parameters and atomic coordinates from single-crystal X-ray diffraction analysis [5,22–29]. The structures were considered as optimized when the stress components were less than 0.01 GPa and the residual forces were less than 0.03 eV/Å. As shown in Figure 1c, the calculated lattice lengths showed satisfactory agreement [30] with the experimental measurements (σ = 0.18 Å; *R*<sup>2</sup> = 0.9992), thereby confirming the reliability of the current method.

Subsequently, we built large-scale supercells of the eleven EM crystals for the shock simulations. As shown in Table 1, the supercells generally contained more than 1000 atoms, with the lattice length in the range of 15~30 Å. There were 24~64 molecules in each supercell, containing 72~192 trigger chemical bonds. The number of the chemical bonds included here was an order of magnitude larger than the traditional *ab initio* MD simulations, and thereby better reflects the randomness and probability characteristics of the chemical reaction kinetics.

**Table 1.** Lattice lengths (Å) of the supercells used for shock simulations. The type of the trigger bonds for each crystal, as well as the number of the composed atoms, molecules and trigger bonds in each supercell, are also given. EMs = energetic materials, BTF = benzotrifuroxan, TATB = insensitive explosive triaminotrinitrobenzene, CL-20 = hexanitrohexaazaisowurtzitane, TNT = trinitrotoluene, HMX = single crystal model of octogen.


#### *2.3. Control Parameter of the Shock Simulation Tests*

All the 11 systems were shocked with the same speed *Vshock* = 9 km/s at a direction perpendicular to the lattice vector. We used this shock wave speed based on a previous classical MD study, in which the breaking of the CL-20 trigger bonds was apparent at this shock speed [31]. The time step for the *ab initio* MD simulation was set to be 0.1 fs. As shown in Table 1, under this shock condition, 10 of the systems (excluding the insensitive explosive TATB) were initiated within 10,000 steps, that is, 1000 fs.

We defined the chemical bonds as broken when they were stretched to a cutoff percentage relative to each equilibrium state. The cutoff criterion was 20% on average, but it varied from 10% to 40% for different types of bonds. The criterion for each type of bond was determined by the statistics of the reaction products of the TATB explosive when the number of product molecules best agreed with the number of the stable clusters with a life span more than hundreds of time steps during the kinetics simulation.

#### **3. Results and Discussion**

#### *3.1. Shock Dynamics of the 11 EM Crystals*

The 11 crystals studied were rapidly compressed under shock, as shown in Figure 2. The temperature and pressure of the systems increased as a function of time during the shock process. The molecules were packed more densely in space and the chemical bonds plastically deformed to break, leading to the decomposition of material. According to our simulation, the N–NO2 bond, N–O bond and C–NO2 bond were the most active chemical bonds to deform and break in the CL-20 molecule, BTF molecule and TATB molecule, respectively, under shock. We therefore denoted these bonds as the "trigger bonds".

Take the shocked ε-CL-20 crystal as an example. The N–NO2 bonds with the exo-spatial orientation with respect to the five-membered imidazolidine ring were the trigger bonds. During the shock process, the molecular conformations at different times are shown in Figure 3a. Both the increase of the trigger bond length and the decrease of the trigger bond strength went in an exponential manner, as shown in Figure 3b,c. When the stretching of the chemical bond reached the cutoff percentage relative to the equilibrium state, we defined it as broken. Therefore, the first breaking of the trigger bond of the shocked ε-CL-20 crystal occurred at *t*<sup>3</sup> = 145.8 fs.

**Figure 2.** Shock dynamics of the (**a**) ε-CL-20 crystal, (**b**) CL-20/HMX cocrystal, (**c**) CL-20/TNT cocrystal, (**d**) CL-20/NMP/H2O cocrystal and (**e**) TATB crystal. In each panel, the temperature and pressure are shown as a function of time. The molecular conformations at the beginning of each decomposition reaction are also plotted.

**Figure 3.** (**a**) Conformation vs time plots for a molecule in the shocked ε-CL-20 crystal, along with the corresponding (**b**) length and (**c**) strength of the trigger bond N–NO2 in each snapshot of (**a**).

The chemical bonds of the eleven crystals studied included N–N, H–C, H–N, H–O, C–C, C–N, C–O, N–O and O–O bonds. For all types of chemical bonds in each crystal, we counted their number as a function of time during the shock process. As shown in Figure 4, the covalent bonds' breaking and recombination of shocked material was a dynamical process. For example, the trigger bonds N–NO2

in CL-20/NMP/H2O started to break at 124.8 fs, but they recombined to the initial state at 189.3 fs. We thereby defined the initiation of the chemical reaction by the time *tinitiation*, from when the breaking of the chemical bonds was always more than their recombination and the number of the trigger bonds decreased continuously. Therefore, the decomposition reaction began at *tinitiation* = 103.6 fs for BTF and at 204.3 fs for CL-20/NMP/H2O.

**Figure 4.** Bond number vs time plots for the shocked ε-CL-20 crystal, BTF crystal, CL-20/NMP/H2O cocrystal and TATB crystal. The chemical bonds counted included N–N, H–C, H–N, H–O, C–C, C–N, C–O, N–O and O–O bonds. The sign of the initiation of the reaction under shock was a continuous reduction in the number of trigger bonds.

Generally speaking, for the energetic crystals containing CL-20 molecules, both the temperature and the pressure increased slowly at the beginning of the shock process, as shown in Figure 2. Then, at ~100 fs, the temperature and the pressure started to drastically increase to higher than 1000 K and higher than 50 GPa. Next, at ~150 fs, the systems reached a gently varied stage, during which the trigger chemical bonds started to break and the decomposition reactions of material began. For the shocked insensitive explosive TATB, both temperature and pressure varied uniformly, and no chemical reactions occurred before 1000 fs, while for the shocked sensitive explosive BTF, the chemical reaction already started at 103.6 fs.

#### *3.2. Theoretical Indicator of Shock Sensitivity: tinitiation*

Because of the lack of experimental shock sensitivity for most of the EMs studied, we proposed using *tinitiation* as a theoretical indicator to characterize the ease of the shock reaction initiation. Because shock sensitivity has been proven to have a satisfactory correlation with the impact sensitivity [1,2], we used the experimental value *h*50% as a reference to compare with *tinitiation*, as shown in Table 2. We note comparing *h*50% values from the same experiments can well reflect the relative sensitivity of different compounds, while comparing those from different experiments was only qualitatively reasonable because of the influence of different experimental conditions used.


**Table 2.** The experimental *h*50% (cm) values and the initiation time of the shock reaction *tinitiation* (fs) of the 11 EMs studied. The strength of the trigger bond *Strigger* (kcal/mol), the temperature *Tinitiation* (K) and the temperature rising rate *TRR* (K/fs) for each crystal are also given for the study of the mechanism of the shock reaction initiation.

As shown in Table 2, there was a satisfactory agreement between the *tinitiation* and the *h*50% values derived from the same experiment. For example, in experiment 1, BTF had the highest sensitivity with *h*50% = 50 cm, and TATB had the lowest sensitivity with *h*50% > 320 cm. Correspondingly, BTF had the shortest delay of shock reaction at *tinitiation* = 103.6 fs, while TATB had the longest delay with *tinitiation* > 1000.0 fs. In experiment 2, the sensitivity order characterized by *h*50% was ε-CL-20 > CL-20/H2O > CL-20/NMP/H2O. Correspondingly, the sensitivity order quantified by *tinitiation* was exactly the same, which was 145.8 fs for ε-CL-20, 174.4 fs for CL-20/H2O and 204.3 fs for CL-20/NMP/H2O. In experiments 3 and 4, the *h*50% for ε-CL-20 was less than those for CL-20/HMX and CL-20/TNT. Consistent with this, *tinitiation* = 145.8 fs for ε-CL-20 was also smaller than *tinitiation* = 156.9 fs for CL-20/HMX and *tinitiation* = 181.8 fs for CL-20/HMX.

According to all the above comparisons between the calculated *tinitiation* and the measured *h*50% values, *tinitiation* is a reproduceable and reliable indicator to calibrate the shock sensitivity.

#### *3.3. Mechanism of Shock Reaction Initiation*

The shock can be simplified into a perfect impulse *f(t)*, which has an infinitely small duration. Its Fourier transform *F*(ω) = +∞ −∞ *f*(*t*)*e*−*j*ω*<sup>t</sup> dt* = *F*<sup>0</sup> implies that the shock causes a constant amplitude response in the entire frequency domain. Therefore, the more characteristic peaks in the vibrational spectra of an energetic crystal, the more modes can be excited and the more heat can be generated under the same shock condition. We plotted the vibrational spectra for the ε, β, γ and ζ polymorphs of CL-20 crystals in Figure 5. The number of characteristic peaks of the vibrational spectra of ζ-CL-20 was 25, which was the least among the four polymorphs. Correspondingly, the generated temperature *Tinitiation* = 1048 K was also the lowest. On the other hand, that peak number was the most for γ-CL-20 (29) and consistent with this, the temperature of *Tinitiation* = 1600 K was the highest. For the other two polymorphs, both the peak number and the temperature fall in the middle.

In order to study the mechanism of the initiation of shock reaction, we calculated for the 11 crystals the strength of the trigger bond *Strigger* and the temperature rising rate (TRR) under shock, as shown in Table 2. In this, the bond strength was quantified by the integrated value of the crystal orbital Hamilton population (COHP) at band energy, and the temperature rising rate was calculated by dividing the temperature (when t = *tinitiation*) by *tinitiation*.

ω

ω

**Figure 5.** Vibrational spectra for ε, β, γ and ζ polymorphs of CL-20 crystals. The characteristic peaks of the vibrational spectra are marked by vertical red lines and are indexed by blue texts. For each of the four polymorphs, the molecular conformation and the temperature (*Tinitiation*) the crystal starts to decay under shock are also shown.

As shown in Table 2, the trigger bond strength of the sensitive explosive BTF was *Strigger* = 42 kcal/mol, while that of the insensitive explosive TATB was three times higher. As well as this, the temperature rising rate of shocked BTF was *TRR* = 21.7 K/fs, while that of TATB was 29 times smaller. For the other EMs containing CL-20 molecules, both the trigger bond strength and the TRR fell in the range between BTF and TATB. In addition, we found that the *tinitiation–TRR* correlation showed a satisfactory power function *<sup>y</sup>* <sup>=</sup> <sup>804</sup> <sup>×</sup> *<sup>x</sup>*<sup>−</sup>0.76, with the coefficient of determination *<sup>R</sup>*<sup>2</sup> <sup>=</sup> 0.995, as shown in Figure 6. Therefore, the ease of the shock reaction initiation was apparently determined by the trigger bond strength and the temperature rising rate under shock. As a simplification, the specific heat capacity of a compound was the amount of heat needed per unit mass in order to raise the temperature by ΔT. Therefore, a compound with a larger specific heat capacity generally has a smaller TRR. We thereby propose that stronger covalent bonds and higher specific heat capacity are beneficial for delaying the time of shock initiation *tinitiation*, that is, reducing the shock sensitivity.

**Figure 6.** (**a**) Correlation between the shock sensitivity and the temperature rising rate under shock; (**b**) is an enlarged plot in a focused range, as marked by gray in plot (**a**).

Among the five CL-20 cocrystals, CL-20/NMP/H2O had the lowest TRR = 5.3 K/fs, while the other EMs had their TRR values in the range of 8.0–8.8 K/fs, as shown in Figure 6b. At the same time, the trigger bond strength of CL-20/NMP/H2O was also the highest, which was *Strigger* = 115 kcal/mol, as shown in Table 2. Therefore, CL-20/NMP/H2O was able to obtain the lowest shock sensitivity among the five cocrystals.

From the explanation above, the initiation of the shock reaction of energetic single crystals is shown to be a process driven by the coupling of heat and pressure. The heat is derived from the mechanical work of the shock compression and is transferred into the vibration of the lattice, the molecules and the chemical bonds of the shocked material. Denser characteristic peaks of vibrational spectrum correspond to a larger amount of heat generated by the shock. Driven by the heat, the temperature of the system quickly increases and the stretch vibrational modes of the chemical bonds are activated. While vibrating, the chemical bonds also endure plastic deformation under the shock compression. When the deformation of the trigger bond is beyond the critical level, the shock reaction begins.

#### *3.4. Shock Sensitivity Bu*ff*er: Intermolecular Hydrogen Bond*

On the basis of the calculated *tinitiation*, we were able to predict the relative sensitivity of all the 11 EMs studied, which was shown to be BTF > ζ-CL-20 > γ-CL-20 > β-CL-20 > ε-CL-20 > CL-20/HMX > CL-20/H2O > CL-20/DNB > CL-20/TNT > CL-20/NMP/H2O > TATB. The predicted order shows a close relationship between the shock sensitivity and the hydrogen bonding amount. For example, BTF contains no hydrogen and it owns the highest sensitivity, whereas TATB contains the most hydrogen and it has the lowest sensitivity.

In Figure 7 we show quantitatively the relationship between the shock sensitivity and the hydrogen bonding amount for the EMs containing CL-20 molecules, wherein the hydrogen bonding amount is represented by the occupied percentage in the Hirshfeld surface of the CL-20 molecules. The correlation is a satisfactory exponential function, in which *tinitiation* ∝ 1/ <sup>1</sup> + *<sup>e</sup>*−*k*(*x*−*x*0) , with the coefficient of determination *R*<sup>2</sup> = 0.9998. The correlation implies that the more hydrogen bonding occurs, the lower the shock sensitivity.

**Figure 7.** Shock sensitivity vs amount of hydrogen bonding for the CL-20 composed EMs studied.

The intermolecular hydrogen bond A:H–D (with ":" representing the electron lone pair, A for acceptor and D for donor) integrates the H–D polar-covalent bond, the A:H nonbond, and the A–D repulsive coupling interaction. Under shock, the hydrogen bonds show their elasticity—the covalent bond segment contracts and the nonbond elongates [33,34]. The special elasticity allows hydrogen bonds to vibrate in a continuous frequency region (<200 cm−1) so that the crystal can absorb more energy from the shock before reaching a temperature that is too high. This is analogous to the function of hydrogen bonds on improving the specific heat capacity of liquid H2O [35]. In order to confirm our hypothesis, we show the relationship between the TRR and the hydrogen bonding amount in Figure 7, which is roughly a power function. This relationship suggests that the hydrogen bonding has a buffering effect on the heating of the system under shock, thereby delaying the initiation time of the chemical reaction *tinitiation*. This is the fundamental reason why cocrystallization with low-sensitive EM components can effectively reduce the sensitivity of CL-20 crystals.

#### **4. Conclusions**

To conclude, we have developed an *ab initio* molecular dynamics method based on the multiscale shock technique and performed shock wave simulation tests for the sensitive explosive BTF, insensitive explosive TATB, four polymorphs of CL-20 crystals and five novel CL-20 cocrystals, with each model containing ~1000 atoms. The main conclusion includes:


Our work is expected to provide a theoretical reference for the understanding, certifying and adjusting of the mechanical sensitivity of the single crystals of novel energetic materials.

**Author Contributions:** L.Z. designed and wrote the original manuscript. Y.Y. and M.X. developed and verified the *ab initio* molecular dynamics code.

**Funding:** This research was funded the National Natural Science Foundation of China [grant number 11604017]. **Acknowledgments:** L.Z. thanks Hui Huang from the China Academy of Engineering Physics for fruitful discussion.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References and Note**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Nanomaterials* Editorial Office E-mail: nanomaterials@mdpi.com www.mdpi.com/journal/nanomaterials

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18