Interaction, Insensitivity and Thermal Conductivity of CL-20/TNT-Based Polymer-Bonded Explosives through Molecular Dynamics Simulation

Binders mixed with explosives to form polymer-bonded explosives (PBXs) can reduce the sensitivity of the base explosive by improving interfacial interactions. The interface formed between the binder and matrix explosive also affects the thermal conductivity. Low thermal conductivity may result in localized heat concentration inside the PBXs, causing the detonation of the explosive. To investigate the binder–explosive interfacial interactions and thermal conductivity, PBXs with polyurethane as the binder and 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane/2,4,6-trinitrotoluene (CL-20/TNT) co-crystal as the matrix explosive were investigated through molecular dynamics (MD) simulations and reverse non-equilibrium molecular dynamics (rNEMD) simulation. The analysis of the pair correlation function revealed that there are hydrogen bonding interactions between Estane5703 and CL-20/TNT. The length of the trigger bonds was adopted as a theoretical criterion of sensitivity, and the effect of polymer binders on the sensibility of PBXs was correlated by analyzing the interfacial trigger bonds and internal trigger bonds of PBXs for the first time. The results indicated that the decrease in sensitivity of CL-20/TNT mainly comes from the CL-20/TNT contact with Estane5703. Therefore, the sensitivity of CL-20/TNT-based PBXs can be further reduced by increasing the contact area between CL-20/TNT and Estane5703. The thermal conductivity of PBXs composed of Estane5703 and CL-20/TNT (0 0 1), (0 1 0) and (1 0 0) crystal planes, respectively, were calculated through rNEMD simulations, and the results showed that only the addition of Estane5703 to the (1 0 0) crystal plane can improve the thermal conductivity of PBX100.


Introduction
In the field of energetic materials, unifying the contradictory properties of high energy and low sensitivity has always been a research hotspot [1]. However, most current singlecompound explosives cannot meet the requirements of high energy and low sensitivity at the same time, which seriously restricts their development and application. For example, 2,4,6,8,10,4,6,8,10, is the most famous highenergy-density compound currently in practical use, featuring high density, high detonation velocity and favorable oxygen balance [2,3]. However, CL-20 cannot meet the safety requirements of modern warfare and new weapons sufficiently because of the relatively high sensitivity [4,5]. With the development of co-crystal explosives in recent years, high energy and low sensitivity have been largely unified in them [6]. A co-crystal is a multicomponent molecular crystal with a fixed stoichiometric ratio formed by two or more neutral molecules under the action of intermolecular non-covalent bonds (hydrogen bonds, van der Waals, π-π conjugation, etc.). By forming a co-crystal, the oxygen balance of some explosives can be effectively improved, the detonation speed and pressure are increased and the safety is also guaranteed [7]. Among them, the more typical CL-20/TNT co-crystal the transition of a polymer between a glassy state and a highly elastic state. As the chain segments move when the glass transition occurs, the degree of change in the volume of the polymer will change. At the same time, due to the movement of the chain segments, the van der Waals force between the chain segments should change considerably around the glass transition temperature. The Estane5703 equilibrium structure obtained through MD simulation was subjected to 200 K to 500 K heating NPT simulation under atmospheric pressure conditions, and the simulation time was 5 ns. Thus, the curve of van der Waals energy with increasing temperature can be obtained, as shown in Figure 1a. In addition, referring to the most commonly used methods in molecular dynamics, the computation of the specific volume-temperature properties of Estane5703 is shown in Figure 1b. The Estane5703 cells were cooled from 450 K to 50 k in 50 k increments, and run equilibration for 5 ns at each temperature. All processes for volume equilibration were performed under NPT conditions. Waals energy with increasing temperature can be obtained, as shown in Figure 1a. In addition, referring to the most commonly used methods in molecular dynamics, the computation of the specific volume-temperature properties of Estane5703 is shown in Figure 1b. The Estane5703 cells were cooled from 450 K to 50 k in 50 k increments, and run equilibration for 5 ns at each temperature. All processes for volume equilibration were performed under NPT conditions. The glass transition temperatures obtained from the van der Waals-temperature curve and volume-temperature curve were 255.7 K and 257.3 K, respectively, while its experimental value was 250 K [29][30][31][32], indicating that the selected force field is applicable. The effectiveness described by van der Waals energy had a direct impact on the investigation of the interaction between the binder and co-crystal. In addition, the thermal conductivity of Estane5703 calculated under the COMPASS force field was 0.196 W·K −1 ·m −1 , which was close to the experimental value of 0.21 W·K −1 ·m −1 . The accuracy of the thermal conductivity of the binder was also directly related to the effectiveness of the thermal conductivity evaluation of PBXs. The molecular dynamics simulation results of the density, thermal conductivity and glass transition temperature of Estane5703 were highly consistent with the experimental values, as shown in Table 1, demonstrating the applicability of the simulation method and force field to Estane5703.  The glass transition temperatures obtained from the van der Waals-temperature curve and volume-temperature curve were 255.7 K and 257.3 K, respectively, while its experimental value was 250 K [29][30][31][32], indicating that the selected force field is applicable. The effectiveness described by van der Waals energy had a direct impact on the investigation of the interaction between the binder and co-crystal. In addition, the thermal conductivity of Estane5703 calculated under the COMPASS force field was 0.196 W·K −1 ·m −1 , which was close to the experimental value of 0.21 W·K −1 ·m −1 . The accuracy of the thermal conductivity of the binder was also directly related to the effectiveness of the thermal conductivity evaluation of PBXs. The molecular dynamics simulation results of the density, thermal conductivity and glass transition temperature of Estane5703 were highly consistent with the experimental values, as shown in Table 1, demonstrating the applicability of the simulation method and force field to Estane5703.

