Abstract:To study the fundamental physical properties and intermolecular interaction of energetic materials under the loading temperature， the first-principle calculation was performed combined with zero-point energy and temperature effect corrections. The accuracy of lattice parameters at experiment temperature （173 K） can be significantly improved， and the deviations between the calculated lattice parameters and available experimental data are within 1%. The unit cell volume change with temperature is relatively reasonable compared with the experimental value at 0-500 K， and their deviation is mainly from the lack of interactions between phonons. Furthermore， the basic thermodynamic properties such as heat capacity， entropy and bulk modulus were predicted， and the results indicate that the lattice parameters and thermal expansion coefficient of FOX-7 have strong anisotropy in 0-500 K. Especially， the thermal expansion coefficient of interlayer direction is higher than that of inner layer direction， which is closely related to the molecular configuration and stacking. Importantly， when the temperature reaches 200 K， the shrinkage of thermal expansion coefficient of FOX-7 is related to the rotation of NO₂ group. The NO₂ group would regulate the intermolecular interaction by changing the dihedral angle with the molecular plane， thereby triggering potential phase transformation of FOX-7. In addition， the bulk modulus under the adiabatic conditions is consistent with the experimental values reasonably， and the evolution of adiabatic bulk modulus with temperature reflects the softening behavior of FOX-7 at the finite temperature. With the increase of temperature， the calculated heat capacity and entropy increase gradually， showing obvious numerical differences under the constant volume and pressure due to the anharmonic effect.

Abstract:The design efficiency of energetic compounds depends on many factors， such as the proportion of potential high performance samples in the screening space and the accurate prediction method of key properties. In this study， we proposed a scheme to improve the overall performance of virtual screening space by pre-screening molecular skeletons， and a method combining high-throughput computing and deep learning is applied to the design of energetic compounds. It was found that there is a moderate positive correlation between crystal density molecular skeleton density of energetic molecules， and the overall density of virtual screening space can be effectively improved by pre-screening high-density molecular skeletons. Based on the density data-set of energetic crystals collected from the crystallography database CCDC， a new density prediction model of energetic crystals was obtained via deep learning， with reliable accuracy and generalization. We took fused-ring energetic molecules as the research object， obtained high-density fused-ring skeletons through skeleton pre-screening， and then the virtual screening space composed of potential high-density molecules was constructed through fragment docking. The formation enthalpy， detonation performance and chemical stability were predicted by quantum chemical calculation and the equation of state of detonation products. Finally， 6 novel energetic molecules with energy level better than RDX and stability better than TNT were selected by performance ranking. This study shows that the overall performance of virtual screening space can be effectively improved by pre-screening molecular skeletons， and on this basis， high-throughput computing and deep learning can be used to achieve efficient design of energetic molecules.

Abstract:The preparations of high-energy pentazolate salts is a research hotspot in the field of energetic materials. The preparation of pentazolate anion is a key step in the preparations of high-energy pentazolate salts. However， as the important precursors for pentazolate anion， the stability of existing arylpentazoles is generally not high. In order to develop new precursors of pentazolate anion with better properties， 18 substituted derivatives of PyN₅ with the electron-withdrawing and electron-donating groups， i.e.， R-PyN₅ （R=－NO₂， －CN， －NF₂， －OH， －OMe， －N（Me）₂）， were designed and studied by using the density functional theory method. The bond dissociation energy（E_BD）， and the activation energy（E_a1） of the bridged C—N bond and the activation energy（E_a2） of the cracking of the N₅ ring were calculated， and the stability of the bridged C—N bond and pentazolate ring were discussed. E_a1 of all molecules is smaller than E_BD， indicating that the breakage of the bridged C—N bond is more likely to follow the path 2 rather than path 1. E_a2 of all molecules is smaller than E_a1， indicating that the stability of the N₅ ring is the key factor to determine the stability of the arylpentazoles. Compared with PhN₅， —N（Me）₂ meta-substituted and bis-substituted compounds have lower E_a1， higher E_a2 and lower ΔE （E_a2－E_a1）. Therefore，—N（Me）₂ meta-substituted and bis-substituted compounds are the most potential precursors of N₅ˉ ion for replacing PhN₅.

