1 Introduction

Energetic materials are essential for both military and civil fields due to their wide applications as armaments, mining, space exploration and fireworks^[1-3]. Sensitivity to detonation (by friction, impact or shock et al.), detonation velocity, thermal stability and crystal morphology are the stringent performance criteria to energetic formulations, so current energetic materials are based on a very small number of compounds. The design and exploitation of advanced energetic materials is an interesting and challenging problem. Theoretical methods play significant parts in the development of advanced energetic materials. Densities and heats of formation (HOFs) are well known to evaluate the explosive performances of energetic materials. Computational approaches have shown their great advantages and been employed to obtain densities and HOFs of demanding materials^[4-8]. Thermal decomposition mechanisms are very important for evaluating the safety and storage reliability of energetic materials. Ab initio MD has obvious advantages in studying thermal decomposition mechanisms of novel materials^[9-16]. Recently, we have studied the thermal decomposition mechanisms of high nitrogen content energetic materials such as RDX^[17-19], tetrazine^[20-24], Fox-7^[25-27] by combining ab initio molecular dynamics (AIMD) method with density functional theory (DFT). It requires no prior experimental knowledge or intuitive assumptions about the decomposition.

The detonation velocity (v_D) and detonation pressure (p_D) are the necessary factors for evaluating the detonation properties of energetic compounds. We have used the VLW equation of state (VLW EOS) to successfully obtain the value of v_D and p_D for energetic compounds, which contain elements such as C, H, N, O, Al, Cu, Zn, by modified VLW code. What is more, properties of some systems are governed to a large extent by the interfaces between unmixed phases, and therefore the interfacial control is very important to tailor basic properties of such systems^[28-32]. Dissipative particle dynamics are capable of providing valuable microscopic and mesoscopic insights into the interfacial behaviors of the immiscible polymer blends, such as the interfacial width and tension^[33-39].

In this paper, the methodologies utilized to deal with theoretical designing of energetic materials in our work were introduced in detail, including evaluating the densities, HOFs, stability and detonation properties of usual energetic compounds. The achievement methods of meso-scale parameters for energetic polymers are also discussed.

2 Densities calculations

Monte-Carlo method is always a good method to calculate the densities of small molecules^[40-42], but this method does not carry conviction at all times, for example energetic metal complexes. Recently, we have employed a new, efficient and convenient method, based on the DFT computations implemented on Material Studio of SGI workstations in the China Academy of Engineering Physics, to predict the densities of a novel environmentally friendly octahedrally coordinated 2D polymeric complexes bis(1, 5-diaminotetrazole)-dichlorozinc(Ⅱ) (Zn(DAT)₂Cl₂)^[43], which was performed by the DMol³ program. The exchange-correlation interaction was treated by functional Perdew, Burke and Ernzerh of generalized gradient approximation (PBE GGA)^[44-49], and applied basis set was double numerical basis set plus d-functions (DND). The density was obtained from the enclosed volume of electron cloud around the molecule divided by the molecular weight (Eq. 1). The enclosed volume (V_en) of each molecule was yielded from on each optimized structure. Enclosed volume shows the volume encloses by the isosurface, within a single repeat unit. The volume is calculated by counting the number of grid points whose value is above (or below) the isovalue, That is to say the enclosed volume is visually that which is on the gray side of the isosurface. Then ρ_en was obtained by the following equation^[50]:

$ {\rho _{{\rm{en}}}}{\rm{ = }}\frac{{{M_{\rm{w}}}}}{{{V_{{\rm{en}}}} \cdot 6.02 \cdot {{10}^{23}} \cdot {{({{10}^{ - 8}})}^3}}} = \frac{{{M_{\rm{w}}}}}{{{V_{{\rm{en}}}} \cdot 0.602{\rm{ }}}}{\rm{g}} \cdot {\rm{c}}{{\rm{m}}^{ - 3}} $

(1)

Where M_w is the molecule weight, g·mol^-1. V_en stands for the enclosed volume of electron cloud around the molecule, ų.

In our work^[43], the isovalue of Cu(DAT)₂Cl₂ is 0.0166. The crystal structures and the shap of Ven of Zn(DAT)₂Cl₂ and Cu(DAT)₂Cl₂ are almost the same, so the isovalue of Zn(DAT)₂Cl₂ is 0.0166, too (Fig. 1). Hence, we got the density of Zn(DAT)₂Cl₂ easily by equation (1), 2.117 g·cm^-3. It is concluded that on the basis of the similar structure and the same isovalue, a new analogues of an exsit compound could be designed and developed. It should be noted that the experimental density of the basic compound must be known, unless we do not need the exact density value.

