Comparative study of melting points of 3,4-bis(3-nitrofurazan-4-yl)furoxan (DNTF)/1,3,3-trinitroazetidine (TNAZ) eutectic compositions using molecular dynamic simulations

Abstract

3,4-Bis(3-nitrofurazan-4-yl)furoxan (DNTF) and 1,3,3-trinitroazetidine (TNAZ) have been widely investigated as important candidate components in melt cast explosives. This paper presents a study of melting point prediction of DNTF/TNAZ eutectic compositions using a molecular dynamics simulations method. The melting points were determined according to the variation of various parameters including specific volume, free volume, diffusion coefficient, specific heat capacity and non-bonded energy, and good agreements were observed from the comparison between calculated values and prevenient experimental data. The binding energies (E_bind) and radial distribution functions (RDFs) were also performed to explore the interactions between DNTF and TNAZ molecules, and results reveal that a weak hydrogen bond from H of TNAZ and O of DNTF vanishes from solid to liquid phase. The detonation performances of DNTF/TNAZ eutectic were also discussed by comparing with other common castable explosives. Results show that the high performance and low melting point make the mixture DNTF/TNAZ (4/6) a good candidate for melt cast explosives. In general, it is worth noting that molecular dynamics simulations provide a useful tool to understand the thermal properties of energetic eutectic compositions.

---

1. Introduction

The melt cast explosive, also known as eutectic explosive, means an explosive mixture with a lower melting point than either component.¹ Compared to mechanical pressing, melt casting is more economical and convenient for large-scale filling in different styles of chamber, and is widely used in military applications such as landmines, high-explosive anti-tank (HEAT) warheads, rocket projectiles and missiles.^2,3 2,4,6-Trinitrotoluene (TNT) is currently the most commonly used casting medium for melt cast explosives because of its low melting point (354 K), high chemical stability, low mechanical sensitivity and capability of binding explosive particles.⁴ Unfortunately, the crystal density (ρ = 1.64 g cm⁻³) and detonation velocity (D = 6940 m s⁻¹) are relative low for TNT, and it is also plagued by shrinking and cracking on cooling.¹ Therefore, explorations of new energetic eutectic systems with higher detonation performance and mechanical properties have attracted widespread attention.^5–8

3,4-Bis(3-nitrofurazan-4-yl)furoxan (DNTF) is known as a typical energetic furoxan compound with high crystal density (ρ = 1.937 g cm⁻³) and detonation velocity (D = 9250 m s⁻¹).^9,10 Due to the low melting point (383 K), DNTF has the capability to be candidate component in melt cast explosives.¹¹ However, the high impact sensitivity makes it fail to meet the strong requirements of safety.¹² Another promising explosive in this class is 1,3,3-trinitroazetidine (TNAZ), which derives its virtues from low melting point (373 K), relative low sensitivity, good performance (30% greater than TNT) and compatibility.^13,14 Some binary and ternary eutectic explosive systems based on TNAZ and DNTF have been studied for thermal and explosive properties. For example, Jiang et al. prepared the eutectic mixture of TNAZ and DNTF, which exhibits lower melting point than each single component and moderate sensitivity closed to TNAZ.¹⁵ The eutectic ternary phase diagrams of TNAZ/DNTF/RDX (TNT) were investigated by Shao et al. with high pressure differential scanning calorimeter.¹⁶ In general, constrained by time-consuming and expense, the efficiency of experimental testing is relative low.

As an alternative approach, molecular dynamic (MD) simulations can reveal both physical and chemical properties at atomistic and molecular scales. Recently, MD simulations have been widely applied to study the structures and properties of polymer-bonded explosives (PBXs),¹⁷ propellants¹⁸ and cocrystal explosives.¹⁹ Melting point (T_m) is the crucial physical properties in production processes and application for eutectic explosives. Some methodologies based on group contribution or group additivity can be applied to predict the melting points of pure energetic materials such as nitramines, nitrate esters and nitroaromatic compounds, but there is no indication that they are suitable to estimate the T_m of binary system.^20–23 By using MD simulations, T_m (or glass transition temperature) can be predicted precisely for organic compound (pure or binary), inorganic compound, polymer and ionic liquid.^24–28 For instance, Yang et al. calculated the glass transition temperature of an epoxy molding compound by conducting the MD simulations, and the results are in good agreement with experimental data.²⁹ Based upon the research of Agrawal et al., it is proved that AMBER-SRT force field is more suited to compute the melting point of TNAZ compared with AMBER force field.³⁰ Li et al. investigated the T_m of trans-1,4,5,8-tetranitro-1,4,5,8-tetraazadacalin (TNAD)-based binary energetic systems, the force field energy and intermolecular interactions were also discussed.³¹ However, theoretical study of the structural–property relationship of DNTF/TNAZ eutectic compositions have not yet been explored.