Abstract:The performance prediction methods for high-throughput energetic molecule design and screening are required to balance accuracy and efficiency. In the present work， the suitability of three common theoretical methods on different levels， including semi-empirical PM6 method， density functional theory method B3LYP/6-31G（d，p）， and high-precision complete basis set CBS-4 method， for heat of formation （HOF） prediction under the atomization scheme for high-throughput energetic molecule design and screening were evaluated. The solid HOF of twenty energetic molecules were compared， and the results of different theoretical levels are found to differ greatly. Based on the predicted HOF， experimental density， and three models （K-J， BKW and VLW）， the detonation performance of ten common energetic molecules were predicted. The results show that B3LYP method possesses the best suitability and efficiency， and the predicted detonation performance is closed to that obtained by CBS method. For example， the average relative deviation of the detonation velocity and detonation pressure predicted by BKW are only 0.4% and 1.2%， respectively. However， both the low-precision PM6 method and the time-consuming CBS method are difficult to balance the requirements of precision and efficiency in high-throughput energetic molecule design and screening. It suggests that， for the high-throughput design of energetic molecules， a medium-precision method is adequate for rapid HOF prediction.

Abstract:ReaxFF molecular dynamics （ReaxFF MD） simulations were adopted to identify the main intermediate products， final products and chemical reactions during 3，4-dinitro-1H-pyrazole （DNP） thermal decomposition. Accelerating rate calorimeter （ARC）-mass spectrometer （MS） technique was adopted to study DNP thermal decomposition properties and identify the gaseous products. The simulated results illustrate that C₃HO₄N₄， C₃HO₃N₄， C₃HO₂N₃， C₃HNO₂， NO₂ are the main intermediate products， and H₂O， CO₂， N₂ are the main final products. MS detected main gaseous products are H₂O， CO₂， N₂ as well. According to the simulation results， the produced time and abundance of the products are obtained as well. Among which C₃HO₃N₄ is the first generated intermediate product， and H₂O is the first generated final product. C₃HO₃N₄ and N₂ are the intermediate and final products with the largest amount， respectively. Additionally， the main chemical reactions in DNP thermal decomposition process are also acquired by molecular dynamics simulations. According to the generation time and abundance of products， the decomposition path of DNP was obtained.

Abstract:Nano-thermites， as one kind of energetic composites， have wide applications. A systematic study on the relationship between the interface structures and properties has great significance for the preparation of the new nano-thermites with excellent performance. The structures and energies of Fe₂O₃（104） and Fe₂O₃（110） surfaces and the structures， bonding properties， and adhesion work of Al（111）/Fe₂O₃（104） and Al（111）/Fe₂O₃（110） interfaces （AFS1， AFS2， AFS3， AFS4 and AFS5） were studied with the periodic density functional theory in this work. Results show that O-terminated Fe₂O₃（104） and Fe₂O₃（110） surfaces and the interfaces formed by these surfaces with Al（111） are more stable than those of the （104） and （110） surfaces of Fe₂O₃ respectively. Among 5 of the Al/Fe₂O₃ interfaces， the interfaces composed by the O-terminated Fe₂O₃（104） and Fe₂O₃（110） surfaces with Al（111）， i.e.， AFS1 and AFS5， have the maximum adhesion work （3.92 J·m^-2 and 3.02 J·m^-2， respectively）， and AFS1 is more stable than AFS5. In these two most stable interfaces， the Al atoms stack on the top position of the O atoms of the Fe₂O₃ surfaces and the binding of Al and Fe₂O₃ surfaces is mainly through the Al-O ionic bonds.

Abstract:In order to realize the non-destructive measurement when predicting the storage life of nitrate ester plasticized polyether （NEPE） propellant， the high temperature accelerated aging， gas content monitoring and uniaxial tensile mechanical property experiment were carried out on NEPE propellant with 10% constant compression strain . The non-destructive storage life prediction model based on characteristic gas contents was proposed through correlation analysis and remaining life prediction model. The results show that during the storage and aging processes， the total amount of CO gas is the largest， reaching more than 1300 mg at different temperatures. The generating rates of NO and CO are growing slowly in the early aging period， and growing faster in the late period. The generating rate of HCl increases rapidly during the eraly and late aging period and slowly in the middle. Maximum tensile strength σ_m and maximum elongation ε_m increase slightly in the early aging period， the former oscillates slightly and the latter gradually increases in the middle period， and both of them decrease sharply in the late period. The correlation between the contents of CO and the maximum tensile strength is largest and there is a single correlation between them at different temperatures. The maximum correlation value reaches about 0.93-0.95. Four life prediction methods of NEPE propellant are established based on traditional and improved aging life prediction model， tensile strength and CO content. The maximum correlation coefficient and estimation results show that the improved prediction model based on the content of CO gas release is most effective.

