Exploring aliphatic nitro azides for plasticizers: a combined DFT and MD investigation

Abstract

Nitro azide plasticizers have been attracting more and more attention due to their excellent performance. To search for new and promising nitro azides for plasticizers, in this work, a systematic theoretical investigation was performed using the density functional theory (DFT) and molecular dynamics (MD) methods. In the first part, a series of aliphatic nitro azides (M1–M8) were designed and studied using the DFT method. The results show M2–M8 all possess high chemical and thermal stabilities. As a plasticizer, it should have not only a good stability but also the ability to improve the mechanical properties of composites. Therefore, in the second part, M2 was taken as an example to explore the plasticizing effect of nitro azides on GAP, an attractive azide binder. GAP/M2 composites with the mass ratios of 77.5/22.5 (I), 56.4/43.6 (II), and 36.5/63.5 (III) were constructed and simulated using the MD method. Results show that M2 has a good compatibility with GAP and can effectively improve the mechanical properties of GAP, which suggests M2 is a promising plasticizer of GAP. Based on the similar structures of M2–M8, M3–M8 may also be a promising plasticizers of GAP and are worth further experimental investigations.

---

1. Introduction

Plasticizers and binders are two important components of propellants and significant studies have been reported on their applications in gun and rocket propellants.¹ One of the earliest binders in energetic materials was a mixture of nitrocellulose and nitroglycerine,² later, various binders like hydroxy-terminated polybutadiene (HTPB), carboxy-terminated polybutadiene (CTPB), hydroxy-terminated polyethers (HTPE), and glycidyl azide polymer (GAP), were developed. GAP has a low glass transition temperature (−45 °C) and a low weight percentage of polymer weight-bearing chains, which results in an energetically favorable binder system.³ However, at low temperatures, GAP is hard and brittle, i.e., it suffers from poor mechanical properties under low temperatures. To overcome this problem, the use of nitro azide (–NO₂ and –N₃) plasticizers in GAP binder is preferred.⁴ So far, the number of synthesized nitro azide plasticizers is limited. To search for new and promising nitro azide plasticizer, in this work, a series of aliphatic nitro azide compounds (M1–M8, Fig. 1) were designed and studied using the density functional theory (DFT) and molecular dynamics (MD) methods. Theoretical studies based on the quantum chemical calculations have provided an effective way to screen the candidates for energetic compounds for decades. The thoughts of this study are as follows:

Fig. 1 Structures of studied compounds.

In the first section, the chemical and thermal stability of all designed compounds were investigated using the DFT method to screen the molecule with the high stability. The compounds, M2–M8, were found. As a plasticizer, it should not only have a good stability, but also can improve the mechanical properties, thus, in the second section, we chose M2 as an example to explore the plasticizing effect of nitro azides on the GAP. Various composites of M2 and GAP were constructed and simulated using the MD method. To better characterize M2, the molecular packing, IR, and NMR spectra of this compound were analyzed in the third section.

This work provides a train of thought to design and screen the compounds as plasticizers for propellants and gives a practical application for GAP/M2 composite. In addition, the analyses of molecular packing and spectra are helpful to identify M2 for the experimental researchers who are interested in this new compound.

2. Computational details

2.1 DFT calculation

Utilizing the Gaussian program package,⁵ all designed compounds were optimized with the DFT-B3LYP^6,7 method and the 6-311++G**^8,9 basis set. The optimized structures were characterized to be the energy minima on the potential energy surface by vibrational analysis and the Cartesian coordinates of all optimized geometries are presented in the ESI.† Chemical stability was predicted by the energy gap (E_g) between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). Thermal stability was examined by calculation of the bond dissociation energies (E_BD) and activation energy (E_a) of all possible pyrolysis processes. The formulas are as follows:

E_g = E_LUMO − E_HOMO (1)

E_BD(A–B) = E_A˙ + E_B˙ − E_A–B (2)

E_a = E_TS − E_R (3)

where E_LUMO and E_HOMO are the energy of HOMO and LUMO; A–B stands for the neutral molecule and A˙ and B˙ for the corresponding radical products after the dissociation of A–B bond; E_A–B, E_A˙ and E_B˙ are their corresponding total energies after the correction of the zero-point energy (ZPE). E_TS and E_R are the total energies after the correction of ZPE for the transition state and reactant, respectively. The reliability of the TS was verified by analyzing the optimized structures, frequencies and IRC (intrinsic reaction coordinate).