In the present work, the melting points of DNTF, TNAZ and DNTF/TNAZ eutectic compositions are predicted by molecular dynamics simulations. By observing the inflection point on different curves of specific volume vs. temperature, free volume vs. temperature, diffusion coefficient vs. temperature, specific heat capacity vs. temperature and non-bonded energy vs. temperature, the T_m of mixture 3 (DNTF/TNAZ with mass ratio of 4/6) are calculated and compared with prevenient experimental results. The analyses of the interactions between DNTF and TNAZ molecules are conducted by studying the binding energies and radial distribution functions, as well as related discussions. In addition, the energetic performances of mixture 3 are also calculated by empirical method.

2. Modeling and computational methods

2.1 Choice of force field

The COMPASS force field,³² which is applicable for simulation of the condensed-phase materials, is chosen to calculate the melting point and intermolecular interactions of the eutectic compositions. We have given accurate prediction of crystal morphologies of DNTF in different solvent system under COMPASS force field.¹² On the other hand, COMPASS force field also has been proven to be able to give accurate simulation of structures and properties of nitramine explosives such as 1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX) and 2,4,6,8,10,12-hexanitro-hexaazaisowurtzitane (CL-20).^33,34 Here, the validation of COMPASS force field was evaluated by comparing with experimental results of lattice parameters, density and lattice energy of TNAZ crystal from the optimized values calculated by COMPASS force field. The lattice energy (E_latt) is calculated from the difference between the total internal energy of the molecule in crystal and the corresponding energy in gas phase. The experimental lattice energy (E_latt) can be evaluated using the following formula:

E_latt = −ΔH_sub − 2RT (1)

where ΔH_sub is experimental sublimation enthalpy, R is the gas constant and T is the temperature.³⁵ From Table 1, it can be noted that all of the parameters obtained from calculation are in good agreement with the experimental values (within 5.0% deviation). Results show that the COMPASS force field is adapted for this simulation model.

Table 1 The comparison of the experimental and optimized lattice parameters, density and lattice energy of TNAZ

Lattice parameter	a/Å	b/Å	c/Å	α/°	β/°	γ/°	ρ/g cm^−3	Elatt/kcal mol^−1
Exp.^36,37	5.733	11.127	21.496	90.00	90.00	90.00	1.86	−24.14
COMPASS	5.451	11.407	21.020	90.00	90.00	90.00	1.95	−24.74
Relative error/%	−4.92	2.52	−2.21	0	0	0	4.84	2.49

---

2.2 Construction and simulation details of models

All simulations were carried out with program Materials Studio 5.5 (Accelrys Inc., USA) in this paper.³⁸ The initial DNTF and TNAZ structures used in the simulation are derived respectively from Sheremetev³⁹ and Archibald,³⁶ as shown in Fig. 1. The cubic simulation boxes of pure DNTF, TNAZ and five DNTF/TNAZ mixtures were constructed by Amorphous Cell module with periodic boundary conditions. The initial densities of each system were achieved according to the additivity of volume ratio for different blends. Table 2 presents the details for the model building, where mixture 1, mixture 2, mixture 3, mixture 4 and mixture 5 correspond to the blend of DNTF/TNAZ with mass ratio of 2/8, 3/7, 4/6, 5/5 and 6/4.

Fig. 1 The molecular structures of DNTF and TNAZ.

Table 2 The modeling details of pure DNTF, TNAZ and DNTF/TNAZ mixtures

System	Number of DNTF	Number of TNAZ	Number of atom	Initial density/g cm^−3
DNTF	50	—	1100	1.93
TNAZ	—	50	850	1.86
Mixture 1	10	65	1325	1.87
Mixture 2	15	57	1299	1.88
Mixture 3	21	51	1329	1.89
Mixture 4	27	44	1342	1.89
Mixture 5	34	37	1377	1.90

---