Interface Interaction
The binding energy was defined as the negative value of the intermolecular interaction energy (E inter ). The interaction energy between the explosive and the polymer binder was calculated by subtracting the individual component energy in the system from the total energy of the whole system. Thus, the binding energy can be expressed as follows [2]: where E bind and E inter represent the binding energy and interaction energy, respectively; E total represents the total energy of PBXs; E CL-20/TNT represents the energy of CL-20/TNT matrix explosives and E Estane5703 represents the energy of Estane5703. For PBXs, the binding energy indicated the strength of the interaction between the polymer binder and the base explosive. In this paper, binding energy was used to reflect the thermal stability of energetic systems. Table 2 shows the binding energy and the energy of the corresponding components between the crystal faces of three different co-crystal explosives and Estane5703. From Table 2, it can be seen that PBX100 has the largest binding energy value and PBX001 has the smallest binding energy value. Due to differences in interfacial area, E bind was used to represent the binding energy per unit interfacial area, where E bind is the E bind of the unit interfacial area. By comparing E bind , it was discovered that the thermal stability follows the order PBX100 > PBX010 > PBX001. The above showed that the interfacial interactions between Estane5703 and the CL-20/TNT crystal planes differed significantly, and the largest value of E bind and E bind indicated that the interactions between Estane5703 and the (1 0 0) crystal plane were the strongest.