Abstract:In order to determine the parameters of the JWL equation of state for unreacted explosives， a method to determine the JWL parameters by using the BP neural network-Genetic Algorithm （BP-GA Algorithm） and the shock Hugoniot was proposed. Firstly， BP neural network is trained to fit the nonlinear system composed of different JWL parameters， and then Genetic Algorithm is used to search the set of JWL parameters with the largest fitness value. The results show that the JWL parameters can be determined by the BP-GA Algorithm when the initial density， detonation velocity， Hugoniot parameters C₀ and S are known. The p-v curves of eight kinds of unreacted explosives determined by BP-GA Algorithm are consistent with those determined by test data， and the R² of eight p-v curves are not less than 0.9995， which proves the high accuracy of BP-GA algorithm.

Abstract:In order to investigate the thermal decomposition mechanism of the energetic cocrystal TKX-55 and the effect of solvent component dioxane （1，4-dioxane， DIO） on the decay of the energetic component 5，5´-bis（2，4，6-trinitrophenyl）-2，2´-bis（1，3，4-oxadiazole） （BTNPBO）， the molecular dynamics simulations on TKX-55 and pure solvent component DIO were carried out with the ReaxFF-lg （Reactive Force Field-Low Gradients） force field. The results show that the initial decomposition reaction of TKX-55 includes the dimerization of energetic molecules， the hydrogen transfer between energetic and solvent components， the ring-opening reaction of 1，3，4-oxadiazole in energetic components， and the dissociation of nitro group. The dimerization reaction facilitates the rapid growth of the subsequent clusters， and the release of the heat and the stable small molecule products are restricted by the formation of a large number of clusters. It is one essential reason for the high heat resistance of TKX-55. For the pure solvent， the heat release and clustering are constrained at low temperatures； while enhanced at elevated temperatures. The main role of DIO molecules in TKX-55 is thought-to adsorb small reactive intermediates （such as OH， NO， NO₂， etc.） and thereby inhibit the decomposition of BTNPBO.

Abstract:Pentazoles are currently a research hotspot in the field of energetic materials， however， the stability of existing pentazoles is generally not high. In order to develop new pentazoles with better properties， based on the analyses of the available structures， 20 substituted derivatives of HN₅ with the electron-withdrawing groups， i.e.，N₅（CH₂）_x-1R（R=—NO₂，—CF₃，—CN， —CHO， —COOH； x=1， 2， 3， 4）， were designed and studied by using the density functional theory method. The bond dissociation energy （E_BD） of the bonds linked with the N₅ ring and the activation energy （E_a） of the cracking of the N₅ ring were calculated and compared with that of some pentazoles substituted by the electron-donating groups， and the effects of substituents on E_BD and E_a were discussed. Results show that the E_a of all molecules is much smaller than E_BD， indicating that the stability of the N₅ ring is the key factor to determine the stability of the pentazoles. When R is directly connected to the N₅ ring （N₅R）， the E_a of N₅R with R being an electron-withdrawing group is smaller than that of N₅R with R being an H or an electron-donating group. The N₅ ring is a strong electron-withdrawing group， and bearing too much or too little negative charges is not conducive to the stability of the N₅ ring.

Abstract:The attachment energy （AE） model and molecular dynamics （MD） methods were used to predict the crystal morphology of 1，1-diamino-2，2-dinitroethylene （FOX-7） under vacuum condition and in eight solvent systems including dimethylsulfoxide （DMSO）， acetone， methanol， N-methylpyrrolidone （NMP）， N，N-Dimethylacetamide （DMAC）， ethylacetate （EA）， water （H₂O）， and DMSO/H₂O （V/V=2/1）. By calculating the interaction energies between the solvent and crystal plane， and the attachment energies under the influence of the solvent， the simulated crystal habit and its aspect ratio were obtained. The results show that FOX-7 crystal has six important growth planes under vacuum： （1 0 ）， （1 0 1）， （0 1 1）， （0 0 2）， （1 1 0）， （1 1 ）. Among them， the area of （0 1 1） plane accounts for the largest proportion， which is the most important crystal plane affecting the crystal morphology of FOX-7. The influence degree of solvent on the aspect ratio of crystal is in the following order： DMSO2O2O. By recrystallization experiments， FOX-7 crystals have a bulk-like shape in DMSO， methanol， and DMSO/H₂O； a rod-like shape in acetone and NMP； a needle-like shape in DMAC and H₂O； a flake-like shape in EA. The theoretical prediction results are in good agreement with the experimental results， which proves that the simulation of the crystal habit of FOX-7 based on the AE model can provide better guidance for the crystallization experiment. Results in thermal properties show that the crystal surface morphology and internal defects affect the phase transition temperature and thermal decomposition temperature of FOX-7. The fewer the crystal defects， the higher the α→β transformation temperature. The larger the crystal aspect ratio， the smaller the particle size， and the lower the first decomposition temperature.