The geometric optimization was implemented for the amorphous molecular models with the Forcite module until the energy of generated amorphous cell is minimized to convergence value of 0.001 kcal mol⁻¹. In order to get a fully relaxed system, an annealing treatment was performed for the above models under the canonical (NVT) ensemble from 480 K to 240 K with an interval of 20 K. This process was run at the cooling rate of 10 K/10 ps under Nosé–Hoover thermostat⁴⁰ to make system uniformly distribute. Afterwards, following MD simulation was carried out with the NPT (isothermal–isobaric) ensemble for another 1 ns at each temperature under atmospheric pressure. The first 500 ps run for equilibrium, and the following 500 ps run for production. The final equilibrium structure from previous temperature was taken as the initial structure for the next dynamic step. During the simulation, the electrostatic interactions were calculated by using the standard Ewald method with accuracy of 0.0001 kcal mol⁻¹, and the van der Waals (vdW) forces were corrected by atom-based summation with a cutoff radius of 15.5 Å. “Nosé–Hoover” and “Berendsen”⁴¹ were selected as the thermostat and barostat method to control the temperature and pressure of systems. The data were recorded with 1 fs sampling interval. For illustration, the equilibrium structure of mixture 3 at 320 K is shown in Fig. 2 as an example.

Fig. 2 The equilibrium structure of the mixture 3 at 320 K.

3. Results and discussion

3.1 The melting points of DNTF and TNAZ

Generally, the density of most substances will decrease from solid to liquid state. That is because the packing structure of solid state is closer than that of liquid state. The specific volume of a substance is defined as the ratio of the substance's volume to its mass. That is to say, in the course of melting phase transition, the specific volume will undergo an abrupt change. Hence, by observing the point of inflection from the curve of specific volume versus temperature, the melting point of desire substance could be achieved. The applicability of this method is appreciated by former investigations.^19,25,31

In Fig. 3 and 4, we present the curve of specific volume vs. temperature for DNTF and TNAZ, respectively. It is observed that the specific volume increases slightly with temperature below 400 K for DNTF, afterwards the change rate of specific volume rises sharply. The calculated T_m for DNTF (384.8 K) is represented by the intersection of two lines obtained by least square fitting, which is only 1.8 K higher than the experimental value of 383 K (Fig. 3).⁹ Similarly, the predicted T_m of TNAZ (373.7 K) provides a reasonable match to experimental value of 374 K (Fig. 4).¹³ The above results demonstrate that COMPASS force field can well describe thermodynamics properties of DNTF and TNAZ, and this method is reliable for the melting point prediction of single energetic compound.

Fig. 3 Specific volume versus temperature relationship for DNTF.

Fig. 4 Specific volume versus temperature relationship for TNAZ.

3.2 The melting points of DNTF/TNAZ eutectic compositions

After satisfying results of pure component are derived, the melting points of DNTF/TNAZ eutectic compositions are also simulated by calculating the specific volume as a function of temperature. Corresponding graphs of specific volume vs. temperature data of five types of mixture with different ratio are shown in Fig. 5. For each system, it is clear that two straight lines are obtained and the specific volume is increased with increasing temperature accompanied by a change in the slope. In addition, the specific volume increases more rapidly with the increase of the temperature above the inflection point of each system, which may differentiate between solid and liquid state. It is also noted that the T_m values of all binary compositions are lower than that of pure DNTF and TNAZ, which indicated the formation of eutectic composition. Because of the mutual diffusion, DNTF and TNAZ are completely soluble in each other and the T_m value is changed due to the altered molecular interactions.

Fig. 5 The specific volume versus temperature relationship for DNTF/TNAZ eutectic compositions. (a), (b), (c), (d) and (e) correspond to the results of mixture 1, mixture 2, mixture 3, mixture 4 and mixture 5, respectively.

Table 3 presents the comparison of the melting points of DNTF/TNAZ eutectic compositions acquired from experiment¹⁵ and MD simulation. Results show that the data getting from both the two methods exhibit same changing tendency as the variation in proportion of DNTF and TNAZ. Each of the predictional melting point gives a reasonable match to the experimental value and the relative errors are all less than 1.0%. It is worth noting that the calculated T_m values for most cases are somewhat higher than the values obtained by DSC studies. That could be ascribed to the higher cooling rate used in MD simulations and also observed by other groups previously.^25,29 According to Table 3, among all five compositions, mixture 3 presents the lowest melting point and it is close to that of TNT (353.6 K). This reveals that mixture 3 meets the requirement of casting procedure.