Pair Correlation Function
The pair correlation function (PCF) describes the probability number density g(r) with the distance of the reference particle, representing the possibility of finding a certain particle at the distance r from the reference particle. Therefore, by analyzing the g(r) of different atom pairs at the Estane5703 and CL-20/TNT interface, an insight into the material structure could be gained through the revealed local spatial ordering. Because the interaction between particles decreased sharply when the distance between particles was longer, only the atoms closer to the interface between Estane5703 and CL-20/TNT needed to be considered. The pair correlation function of the interfacial structure atom pair of PBX001, PBX010 and PBX100 are shown in Figure 2. For convenience, the hydrogen atom, oxygen atom and nitrogen atom in CL-20/TNT are denoted as H 1 , O 1 and N 1 , respectively. in PBX100 as an example. Although there were hydrogen bond interactions in the in structure composed of the three crystal planes of Estane5703 and CL-20/TNT co-c the peaks in the hydrogen bond range varied in intensity. For PBX001, in the hyd bond range, it could be found that the g(r) value of O1-H2 was significantly less tha of H1-O2. Similarly, the g(r) value of O1-H2 was larger than that of O1-H2 in PBX01 the g(r) values of O1-H2 and H1-O2 were almost the same in PBX100. For N1−H2 in 2, the curve had no peak for hydrogen bond interactions, implying that only vdW electrostatic interactions existed between N1−H2 pairs. Non-bonding interactions between atoms include hydrogen bonds, van der Waals (vdW) and electrostatic interactions. In general, the distance range for hydrogen bonding interactions was 2.0-3.1 Å and that of the vdW and electrostatic interactions was usually 3.1-5.0 Å. When the distance between two atoms was farther than 5.0 Å, the interaction was quite weak. It can be seen from Figure 2 that in the hydrogen bond range, the PCF curves of O 1 -H 2 and H 1 -O 2 all displayed comparatively high peaks, indicating that the hydrogen bonds exist in O 1 -H 2 pairs and H 1 -O 2 pairs for PBX001, PBX010 and PBX100. Figure 3 displays a visual representation of the hydrogen bonding interactions of O 1 -H 2 in PBX100 as an example. Although there were hydrogen bond interactions in the interface structure composed of the three crystal planes of Estane5703 and CL-20/TNT co-crystal, the peaks in the hydrogen bond range varied in intensity. For PBX001, in the hydrogen bond range, it could be found that the g(r) value of O 1 -H 2 was significantly less than that of H 1 -O 2 . Similarly, the g(r) value of O 1 -H 2 was larger than that of O 1 -H 2 in PBX010, but the g(r) values of O 1 -H 2 and H 1 -O 2 were almost the same in PBX100. For N 1 −H 2 in Figure 2, the curve had no peak for hydrogen bond interactions, implying that only vdW and electrostatic int For details, the integral areas of the PCF curves were calculated. The difference in the interaction between the three atoms pairs can be clearly compared by integrating the PCF curve, as shown in Table 3. The value of the integrated area of the correlation function can reflect the magnitude of the interaction strength [35]. It can be concluded from Table 3 that the strength of the hydrogen bonds in the three interface structures is less than that of the van der Waals interaction, which indicates that there are hydrogen bond interactions in the interface structure but the van der Waals interaction is greater than the hydrogen bond. From Table 3, it can also be seen that the integral areas for the hydrogen bonds range in PBX001 were the smallest. The integration area of O 1 -H 2 was the largest in PBX010 and the integration area of H 1 -O 2 was the largest in PBX100. For distance in the hydrogen bond range, PBX010 had the smallest integral area, which also coincided with the smallest value of its binding energy. eractions existed between N 1 −H 2 pairs. For details, the integral areas of the PCF curves were calculated. The differen interaction between the three atoms pairs can be clearly compared by integrating curve, as shown in Table 3. The value of the integrated area of the correlation fun reflect the magnitude of the interaction strength [35]. It can be concluded from Ta the strength of the hydrogen bonds in the three interface structures is less than th van der Waals interaction, which indicates that there are hydrogen bond intera the interface structure but the van der Waals interaction is greater than the hydrog From Table 3, it can also be seen that the integral areas for the hydrogen bonds PBX001 were the smallest. The integration area of O1-H2 was the largest in PBX the integration area of H1-O2 was the largest in PBX100. For distance in the hydrog range, PBX010 had the smallest integral area, which also coincided with the small of its binding energy.

