Compression Behavior and Vibrational Properties of New Energetic Material LLM-105 Analyzed Using the Dispersion-Corrected Density Functional Theory

The 2,6-diamino-3,5-dinitropyrazine-1-oxide (LLM-105) is a newly energetic material with an excellent performance and low sensitivity and has attracted considerable attention. On the basis of the dispersion-corrected density functional theory (DFT-D), the high-pressure responses of vibrational properties, in conjunction with structural properties, are used to understand its intermolecular interactions and anisotropic properties under hydrostatic and uniaxial compressions. At ambient and pressure conditions, the DFT-D scheme could reasonably describe the structural parameters of LLM-105. The hydrogen bond network, resembling a parallelogram shape, links two adjacent molecules and contributes to the structure stability under hydrostatic compression. The anisotropy of LLM-105 is pronounced, especially for Raman spectra under uniaxial compression. Specifically, the red-shifts of modes are obtained for [100] and [010] compressions, which are caused by the pressure-induced enhance of the strength of the hydrogen bonds. Importantly, coupling modes and discontinuous Raman shifts are observed along [010] and [001] compressions, which are related to the intramolecular vibrational redistribution and possible structural transformations under uniaxial compressions. Overall, the detailed knowledge of the high-pressure responses of LLM-105 is established from the atomistic level. Uniaxial compression responses provide useful insights for realistic shock conditions.


Introduction
Energetic materials (EMs), such as explosives, propellants, and pyrotechnics, are widely used for civilian and military purposes [1,2]. EMs can undergo multiple phase transitions during detonation, in which the molecules are modified and rearranged under high-temperature and high-pressure conditions. Furthermore, the evolutions of the pressure-dependent physical and chemical properties of EMs are important in the control of their sensitivity, performance, and safety [3,4]. Recently, 2,6-diamino-3,5-dinitropyrazine-1oxide (LLM-105), a new generation of high-energy and low-sensitivity EM, was synthesized by Lawrence Livermore Laboratories in the United States, and its foundational properties attracted attention.
Single-crystal X-ray diffraction experiment [5] showed that the LLM-105 crystal contains C, H, N, and O with a space group of P2 1 /n in the monoclinic crystal, in which every H atom is involved in intramolecular hydrogen-bonding interactions with its neighboring O atom and in intermolecular interactions with adjoining molecules [6]. At ambient conditions, LLM-105 with a wave-like π-π stacking arrangement possesses rich intra-and Molecules 2021, 26, 6831 2 of 12 intermolecular hydrogen-bonded networks, which buffer the perturbations of the external environment on the LLM-105 crystal [7,8].
Several experimental studies showed that the LLM-105 crystal exhibits a strong mechanical and thermal stability under external temperature and pressure loading. X-ray diffraction, Raman spectroscopy, and infrared spectroscopy are powerful experimental tools to elucidate the changes in the crystalline structure under extreme conditions. Gump et al. [9]. and Stavrou et al. [10]. used X-ray diffraction technology to determine the isothermal equation of the state of unreacted LLM-105 crystal. Their results demonstrated that the phase structure of LLM-105 remained stable at a pressure up to 20 GPa and temperature up to 513 K. Recently, Zhang et al. [11]. confirmed that the strong intermolecular and intramolecular hydrogen bond network of LLM-105 crystal contributes to the stability of the structure through high-pressure Raman and infrared spectroscopy. Importantly, the first phase transition of LLM-105 is confirmed near the hydrostatic pressure of 30 GPa, in which several Raman peaks disappear, and vibrations show an abrupt change at about 30 GPa. In addition, by analyzing the photoluminescence and absorption spectra [12], the LLM-105 crystal has an adjustable indirect band gap and a phase change in the electronic structure around 10 GPa.
Alternatively, theoretical calculations can give several important insights at atomic and molecular levels on the structural evolution and anisotropy properties of LLM-105 crystals under hydrostatic compression. Wu et al. [13]. used the standard density functional theory (DFT) to simulate the structural changes in the LLM-105 crystal at a hydrostatic pressure of 0-50 GPa. They pointed out that the evident irregular changes in lattice parameters, including unit-cell angles, bond lengths, bond angles, and band gaps, were closely related with the structural transformations of LLM-105. However, Stavrou et al. [10] conducted first-principles molecular dynamics simulations to verify that the ambient pressure phase of LLM-105 could remain stable up to 20 GPa. On the basis of the dispersion-corrected DFT calculations, Manaa et al. [14]. found no evidence for structural phase transitions of up to 45 GPa by calculating the equation of state (EOS) of the LLM-105 crystal under hydrostatic compression. Zong et al. [15]. used the generalized gradient approximation parameterized by Perdew-Burke-Ernzerhof (GGA-PBE method) to simulate the properties of LLM-105 at a pressure of 0-100 GPa and found that the structural phase change occurs at 35 GPa. In addition, the decomposition mechanism of LLM-105 is studied on the basis of DFT. Zhang et al. [16]. showed that the intramolecular hydrogen transfer is reversible during the decomposition process of LLM-105, and the reversible hydrogen transfer can buffer the external stimuli caused by energy transfer and slight structural changes.
Although considerable experimental and theoretical efforts focused on the compression behaviors and thermal decomposition mechanism [17][18][19][20][21][22][23], the anisotropic vibrational property and structural modification of LLM-105 under different pressure loading conditions require further clarification. In this work, we performed systematic first-principles calculations to elucidate the crystalline structure and vibrational properties under hydrostatic and uniaxial compressions. Under high-pressure conditions, the hydrogen bond network between adjacent molecules maintained the stability of the crystal structure. The red-shift of vibrational modes was observed under hydrostatic and uniaxial compressions. The discontinuous change in high-wavenumber modes along a [001] orientation compression indicates the onset of possible structure transformations. The in-depth research in this area aids in the understanding of the structural response of the LLM-105 crystal under a high pressure, further providing theoretical support for understanding the stability of EMs.