Table 3 The melting points of DNTF/TNAZ eutectic compositions

Composition	Tm/K (MD simulation)	Tm/K (ref. 15)	Relative error/%
Mixture 1	364.5	365.8	−0.4
Mixture 2	362.2	359.9	0.6
Mixture 3	355.4	352.2	0.9
Mixture 4	359.7	358.2	0.4
Mixture 5	366.0	364.0	0.5

---

To explore further possibilities of melting point simulation, the different methods are examined by taking mixture 3 as an example. The free volume (V_f) means the unoccupied space in the molecular matrix, which can be evaluated using the following equation:

V_f = V − V_occ (2)

where V is specific total volume, and V_occ is occupied volume. V_occ can be identified with van der Waals volume of the molecules. The specific total volume has an abrupt change during the transformation from solid to liquid phase. Accordingly, the free volume will come through a discontinuous change when the material undergo melting process or glass transition.^42–44 In Fig. 6, the free volumes of mixture 3 are reported with respect to temperature by using 1.0 Å as the radius of the probe sphere. It can be seen that there is only a slight increase of free volume when the temperature is lower than 340 K, afterwards the rate of increasing of free volume is rising nine times faster. The obvious change in slope signifies the melting point of the mixture as 356.0 K, which is also found to be comparable with the experimental results.

Fig. 6 The free volume versus temperature relationship for mixture 3.

The mean square displacement (MSD) analysis is usually used to study the mode of displacement of a given particle followed over time. The diffusion coefficient (D), which is evaluated as the derivative of MSD with respect to time, can be measured by the Einstein relation as follows:

Generally, the diffusion coefficient proceeds at a higher rate in liquids than that of solids. Fig. 7 shows all atoms' MSD versus time curves at temperatures of 240 K, 280 K, 320 K, 360 K, 400 K and 440 K, for the first 400 ps. The results of MSD for all temperatures are depicted in Fig. S1 of the ESI.† It is noted that all MSD curves exhibit good linearity as a function of time, and the diffusional behavior continuously increases with temperature. The diffusion coefficient as a function of temperature can be seen in Fig. 8 by calculating from the slope of the MSD versus time during the 100–300 ps time interval. An obvious discontinuity is shown at 351.8 K for mixture 3 as the transition point from solid to liquid. The difference between calculation of diffusion coefficient and experimental result of the melting temperature is only 0.4 K, which provides positive evidence to validate the reliability of this method.

Fig. 7 The mean square displacement versus time curves for mixture 3 at different temperatures.

Fig. 8 The diffusion coefficient versus temperature relationship for mixture 3.

The specific heat capacity C_p of an object is the quantity of heat demanded to raise the temperature of unit mass by 1 K at constant pressure. From solid to liquid state, the value of C_p will show the different change rate due to the different molecular packing. From MD simulations, C_p can be calculated by the following equation:^25,45

where E_k is the kinetic energy, U is the potential energy, k_B is Boltzmann constant, T is the temperature, P is the pressure and V is the volume. The simulated specific heat capacity with respect to the temperature for mixture 3 is shown in Fig. 9. Although no experimental data of heat capacity of DNTF/TNAZ eutectic compositions are available, the predictional C_p remains credible because that it (1.06 J g⁻¹ K⁻¹) is close to that of pure DNTF (0.94 J g⁻¹ K⁻¹)⁴⁶ and TNAZ (1.14 J g⁻¹ K⁻¹)⁴⁷ at 300 K. By comparison with experimental melting point, the simulation value is slightly lower with relative error of 1.1%. Results show that the melting behavior of eutectic compositions also can be obtained by studying the variation of specific heat capacity.

Fig. 9 The specific heat capacity versus temperature relationship for mixture 3.

For a certain molecular system, the total potential energy U_total can be obtained by:

U_total = U_stretch + U_bend + U_torsion + U_non-bond (5)

