Molecular dynamic simulation of performance of modified BAMO/AMMO copolymers and their effects on mechanical properties of energetic materials

Based on molecular dynamic method, densities, mechanical behavior and mechanical performance of P(BAMO/ AMMO) (Polymer 1) and two novel modified P(BAMO/AMMO) (Polymer 2: containing amino group, Polymer 3: containing nitro group), and their effects on mechanical properties of four energetic materials are investigated, the main results are as follow: Polymer 2 (1.235 g/cm3, 240 ± 5 K) and Polymer 3: 1.281 g/cm3, 181 ± 3 K) possess higher densities and lower glass transition temperatures than Polymer 1 (1.229 g/cm3, 247 ± 4 K). The modification makes Polymer 1 difficult to expand, improves its mechanical properties, but has few effect on its diffusion coefficient at same temperature and state. In addition, three binders are compatible with TNT, HMX and CL-20, and may react with DNTF. All polymers particularly improve rigidity of four energetic materials, and enhance their ductility except Polymer 2 on TNT. The ability of Polymer 2 and Polymer 3 improving rigidity (except Polymer 3 on HMX) and ductility of TNT and HMX is inferior to that of Polymer 1, but it is contrary for CL-20 and DNTF (except Polymer 2 on rigidity of DNTF). Moreover, Polymer 2-based interfacial crystals exhibit higher rigidity than Polymer 3-based interfacial crystals.

Constructions of models and simulation methods for polymers. The Fig. 2). After structure optimization, the cell of P(BAMO/ AMMO) and modified P(BAMO/AMMO) were relaxed by 50 ps constant particle number, pressure, and temperature (NPT) ensemble, 50 ps constant particle number, volume, and temperature (NVT) ensemble, and annealing simulation (temperature from 300 to 500 K). Then, the structure of minimum energy was used to perform NPT-MD simulation. The total calculated time was 1 ns, and the step time was 1 fs. The given temperature and pressure were 298 K and 100 kPa, respectively. Anderson method and Berendsen method were used to control temperature and pressure, respectively 29,30 . Initial velocity was sampled by Maxwell distribution, and velocity Verlet arithmetic was utilized 31 . Van der Waals force and Electrostatic interaction were calculated by atom-based method and Ewald method, respectively. The balanced structure obtained from 1 ns NPT-MD simulation was repeatedly performed 500 ps NPT-MD simulation, where the temperature was successively set from 513 K to 33 K, and the interval of temperature was 20 K. After the simulation finished, the balanced structures were performed 100 ps NVT-MD simulation for every temperature. Then, every system in different temperatures was relaxed again. Ultimately, the relaxed structures in different temperatures were performed 1 ns NPT-MD simulation in their corresponding temperatures. The final results were applied to calculate volume (V), mean square displacement (MSD), non-bond energy (NBE) and mechanical properties.