Abstract:In order to analyze the factors affecting the ignition delay time in the gas phase reaction process between hydrazine fuel and NO₂， the reaction process was simulated by using density functional theory. The reaction activities， active sites， potential energy surfaces and reaction rate constants in hydrogen extraction reaction process of hydrazine （N₂H₄）， methyl hydrazine （MMH） and unsymmetrical dimethyl hydrazine （UDMH） were calculated. The results showed that the energy difference between the highest occupied orbital and the lowest vacant orbital of UDMH was the smallest among the three hydrazine fuels， which was 0.20522 eV， indicating that UDMH had the highest activity， so it has the fastest reaction rate with NO₂， which in line with the characteristic of the shortest ignition delay time. The active sites of three hydrazine fuels were identified， N（1） or N（4） for N₂H₄， N（1） for MMH and N（1） for UDMH. The active sites of hydrogen extraction reaction of three hydrazine fuels were calculated， it was found that the reaction barrier of UDMH is the smallest， which is 3.589 kJ ·mol^-1， and the reaction rate constant is the largest， which is 9.81×10⁵ L·s^-1·mol^-1， which is consistent with the shortest ignition delay time， it is concluded that in hydrazine fuel， the smaller the hydrogen extraction reaction barrier with NO₂， the larger the reaction rate constant， and the shorter the ignition delay time.

Abstract:In order to obtain high energy density liquid propellant fuels and increase the payload of the launch vehicle， 20 different methyl-substituted bicyclobutyl derivatives were designed， and the influence of the structure of bicyclobutyl derivatives on performance was studied through theoretical calculations. Results show that with the increasing number of methyl substituents， the heat of formation and specific impulse of bicyclobutyl derivatives show a decreasing trend. When the substituent is para-substituted， its molecular stability is the best， and the heat of formation and specific impulse are larger， while the ortho-position substitution has a weakening effect on the heat of formation and specific impulse of bicyclobutyl derivatives. Among the designed compounds， the specific impulse of bicyclobutyl is the highest. When the mixing ratio of bicyclobutyl and liquid oxygen is 28.5∶71.5， the specific impulse can reach 304.52 s， and the main combustion products are CO（34.64%）， CO₂（13.89%） and H₂O （29.54%）. The comprehensive performance of all designed products is better than that of rocket kerosene. This study provides theoretical support for the design and synthesis of high-energy fuels.

Abstract:To study the influence of crystal defects on the initial reaction of hexanitrohexaazaisowurtzitane （CL-20）， molecular dynamic simulation and ReaxFF-lg reactive force field are used to study the initial reaction path， thermal decomposition products， and reaction kinetics of CL-20 with vacancy defects at high temperature （1500-3500 K）. The results show that the initial decomposition path of CL-20 with vacancy is the breaking of N─NO₂ bond， the same as that of perfect crystal. The vacancy defects prove to increase the frequency of ring-opening reactions and the production of NO₂. Compared with perfect CL-20， it can be seen that the vacancy defects would reduce the CL-20 activation energy barrier and accelerate its thermal decomposition process. The reaction rate constants of CL-20 with 16.7% vacancies are 1.7 and 1.4 times higher than that of perfect CL-20 at 2000 K and 3000 K， respectively. The CL-20 molecules around the vacancy are easier to decompose， leading to the increase of the sensitivity of CL-20.

Abstract:The study of polynitrogen pentazolate salts， which are usually achieved from the precursors of N₅ covalent compounds， is the hotspot of research in the field of new energetic materials. In most cases， the stability of N₅ covalent compounds will significantly affect the possibility of successful preparations of pentazolate salts. Herein， the calculations of dissociation energy （E_BD） of the bonds in the straight side chain and activation energy （E_a） of the N₅ ring were carried out for the selected eighteen non-aryl substituted N₅ compounds （R—N₅ or N₅—R—N₅） by using the B3LYP/6-31G** method of density functional theory， and meanwhile the influence of the side chain on their stabilities and pyrolysis mechanism were investigated. When R is the hydroxyl or amino group， the side chain and the N₅ ring are more prone to break， making it difficult to obtain the N₅^- ring. When R is alkyl， the E_a of the N₅ ring cleavage is relatively larger， making it more likely to produce the N₅^- ring， and the stability of the side chain′s C—N bond as well as the N₅ ring will be little affected by the length of the alkyl chain. The sequential cleavage of two N₅ rings occurs in the bicyclic molecular structures and the energy barrier of the second ring is higher than that of the first one， resulting in the formation of N₂ and azide. The C—C bond on the side chain of the molecule will be broken before the break of C—N bond， which may significantly reduce the E_BD of the C—N bond but have little effect on the E_a of the N₅ ring. Therefore， for the preparations of pentazolate salts from covalent pentazoles， cutting off the C—C bond first may be more conducive to obtaining N₅^- ring.