3 Heats of formation (HOFs)

HOF is one of the most important quantities used to assess the energetic properties of high energy high density materials, because the heat release upon decomposition or combustion is a radical factor to determine detonation or propellant performance^[51-52]. However, it is often difficult to measure HOF via experiments due to the danger and difficulty. Quantum chemical calculation has been widely used to solve this problem. The theoretically predicted values of HOFs can be in good agreement with experiments.

The method of atomization scheme^[53-54] (scheme 1) has been employed very successfully to calculate HOFs of energetic compounds. We have efficiently obtained the HOFs of high-nitrogen energetic substituted s-tetrazine compounds^[40], Cu(DAT)₂Cl₂, Zn(DAT)₂Cl₂ et al.^[43], by simplifying the HOFs equations of molecule M at 298 K (Eq.2). One strength of Eq.2 with respect to traditional atomization scheme methods is that ε₀(M) is the only one parameter that needs to be calculated, and it can be performed by soft packages such as Gaussian.

In the present equation 2, both the HOF of atom X and the corrected enthalpy of molecule M from 0 K (the standard state) to 298 K (gaseous state) needed not to be calculated.

HOF of M at 298 K can be written as Equations:

$ \begin{array}{l} \Delta H_{\rm{f}}^0\left( {{\rm{M}}, 298\;{\rm{K}}} \right) = \Delta H_{\rm{f}}^0\left( {{\rm{M}}, 0\;{\rm{K}}} \right) + \sum\limits_{{\rm{atom}}} {x(H_{\rm{x}}^0\left( {0\;{\rm{ K}}} \right)}-H_{\rm{x}}^0\left( {298\;{\rm{K}}} \right)) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }}(H_{\rm{M}}^0\left( {298\;{\rm{ K}}} \right)-H_{\rm{x}}^0\left( {0\;{\rm{K}}} \right))\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }} = \Delta H_{\rm{f}}^0\left( {{\rm{M}}, 0\;{\rm{K}}} \right) + \sum\limits_{{\rm{atom}}} {x(H_{\rm{x}}^0\left( {0\;{\rm{K}}} \right)}-H_{\rm{x}}^0\left( {298\;{\rm{K}}} \right)) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }}({H_{{\rm{corr}}}} - {\varepsilon _{{\rm{ZPE}}}}\left( {\rm{M}} \right))\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{{\rm{atom}}} x \Delta H_{\rm{f}}^0\left( {{\rm{X}}, 0\;{\rm{ K}}} \right) - \sum\limits_{{\rm{atom}}} {x{\varepsilon _0}} \left( {\rm{X}} \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }}\sum\limits_{{\rm{atom}}} {x(H_{\rm{x}}^0} \left( {298\;{\rm{K}}} \right) - H_{\rm{x}}^0\left( {0\;{\rm{ K}}} \right)) + ({\varepsilon _0}\left( {\rm{M}} \right) + {H_{{\rm{corr}}}}) \end{array} $

(2)

Where, ε₀(M) and ε₀(X) are the total energies of molecule M and each element that makes up M at 0 K, and x stands for the number of atoms of X in M. ΔH_f⁰(M, 298 K) is the enthalpies of formation of the radicals. ΔH_f⁰(X, 0 K) stands for the HOFs of atoms X at 0K which can be found from ref.^[53]. And H_corr is thermal correction to enthalpy from 0 to 298 K.Finally, the HOFs values can be calculated using our own computer code.

4 Thermal decomposition mechanisms

The thermal decomposition mechanisms of high nitrogen energetic compounds are of great interest. Recently, we have studied the thermal decomposition mechanisms of RDX^[17-19], tetrazine^[20-22] and Fox-7^[23-25]. The reaction channels were simulated theoretically with a molecule at high temperature in a number of trajectories, using the DFT method in combination with AIMD method with a plane wave basis set and pseudo potentials. This method requires no prior experimental knowledge or intuitive assumptions about the decomposition. Moreover, various decomposition channels of the compounds studied can be decided according to their occurrence frequency^[20-24].