Constructions of models and simulation methods for PBXs. The Morphology Module was applied
to evaluate the growth face of TNT, HMX, DNTF and CL-20 crystals in vacuum, and the original cell of all crystals are obtained from Cambridge Crystallographic Data Centre (CCDC). The growth faces of four crystals are presented in Fig. 3, and their first two growth faces and corresponding percentage are listed in Table 1.  www.nature.com/scientificreports/ The construction of polymer bonded explosives models is as follow: taking CL-20 and Polymer 1 as an instance (Fig. 4): firstly, the 4*3*3 supercell of CL-20 was built, and was cleaved along the main growth face (0 1 1) ; secondly, the surface (0 1 1) was built as a crystal without vacuum; thirdly, the Polymer 1 chain was constructed as amorphous cell with same size as underside of crustal (0 1 1); finally, the crystal (0 1 1) and amorphous cell of Polymer 1 were built layer and their interfacial crystal was obtained. Repeating that method, the other interfacial crystals were acquired. It is worth noting that all steps were needed to be performed structure optimization to make energies of systems exhibit minimums after every step, because the system were deprived initial structure and possessed higher energy than its most stable structure after every step. As listed in Table 2, it is the number of total atoms of interfacial crystals and the mass ratio of polymer binder in every system. And the mass ratio of polymer binders in the systems was about 7.4%-8.4%.
The optimized structures of interfacial crystal and four supercell were performed 1 ns NPT-MD simulations, and the calculated method and set was same as chapter "2.2 Constructions of models and simulation methods for polymers". The final 300 ps were used to predicted mechanical properties.  where σ i is the stress tensor(GPa), ε j is the strain tensor(GPa), and C ij is the 6 × 6 stiffness matrix of elastic constants. When the material is regarded as an isotropic material, the stiffness matrix of the stress-strain behavior can be fully expressed by Lamé coefficients (λ and μ), as follow Scheme 1. Then, the Young's modulus (E, GPa), Bulk modulus (K, GPa), Shear modulus (G, GPa) and Poisson's ratio (γ) can be described by Lamé coefficients as follows 33-35 : Results and discussions Density and glass transition temperature of polymers. The densities of Polymer 1, Polymer 2 and Polymer 3 are 1.229 g/cm 3 , 1. 235 g/cm 3 and 1.281 g/cm 3 , respectively. Therein, the density of Polymer 1 is close to its the reference value (1.25 g/cm 3 ) and the relatively error is 1.7%, which indicates that the method applied in this work is reliable 14 . Otherwise, it is found that the -NH 2 little improve the density of Polymer 1, and the -NO 2 well increase the of Polymer 1.
The glass transition temperature (T g ) is probable one of the most important properties of polymers because it determines the processing and working temperature range. A large number of methods have been applied to evaluate the T g of polymer based on molecular dynamic method 20,21 . For each temperature, the volume of the polymer model was obtained from the whole duration of MD simulation and averaged. It is worth noting that every 1 ns NPT-MD simulation in different temperature is performed based on 500 ps NPT-MD and 100 ps NVT-MD relaxation in corresponding temperature, therefore, the errors of averaged volume at each temperature can be negligible. The volume of three polymers at different temperature is listed in Table 3, and the curves of their volume versus temperature are presented in Fig. 5. It is easy to see, there are obvious discontinuities in the slope of the curves, which means the polymers exist the glass transition, namely, the polymer transforms from rubbery state to glassy state. The T g is predicted by performing the segmental linear regression of the data. Eventually, the T g of Polymer 1 Polymer 2 and Polymer 3 are respectively 243.7 K, 235.9 K and 181.1 K. The T g of Polymer 1 is close to its literature value (244.7 K), and the relative error is 0.4%, which also shows that this work is reliable 16 . Meanwhile, it is obvious that the introduced functional groups all make the T g of P(BAMO-AMMO) decrease, and nitro group particularly reduces the T g .
Otherwise, the volumetric coefficient of thermal expansion (α) of three polymer binders can be also evaluated by the curves of volume vs temperature, and the α is defined by 36 : where V 0 is the original volume of equilibrium system before cooling.
The (∂V/∂T) values of three systems in rubbery state and glassy state and α values in three temperatures are shown in Table 4. When the polymers exhibit rubbery state in same temperature, it is easy to find the α values of (1) www.nature.com/scientificreports/ Polymer 3 and Polymer 2 are close to each other and are lower than that of Polymer 1. It indicates that Polymer 1 is easiest to expand at room and high temperature, which may be easy to introduce damage of the propellants or PBXs systems. In glassy state, the Polymer 1 is easiest to expand and Polymer 2 is much difficult to expand. It means that Polymer 2 and Polymer 3 are more difficult to expand than Polymer 1 in the same state and temperature. Otherwise, when temperature is 181.1 K < T < 235.9 K, Polymer 3 significantly expand than Polymer 1 and Polymer 2, when temperature is 235.9 K < T < 243.7 K, Polymer 3 and Polymer 2 particularly expand than Polymer 1. Because the polymers continue to move in system and its movement is constrained by space and temperature, the T g can also be predicted by MSD vs temperature curve. As is shown in Fig. 6, there are the MSD vs temperature curves of P(BAMO/AMMO) and modified P(BAMO/AMMO). Meanwhile, the specific values of their MSD are listed in Table 5. Ultimately, the T g values of Polymer 1, Polymer 2 and Polymer 3 evaluated by the fitted curves are 245.4 K, 238.9 K and 182.6 K, respectively. The acquired variation trend is agreed with that obtained by volume-temperature data, namely, nitro group significantly reduces the T g of P(BAMO/AMMO). In addition, it is found that the T g values of three polymers obtained by MSD-temperature data are all slightly higher than those acquired by volume-temperature data.
The migration of plasticizer is related to its diffusion coefficient (D), which depends on the MSD and time (t) of MD simulation. Therefore, D is described in Eq. (6) based on the relation of Einstein 37,38 : where r(t) is coordinate of plasticize in t, and r(0) is the original coordinate.
The relationship of MSD is calculated by Eq. (7) 38 : The D is ultimately estimated by Eqs. (6) and (7) 38 : The migration of three polymers almost in low temperature (233 K), room temperature (293 K) and high temperature (333 K) are compared. The Fig. 7 plots the MSD-t curves and their fitting curves of Polymer 3 in   Non-bond energy (NBE) also possesses a discontinuity near the glass transition temperature 21 . Therefore, the T g can be estimated by fitting NBE-temperature curves (as shown in Fig. 8). And the related data are presented  In a conclusion, based on calculating V, MSD and NBE, the T g values of three polymer are respectively 247 ± 4 K, 240 ± 5 K and 181 ± 3 K by assessing the range of every kind of T g . The -NH 2 and -NO 2 group all decrease the glass transition temperature of P(BAMO/AMMO), and the reducing effect of -NO 2 group is better.     Table 8. It is easy seeing that the E and G values of Polymer 1 have a little diminution when introducing -NH 2 and -NO 2 group, however, its K value slightly increase.  www.nature.com/scientificreports/ Therefore, there is a conclusion that the Polymer 2 and Polymer 3 become more plastic than Polymer 1, namely, the former is more effortless to deformed, which may make the propellants or PBXs systems easier to process. Ratio of K/G can be also used to evaluate the ductility of materials and with higher value expected. Therefore, the Polymer 2 and Polymer 3 possess better ductility than the Polymer 1 by comparing their K/G values. Miscibility means whether the polymer and energetic materials can mix in all proportions and ultimately form a homogeneous state. And cohesive energy density and solubility parameter are important index to estimate the miscibility of two components, which can be calculated by MD methods. If the solubility parameters of two components are more close, their compatibility are better, namely, they are miscible. As presented in Table 9, there are the cohesive energy density and solubility parameter of P(BAMO/AMMO), modified P(BAMO/AMMO) and four energetic materials. The cohesive energy density and solubility parameters of three polymer are close to each other despite of introducing functional groups, which indicates the functional groups have little influence on the solubility parameter of P(BAMO/AMMO). Otherwise, the solubility parameters of four energetic materials and three binders are 28.00-31.70 (J/cm 3 ) 1/2 and 19.29-19.92 (J/cm 3 ) 1/2 , respectively. It is obtained the number of their difference value (|Δδ|) is 8.08 (J/cm 3 ) 1/2 <|Δδ|< 12.41 (J/cm 3 ) 1/2 . Therein, the |Δδ| of DNTF is largest than the other when corresponding to same polymer. Otherwise, previous reports 16 , as listed in Table 10, proved that Polymer 1 was compatible with HMX and CL-20, but medium reacted with DNTF. Therefore, it can be concluded that Polymer 2 and Polymer 3 are compatible with TNT, HMX and CL-20, but may be incompatible with DNTF as same as Polymer 1.