N-NO 2 Trigger Bond Length
Compared with C-NO 2 , N-NO 2 is a more easily broken trigger bond. Additionally, it is well known that CL-20 is more sensitive than TNT, and the CL-20 component in CL-20/TNT decomposes first in detonation. Therefore, the N-NO 2 bond in CL-20 was chosen as the trigger bond in this work. The maximum trigger bond length (L max ) and the average trigger bond length (L ave ) were used to characterize the sensitivity of CL-20/TNT co-crystal-based PBXs. In addition, since the interaction of the polymer binder at the CL-20/TNT interface will cause the sensitivity of CL-20/TNT at the interface to be different from that of CL-20/TNT that is not in contact with the binder, the trigger bond of PBXs was divided into both interface and internal parts for analysis. Table 4 lists the maximum bond length L max and the average bond length L ave of N-NO 2 . It can been seen from Table 4 that the bond lengths of N-NO 2 are all smaller than those of CL-20/TNT, for PBX001, PBX010 and PBX100, which indicates that the addition of the binder Estane5703 can make the sensitivity decrease. The L ave at the interface of PBX001 and the L max of PBX001 are essentially unchanged compared to CL-20/TNT, which shows that placing Estane5703 on the (0 0 1) crystal plane of CL-20/TNT has less effect on the sensitivity. After that, the L max of PBX010 and PBX100 are smaller than the L max of CL-20/TNT, and the L ave of PBX010 and PBX100 are markedly smaller than the L ave of CL-20/TNT, which indicates that the addition of the binder Estane5703 made the sensitivity decrease of PBX010 and PBX100 more obvious. In addition, it can be found that the L max and L ave at the interface are smaller than the internal L max and L ave , which indicates that Estane5703 decreased the sensitivity of CL-20/TNT, mainly due to the decreased sensitivity of CL-20/TNT in contact with Estane5703. At the same time, this also shows that the sensitivity of CL-20/TNT can be reduced by increasing the contact area between the polymer binder Estane5703 and CL-20/TNT.

Thermal Conductivity
For energetic materials, if their thermal conductivity is low, the local temperature may rise quickly because of the localized heat concentration inside the PBXs, which may detonate the energetic materials. Figure 4 shows the thermal conductivity along the c-axis of PBXs with three different crystal planes and CL-20/TNT with three different crystal planes. Due to the different atomic arrangements in the three crystal planes (0 0 1), (0 1 0) and (1 0 0), the thermal conductivity of the CL-20/TNT co-crystal was also different. Among them, the thermal conductivity of the (0 0 1) crystal plane along the c-axis was the best. The thermal conductivity of PBX001 and PBX010 along the c-axis was lower than their corresponding CL-20/TNT co-crystals. The decreased thermal conductivity of PBX001 formed by adding the adhesive Estane5703 in (0 0 1) might have increased the sensitivity of CL-20/TNT, which is consistent with the conclusion from the PBX001 trigger bonds. The thermal conductivity of PBX010 formed by adding the adhesive Estane5703 in the (0 1 0) plane was almost unchanged compared with that of CL-20/TNT. The thermal conductivity of PBX100 was greater than CL-20/TNT, indicating that adding Estane5703 in the (1 0 0) crystal plane can improve the thermal conductivity of PBX100.

Amorphous Polymer Chain
The molecular structure of polyurethane (Estane5703) is shown in Figure 5a. This polymer is a copolymer composed of hard segments and soft segments. In this paper, the number of repeating units of hard segments was m = 1, the number of repeating units of soft segments was n = 4 and the end groups were saturated with H and CH3 to form the initial stage of Estane5703. The Monte Carlo method was used to generate twist angle in the chain to obtain the random conformation of the Estane5703 molecular chain. Then, six Estane5703 molecular chains were put into a periodic box with a low density, and a density-convergent Estane5703 amorphous model was obtained by using the high-low pressure dynamic simulation method [36]. This was the initial model of the polymer binder in this paper.