where U_stretch, U_bend and U_torsion respectively represent the bond-stretching, bond-bending and torsion potentials, which are classified as bonded terms. In the COMPASS force field, the non-bonded terms (U_non-bond) include van der Waals interaction (represented by LJ-9-6 function,⁴⁸ U_LJ) and Coulomb interaction (U_Coulomb). For two interaction sites i and j separated by a distance of r_ij, the above-mentioned item is given by:where ε_ij denotes the Lennard-Jones energy, the q_i and q_j denotes partial charges on sites i and j. The non-bonding attractions enable this composition to exist as solid state at low temperature, but cannot withstand disruptions at higher temperature caused by thermal energy. Therefore, it is feasible to obtain the indication of melting by observing the variation of non-bonded energy with temperature. The simulated melting point or glass transition temperature reported by previous studies according to this method match well with experimental data.^29,31,49 From the relationship plotted in Fig. 10, it can be noted that the general trend of non-bonded energy is descended with the temperature increase, and the break in its slope means the melting point. The T_m calculated from non-bonded energy vs. temperature is a bit higher than the experimental value, while it is practically identical to the result from specific volume–temperature relationship.

Fig. 10 The non-bonded energy versus temperature relationship for mixture 3.

Comparatively speaking, out of these five options, the methods based on specific volume and diffusion coefficient perform better in melting point prediction of DNTF/TNAZ eutectic composition for the simple calculation process and precise results. On the other hand, the results from free volume, specific heat capacity and non-bonded energy suffer from relative large error and nonlinear data, therefore they are less successful than the first two methods.

3.3 Binding energy and radial distribution function analysis

To better understand the chemical–physical properties of DNTF/TNAZ eutectic compositions, the interactions between DNTF and TNAZ molecules are investigated based on MD simulations in this section. The binding energy (E_bind) means the amount of energy required to separate a group of particles from the whole system, which is defined as the negative value of the intermolecular interaction energy (E_inter). For DNTF/TNAZ eutectic system, the molecular interactions can be calculated by the total energy and individual component energy, and E_bind between DNTF and TNAZ can be estimated as follows:

E_bind = −E_inter = −(E_bind − E_DNTF − E_TNAZ) (7)

where E_total is the average total energy of the DNTF/TNAZ blend, E_DNTF and E_TNAZ are the average energies of DNTF and TNAZ, respectively. E_bind represents the strength of intermolecular interactions among different components, which is used to evaluate the compatibility and thermal stability of the eutectic systems. Table 4 lists the binding energies (E_bind) of mixture 2, 3 and 4 at different temperatures. For all three mixtures, it is obviously that the binding energies decrease gradually with the increasing of temperature. which indicates the decline of thermodynamic stability and also can be found in other systems such as cocrystal.⁵⁰ Furthermore, it is found that the order of E_bind of different compositions is written as follows: mixture 4 > mixture 3 > mixture 2, which suggests that the compatibility is improved with the increasing of the proportion of DNTF. However, this conclusion has not been confirmed due to the lack of support from experimental data.

Table 4 The binding energies of mixture 2, 3 and 4^a

System	T/K	Etotal	EDNTF	ETNAZ	Ebind
a All energies are in kcal mol^−1.
Mixture 2	280	−4098.49	−11.26	−3311.79	775.44
	320	−3910.59	31.56	−3203.63	738.52
	360	−3695.23	65.62	−3053.48	707.37
	400	−3469.95	118.55	−2907.06	681.44
	440	−3310.32	143.24	−2826.67	626.89
Mixture 3	280	−3853.25	−68.81	−2858.67	924.77
	320	−3676.50	−31.8	−2758.96	885.74
	360	−3411.78	35.22	−2589.13	857.87
	400	−3262.99	62.70	−2542.19	783.50
	440	−3055.23	112.42	−2407.53	760.12
Mixture 4	280	−3551.91	−189.61	−2369.98	992.32
	320	−3356.52	−121.05	−2246.66	988.81
	360	−3148.78	−37.35	−2159.75	951.68
	400	−2982.41	−14.31	−2057.55	910.55
	440	−2835.18	64.72	−2020.45	879.45

---