Mechanical properties of different PBXs. Energetic materials used as oxidation are important solid
components in solid rocket propellant and polymer bonded explosives, whose mechanical properties significantly depend on energetic materials. Meanwhile, binders which are applied to bond solid components also play an important role in improving the mechanical properties of the propellants of PBXs system. The mechanical parameters of different energetic material/binder systems calculated by MD method are presented in Table 11-14. Moreover, their mechanical modulus are plotted in Fig. 9.
Binder/TNT interfacial crystals. The percentage of (0 0 2) in total growth face of TNT is regarded as 38.6% in this paper, which consists of (0 0 2) and its symmetry plane (0 0 − 2) possessing same proportion. The (2 0 0) is in the same case with (0 0 2). Therefore, (0 0 2) and (2 0 0) are regarded as the first two growth faces of TNT. As listed in Table 11, there are mechanical parameters of polymer binder/TNT interfacial crystals and pure TNT. It can be found that the mechanical modulus of all interfacial crystals are significantly lower than those of pure TNT, which suggests that three polymer binders can well decrease the rigidity of TNT. Meanwhile, when the binder is Polymer 1, the Young's modulus (E), bulk modulus (K) and shear modulus (G) values of corresponding interfacial crystals are smallest in spite of they obtained from TNT (0 0 2) or TNT (2 0 0) (except the K value of Polymer 1/TNT (0 0 2) crystal), which indicates that modified P(BAMO/AMMO) do not well enhance the ductility of TNT compared with P(BAMO/ AMMO). Otherwise, the mechanical modulus of Polymer 3/TNT  Binder/HMX interfacial crystals. Compared with pure HMX crystal, the mechanical modulus of binder/HMX interfacial crystals all significantly reduce, which indicates that binders all effectively improve the isotropy of HMX crystal (as shown in Table 12). The E, K and G values of polymer binder/HMX (0 1 1) (the first growth face) crystal are particularly lower than those of polymer binder/HMX (1 1 0) (the second growth face) crystal when binder are Polymer 1 and Polymer 3. However, while binder is Polymer 2, the mechanical modulus of the E and G values of polymer binder/HMX (0 1 1) crystal are even little larger than those of binder/HMX (1 1 0) crystal, which notes the effects of Polymer 2 on the two main growth faces of HMX have no obvious distinction. Otherwise, as presented in Fig. 8b, by comparing the elastic constants of binder/HMX interfacial crystal, when interface is (0 1 1) face of HMX crystal, the binders increase the plasticity of interfacial crystals in the order as follow: Polymer 3 > Polymer 1 > Polymer 2. Nevertheless, while being (1 1 0) face, the order is Polymer Table 11. The mechanical parameters of polymer/TNT interfacial crystal.    Binder/DNTF interfacial crystals. The data of mechanical parameters of binder/DNTF interfacial crystals and pure DNTF are listed in Table 13. It is easy to see that three polymer binders also effectively enhance the isotropy of DNTF crystal. As the first growth face, the mechanical modulus of interfacial crystals corresponding to DNTF (0 1 1) face are lower than those corresponding to DNTF ( Binder/CL-20 interfacial crystals. It can be found in Table 14 1 1) (1 0 1) (0 1 1) (1 0 1) (0 1 1) (1 0 1

