ReaxFF Molecular Dynamics Simulation of Hydrostatic and Uniaxial Compression of Nitrate Energetic Materials

The physical and chemical properties of typical nitrate energetic materials under hydrostatic compression and uniaxial compression were studied using the ReaxFF/lg force field combined with the molecular dynamics simulation method. Under hydrostatic compression, the P–V curve and the bulk modulus (B0) obtained using the VFRS equation of state show that the compressibility of the three crystals is nitroglycerine (NG) > erythritol tetranitrate (ETN) > 2,3-bis-hydroxymethyl-2,3-dinitro-1,4-butanediol tetranitrate (NEST-1). The a- and c-axis of ETN are easy to compress under the action of hydrostatic pressure, but the b-axis is not easy to compress. The b-axis of NEST-1 is the most compressible, while the a- and c-axis can be compressed slightly when the initial pressure increases and then remains unchanged afterward. The a-, b-, and c-axes of NG all have similar compressibilities. By analyzing the change trend of the main bond lengths of the crystals, it can be seen that the most stable of the three crystals is the N–O bond and the largest change is in the O–NO2 bond. The stability of the C–O bond shows that the NO3 produced by nitrates is not from the C–O bond fracture. Under uniaxial compression, the stress tensor component, the average principal stress, and the hydrostatic pressure have similar trends and amplitudes, indicating that the anisotropy behaviors of the three crystals ETN, NEST-1, and NG are weak. There is no significant correlation between maximum shear stress and sensitivity. The maximum shear stresses τxyand τyz of the ETN in the [010] direction are 1.5 GPa higher than τxz. However, the maximum shear stress of NG shows irregularity in different compression directions, indicating that there is no obvious correlation between the maximum shear stress and sensitivity.


INTRODUCTION
In recent years, the main research on energetic materials has focused on the thermal decomposition mechanism, 1−4 "hot spot" hypothesis, 5 impact explosion, 6,7 laser ignition, 8−10 and so on. Pentaerythritol tetranitrate (PETN, C(CH 2 ONO 2 ) 4 ) is widely used in military and civilian applications as a typical nitrate energetic material. 11 The impact of PETN on anisotropy has attracted people's interest and achieved many research results. 12−18 Dick 12 first found that the impact initiation of PETN strongly depends on the orientation of the crystal axis relative to the direction of impact, that is, PETN is sensitive to the impact in the [110] and [001] directions rather than in other directions. Dick et al. 13−16 proposed the steric hindrance model, through experiments and calculations, to explain the impact of PETN on anisotropy. This model pointed out that there were available sliding systems in the direction of impact that were insensitive to the initiation of impact. These systems could reduce the shear stress of shock waves due to transmission. The lack of available sliding systems causes the block in the insensitive direction. However, subsequent studies confirmed that the model had failed to expand to other energetic materials. 17,18 With the continuous update of the synthesis technology of energetic materials, 19 Klapoẗke et al. 20 synthesized PETN analogs by the nitrification reaction of tetramethylsilane in nitric acid, Si− PETN(Si(CH 2 ONO 2 ) 4 ), that is, substituting the C-atom located at the center of PETN with the Si atom. Si−PETN has a much higher sensitivity than PETN, slight friction can lead to an explosion, and the extremely high sensitivity of Si− PETN makes it difficult to study its sensitivity and other properties by experimental methods. 21−23 Liu et al., 24 through the method of density functional theory, predicted the possible chemical reaction path of the single-molecule Si−PETN and concluded that the cause of the high sensitivity of Si−PETN was the generation of low-energy barriers and highly heating Si−O bonds.
It can be seen from the very high sensitivity of Si−PETN that changing the molecular structure will have a great influence on the molecular properties. At present, the unified relationship between the crystal structure, molecular frame, and sensitivity has not yet come into consistent conclusion. 25,26 Kamlet et al. 27 have put forward a unique theory in terms of the structure, sensitivity, and thermal stability of organic highly explosive materials. However, the latest studies show that there is no single universal relationship between the molecular structure and initiative reactivity. 28 Common nitrate energetic materials similar to PETN structures are erythritol t e t r a n i t r a t e ( E T N , c h e m i c a l f o r m u l a (CH 2 ONO 2 ) 2 (CHONO 2 ) 2 ), 29 NEST-1 (2,3-bis-hydroxymethyl-2,3-dinitro-1,4-butanediol tetranitrate, chemical formula (CH 2 ONO 2 ) 4 (CNO 2 ) 2 ), 30 nitroglycerine (NG, chemical formula H(CONO 2 )(CH 2 ONO 2 ) 2 ), 31 and so on. The study on ETN began to increase due to the large-scale industrial production of meso-erythritol, 32 ETN's precursor, as a sweetening agent. ETN can be made by the nitrification of acetyl nitrate. 33,34 Oxley et al. 35 thought that ETN and PETN have the same number of nitro groups and oxygen balance; PETN can be replaced in some formulations with ETN. Furman et al. 36 used the symmetric plate impact method to study the physicochemical phenomena under the impact of ETN with nanoholes and pointed out that the decomposition mechanism of the condensed ETN followed the monomolecular path, which was consistent with studies from Oxley et al. 37 NEST-1 has similar sensitivity 30 to PETN; Landerville et al. 38 predicted the isothermal fluid static state equation of NEST-1 by DFT, which was well consistent with the experiment. NG has been used especially as an active ingredient in the manufacture of explosives, commonly used in construction and demolition industries. Waring et al. 39 studied the thermal decomposition of vapor and liquid NG in a certain temperature range and pointed out that the reaction rate significantly changed with the ratio of nitroglycerin mass to reactor volume, and the elimination reaction of NO 2 was the initial reaction path.
The above studies showed that ETN, NEST-1, and NG have similar structures to PETN. Thermodynamic and chemical properties similar to PETN were also obtained, but the relationship between the molecular structure and the physicochemical properties of molecules was still unclear. This paper, using ETN, NEST-1, and NG as the research object, simulated the physicochemical properties presented under loading conditions of hydrostatic compression and uniaxial compression using the ReaxFF/lg force field combined with the molecular dynamics method to provide references for problems of the synthesis of energetic materials.