The radial distribution functions (RDFs) g(r) is often used to investigate the local structure evolution during the phase transition process by quantifying the distributions of one atom around the reference atom. Generally, the distance range (r) of hydrogen bonding interaction is 2.0–3.1 Å and for vdW interaction is 3.1–5.0 Å. When distance range is larger than 5.0 Å, the Coulomb interaction plays a major role. In the following, two kinds of atom pairs (H–O and H–N) are studied separately. The H atoms are derived from TNAZ molecules, while the O and N atoms are taken from DNTF molecules. Fig. 11 shows the RDFs for H–O and H–N atomic pairs in DNTF/TNAZ eutectic composition at different temperatures of 320 K, 340 K, 360 K and 380 K by taking mixture 3 as example. In Fig. 11(a), the first peaks appear at 2.7 Å for r(H–O) at 320 K and 340 K, which corresponds to the formation probability of the weak hydrogen bond. On the other hand, for r(H–O) at 360 K and 380 K, the first peaks at 3.1 Å are observed, representing the existence of vdW interactions. That is to say, the formation probability of the hydrogen bond between H atoms and O atoms decreases with the increasing of temperature, and simultaneously the packing of molecules becomes more and more loose, which results in the appearance of phase transition. The fact that the experimental T_m locates in the interval of 340 K and 360 K may verify the validity of RDFs analysis. In addition, based on the concept presented by Allen et al., the peaks of RDFs in liquid phase are usually blunt and low due to the unpredictable fluidity.⁵¹ It also can be confirmed by comparison of curves of 320 K and 380 K in Fig. 11(a). On the contrary, the first peaks of r(H–N) occur in the range of 3.3–3.5 Å from Fig. 11(b), which suggests that the interactions between H atoms and N atoms are mainly vdW and Coulomb interactions.

Fig. 11 RDFs analysis for mixture 3 at different temperatures. (a) Corresponds to H–O atomic pairs, (b) corresponds to H–N atomic pairs.

3.4 Detonation performance

Detonation represents a form of reaction in which the explosive substance undergoes chemical reactions at high speed and produces a shock wave. The detonation performance of an explosive principally determined by the velocity of detonation (D) and detonation pressure (P). For a C, H, O, N type explosive, there is a simple linear relationship between the velocity of detonation and the detonation factor F, which can be derived from the chemical constitution and molecular structure as follows:⁴where nO, nN, nH are the numbers of oxygen, nitrogen and hydrogen atoms in a molecule, nB is the number of oxygen atoms in excess of those already available to form CO₂ and H₂O, nC is number of oxygen atoms double bonded directly to carbon, nD is the number of oxygen atoms singly bonded directly to carbon, nE is the number of nitrate groups either as nitrate-esters or nitrate salts. A = 1 if the compound is aromatic, otherwise A = 0; G = 0.4 for a liquid explosives and G = 0 for a solid explosives. The D values of some mixed explosives calculated with these equations have been proven to be closed to the experimental data.⁵² The detonation pressure fits the following relationship:⁴

P = 93.3D − 456 (10)

Table 5 lists the predicted detonation velocity and pressure of mixture 3 and the other three types of common castable explosives. Inevitably, the D and P of COMP B-3, CYCLOTOL and OCTOL exhibit relative low value, which are restricted by the poor performance of TNT. On the other hand, benefit from the high detonation properties of both DNTF and TNAZ, the D and P of mixture 3 are obviously higher than that of TNT-based explosives. Combined with the fact of low melting point, it is suggested that mixture 3 has considerable potential to be applicated as melt cast explosive.

Table 5 The detonation performances for mixture 3, COMP B-3, CYCLOTOL and OCTOL^a

Explosive	Component	D/km s^−1	P/GPa
a RDX is 1,3,5-trinitroperhydro-1,3,5-triazine, and HMX is octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine.
Mixture 3	DNTF/TNAZ = 40/60	8.79	36.4
COMP B-3	RDX/TNT = 60/40	7.93	28.4
CYCLOTOL	RDX/TNT = 70/30	8.25	31.4
OCTOL	HMX/TNT = 75/25	8.56	34.3

---

4. Conclusions

In this paper, molecular dynamics simulations were conducted to investigate the melting points of DNTF/TNAZ eutectic compositions on molecular point of view. Five DNTF/TNAZ systems were constructed with amorphous structure by performing NPT-MD calculations. We firstly studied the T_m of DNTF, TNAZ and DNTF/TNAZ eutectics based on the variation of specific volume as a function of temperature. Results show that the predicted T_m for different systems are all well correlated with the experimental data and the relative errors are within 1.0%. Take mixture 3 for example, the free volume, diffusion coefficient, specific heat capacity and non-bonded energy are all proved to be effective parameters in justifying the phase transition temperature. The inflection point observed from the curves of diffusion coefficient versus temperature provides the most accurate T_m value of mixture 3. The binding energy values between DNTF and TNAZ decrease gradually as the temperature rises, indicating the decline of thermodynamic stability. The radial distribution function analysis reveals that a weak hydrogen bond from H of TNAZ and O of DNTF exists in the solid phase, and vanishes after the phase transition. Finally, compared with the current TNT-based castable explosives, the detonation properties of DNTF/TNAZ eutectic are obviously improved. These simulations are helpful to offer the variation of related parameters during melting process at the molecular level, and may also provide an effective theoretical method for exploring new eutectic explosives with improved performances.