Lattice Parameters of LLM-105 at Ambient and Pressure Conditions
Here, the lattice parameters of LLM-105 and the internal coordinates of the atoms are fully relaxed on the basis of DFT-D calculations at zero pressure. In Table 1, the calculated results are compared with previous experimental and theoretical data [5,9,10,13,15,18]. The DFT calculations with DFT-D2 schemes significantly improved the description of the crystal structure at ambient conditions. The calculated lattice parameters are in agreement with experimental values [5] and superior to previous first-principles calculations based on the force field [18] and parameterized CA-PZ and PW91 scheme [13], and close to the theoretical result by Zong et al. [15]. Specifically, the calculated lattice parameters are consistent with the experimental measurements. The relative errors of lattice parameters (including a, b, c, V, and β) are below 1.1%, indicating that the current PBE-D2 scheme gives a reasonable description of intermolecular interactions, including intermolecular hydrogen bonding.  Figure 1 shows the evolutions of cell volume and lattice parameters of LLM-105 under pressure up to 10 GPa. Overall, the cell volume and lattice parameters of LLM-105 from the DFT-D calculations are reasonably consistent with experimental XRD data [10,11]. Figure 1a shows the variation in the cell volume of LLM-105 from our DFT-D calculations and previous experimental results as the pressure increases. Current DFT-D calculations reproduce the trend of pressure-induced volume compression from experiments, in which the volume compression at 10 GPa is up to 22%. The pressure-induced reductions in lattice parameters are displayed in Figure 1b. The lattice parameter b exhibits a significant reduction compared to lattice parameters a and c as the pressure increases, indicating that the LLM-105 crystal is easily compressible along the b-axis. This phenomenon occurs because the weak interlayer vdW interaction along the b-axis is more compressible than the relatively strong hydrogen bonds in the molecules along the aand c-axes. The anisotropic compressibility is common feature for EMs [24,25] and reflects the strength of intermolecular interactions in different directions, which is related to the anisotropic molecular responses of LLM-105 under pressure loading. In addition, third-order Birch-Murnaghan EOS are applied to fit the bulk modulus (B0) of LLM-105 to evaluate its stiffness on the basis of the calculated P-V relationship. The third-order Birch-Murnaghan EOS are given as follows: [26]  In addition, third-order Birch-Murnaghan EOS are applied to fit the bulk modulus (B 0 ) of LLM-105 to evaluate its stiffness on the basis of the calculated P-V relationship. The third-order Birch-Murnaghan EOS are given as follows [26]: where P represents the applied pressure, and V 0 and V represent volumes at ambient and pressure conditions, respectively. B 0 represents the bulk modulus and B 0 represents its first pressure derivative. Table 2 summarizes the bulk modulus (B 0 ) and its derivative (B 0 ) of LLM-105, as well as previous experimental and theoretical results [9][10][11]14,15,18]. Our calculated B 0 and B 0 are 16.03 and 8.67 GPa, respectively, which are in accordance with the reported experimental and theoretical data.  [11] 19.23 6.70 Calc. [15] 16.5 9.4 Calc. [18] 13.8 11.7