Amorphous Polymer Chain
The molecular structure of polyurethane (Estane5703) is shown in Figure 5a. This polymer is a copolymer composed of hard segments and soft segments. In this paper, the number of repeating units of hard segments was m = 1, the number of repeating units of soft segments was n = 4 and the end groups were saturated with H and CH 3 to form the initial stage of Estane5703. The Monte Carlo method was used to generate twist angle in the chain to obtain the random conformation of the Estane5703 molecular chain. Then, six Estane5703 molecular chains were put into a periodic box with a low density, and a density-convergent Estane5703 amorphous model was obtained by using the high-low pressure dynamic simulation method [36]. This was the initial model of the polymer binder in this paper.

Amorphous Polymer Chain
The molecular structure of polyurethane (Estane5703) is shown in Figure 5a. This polymer is a copolymer composed of hard segments and soft segments. In this paper, the number of repeating units of hard segments was m = 1, the number of repeating units of soft segments was n = 4 and the end groups were saturated with H and CH3 to form the initial stage of Estane5703. The Monte Carlo method was used to generate twist angle in the chain to obtain the random conformation of the Estane5703 molecular chain. Then, six Estane5703 molecular chains were put into a periodic box with a low density, and a density-convergent Estane5703 amorphous model was obtained by using the high-low pressure dynamic simulation method [36]. This was the initial model of the polymer binder in this paper.

Models of CL-20/TNT and PBXs
Based on the data of CL-20/TNT co-crystal obtained using X-ray analysis [37], a (2 × 1 × 1) supercell model was built. Since the molecular arrangements of CL-20/TNT in the three crystal plane directions (0 0 1), (0 1 0) and (1 0 0) were quite different, these three crystal planes were selected to study the interfacial interaction and thermal conductivity of CL-20/TNT-based PBX. First, the CL-20/TNT supercell model was cut along the three crystal plane directions (0 0 1), (0 1 0) and (1 0 0), respectively. The c lattice heights of the obtained (0 0 1), (0 1 0) and (1 0 0)  Taking the (0 0 1) crystal plane model as an example, the steps to establish a PBX model for the corresponding crystal plane were as follows. An Estane5703 molecular chain was chosen from the initial model of the polymer binder and put in a period box with the same size as the (0 0 1) crystal plane model. Then, keeping the surface (a × b) unchanged, the height of the c-axis was adjusted gradually until the theoretical density of Estane5703 was reached; each adjustment needed the molecular dynamics simulation running to equilibrium state. After that, the (0 0 1) crystal plane model was merged at the two ends of the Estane5703 period box perpendicular to the c-axis, respectively. Therefore, the PBX001 structure was obtained, containing 1994 atoms, as shown in Figure 6. The (0 1 0) and (1 0 0) crystal plane structures obtained through the same method are denoted as PBX010 and PBX100, respectively. In this way, the weight percentages of the binders in the PBXs are all about 5.2%, within the general control range of polymer binder content in explosives [21].