2.2 MD simulation

Various GAP/M2 composites (I, II, and III) were constructed by addition of different ratio of M2 into GAP using the Materials Studio program package.¹⁰ The initial densities, the number of chains and atoms, and the weights of all constructed amorphous models, as well as the mass ratio of GAP/M2 in I, II, and III, are listed in Table 1. Considering the comparability of the models and the computer resources, the atom number and molecular weight of all involved models were taken to be about 1200 and 10 000–12 000, respectively. The mass ratios of I, II, and III were obtained according to the number of molecules of GAP and M2 in the blends. Here, GAP single chain is constructed by 20 repeat units (n). It has been proved that n = 20 is the most appropriate chain length to represent a GAP chain.¹¹ The crystal density of M2 was predicted using the molecular mechanics method with the Dreiding force field¹² and that of GAP came from ref. 13. The densities of the composites were calculated based on the weight percentage of each component in the blending systems. Energy minimization first and then the MD simulation were performed on these models. The details are similar to that presented in our previous work,¹¹ so they are not stated here. When the systems reach the equilibrium, the compatibility between GAP and M2, and the mechanical properties of various composites were studied.

Table 1 Parameters for all simulated amorphous units

	Initial density (g cm^−3)	Number of chains	Mass ratios	Number of atoms	Weight of molecule
M2	1.438	100		1200	11 800
GAP	1.300 (ref. 13)	5		1216	10 000
I	1.329	4/20	77.5/22.5	1212	10 312
II	1.357	3/40	56.4/43.6	1209	10 634
III	1.384	2/60	36.5/63.5	1206	10 956

---

The compatibility was evaluated by the solubility parameter (δ), which was calculated as follows:

CED = E_coh/V (5)

E_coh = E_vac − E_bulk (6)

where CED is the cohesive energy density, E_coh is the cohesive energy, E_vac and E_bulk are the energies in the vacuum state and amorphous state, respectively, and V is the molar volume. MD calculation provides an easy way to evaluate the CED of polymers and has been used in many studies.^14–16