SIMULATION DETAILS
The ReaxFF 40 reaction force field based on the first principle can be used under extreme conditions to describe the formation and fracture of the intermolecular bonds of energetic    materials 41,42 under near-real conditions through the relationship between bond length, bond level, and bond energy as well as the dynamic transfer of charges between atoms. ReaxFF/lg 43 added a description of the London dispersion force on the basis of ReaxFF, that is, amending the long-range interaction force between molecules to more accurately describe the structure and density of molecular crystals. The calculated crystal structure and equation of states showed good consistency with the experimental results.
To obtain a higher computational accuracy 44 and fully observe the physicochemical properties of the three molecular structures under the loading conditions of hydrostatic compression and uniaxial compression, 3 × 8 × 3, 5 × 2 × 5, and 5 × 3 × 6 supercells were constructed for ETN, NEST-1, and NG, respectively, to obtain the approximate atomic number and volume of the three structures (see details in Supporting Information Figure S1). Table 1 shows the specific parameters of the model, including lattice parameters and amounts of molecules and atoms. Figure 1 shows a monomolecular model of the three structures. The training set of the ReaxFF/lg force field did not contain the three structures of ETN, NEST-1, and NG, so the applicability of the ReaxFF/lg force field needed to be verified. The super crystal cell was structurally optimized by energy minimization; then the system was placed in the NVT ensemble (300 K) running 10 ps, consequently running 40 ps in the NPT ensemble (300 K, 0 GPa) to optimize atomic coordinates. After optimization, the densities of ETN, NEST-1, and NG structures were 1.786, 1.872, and 1.743 g cm −3 , respectively, which were close to experimentally measured densities 29−31 of 1.773, 1.917, and 1.846 g cm −3 , respectively. The measured lattice constant was consistent with the experimental value (see details in Supporting Information Table S1), proving that the ReaxFF/lg force field had good mobility and could be used to research on ETN, NEST-1, and NG.
Simulate the hydrostatic compression of ETN, NEST-1, and NG, compressing their super crystal cells from V/V 0 = 1.00− 0.60, where V 0 is the predicted volume below zero pressure. In each compression step, the shape and atomic coordinate of the unit crystal cell were allowed to relax without symmetry constraints at a constant volume. Uniaxial compression was done in three directions: [100], [010], and [001]. For each compression direction, the components of the lattice vector were compressed from V/V 0 = 1.00−0.70. During each compression process, only an atomic coordinate offset was allowed, while the crystal cells were fixed. The simulation process uses conditions of the periodic boundary. The simulations were performed using ReaxFF/lg potential in a large-scale atomic/molecular massively parallel simulator 45,46 with a time step of 0.1 fs.