The Vienna Ab Initio Simulation Package (VASP) program^[56-59] is based on AIMD method. It can be used to study the thermal decomposition trajectories of energetic materials, since it utilizes a traditional self-consistency scheme to evaluate the instantaneous electronic ground state at each MD step by conjugate gradient minimization of the total electronic energy. The detail calculation steps are described as follows: Firstly, the molecule studied is put in a cubic box (with a length of such as 8Å) to imitate gas phase condition. Secondly, the Perdew-Wang^[58] gradient correction will be added to the exchange-correlation functional^[59]. A plane-wave basis set with the cut-off energy of such as 348.1 eV should be used for the electron wave function. At last, the optimized Vanderbilt ultra-soft pseudo potentials^[60] supplied with the VASP package^{[61, 62]} are used for C, N and H atoms for the core region.

Here, we would like to review some examples where the applicability of VASP package is illustrated for a range of systems.

By VASP package, the thermal decomposition trajectories of five simple hydronitrogen compounds (N₂H₂, N₂H₄, N₃H, N₄H₂ and N₄H₄) and six tetrazine compouds (s-tetrazine, 3, 6-diamino-1, 2, 4, 5-tetrazine(DAT), 3, 6-dihydrazino-1, 2, 4, 5-tetrazine (DHT), 3, 6-diazino-1, 2, 4, 5-tetrazine (DiAT), 3, 6-bis(1H-tetrazol-5-amino)-1, 2, 4, 5-tetrazine (BTATz), and 3, 3′-Azobis(6-amino-1, 2, 4, 5-tetrazine) (DAAT)) were studied. The reaction channels were studied by Gaussian03^[63]at B3LYP/6-311G(d, p) level to locate the local minimum points and the transition structures. Vibrational frequencies were calculated at the same level to take account of the zero point energy and to identify the transition structures. To obtain more accurate and reliable reaction information, the high accuracy single point calculations were further performed at CCSD(T)=full/6-311+G(3df, 2p) level for these simple hydronitrogen compounds. For s-tetrazine, the single point energy calculations at CCSD(T)/6-311G(d, p), B3LYP/6-311+G(2df, 2p), G3MP2B3, G3B3 and CCSD(T)/6-311+G(2df, 2p) levels were performed and compared, and G3MP2B3 was found both accurate and efficient. So, G3MP3B3 was selected to calculate the larger tetrazine derivatives, but the much larger BTATz and DAAT. For the reaction pathways with similar energy barriers, the rate constants were calculated with chemical reaction kinetics method to verify the main thermal decomposition pathway.

Based on the results of the simple hydronitrogen compounds and s-tetrazine without substituents, the thermal decomposition mechanisms and rate-determining steps of five tetrazine derivatives DAT, DHT, DiAT, BTATz, DAAT (Fig. 2) were studied in detail to illuminate the influence of substituents and bridge connecting units on stability of tetrazine rings.

The conclusions were drawn^[22]: (1) The tetrazine ring in the tetrazine derivatives may be broken by three modes: concerted triple dissociation, concerted double dissociation or single dissociation, however the concerted triple dissociation is generally dominant. (2) The stability of the tetrazine ring can evidently be strengthened by delocalization of these substituents with nitrogen. (3) The reaction between the substituents and tetrazine ring (such as H transfer) is not main thermal decomposition pathway and the thermal decomposition mechanisms of the substituent R are similar to that of its simple compound HR. (4) Tetrazines exhibit two principal modes of decomposition, which are the ring dissociation and the reaction of substituent groups. If the stability of the substituent is better than that of the tetrazine ring, decomposition occurs first through breaking of the ring. Moreover, some new conclusions about the effect of intermolecular reaction on the thermal decomposition mechanisms of tetrazines were achieved. The intermolecnlar reaction can make the decomposition pathway of molecular in crystal different from that of unimolecular.

On the other hand, some deficiencies for such a study still exit. On the DFT part, a larger box and higher cutoff energy for the plane-wave basis set would improve the accuracy. On the MD trajectory part, longer simulation time, both for the sampling of starting geometries and for the duration of each trajectory, would be essential if accurate branching ratio is needed. For the purpose of understanding the decomposition mechanism, these deficiencies are compensated by a separate set of calculations, using the conventional Gaussian-based molecular orbital method to locate the reaction barrier and transition structure for each type of the reaction channels observed in the trajectory study. In essence, the trajectory study based on ab initio MD method provides a lead from first principles for the elucidation of the reaction mechanisms by ab initio MO method.

5 Detonation properties calculations