Raman Spectra of LLM-105 at Ambient and Hydrostatic Conditions
Raman spectroscopy is used to determine the crystal structure and chemical bonding of materials by identifying the vibrational modes of molecules in experimental and theoretical studies [27]. An atomistic level spectroscopic simulation can demonstrate experimental findings, thereby providing deep insights into structural modification and noncovalent interactions. For the space group P2 1 /n (2/m), a group theoretical analysis gives 228 vibrational modes (Γ = 57A g + 57B g + 57A u + 57B u ) for the LLM-105 crystal at ambient conditions. Among them, the calculated vibrations of even parity (g) are Raman-active. Thus, 57 Raman vibrations of each symmetry are present.
In Figure 2, the calculated and experimental Raman spectra [11] of LLM-105 are present over a broad wavenumber range. The relative wavenumber and intensity of characteristic peaks are essentially consistent between calculated and experimental spectra. The high-wavenumber C-H and O-N-O stretch patterns have a larger deviation than other internal modes due to the lack of anharmonic effects and insufficient description of intermolecular interactions. In addition, the characteristics of the calculated and experimental internal/lattice vibrational modes, including wavenumber, assignment, and pressuredependence, are presented in Table S1. Several coefficients of pressure dependence from theoretical calculations deviate from the experiments due to a lack of temperature effects under harmonic approximation.
internal modes due to the lack of anharmonic effects and insufficient description of intermolecular interactions. In addition, the characteristics of the calculated and experimental internal/lattice vibrational modes, including wavenumber, assignment, and pressure-dependence, are presented in Table S1. Several coefficients of pressure dependence from theoretical calculations deviate from the experiments due to a lack of temperature effects under harmonic approximation. The pressure dependence of Raman spectra under hydrostatic compressions is shown in Figure 3a, and detailed comparisons of the selected Raman shift between calculated and experimental results are presented in Figure 3b-d. The Raman spectra in the experiment are concentrated in the middle-and low-wavenumber regions (<1700 cm −1 ). The trends of the calculated Raman shift of most modes are in agreement with the experimental measurements, which indicated that the current DFT-D scheme is relatively reliable. Additionally, all of these vibrational modes are blue-shifted due to the reduced intermolecular distance as a result of the enhanced intermolecular interactions under pressure. Importantly, the high-wavenumber modes (v1-v4) associated with the motion of the amino group are predicted at 3200−3500 cm −1 . The asymmetric NH2 stretching (v1 and v2) has higher wavenumber than symmetric ones (v3 and v4). The NH2 stretching modes, including asymmetric v2, symmetric v3, and symmetric v4 vibrations, present a red-shift at 0-10 GPa, indicating that the LLM-105 crystal tends to be unstable with an elevated pressure. However, the wavenumber of the asymmetric v1 vibrational mode is nearly identical in the entire pressure range. From the vibrational pattern of LLM-105 in Figure 4a,b, modes v1 and v3 are from the motion of the amino group at site A. The N1-H2 bond participates in the v1 vibrational mode with a relatively strong contribution, whereas the N1-H3 bond participates in the v3 vibrational mode. In Figures 4c and S1, the configuration of the hydrogen bond N1-H2···O1 is stable with an initial bond length and angle at 0-10 GPa; this is why the wavenumber of the v1 vibrational mode remains almost constant. The length of the N1-H3 bond does not change, whereas the bond angle of H2-N1-H3 continues to decrease, which can cause the red-shift of the v3 vibration mode. In addition, the hydrogen The pressure dependence of Raman spectra under hydrostatic compressions is shown in Figure 3a, and detailed comparisons of the selected Raman shift between calculated and experimental results are presented in Figure 3b-d. The Raman spectra in the experiment are concentrated in the middle-and low-wavenumber regions (<1700 cm −1 ). The trends of the calculated Raman shift of most modes are in agreement with the experimental measurements, which indicated that the current DFT-D scheme is relatively reliable. Additionally, all of these vibrational modes are blue-shifted due to the reduced intermolecular distance as a result of the enhanced intermolecular interactions under pressure. Importantly, the highwavenumber modes (v 1 -v 4 ) associated with the motion of the amino group are predicted at 3200−3500 cm −1 . The asymmetric NH 2 stretching (v 1 and v 2 ) has higher wavenumber than symmetric ones (v 3 and v 4 ). The NH 2 stretching modes, including asymmetric v 2 , symmetric v 3 , and symmetric v 4 vibrations, present a red-shift at 0-10 GPa, indicating that the LLM-105 crystal tends to be unstable with an elevated pressure. However, the wavenumber of the asymmetric v 1 vibrational mode is nearly identical in the entire pressure range. From the vibrational pattern of LLM-105 in Figure 4a,b, modes v 1 and v 3 are from the motion of the amino group at site A. The N 1 -H 2 bond participates in the v 1 vibrational mode with a relatively strong contribution, whereas the N 1 -H 3 bond participates in the v 3 vibrational mode. In Figure 4c and Figure S1, the configuration of the hydrogen bond N 1 -H 2 ···O 1 is stable with an initial bond length and angle at 0-10 GPa; this is why the wavenumber of the v 1 vibrational mode remains almost constant. The length of the N 1 -H 3 bond does not change, whereas the bond angle of H 2 -N 1 -H 3 continues to decrease, which can cause the red-shift of the v 3 vibration mode. In addition, the hydrogen bond network between the adjacent molecules always maintains a parallelogram-like configuration (O 1 ···H 1 ···O 2 ···H 2 , Figure 4b,c. When the pressure increases, the O···H-bonded molecular pairs, via nearby molecules and the cross-stacking of molecular layers for LLM-105, seem to be tightly confined, in which such special configurations can contribute to maintain the stability of the crystal structure under high-pressure conditions. bond network between the adjacent molecules always maintains a parallelogram-like configuration (O1···H1···O2···H2, Figure 4b,c. When the pressure increases, the O···H-bonded molecular pairs, via nearby molecules and the cross-stacking of molecular layers for LLM-105, seem to be tightly confined, in which such special configurations can contribute to maintain the stability of the crystal structure under high-pressure conditions.   bond network between the adjacent molecules always maintains a parallelogram-like configuration (O1···H1···O2···H2, Figure 4b,c. When the pressure increases, the O···H-bonded molecular pairs, via nearby molecules and the cross-stacking of molecular layers for LLM-105, seem to be tightly confined, in which such special configurations can contribute to maintain the stability of the crystal structure under high-pressure conditions.    Table S1 shows that most vibrational modes of LLM-105 are combinations of cation and anion vibrations at ambient conditions. The Raman frequencies of these modes increase with the elevated pressure in most cases below 1600 cm −1 . The highest coefficient (8.13 cm −1 ·GPa −1 ) is observed for the low-frequency internal mode v 29 , whereas lattice modes have coefficients less than 2.99 cm −1 ·GPa −1 .

Stress Tensor and Shear Stress
The high-pressure behavior of EMs under uniaxial compressions can be associated with their anisotropic shock sensitivity. The stress tensor and its derived shear stress are calculated to evaluate the anisotropic response under uniaxial loading and to quantify this sensitivity. On the basis of the calculated diagonal elements of the stress tensor, the principal stresses are obtained under uniaxial compressions in Figure S2. Three principal stresses, σ xx , σ yy , and σ zz , along the [100] and [010] orientations appear to increase under the uniaxial compression ratio (V/V 0 ). However, the difference in variation trends and amplitudes clearly indicates the significant anisotropic behavior of LLM-105. For example, the largest principal stress σ zz along the [001] orientation is 17.06 GPa at V/V 0 = 0.78, which is larger than the σ xx (12.07 GPa) along the [100] orientation and the σ yy (11.94 GPa) along the [010] orientation. Thus, we suggest that the [001] orientation can be more sensitive than the [100] and [010] orientations. Furthermore, the average principal stress σ, defined as σ = σ xx + σ yy + σ zz /3, is calculated as the function of V/V 0 in Figure S3, which demonstrates the anisotropic behavior of LLM-105 under compression.
Shear stresses, which are defined as τ xy = σ xx − σ yy /2, τ xz = (σ xx − σ zz )/2 and τ yz = σ yy − σ zz /2, are the major parameters used to calculate the mechanical anisotropy. The calculated shear stresses as a function of V/V 0 under two uniaxial compressions are displayed in Figure 5 Table S1 shows that most vibrational modes of LLM-105 are combinations of cation and anion vibrations at ambient conditions. The Raman frequencies of these modes increase with the elevated pressure in most cases below 1600 cm −1 . The highest coefficient (8.13 cm −1 ·GPa −1 ) is observed for the low-frequency internal mode v29, whereas lattice modes have coefficients less than 2.99 cm −1 ·GPa −1 .

Stress Tensor and Shear Stress
The high-pressure behavior of EMs under uniaxial compressions can be associated with their anisotropic shock sensitivity. The stress tensor and its derived shear stress are calculated to evaluate the anisotropic response under uniaxial loading and to quantify this sensitivity. On the basis of the calculated diagonal elements of the stress tensor, the principal stresses are obtained under uniaxial compressions in Figure S2      Most vibrational frequencies shift gradually toward higher frequencies (i.e., blue-shift) with increasing pressure due to reduced interatomic distances. This behavior is consistent with the previous, high-pressure Raman experiment of the LLM-105 crystal [20]. The observed blue-shift in most vibrational frequencies is caused by the strengthening of chemical bonds. The inter-and intramolecular hydrogen bonds are enhanced under compression. Several vibrational modes show red-shift for three uniaxial compressions. The discontinuity of the Raman shift is observed along the [001] compression, which indicates the onset of possible phase transitions under uniaxial compressions. To understand such behavior, we examined the changes in the amino group of LLM-105, since it could provide insight into the stability of this crystal at extreme conditions in Figure 7.   For the [100] orientation, the vibrational frequencies of the v 1 and v 2 modes are almost identical in the studied pressure range, which is similar to the case under hydrostatic conditions. The chemical bonds (N 1 H 1 and N 2 H 3 ) related to the v 1 and v 2 modes are insensitive at 0-10 GPa in Figure 7a. The compressibility values of N 1 H 1 and N 2 H 3 bonds are only 0.29% and 0.19%, respectively. By contrast, the N 1 H 2 and N 2 H 4 bonds lengthen with the increased pressure and are the main cause of the red-shifts of high-wavenumber modes (v 3 and v 4 ). Similarly, for [010] orientation, the v 1 mode shows a remarkable red-shift, which is due to the lengthening of N 1 H 1 covalent bonds.
In addition, coupling modes are observed in the studied pressure range, which can be related to the intramolecular vibrational redistribution of Ems [28]. In Figure 7b, internal modes v 3 and v 4 have the same A g symmetry and vibrational patterns (NH 2 symmetrical stretching). The initially low-intensity ν 3 mode is enhanced and shifts with the increasing pressure to a higher frequency. Above V/V 0 = 0.92, the intensity exchange and the avoided frequency crossing between the two modes are clearly observed. Such coupling mode contributes to the intramolecular vibrational energy transfer by the intensity exchange and avoided frequency crossing between different modes.
For the [001] orientation, all NH bonds show discontinuous changes at V/V 0 of 0.82. The shortening lengths of the N 1 -H 2 bonds result in the blue-shift of the v 3 vibrational mode, whereas the lengthening of other N-H covalent bonds can be treated as the main reason that leads to the red-shifts of high-wavenumber modes.
Overall, the vibrational properties of LLM-105 under hydrostatic and uniaxial conditions are simulated to examine its pressure effects and anisotropy. Based these results, the intermolecular hydrogen-bonded networks shrink and/or rearrange when the lattice space is rapidly compressed, which results in the redistribution of the electronic density at the atoms. As a result, the intramolecular geometry and the intermolecular interactions of LLM-105 are further modified and presented as distinct anisotropic responses.

Computational Methods
The first-principle calculations based on DFT were performed using the Cambridge Sequential Total Energy Package [29]. The exchange-correlation interaction was treated within the GGA-PBE function [30]. The semiempirical dispersion correction by Grimme (DFT-D2) [31] was used to account for the intermolecular noncovalent interactions of LLM-105 crystals. The Brillouin zone sample of Monkhorst-Pack [32] k grid of 0.05 Å −1 was adopted, and a kinetic energy cutoff of 1000 eV for the plane wave basis was used. The periodic boundary conditions of LLM-105 was employed and the unit cell parameters and atomic coordinates were fully relaxed by the Broyden, Fletcher, Goldfrab, and Shannon algorithms to obtain the equilibrium crystal structure of LLM-105 [33]. Among them, the convergence criteria of the maximum total energy of 5 × 10 −6 eV per atom, with a maximum stress of 0.02 GPa, maximum force of 0.01 eV·Å −1 , and maximum displacement of 5 × 10 −4 Å were adopted for all calculations. Hydrostatic compressions were applied from 0 GPa to 10 GPa at a step of 0.5 GPa. Uniaxial compression was applied up to 78% of the equilibrium cell volume in steps of 2% along lattice directions. At each compression step, only atomic coordinates were allowed to relax, and the lattice was fixed. The maximum stress under uniaxial compression was less than 16 GPa in each lattice direction for improved contrast with hydrostatic compression.
On the basis of the optimized structure, the Raman spectrum was obtained within the density functional perturbation theory [34]. The wavenumbers of vibrational modes were obtained from the diagonalizable dynamical matrix. The Raman intensity was calculated using the linear response formalism combined with the dispersion correction method. The experimental factors of temperature (298 K) and incident light wavelength (514.5 nm) were considered for real Raman intensity. Then, the spectrum band was broadened by the Lorentzian function with 15 cm −1 width.