RESULTS AND DISCUSSION
3.1. Hydrostatic Compression. At 300 K, under hydrostatic compression, the P−V curve of ETN, NEST-1, and NG is shown in Figure 2. When V/V 0 = 0.6, the hydrostatic pressure of NG was up to 30 GPa, which is similar to the pressure of PETN blast at the CJ point. 47 For the same amount of compression, the hydrostatic pressure of ETN was more than 30 GPa. Fedorov et al. 48 obtained similar results when calculating the hydrostatic compression of ETN and PETN using the first principle. The hydrostatic pressure of NEST-1 was always very high; when V/V 0 = 0.6, the hydrostatic pressure was up to 44 GPa, indicating that NEST-1 was the most difficult to compress in the three crystals.
The volume modulus (B 0 ) of the crystal and the first derivative (B′ 0 ) of the volume modulus to the pressure could be obtained by fitting of the P−V curve. The bulk modulus is used to characterize the ability of a material to resist deformation under external pressure. 49 For energetic materials, volume modulus describes the compressibility of materials, which is related to the impact sensitivity of materials. There are two common ways to fit the P−V curve: the third-order Birch−Murnaghan (BM) 50 where B 0 is the volume modulus, V is the compressed volume, V 0 is the initial volume, and B′ 0 is the first derivative of B 0 to the pressure P. This paper uses the VFRS equation to fit the P−V curves of the three crystal forms. The fitted curve is the solid line in Figure 2. B 0 and B′ 0 obtained from the fitting are shown in Table 2 (this work). It can be known from the PETN data in Table 2 that for the same state equation (BM), the obtained B 0 and B′ 0 are also quite different. This is because the result of the fitting depends on the range of the fitted data and the number of data points. 54 The results obtained will vary greatly depending on the value range. This is because the result of fitting depends on the range of the fitted data and the amount of the data points, so the results from different ranges of values can vary greatly. Because B 0 is the bulk modulus at atmospheric pressure, the pressure range for fitting chosen in ACS Omega http://pubs.acs.org/journal/acsodf Article this paper is (0 GPa, 10 GPa). Different numbers of data points are selected for multiple fittings, and the average values of B 0 and B′ 0 are obtained. The order from small to large of volume modulus of the three crystals is B 0 NEST-1 > B 0 ETN > B 0 NG , indicating that the compressibility of the three crystals is NG > ETN > NEST-1, that is, NG can be most easily compressed and NEST-1 is the most stable, which is also consistent with the fact that the previous P−V curve in the NEST-1 hydrostatic pressure is always high. Energetic materials and materials with a higher volume modulus are not easy to deform and react chemically under pressure. The three crystals have different compressibilities, which is mainly related to the molecular configuration and arrangement of crystals. Figure 3 shows the changing curves of the lattice constant of the three crystals under hydrostatic pressure. The a-and c-axis of ETN are easy to be compressed under hydrostatic pressure, while the b-axis is not easy to be compressed. The b-axis of NEST-1 is the easiest to be compressed, while the a-and c-axis can be compressed slightly when the initial pressure increases and then remains unchanged. The a-, b-, c-axes of NG have similar compressibilities. Figure 4 shows the curve of the key bond length with the compression ratio change in the three crystals, and the corresponding bond length can be found in Figure 1. To get more general results, the similar bonds in the crystal were averaged. It shows that the most stable bond in the three crystals is the N−O bond and the O−NO 2 bond has the largest amplitude of variation. Hiskey et al. 55 pointed out that the first step of thermal decomposition of the nitric ester is the cleavage of the O−NO 2 bond to form NO 2 . The C−C bond corresponding to the red icon is the C−CH 2 ONO 2 bond, which also shows a larger change than the other bonds. The unique bond of NEST-1 is C−NO 2 , presenting pretty high compressibility (pink line in Figure 4b). Comparing Figure  4a,c, we find that the C−H bond located on the main chain is more active than the C−H bond in CH 2 . The C−O bonds in the three crystals show good stability, which also shows that NO 3 generated by the nitrate energetic materials is not from C−O bond breakage.
3.2. Uniaxial Compression. The properties of energetic materials in uniaxial compression are particularly important for the plastic deformation and mechanical chemical reactions of energetic materials. Stress tensor is used to describe anisotropic response under the impact load. Figures 5−7 show the change curves of the stress tensor components σ xx , σ yy , and σ zz with the volume compression ratio of the three crystals under uniaxial compression. To make a comparison with hydrostatic pressure, the corresponding stress during hydrostatic compression is increased in each figure. In general, as V/V 0 decreases, the principal stresses σ xx , σ yy , and σ zz along the