All detonation properties were calculated by modified VLW code basing on the VLW equation of state (VLW EOS) (Scheme 2)^[64-66], law of mass conservation, law of momentum conservation, law of energy conservation and C-J detonation conditions^[67]. In the scheme 2, N, k and T^* is Avogadro Constant, Boltzmann constant and the nondimensional temperature, respectively. σ and ε are the Lennard-Jones potential parameters (Table 1)^[68-70]. The VLW code only needs the chemical compositions, heats of formation, and densities of compounds as input, and can calculate the detonation velocities, pressures and heats of the solid, liquid and gaseous explosives and propellants.

Where,

$ \begin{array}{l} {B^*}({T^*}) = \sum\limits_{j = 0}^\infty {{b^{(j)}}{T^{*-\left( {2j + 1} \right)/4}}} ;{b^{(j)}} =-\frac{{{2^{j + 1/2}}}}{{4j!}}\Pi \left( {\frac{{2j-1}}{4}} \right);\\ {b_0} = \frac{2}{3}{\rm{ \mathsf{ π} }}N{\sigma ^3};{T^*} = \frac{{kT}}{\varepsilon };w = \frac{{{b_0}}}{V} \end{array} $

Then with the calculated values of HOFs (Eq.2) and ρ_en obtained by enclosed volume of electron cloud (Eq.1), the detonation velocity (v_D) and pressure (p_D) of tetrazole energetic metal complexes (TEMCs) that contain Co, Cu and Zn were evaluated with the modified VLW code (Table 2)^[43]. It should be noted that the original VLW code can only deal with the compounds that compose with C, H, N, O, F and Al.

The calculated v_D of BNCP is closed to that of ref.[73], which indicates the results calculated by the modified VLW code to be credible. Accordingly, Cu(DAT)₂Cl₂ and Zn(DAT)₂Cl₂ are suggested to be energetic material candidates.

6 Interface properties calculation

Dissipative particle dynamics (DPD) method is a coarse-grained particle-based mesoscale dynamics simulation^[74]. It is useful for investigating the interface properties of energetic materials, especially for the energetic polymers. In a word, this method can give physical insight into the problem. Zhou^[75] used DPD method to study the interface properties of the immiscible A/B homopolymer blend systems in the presence of nanorods. All the calculations were based on the assumption that the nanorods were constructed by a string of beads connected through the spring forces. Then their rigidity was achieved depending on the angle forces.

The interaction between non-bonded DPD particles can be expressed by a conservative force FC, a dissipative force FD and a random force FR, respectively (Eq.4). Additionally, the interaction between bonded DPD particles can be described by a harmonic spring force FS and a angle force FA. The total force and three non-bonded interactions are as follows:

$ {f_i} = \sum\limits_{i \ne j} {F_{ij}^C + F_{ij}^D + F_{ij}^R + F_{i, i + 1}^S + F_{i\;1, i + 1}^A} $

(3)

$ \begin{array}{l} F_{ij}^C =-{\alpha _{ij}}{w^C}\left( {{r_{ij}}} \right){e_{ij}}\\ F_{\mathit{ij}}^D =-\gamma {w^D}\left( {{r_{ij}}} \right)\left( {{e_{ij}} \cdot {v_{ij}}} \right){e_{ij}}\\ F_{ij}^R =-\sigma {w^R}\left( {{r_{ij}}} \right){\xi _{ij}}\Delta {t^{ - 0.5}}{e_{ij}} \end{array} $

(4)

Where r_ij=r_i-r_j, r_ij= |r_ij| and v_ij=v_i-v_j. ξ_ij is a random number with zero mean and unit variance. α_ij stands for the maximum repulsion which reflects the chemical characteristics of interacting particles. w^C, w^D, and w^R are three weight functions. For w^C, a simple form is chosen as w^C(r_ij)= 1-r_ij when r_ij < 1, and w^C(r_ij)= 0 when r_ij≥1. On the other hand, to satisfy the fluctuation-dissipation theorem, w^D and w^R have a certain relationship^[76]:

$ {w^D}\left( {{r_{ij}}} \right) = {\left[{{w^R}\left( {{r_{ij}}} \right)} \right]^2};{\rm{ }}{\sigma ^2} = 2\gamma {k_B}T $

(5)

Here they also utilized the same expressions as w^C to describe w^D and w^R, and this method has been proved correct by Groot and Warren^[77]. In addition, the forces describing the connected particles were obtained by the differential of spring and angle-potential:

$ \begin{array}{l} F_{(i, i + 1)}^S =- \nabla U_{(i, i + 1)}^S;U_{(i, i + 1)}^S = \sum\limits_i {0.5{k_S}{{({l_{(i, i + 1)}}- {l_0})}^2}} ;\\ F_{(i- 1, i + 1)}^A = - \nabla U_{(i - 1, i, i + 1)}^A;U_{(i - 1, i, i + 1)}^A\\ \;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_i {{k_A}[1-{\rm{cos}}({\varphi _{i-1, i, i + 1}}-{\varphi _0})]} \end{array} $

(6)

Where l_{(i, i+1)} is the bond length between the connected two particle i and i+1, $φ$_{(i-1, i, i+1)} is the bond-angle of the adjacent three particle i-1, i and i+1.

In DPD, the polymers can be represented by the particles connected with the spring force. If the angle force is introduced additionally, those polymers containing rigidity or semi-rigidity segments can be also described successfully. Recently, we have used the two bonded forces to reinforce the stiffness of nanorods. Firstly, the equilibrium bond length l₀=0.5 was fixed, and the fluctuation of bond lengths was confined by a large spring coefficient k_S=50. Secondly, we set the equilibrium angle $φ$_{(i-1, i, i+1)}=π and chose a larger bending coefficient k_A=100. The angle, driven by the angle force between two consecutive bonds, is hoped to be close to the value of π providing satisfactory rod-like formation. Actually, the method we used is similar to that Chou used^[78], except that they mainly emphasized the function of the spring force. Therefore, they utilized a larger spring coefficients (k_S=100) and a comparative small bending constant (k_S=20).

In the present system, three DPD particles (A, B, and R) were used to describe the homopolymer and nanorods. For the two immiscible homopolymers (A_x and B_y), x=y=20, k_S=4 and l₀=0 were fixed. Accordinglly, the unstretched bond length of the nanorod was l₀, and its aspect ratio was approximately equal to (n-1)× l₀~f(n). Finally, the aspect ratio of nanorods can be described by a simple parameter, n.

It should be noted that all the DPD simulations were performed on the commercial molecular modeling software package Materials Studio (MS) program. The radius of interaction, the particle mass and the temperature as r_c=m=k_BT=1, σ=3 were all on the basis of the defaults value of MS program. We mainly focus on the effect of nanorods′ length and volume fraction on the interface of the immiscible homopolymer blends. Therefore, only one type of parameters were used as x=y=20 and α_AB=50 according to the work of Qian et al.^[79]. More details relationship between a_ij or l and the interfacial characteristics has been reported by Zaman^[80]. We set the repulsive parameters α_AR=α_BR=25. Then the modified nanorods have the function of block copolymer, which means that nanorods can be immersed in both A and B phase at the same time. Since the nanorod can be modified by different functional groups, the interaction between nanorods and blends can effectively be adjusted, and the hypothesis should come out true^[81-82].

A recent experimental work by Composto′s group^[83] on gold nanorods implies that an experiment of surface-modified nanorods in binary polymer blend may be devised very well. The interaction between the nanorods and the two polymers can be achieved through the surface modification. Then the theoretical work of Hore et al.^[84-85] also point out that the nanorods should be considered as a viable emulsifying agent for immiscible polymer blends.

7 Conclusion

Improved theoretical methods are of significant importance for evaluating the properties of energetic materials. In this review, the main methodologies, that our team members usually utilized to perform several indispensable properties, were introduced in detail.

(1) Based on each optimized structure, the enclosed volume (V_en) of electron cloud around the molecule was calculated. ρ_en, which can be used as the crystal density of a energetic metal complex, was obtained by equation (1);

(2) The HOF of energetic compounds can be obtained by atomization scheme 1;

(3) The combined ab initio MD and ab intio MO study can successfully revealed the dissociation mechanisms for some simple energetic materials. The deficiencies for decomposition mechanism study are compensated by a separate set of calculations, using the conventional Gaussian-based molecular orbital method to locate the reaction barrier and transition structure for each type of the reaction channels observed in the trajectory study. In essence, the trajectory study based on ab initio MD method may provide the leads from first principles for the elucidation of the reaction mechanisms by ab initio MO method;

(4) By the modified VLW code, the detonation properties can be effectively calculated;

(5) Dissipative particle dynamics (DPD) are useful analysis tools for studying the self-assembly of energetic polymers and can give physical insight into the problem. The nanorods are mainly described by the angle force and interact favorably with the two homopolymer through the three forces (F^C, F^D and F^R). The comparisons with the experimental and theoretical studies proved that DPD is intrinsically promising in the simulations.