Acknowledgements

We gratefully acknowledge the support for this work by the National Natural Science Foundation of China (Grant No. 21373157 and Grant No. 51511130036).

References

1. P. Ravi, D. M. Badgujar, G. M. Gore, S. P. Tewari and A. K. Sikder, Propellants, Explos., Pyrotech., 2011, 36, 393–403.
2. D. Sun and S. V. Garimella, Propellants, Explos., Pyrotech., 2005, 30, 369–380.
3. D. Cao, Y. Li, Y. Du, J. Wang and Y. Li, Chin. J. Energ. Mater., 2013, 21, 157–165.
4. J. P. Agrawal, High Energy Materials: Propellants, Explosives and Pyrotechnics, Wiley-VCH & Co. KGaA, Weinheim, 2011.
5. R. Vijayalakshmi, S. Radhakrishnan, P. Shitole, S. J. Pawar, V. S. Mishra, R. K. Garg, M. B. Talawar and A. K. Sikder, RSC Adv., 2015, 5, 101647–101655.
6. L. Chen, H.-R. Li, R.-J. Xu, Y. Xiong, T. Xu, B.-Y. Huang and Y.-J. Shu, Propellants, Explos., Pyrotech., 2014, 39, 217–223.
7. A. Elbeih and S. Zeman, Cent. Eur. J. Energ. Mater., 2014, 11, 501–513.
8. Y.-H. Shao, X.-N. Ren and Z.-R. Liu, J. Therm. Anal. Calorim., 2010, 101, 1135–1141.
9. Y. Zhou, B. Wang, J. Li, C. Zhou, L. Hu, Z. Chen and Z. Zhang, Acta Chim. Sin., 2011, 69, 1673–1680.
10. R. Tsyshevsky, P. Pagoria, M. Zhang, A. Racoveanu, A. DeHope, D. Parrish and M. M. Kuklja, J. Phys. Chem. C, 2015, 119, 3509–3521.
11. A. I. Stepanov, D. V. Dashko and A. A. Astrat'ev, J. Mol. Graphics Modell., 2012, 9, 329–342.
12. N. Liu, Y. Li, S. Zeman, Y. Shu, B. Wang, Y. Zhou, Q. Zhao and W. Wang, CrystEngComm, 2016, 18, 2843–2851.
13. N. Sikder, A. K. Sikder, N. R. Bulakh and B. R. Gandhe, J. Hazard. Mater., 2004, 113, 35–43.
14. M. B. Talawar, A. Singh, N. H. Naik, B. G. Polke, G. M. Gore, S. N. Asthana and B. R. Gandhe, J. Hazard. Mater., 2006, A134, 8–18.
15. Q. Jiang, H. Wang, Y. Luo, W. Wang, H. Wang, J. Gao and K. Zhao, Initiators Pyrotech., 2014, 3, 42–45.
16. Y.-H. Shao, X.-N. Ren, Z.-R. Liu and X. Zhang, J. Therm. Anal. Calorim., 2011, 103, 617–623.
17. Y. Yu, S. Chen, X. Li, J. Zhu, H. Liang, X. Zhang and Q. Shu, RSC Adv., 2016, 6, 20034–20041.
18. Y. Zhao, X. Zhang, W. Zhang, H. Xu, W. Xie, J. Du and Y. Liu, J. Phys. Chem. A, 2016, 120, 765–770.
19. R. Sarangapania, S. T. Reddy and A. K. Sikder, J. Mol. Graphics Modell., 2015, 57, 114–121.
20. K. G. Joback and R. C. Reid, Chem. Eng. Commun., 1987, 57, 233–243.
21. M. H. Keshavarz, J. Hazard. Mater., 2006, A138, 448–451.
22. M. H. Keshavarz and H. R. Pouretedal, J. Hazard. Mater., 2007, 148, 592–598.
23. M. H. Keshavarz, J. Hazard. Mater., 2009, 171, 786–796.
24. S. W. Watt, J. A. Chisholm, W. Jones and S. Motherwell, J. Chem. Phys., 2004, 121, 9565–9573.
25. D.-X. Li, B.-L. Liu, Y.-S. Liu and C.-L. Chen, Cryobiology, 2008, 56, 114–119.
26. T. Arima, K. Idemitsu, Y. Inagaki, Y. Tsujita, M. Kinoshita and E. Yakub, J. Nucl. Mater., 2009, 389, 149–154.
27. Y. Zhang and E. J. Maginn, Phys. Chem. Chem. Phys., 2012, 14, 12157–12164.
28. Y. Zhang and E. J. Maginn, J. Chem. Phys., 2012, 136, 144116.
29. S. Yang and J. Qu, Polymer, 2012, 53, 4806–4817.
30. P. M. Agrawal, B. M. Rice, L. Zheng, G. F. Velardez and D. L. Thompson, J. Phys. Chem. B, 2006, 110, 5721–5726.
31. X. Li and X. Ju, Chin. J. Chem. Phys., 2014, 27, 412–418.
32. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
33. J.-J. Xiao, W.-R. Wang, J. Chen, G.-F. Ji, W. Zhu and H.-M. Xiao, Comput. Theor. Chem., 2012, 999, 21–27.
34. X.-J. Xu, H.-M. Xiao, J.-J. Xiao, W. Zhu, H. Huang and J.-S. Li, J. Phys. Chem. B, 2006, 110, 7203–7207.
35. G. Chen, C. Chen, M. Xia, W. Lei, F. Wang and X. Gong, RSC Adv., 2015, 5, 25581–25589.
36. T. G. Archibald, R. Gilardi, K. Baum and C. George, J. Org. Chem., 1990, 55, 2920–2924.
37. M. Sućeska, M. Rajić, S. Matečić-Mušanić, S. Zeman and Z. Jalový, J. Therm. Anal. Calorim., 2003, 74, 853–866.
38. Material Studio 5.5, Acceryls Inc., San Diego, 2010.
39. A. B. Sheremetev, E. A. Ivanova, N. P. Spiridonova, S. F. Melnikova, I. V. Tselinsky, K. Yu. Suponitsky and M. Yu. Antipin, J. Heterocycl. Chem., 2005, 42, 1237–1242.
40. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069–4074.
41. H. J. C. Berendsen, J. P. M. Postma, W. F. Gunsteren, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
42. G. Dlubek, D. Kilburn, V. Bondarenko, J. Pionteck, R. Krause-Rehberg and M. A. Alam, 24th Arbeitskreistagung Nichtkristalline Strukturen of DGK, Jena, September, 2003.
43. S. Palzer, Agglomeration of Dehydrated Consumer Foods, in Granulation, ed. A. D. Salman, M. J. Hounslow and J. P. K. Seville, Elsevier Science, Amsterdam, 2006.
44. T. G. Fox and S. Loshaek, J. Polym. Sci., 1955, 15, 371–390.
45. Z. Rao, S. Wang and F. Peng, Int. J. Heat Mass Transfer, 2013, 64, 581–589.
46. H.-X. Gao, F.-Q. Zhao, R.-Z. Hu, K.-Z. Xu, H. Zhang, P. Wang, Z.-M. Du, S.-Y. Xu, J.-H. Yi, H.-X. Ma, C.-R. Chang and J.-R. Song, Chem. Res. Chin. Univ., 2008, 29, 981–986.
47. H. Zhao, R. Hu, X. Wang, F. Zhao, H. Gao, H. Zhang, X. Zhang, Y. Feng and H. Ma, Acta Chim. Sin., 2009, 67, 2536–2540.
48. S. Lifson, A. T. Hagler and P. Dauber, J. Am. Chem. Soc., 1979, 101, 5111–5121.
49. H. Feng, J. Zhou and Y. Qian, J. Chem. Phys., 2011, 135, 144501.
50. T. Sun, J. J. Xiao, Q. Liu, F. Zhao and H. M. Xiao, J. Mater. Chem. A, 2014, 2, 13898–13904.
51. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, New York, 1987.
52. H. Muthurajan, R. Sivabalan, M. B. Talawar and S. N. Asthana, J. Hazard. Mater., 2004, A112, 17–33.

---

Footnote

† Electronic supplementary information (ESI) available: All of the data used in curve fitting are listed in Tables S1–S6. The results of MSD for all temperatures are shown in Fig. S1. See DOI: 10.1039/c6ra12041e

---