ACS Omega
http://pubs.acs.org/journal/acsodf Article trend and amplitude of the stress tensor components σ xx and σ zz are close, but thare is a difference that the value of σ xx along the [100] direction is 3.2 GPa greater than the value of σ xx along the [010] direction, which is slightly smaller than the changing value of σ zz . In terms of NEST-1, the stress tensor component σ yy shows a significant change under uniaxial compression. When the volume is compressed to 0.2, the σ yy value in the [100] direction is greater than the hydrostatic pressure and the σ yy value in the other two directions. When V/V 0 is 0.7, the σ yy value in the [010] direction is 4.3 GPa greater than the value of σ yy . However, the change in stress tensor component σ zz is not obvious; the σ zz along the [100] and [001] directions and hydrostatic pressure show almost the same curve. The change trend and amplitude of the stress tensor components σ xx and σ yy are approximate. The difference is that the value of σ xx along the [100] direction is 3.1 GPa greater than the value of σ xx along the [001] direction, which is lower than the value of σ yy .
In terms of NG, the stress tensor component σ yy shows a large change under uniaxial compression. When the volume is compressed, σ yy along the [010] direction and the hydrostatic pressure remain highly consistent, and the σ yy values are higher than the σ yy values along the [100] and [001] directions. When V/V 0 is 0.7, the σ yy value in the [010] direction is 3.4 GPa greater than the σ yy value in the [100] direction. However, the change in stress tensor component σ xx is not obvious; the change trend and amplitude of the stress tensor components σ zz and σ yy are approximate. The difference is that the value of σ zz in the [001] direction is 2.9 GPa greater than the value of

ACS Omega
http://pubs.acs.org/journal/acsodf Article σ zz in the [010] direction, which slightly less than the change in σ yy . Figure 8 shows the average principal stress of the three crystals under uniaxial compression σ ̅ (σ ̅ = (σ xx + σ yy + σ zz )/3) and the curve of hydrostatic pressure with volume compression ratio. Obviously, the average principal stress σ ̅ and hydrostatic pressure have similar changing trends and amplitudes, suggesting that the anisotropic behavior of the three crystals, ETN, NEST-1, and NG, is weaker than that exhibited from PETN. This may be due to the high symmetry of the crystal and molecule.
3.3. Shear Stress. Shear stress is an important driving force for plastic deformation of EM crystals and is the main physical variable to characterize the anisotropism of mechanical machinery, so it is necessary to calculate the maximum shear stress for each uniaxial compression. Under the uniaxial compression, the shear stress is maximum when the angle between crystal and compression reaches 45°, where τ xy = (σ xx − σ yy )/2, τ yz = (σ yy − σ zz )/2, and τ xz = (σ xx − σ zz )/2. Figure 9 shows the curves of maximum shear stresses τ xy , τ yz , and τ xz of three kinds of crystals along the volume compression ratio V/V 0 change. It can be seen from the comparison that the change trends of the maximum shear stress of ETN and NEST-1 are similar, while the change of the maximum shear stress of NG is more complicated. It can be seen from Figure 9a Figure 9a,d,e, it can be seen that τ yz and τ xy of ETN in the [010] direction are higher than τ xz . Also, τ yz and τ xy increase exponentially in the [010] direction. It can be observed from Figure 9b,g,h that both τ xy and τ xz of NEST-1 in the [100] direction are higher than τ yz . The maximum shear stress of NG shows no obvious rule in different compression directions, but its value is below 1.5 GPa, which is the smallest of the three kinds of crystals.

CONCLUSIONS
This paper studied the physicochemical properties of typical nitrate energetic materials under hydrostatic and uniaxial compression by combining the ReaxFF/lg field with molecular dynamics simulation. Under hydrostatic compression, the P−V curve and the volume modulus (B 0 ) obtained by fitting the VFRS state equation indicated that the compressibility of the three kinds of crystals is NG > ETN > NEST-1. Due to the different molecular configuration and crystal arrangement, the three kinds of crystals have different change trends in the lattice constants of hydrostatic pressure. The a-and c-axis of ETN can be compressed easily under hydrostatic pressure, while the b-axis cannot be. The b-axis of NEST-1 is the easiest to be compressed, while the a-and c-axis can be compressed slightly when the initial pressure increases. The a-, b-, and caxes of NG have similar compressibilities. By analyzing the trend of the change of main bond length of the crystal, it can be seen that the N−O bond is the most stable in the three kinds of crystals, and the O−NO 2 bond has the largest change. The stability of the C−O bond shows that the NO 3 generated by the nitrite ester containing energetic materials is not from the C−O bond breakage.
The average principal stress σ ̅ and the hydrostatic pressure of the three kinds of crystals under uniaxial compression have similar trends and amplitudes, indicating that the anisotropic behavior of the three kinds of crystals, ETN, NEST-1, and NG, is weak. The maximum shear stresses τ xy and τ yz of ETN in the [010] direction are 1.5 GPa higher than τ xz . However, the maximum shear stress of NG shows irregularity in different compression directions, implying that there is no obvious correlation between the maximum shear stress and sensitivity. Future work may study the effects of maximum shear stress on crystal sliding systems, which will provide references for understanding the physicochemical properties of energetic materials and the problems of synthesis of energetic materials.
Model of ETN, NEST-1, and NG supercell ( Figure S1); unit cell parameters and the volume of ETN, NEST-1, and NG crystal under normal temperature and pressure (Table S1) (PDF)