Mechanical parameters, like tensile (Young's) modulus (E), bulk modulus (K), shear modulus (G), Poisson's ratio (γ), and Cauchy pressure (C), were predicted using eqn (7).

where μ and λ are the Lamé coefficients. From the statistical mechanics theory of elasticity,¹⁷ the most general relationship between stress and strain can be stated by the generalized Hooke's law: σ_i = C_ijε_j (i, j = 1–6), where C_ij are the elements of the elastic constant matrix, i.e., elastic coefficients. They manifest that there are different elastic effects everywhere in materials. Because of the existence of the strain, the elastic coefficient matrix of a material should satisfy the formula: C_ij = C_ji, even for an extremely anisotropic body and there are 21 independent elastic coefficients. For an isotropic solid, there are only two independent elastic coefficients (C₁₁ and C₁₂). According to the C₁₁ and C₁₂, the Lamé coefficients (μ and λ) can be calculated by the formula of μ = (C₁₁ − C₁₂)/2 and λ = C₁₂.

3. Results and discussion

3.1 Stability of all designed compounds

3.1.1 Chemical stability. Based on the frontier orbital theory,¹⁸ the energy gap E_g between HOMO and LUMO determines the molecular reactivity, such as the ability to absorb light and to react with other species. The molecule with a smaller E_g is expected to have a higher reactivity and a lower stability in the chemical or photochemical processes with electron transfer or leap.¹⁹ Table 2 summarizes E_HOMO, E_LUMO, and E_g of all designed compounds. E_g has the order of M1 > M2 > M3 > M4 > M5 > M6 > M7 > M8, suggesting the chemical or photochemical stability have the decreasing trends from M1 to M8. It is worth noting that all designed compounds have higher chemical stability (E_g = 4.89–5.44 eV) than the synthesized nitro azide compound 1,3-diazido-2-methyl-2-nitropropane (E_g = 4.85 eV (ref. 20)).

Table 2 Energies of the frontier orbits and their gaps

	M1	M2	M3	M4	M5	M6	M7	M8
EHOMO/a.u.	−0.30062	−0.29052	−0.28096	−0.27608	−0.27187	−0.26951	−0.26741	−0.26616
ELUMO/a.u.	−0.10056	−0.10096	−0.09511	−0.09233	−0.08981	−0.08833	−0.08708	−0.08630
Eg/eV	5.44	5.16	5.06	5.00	4.95	4.93	4.91	4.89

---

3.1.2 Thermal stability. Two possible decomposition reactions, i.e., homolysis of C–NO₂ and breakage of N–N₂, were considered for the pyrolysis of M1–M8. It is found that C–NO₂ ruptures directly with formation of two radical products, while the breakage of N–N₂ needs to pass through a TS with release of N₂. The required energies of the two processes are marked as E_BD(C–NO₂) and E_a(N–N₂), respectively (Table 3). Generally, the bond that requires the minimum energy to break is the weakest and is most likely to be the trigger bond. Obviously, for each compound, E_a(N–N₂) is much lower than E_BD(C–NO₂), suggesting that the pyrolysis of M1–M8 are all initiated from the rupture of N–N₂. It is observed that in these processes, the H transfer happens. The energy diagram for the pyrolysis of M2 is plotted in Fig. 2 as an example.

Table 3 Bond dissociation energy of C–NO₂ (kJ mol⁻¹) and activation energy of N–N₂ (kJ mol⁻¹)

	M1	M2	M3	M4	M5	M6	M7	M8
EBD(C–NO2)	166.21	215.29	213.46	218.14	220.37	221.93	222.81	222.38
Ea(N–N2)	139.19	152.17	156.91	158.74	159.59	159.97	160.21	160.38

---

Fig. 2 Potential energy diagrams for the pyrolysis of M2.

From Table 3, we can also see that with the increasing number of –CH₂– group, the E_a(N–N₂) increases fast from M1 to M2, and then increase slowly from M2 to M8. Except for M1, other compounds, i.e., M2–M8, have higher thermal stability than previously researched azide compounds.^21,22 That is to say, M2–M8 all have the potential as the energetic plasticizers. To study their plasticizing effects on GAP, the one with the lowest E_a(N–N₂), i.e., M2, was chosen as the target to do the following researches. Since M2–M8 have similar structures, if M2 has a good plasticizing effect on GAP, M3–M8 with higher stability will also be the promising plasticizers of GAP.

3.2 Plasticizing effect of M2 on GAP

Plasticizing is the process that polymer and plasticizer dissolve with each other, therefore, they should have a good compatibility first. The compatibility between GAP and M2 is predicted by the solubility parameter δ. The predicted δs of GAP and M2 are 18.42 and 24.28 MPa^0.5, respectively, with a difference of 5.86 MPa^0.5. Since the materials with similar δs are thermodynamically compatible²³ and according to the principle proposed by Greenhalgh et al.,²⁴ the system is miscible with Δδ < 7.0 MPa^0.5, whereas immiscible with Δδ > 10.0 MPa^0.5, it can be seen that GAP and M2 are compatible. The simulated δ of GAP gives a good consistent result with its experimental value (18.0–18.8 MPa^0.5 (ref. 25)), suggesting the constructed models are reliable and the predicted δ of M2 is believable.

Mechanical parameters, i.e., E, K, G, γ, K/G and C, of blends I–III, as well as those of GAP, are listed in Table 4. Obviously, I–III have lower E, while higher γ, K/G and C than GAP. Generally, a larger tensile modulus corresponds to a stronger rigidity,²⁶ a greater value of K/G and γ to a better tenacity and malleability,²⁷ a more positive C to a better ductility.¹⁷ That is to say, introduction of M2 into GAP makes rigidity decrease while tenacity, malleability, and ductility increase, i.e., M2 can enhance the mechanical properties of GAP. In addition, it can be found that the mechanical properties of GAP increase with the increasing content of M2. One can choose the optimal mass ratio of GAP/M2 according to the practical requirement. It is worth to mention that the plasticizing effect of M2 on GAP is comparable to that of DIANP (1,5-diazido-3-nitrazapentane),¹⁰ an azide plasticizer having received considerable attentions.^11,28–30 The equilibrated amorphous cells of I–III are displayed in Fig. 3.

Table 4 Elastic mechanical properties of various formulas

	E/GPa	K/GPa	G/GPa	γ/GPa	K/G	C/GPa
GAP	4.49	3.75	1.73	0.30	2.17	0.67
I	3.82	3.86	1.43	0.33	2.69	1.53
II	3.97	4.41	1.47	0.35	3.00	2.41
III	3.74	4.73	1.37	0.37	3.46	2.48

---

Fig. 3 Equilibrated amorphous cells of I–III (yellow shadows for M2).

3.3 Molecular packing and spectral properties of M2

3.3.1 Molecular packing. The polymorph of M2 was predicted by searching the molecular packings among seven most possible space groups (C2/c, P2₁, P2₁/C, P1̄, P2₁2₁2₁, Pbac, and Pna2₁) using the Dreiding force field and the Polymorph module. The reliability of Dreiding force field in predicting the crystal structure has been proved by many studies.^31,32 The molecular packings of M2 with the lowest energy in their respective space group are collected in Table 5. Among them, the structure with the P2₁/C symmetry has the lowest energy (−18.216 kcal per mole per cell) suggesting that the crystal structure of M2 belongs most probably to the P2₁/C space group, since, generally, the most stable polymorph possesses the lowest Gibbs free energy (or total energy at 0 K). The unit cell of the most possible packing of M2 is presented in Fig. 4 and the corresponding cell parameters are a = 8.103 Å, b = 10.634 Å, c = 12.434 Å, α = 90.000°, β = 149.945°, γ = 90.000°, Z = 4, ρ = 1.438 g cm⁻³.

Table 5 Cell parameters predicted with the Dreiding force field

Parameters	C2/c	P21	P21/C	P1̄	P212121	Pbac	Pna21
Z	8	2	4	2	4	8	4
ρ (g cm^−3)	1.436	1.420	1.438	1.399	1.431	1.440	1.416
E (kcal per mole per cell)	−17.490	−17.570	−18.216	−17.509	−17.714	−17.921	−17.491
a (Å)	9.167	4.613	8.103	9.664	8.792	8.012	9.084
b (Å)	12.757	9.201	10.634	5.247	10.087	12.231	12.755
c (Å)	9.379	7.260	12.434	6.241	6.077	10.931	4.699
α (°)	90.000	90.000	90.000	99.977	90.000	90.000	90.000
β (°)	149.945	118.258	149.945	102.802	90.000	90.000	90.000
γ (°)	90.000	90.000	90.000	111.527	90.000	90.000	90.000

---

Fig. 4 Most possible crystal structure of M2.

3.3.2 Spectral properties. IR and NMR (¹³C, ¹H, and ¹⁴N) spectra of M2 were predicted and shown in Fig. 5. The factor 0.96 has been adopted to scale the IR frequencies because the calculated frequencies by DFT are usually overestimated in comparison with those observed experimentally.³³ For the complexity of vibration modes, only some characteristic vibration modes are assigned here. The strongest peak located at 2160 cm⁻¹ belongs to the –N₃ asymmetric stretching vibration; 1560 cm⁻¹ to the –NO₂ asymmetric stretch vibration; 1290 cm⁻¹ to the –N₃ symmetric stretching vibration; and 2900 cm⁻¹ to the –CH₂ symmetric stretch vibration. For NMR spectrum, the chemical shifts of C and H are relative to those of tetramethylsilane and the chemical shifts of N are relative to those of ammonia. The peaks at 3.61, 3.89, 4.21, and 4.31 ppm correspond to the four H atoms, at 5.77, 238.89, 268.79, 413.74 ppm to N2, N4, N3, N1, respectively, at 52.15 and 78.49 ppm to C2 and C1, respectively. The atom numberings can be seen in Fig. 4.

Fig. 5 NMR (top) and IR (bottom) spectra of M2.

4. Conclusions

A series of nitro azides (M1–M8) were designed and studied using the DFT and MD methods. DFT investigation results suggest that M2–M8 have high thermal stability. MD simulations on GAP/M2 deduce that M2–M8 have good compatibility with GAP and can effectively improve the mechanical properties of GAP. That is to say, M2–M8 are promising plasticizers of GAP and worth further experimental researches. More structural information of M2 was predicted by analyzing the molecular packing and IR and NMR spectra, which is helpful to identify M2 for the experimental researchers who are interested in this new compound. The systematic investigation provides a train of thought to design and screen the plasticizers for propellants.

Acknowledgements

We gratefully thank the National Natural Science Foundation of China (no. 21403110), the Research Fund for Natural Science Foundation of Jiangsu Province (no. BK20130755), and the “Excellent Plan and Zijin Star” Research Foundation of NUST for their support to this work. Yang would like to thank the Innovation Project for Postgraduates in Universities of Jiangsu Province (grant no. KYLX_0346) for financial support.

References

1. A. Provatas, Energetic Polymers and Plasticisers for Explosive Formulations – A Review of Recent Advances, DSTO Aeronautical and Maritime Research Laboratory, Melbourne Victoria, 2000.
2. W. Leeming, E. Marshall, H. Bull, M. Rodgers and N. Paul, International Annual Conference of ICT, Karlsruhe, Germany, 1996.
3. I. Mishra, D. Rickenbaugh, C. Ashmore and A. Mehrotra-Dev, JANNAF Propulsion Meeting, 1984.
4. K. Ghosh, J. Athar, S. Pawar, B. Polke and A. Sikder, J. Energ. Mater., 2012, 30, 107–123.
5. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Suzerain, M. A. Robb, J. J. R. Cheeseman, J. A. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, Revision B.05, Wallingford, 2004.
6. A. D. Becke, J. Chem. Phys., 1992, 97, 9173–9177.
7. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
8. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
9. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–655.
10. Materials and Studio, Version 4.4, Accelry software Inc., San Diego, CA, 2008.
11. J. Yang, X. Zhang, P. Gao, X. Gong and G. Wang, RSC Adv., 2014, 4, 41934–41941.
12. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
13. K. Selim, S. Özkar and L. Yilmaz, J. Appl. Polym. Sci., 2000, 77, 538–546.
14. J. Gupta, C. Nunes, S. Vyas and S. Jonnalagadda, J. Phys. Chem. B, 2011, 115, 2014–2023.
15. L. Huynh, J. Grant, J.-C. Leroux, P. Delmas and C. Allen, Pharm. Res., 2008, 25, 147–157.
16. K. Pajula, M. Taskinen, V.-P. Lehto, J. Ketolainen and O. Korhonen, Mol. Pharmaceutics, 2010, 7, 795–804.
17. J. H. Weiner, Statistical Mechanics of Elasticity, DoverPublications. com, New York, 1983.
18. W. J. Hehre, L. Radom, P. v. R. Schleyer and J. A. Pople, Ab Initio Molecular Orbital Theory, Wiley, New York, 1986.
19. H. Dong and F. Zhou, Performance of High Energetic Explosives and Related Compounds, Science Press, Beijing, 1989.
20. J. Yang, H. Yan, X. Zhang, G. Wang and X. Gong, Struct. Chem., 2014, 25, 931–940.
21. J. Yang, H. Yan, G. Wang, X. Zhang, T. Wang and X. Gong, J. Mol. Model., 2014, 20, 2148–2160.
22. J. Yang, F. Wang, J. Zhang, G. Wang and X. Gong, J. Mol. Model., 2013, 19, 5367–5376.
23. J. N. Lee, C. Park and G. M. Whitesides, Anal. Chem., 2003, 75, 6544–6554.
24. D. J. Greenhalgh, A. C. Williams, P. Timmins and P. York, J. Pharm. Sci., 1999, 88, 1182–1190.
25. Z. B. Hu and X. X. Gan, Chin. J. Energ. Mater., 2004, 12, 62.
26. H. C. Andersen, J. Chem. Phys., 2008, 72, 2384–2393.
27. S. Pugh, London, Edinburgh, Dublin Philos. Mag. J. Sci., 1954, 45, 823–843.
28. Y. P. Ji, F. L. Gao, R. Han, B. Chen, Y. L. Wang, W. X. Liu, Y. J. Liu and Y. L. Yao, Chin. J. Energ. Mater., 2013, 21, 612–615.
29. J. L. Wang, Y. P. Ji, F. L. Gao, W. Guo and S. T. Ren, Chin. J. Energ. Mater., 2011, 19, 693–696.
30. J. Yang, X. Gong and G. Wang, Comput. Mater. Sci., 2014, 95, 129–135.
31. V. Ghule, P. Jadhav, R. Patil, S. Radhakrishnan and T. Soman, J. Phys. Chem. A, 2009, 114, 498–503.
32. F. Wang, H. Du, J. Zhang and X. Gong, J. Phys. Chem. A, 2011, 115, 11788–11795.
33. A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502–16513.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra15302b

---