Abstract:In order to analyze high temperature phase transformation of hexanitrohexaazaisowurtzitane （CL-20）， phase transformation temperatures and coefficients of thermal expansion of ε-， β-， and γ-CL-20 were studied via ReaxFF-lg reactive force field molecular dynamics， with modified valence potential intercept. To validate the applicability of selected force field，the density， cell constant， lattice energy， and sublimation enthalpy for three types of CL-20 at room temperature were calculated. The third order Birch-Murnaghan equation of state was used to fit the p-V curve of ε-CL-20， with pressure ranging from 0 to 280 GPa. And the variation of bulk modulus （B₀） and its partial derivative to pressure （B′₀） with the increase of pressure is analyzed. High temperature phase transformation analysis shows that ε-and γ-CL-20 change phases at 398-423 K， of which the ε→γ phase transition occurs at atmospheric pressure， while the γ→ε phase transition needs 0.5 GPa or higher pressure； β-CL-20 transforms to ε crystal form at 448 K. The thermal expansion coefficient analysis shows that there is no obvious anisotropy in the high temperature thermal expansion process of ε-CL-20， while β- and γ-CL-20 show anisotropy in c direction and b direction， respectively. Results show that the modified ReaxFF-lg reactive force field is suitable for the study of phase transition of ε-， β-， and γ-CL-20 at high temperature and high pressure， while the accuracy of thermal expansion of β- and γ-CL-20 needs to be further improved.

Abstract:Molecular dynamics simulation is an important method to predict the shock sensitivity of energetic materials， yet it is computationally expensive and needs to use force fields that may be unavailable. Here， an algorithm was designed and implemented in a computer program in Python for calculating the Steric Hindrance Index （SHI）， which is a descriptor for evaluating shock sensitivity. The algorithm 1） compresses the crystal unit cell of an energetic material keeping the molecular unit rigid to simulate deformation under shock； 2） establishes a new rectangular coordinate system for the specific slip system and rotates the cell to deal with general shock directions and slip systems； 3） assigns molecular units to layers based on the coordinate of their centroid； 4） calculates the overlapped area of each two adjacent layers after projection along the slip direction； and 5） obtains SHI by normalization of overlapped areas. For PETN， BTF， RDX， and TNT at a compression ratio of 0.1， the calculated average SHI are 0.8707， 0.7940， 0.4228， and 0.0924， respectively， which is consistent with the decreasing order of impact sensitivity mentioned in references. SHI classifies the slip systems in line with those based on molecular dynamics simulations， yet with better computing efficiency and methodological applicability.

Abstract:In order to speed up the development of new energetic materials and reduce the time and resource consumption caused by a large number of experiments， a method for predicting enthalpy of formation of energetic materials is proposed based on the theory of material genetic engineering. Firstly， the collected atomic coordinate data representing the molecular structure of energetic materials were converted into a coulomb matrix representing the cartesian coordinate system in the molecule to eliminate the influence of translation， rotation， index order and other operations on the prediction of enthalpy of formation. Then， the enthalpy of formation of energetic materials was predicted according to the proposed fusion model of Convolutional Neural Network （CNN） and Bi-directional Long Short-term Memory Network （Bi-LSTM） based on Attention mechanism. In this way， not only can the characteristics of the data be extracted effectively， but also the correlation between the data and the lack of long-term dependence can be fully considered. Meanwhile， the influence of important characteristics on the prediction results can be highlighted. The comparison of experimental results shows that the proposed method based on deep learning has the lowest experimental error in the prediction of enthalpy of formation. Its Mean Absolute Error （MAE）， Mean Absolute Percentage Error （MAPE）， Root Mean Square Error （RMSE） and Root Mean Squared Logarithmic Error （RMSLE） are 0.0374， 1.32%， 0.0541 and 0.028， respectively. The prediction goal of "structure-performance" is realized， and a new method is provided for the prediction of enthalpy of formation of energetic materials.