Conclusion
Based on molecular dynamic method, two novel modified P(BAMO/AMMO) binders are designed, and their densities, glass transition temperatures (T g ), volumetric coefficients of thermal expansion and mechanical properties are compared with P(BAMO/AMMO). Meanwhile, the effects of three polymer binders on mechanical properties of representative energetic materials (TNT, HMX, CL-20 and DNTF) of three generations are investigated, and the main results are as follow: 1. The densities of Polymer 1, Polymer 2 and Polymer 3 are respectively 1.229 g/cm 3 , 1.235 g/cm 3 and 1.281 g/ cm 3 , which shows introducing functional groups all increase the density of P(BAMO/AMMO). The T g values of Polymer 1, Polymer 2 and Polymer 3 are respectively 247 ± 4 K, 240 ± 5 K and 181 ± 3 by calculating their V, MSD and NBE in different temperature, which suggests the the reducing effect of -NO 2 group is better. 2. Except temperature is greater than T g of Polymer 3 and less than T g of Polymer 1, the Polymer 1 is much easiest to expand in whether rubbery state or glassy state, which may restrict its application. The introduced groups has few effect on diffusion coefficient of Polymer 1. The E and G values of modified P(BAMO/AMMO) are all lower than those of P(BAMO/ AMMO), despite it is contrary for the K value. It may be concluded that modified P(BAMO/AMMO) possess better mechanical properties . 3. The solubility parameters of three polymer binders approximate to each other, and their difference value (|Δδ|) with four energetic materials is 8.08 (J/cm 3 ) 1/2 <|Δδ|< 12.41 (J/cm 3 ) 1/2 , there is a conclusion that three polymers are compatible with TNT, HMX and CL-20, but Polymer 2 and Polymer 3 may react with DNTF as same as Polymer 1. 4. Comparing the mechanical properties of interfacial crystal with corresponding crystal, it can be found that the effects of Polymer 2 and Polymer 3 reducing the rigidity of TNT are inferior to that of Polymer 1, but it is opposite for CL-20. Otherwise, there is only that the ability of Polymer 3 improving the rigidity of HMX and DNTF is superior to that of Polymer 1. Meanwhile, it is worth noting that the rigidity of Polymer 2-based systems are all higher than those of Polymer 3-based systems for four energetic material. 5. All polymers can improve ductility of four energetic materials (except Polymer 2 on TNT and Polymer 1 on CL-20). For TNT and HMX, the ability of Polymer 1 improving their ductility are superior to Polymer 2 and 3. However, it is contrary for DNTF and CL-20.