Models of CL-20/TNT and PBXs
Based on the data of CL-20/TNT co-crystal obtained using X-ray analysis [37], a (2 × 1 × 1) supercell model was built. Since the molecular arrangements of CL-20/TNT in the three crystal plane directions (0 0 1), (0 1 0) and (1 0 0) were quite different, these three crystal planes were selected to study the interfacial interaction and thermal conductivity of CL-20/TNT-based PBX. First, the CL-20/TNT supercell model was cut along the three crystal plane directions (0 0 1), (0 1 0) and (1 0 0), respectively. The c lattice heights of the obtained (0 0 1), (0 1 0) and (1 0 0)  Taking the (0 0 1) crystal plane model as an example, the steps to establish a PBX model for the corresponding crystal plane were as follows. An Estane5703 molecular chain was chosen from the initial model of the polymer binder and put in a period box with the same size as the (0 0 1) crystal plane model. Then, keeping the surface (a × b) unchanged, the height of the c-axis was adjusted gradually until the theoretical density of Estane5703 was reached; each adjustment needed the molecular dynamics simulation running to equilibrium state. After that, the (0 0 1) crystal plane model was merged at the two ends of the Estane5703 period box perpendicular to the c-axis, respectively. Therefore, the PBX001 structure was obtained, containing 1994 atoms, as shown in Figure 6. The (0 1 0) and (1 0 0) crystal plane structures obtained through the same method are denoted as PBX010 and PBX100, respectively. In this way, the weight percentages of the binders in the PBXs are all about 5.2%, within the general control range of polymer binder content in explosives [21]. The above models for CL-20/TNT co-crystal and CL-20/TNT-based PBXs were allowed to evolve dynamically in isothermal-isobaric (NPT) ensemble with Andersen temperature control using the stochastic collision method and Parrinello-Rahman [38] pressure control, fully relaxing all cell parameters at atmospheric pressure. The van der Waals The above models for CL-20/TNT co-crystal and CL-20/TNT-based PBXs were allowed to evolve dynamically in isothermal-isobaric (NPT) ensemble with Andersen temperature control using the stochastic collision method and Parrinello-Rahman [38] pressure control, fully relaxing all cell parameters at atmospheric pressure. The van der Waals interactions were truncated at 12.5 Å with long-range tail correction, and the electrostatic interactions were calculated via the standard Ewald summation [39]. The equations of motion were integrated with a step of 1 fs. An equilibration run was performed for 5 ns. After the equilibration run, production runs of 1 ns were performed, during which data were collected with 10 fs sampling intervals for analysis. In this paper, all simulations were conducted utilizing a COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) force field [40], which was suitable for MD simulation of the condensed phase, especially for nitramine explosives [16,[41][42][43][44].

Models and Methods for Thermal Conductivity
The thermal conductivity of PBXs was calculated through a reverse non-equilibrium MD (rNEMD) [45,46], whose advantage was fast convergence of the temperature gradient in a non-equilibrium stable system. For the rNEMD, we used the Muller-Plathe algorithm to exchange kinetic energy between atoms in the hot region and cold region.
The first step was to establish an appropriate thermal conductivity model, and then apply heat flux to the model in one direction to generate a temperature gradient in the same direction. The thermal conductivity can be further obtained by calculating the heat flux along the temperature gradient. Regarding the thermal conductivity model, for the PBX model, in order to ensure that the heat flux completely passed through the interface between CL-20/TNT and Estane5703, two identical PBX models were merged into their thermal conductivity model along the pseudo-flow direction of heat, that is, the c-axis direction. Figure 7 displays the thermal conductivity model of PBX001 as an example where the heat flux is along the c-axis from the hot region to the cold region. interactions were truncated at 12.5 Å with long-range tail correction, and the electrostatic interactions were calculated via the standard Ewald summation [39]. The equations o motion were integrated with a step of 1 fs. An equilibration run was performed for 5 ns After the equilibration run, production runs of 1 ns were performed, during which data were collected with 10 fs sampling intervals for analysis. In this paper, all simulations were conducted utilizing a COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) force field [40], which was suitable for MD simulation o the condensed phase, especially for nitramine explosives [16,[41][42][43][44].

Models and Methods for Thermal Conductivity
The thermal conductivity of PBXs was calculated through a reverse non-equilibrium MD (rNEMD) [45,46], whose advantage was fast convergence of the temperature gradien in a non-equilibrium stable system. For the rNEMD, we used the Muller-Plathe algorithm to exchange kinetic energy between atoms in the hot region and cold region.
The first step was to establish an appropriate thermal conductivity model, and then apply heat flux to the model in one direction to generate a temperature gradient in the same direction. The thermal conductivity can be further obtained by calculating the hea flux along the temperature gradient. Regarding the thermal conductivity model, for the PBX model, in order to ensure that the heat flux completely passed through the interface between CL-20/TNT and Estane5703, two identical PBX models were merged into their thermal conductivity model along the pseudo-flow direction of heat, that is, the c-axis direction. Figure 7 displays the thermal conductivity model of PBX001 as an example where the heat flux is along the c-axis from the hot region to the cold region. The simulated CL-20/TNT or PBX system was divided into 50 equal slices along the direction of heat flux (c-axis in this paper), where the middle slice is the heat source layer (hot region) and the first slices of the system are the heat sink layer (cold region). The coldest atom in the hot region and the hottest atom in the cold region were paired up and their velocities were exchanged [47]. The exchange could effectively swap kinetic energies of atoms, assuming their masses were the same. If the masses were different, a velocity swap was performed relative to the motion of the center of mass of the atoms to conserve kinetic energy. Thus, the energy could be transferred from the hot region to the cold re gion, and a temperature gradient could be induced in the simulated box. The heat flux can be obtained by calculating the exchanged amount of kinetic energy according to the fol lowing formula: where Nswap is the number of swaps to perform for energy conversion; tswap is the intervals to perform kinetic energy exchange; m is the mass of the exchanged atoms and vh and v are the speeds of the hottest and coldest atoms, respectively. Then, the thermal conductiv ity can be calculated using the Fourier formula: The simulated CL-20/TNT or PBX system was divided into 50 equal slices along the direction of heat flux (c-axis in this paper), where the middle slice is the heat source layer (hot region) and the first slices of the system are the heat sink layer (cold region). The coldest atom in the hot region and the hottest atom in the cold region were paired up and their velocities were exchanged [47]. The exchange could effectively swap kinetic energies of atoms, assuming their masses were the same. If the masses were different, a velocity swap was performed relative to the motion of the center of mass of the atoms to conserve kinetic energy. Thus, the energy could be transferred from the hot region to the cold region, and a temperature gradient could be induced in the simulated box. The heat flux can be obtained by calculating the exchanged amount of kinetic energy according to the following formula: where N swap is the number of swaps to perform for energy conversion; t swap is the intervals to perform kinetic energy exchange; m is the mass of the exchanged atoms and v h and v c are the speeds of the hottest and coldest atoms, respectively. Then, the thermal conductivity can be calculated using the Fourier formula: in which A is the cross-sectional area in the direction of heat flux and ∂T⁄∂z is the derivative of the system temperature with respect to the direction of heat flux. Due to the exchange of the velocity of the related atoms in the heat source layer and the heat sink layer, the temperature gradient in the system could be generated. The temperature gradient of the system was obtained by calculating the temperature of 50 equal parts of the uniform lamella, as shown in Figure 8. The temperature of each layer can be calculated using the following formula: where N is the number of atoms in the slice (on average, each slice contains 80 atoms); k B is Boltzmann's constant; p j is the momentum of the atom and m is the mass of the atom.
in which A is the cross-sectional area in the direction of heat flux and ∂T⁄∂z is the derivative of the system temperature with respect to the direction of heat flux. Due to the exchange of the velocity of the related atoms in the heat source layer and the heat sink layer, the temperature gradient in the system could be generated. The temperature gradient of the system was obtained by calculating the temperature of 50 equal parts of the uniform lamella, as shown in Figure 8. The temperature of each layer can be calculated using the following formula: where N is the number of atoms in the slice (on average, each slice contains 80 atoms); kB is Boltzmann's constant; pj is the momentum of the atom and m is the mass of the atom. In this thermal conductivity simulation, the equations of motion were integrated with a step of 1 fs, and the temperature and pressure were controlled via the Nose-Hoover thermostat and barostat [48,49]. Van der Waals interactions were truncated atom pairs with interatomic distances greater than 12.5 Å, coupled with long-range tail correction. And electrostatic interactions were processed using standard Ewald summation. The simulation systems were first equilibrated under the NPT ensemble at 300 K and atmospheric pressure run for 2 ns. Then, the exchange of kinetic energies was performed for 2 ns. The temperature profile was computed by averaging events in a time interval of 20 ps. Five independent equilibration and NEMD simulations were performed for each crystal orientation. These computations were all carried out using the software LAMMPS (5 June 2019) (Large-scale Atomic/Molecular Massively Parallel Simulator) [50].

Conclusions
In this paper, the interface interaction and sensitivity of PBXs formed by adding the polymer binder Estane5703 in the (0 0 1), (0 1 0) and (1 0 0) crystal planes of CL-20/TNT were studied through molecular dynamics simulations. The simulations involved binding energy calculation, PCF analysis, bond length of the N-NO2 trigger bond and thermal conductivity.
From the calculated binding energies, the values of Ebind and Ebind' between Es-tane5703 and CL-20/TNT for PBX100 were the largest. It was found that PBX100 has the best thermal stability, which means that Estane5703 combined has the best thermal stability in the (1 0 0) crystal plane. Additionally, the PCF analysis of atom pairs in the interfacial structure indicated that there exists hydrogen bond between CL-20/TNT and Estane5703 In this thermal conductivity simulation, the equations of motion were integrated with a step of 1 fs, and the temperature and pressure were controlled via the Nose-Hoover thermostat and barostat [48,49]. Van der Waals interactions were truncated atom pairs with interatomic distances greater than 12.5 Å, coupled with long-range tail correction. And electrostatic interactions were processed using standard Ewald summation. The simulation systems were first equilibrated under the NPT ensemble at 300 K and atmospheric pressure run for 2 ns. Then, the exchange of kinetic energies was performed for 2 ns. The temperature profile was computed by averaging events in a time interval of 20 ps. Five independent equilibration and NEMD simulations were performed for each crystal orientation. These computations were all carried out using the software LAMMPS (5 June 2019) (Large-scale Atomic/Molecular Massively Parallel Simulator) [50].

Conclusions
In this paper, the interface interaction and sensitivity of PBXs formed by adding the polymer binder Estane5703 in the (0 0 1), (0 1 0) and (1 0 0) crystal planes of CL-20/TNT were studied through molecular dynamics simulations. The simulations involved binding energy calculation, PCF analysis, bond length of the N-NO 2 trigger bond and thermal conductivity.
From the calculated binding energies, the values of E bind and E bind' between Es-tane5703 and CL-20/TNT for PBX100 were the largest. It was found that PBX100 has the best thermal stability, which means that Estane5703 combined has the best thermal stability in the (1 0 0) crystal plane. Additionally, the PCF analysis of atom pairs in the interfacial structure indicated that there exists hydrogen bond between CL-20/TNT and Estane5703 in all three crystal planes. The hydrogen bonds mainly came from H 1 −O 2 and O 1 −H 2 , but there were far fewer hydrogen bonds present in the N 1 −H 2 atomic pairs. The PCF also indicated that the number of strong vdW and electrostatic interacting atomic pairs is much greater than the number of hydrogen bond interacting atomic pairs. By analyzing the length of the N-NO 2 trigger bond, it was found that the L max and L ave at the interface are both smaller than the internal L max and L ave , indicating that the main reason for the Estane5703 reducing the sensitivity of CL-20/TNT is that it reduces the contact with Estane5703. Therefore, the sensitivity can be reduced by increasing the contact area of the binder Estane5703 with the CL-20/TNT co-crystal. The thermal conductivity of PBXs showed that merging Estane5703 in the (0 0 1) and (0 1 0) crystal plane can lead to a decrease in the thermal conductivity of CL-20/TNT in the corresponding crystal planes, and only merging Estane5703 in the (1 0 0) crystal plane can improve the thermal conductivity of CL-20/TNT in the corresponding crystal planes. This research work is helpful to deeply understand the interaction between the binder and the explosive in polymer-bonded explosives, and to further reveal the mechanism of the reduction in the sensitivity of the explosive in polymer-bonded explosives.