Stability, electronic properties, and decomposition mechanisms of 2,6-diamino-3,5-dinitropyrazine-1-oxide crystals with different vacancy defects

Abstract

The effects of vacancies on the energies, electronic properties and thermal decomposition mechanisms of the 2,6-diamino-3,5-dinitropyrazine-1-oxide (LLM-105) crystal were investigated by using the density functional tight binding (DFTB) method and DFTB-based molecular dynamics (DFTB-MD). It is discovered that the presence of the vacancies can significantly alter the formation energy and binding energy of the LLM-105 crystal. Additionally, the existence of vacancies can decrease the band gap of the LLM-105 crystal, thereby facilitating its electron transitions. In addition, the dynamics simulations indicate that the introduction of the vacancies does not alter the primary decomposition mechanisms and types of products of the LLM-105 crystal. However, the presence of vacancies can increase the diffusion coefficients of carbon, hydrogen, oxygen, and nitrogen atoms during the decomposition of LLM-105. Furthermore, the presence of the vacancies can change the evolution of the numbers of various chemical bonds over time.

---

1. Introduction

The crystal growth process of explosives is influenced by various factors, including temperature, solution supersaturation,¹ solvent type, and all types of impurity ions in the solution. These factors can lead to the formation of many defects, such as dislocations, twins, doping, vacancies etc., which can significantly alter the physical and chemical properties of the explosives. Furthermore, these defects will also impact the explosion process and the initial decomposition mechanisms of energetic materials. It is commonly understood that the energy distribution in explosives containing vacancies differs significantly from that in a perfect crystal, and their initial decomposition tends to occur primarily in the vicinity of the vacancies. This phenomenon is closely related to the “hot spot” hypothesis, which is based on a large body of experimental evidence: explosives with defects are easily triggered to decompose under mild conditions, whereas ideal explosive crystals usually require a stronger external stimulus for decomposition.²

Currently, various theoretical research methods have been employed to investigate the influence of crystal defects on explosives. Liu et al.³ and He et al.⁴ utilized the density functional theory with dispersion correction (DFT-D) and charge self-consistent tight-binding density functional (SCC-DFTB) in conjunction with the multi-scale shock technique (MSST) method to study the effects of molecular holes on the 3-nitro-1,2,4-triazole-5-one (NTO) crystals, respectively. They found that the presence of defects has a significant effect on the electronic structure of explosive crystals, especially on the band gap values. The vicinity of the void was more prone to the formation of reaction centers, where explosive molecules within these reaction centers were more likely to undergo decomposition. Yan et al.⁵ discovered that the existence of vacancy defects decreases the cohesion energy density (CED) and increases the maximum bond lengths (L_max) of N–NO₂ for 1,3,5-trinitro-1,3,5-triazine (RDX) crystals through molecular dynamics (MD) simulations, suggesting that the presence of defects enhances the sensitivity of the RDX crystals. In addition, they observed that a higher defect concentration and a more uneven distribution of the defects lead to the higher sensitivity of RDX. Li et al.⁶ investigated the effects of void defects on hexanitrohexaazaisowurtzitane (CL-20) co-crystals under impact. Their findings indicate that the CL-20 co-crystals containing void defects generate localized regions of high temperature, known as hot spots, upon impact. These hot spots have a considerable influence on the decomposition process of the molecules surrounding the voids. They also found that the presence of voids in the CL-20 co-crystals also diminishes safety, leading to a significant increase in the sensitivity of high-energy cocrystals such as CL-20/trinitrotoluene (CL-20/TNT) and CL-20/1,3-dinitrobenzene (CL-20/DNB).⁶ Moreover, many studies have demonstrated that the presence of defects can modify the time required for complete decomposition, the reaction rate constants, and the reaction activation barrier of the explosives. Through reactive force field MD (ReaxFF-MD) simulations, Zhou et al.⁷ found that the presence of molecular vacancy defects does not alter the frequency of certain reactions such as N–NO₂ bond cleavage, HONO dissociation, and main ring cleavage of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX), but it increases the reaction rate constant of HMX and decreases its reaction activation energy. Ji et al.⁸ found that the presence of vacancy defects not only affects the thermal decomposition rate of benzotrifuroxan(BTF) molecules, but also shifts the initial decomposition position of the BTF molecules at high temperatures. They also observed that in the early stages of the reaction at 2500 K, the NTO molecules adjacent to the impurity molecules decomposed first. Additionally, the NTO molecules within the impurity-containing system decomposed at a faster rate compared to those in the ideal system.⁹ Furthermore, the crystal defects can have an impact on the sensitivity and mechanical properties of the explosives.^10,11

2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105) is a common ingredient in various explosive formulations and solid propellants due to its high energy content, exceptional insensitivity and thermal stability.^12,13 The energy value of the LLM-105 crystal has been reported to be approximately 85% of that of HMX and 15% higher than that of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) and its insensitivity is very close to that of TATB. It exhibits superior stability compared to most known explosives when subjected to external stimuli such as shock, friction, and sparks.^14,15 Theoretical studies have indicated that the wave-like π–π stacking mode and hydrogen bonding in the LLM-105 crystal are important factors contributing to its insensitivity under external stimuli.^16,17 Hence, LLM-105 has been widely utilized and extensively studied for its structure and properties. Manaa et al.¹⁸ reported the unreacted equation of state (EOS) of the crystal LLM-105 under hydrostatic compression of up to 45 GPa. Ma et al.¹⁹ used DFT combined with thermodynamic theory to predict the detonation performance of LLM-105.

In this work, perfect and defect-containing LLM-105 crystals with different vacancy concentrations were built. Then, their energy, electronic properties, initial decomposition mechanisms, and diffusion behavior at high temperature were studied by DFTB and DFTB-MD methods. The influence of the vacancy defects on the crystal structure and decomposition mechanism of the explosive was discussed. Our work may provide new ideas for understanding the behaviors of explosive crystals with defects.^20–23

2. Computational methods

All calculations were performed using the third-order corrected self-consistent charge DFTB method with first-principles London dispersion correction D3 (SCC-DFTB-D3) as implemented in the CP2K package.^24–27 This approach was based on a second-order expansion of the Kohn–Sham total energy in terms of charge density fluctuations by DFT, which can describe the total energies, atomic forces, and charge transfer self-consistently.²⁸ Recently, DFTB has been proven to possess comparable accuracy to DFT in studying the structure and properties of explosive crystals.^29,30

The SCC tolerance was 10⁻⁶ au and the conjugate gradient (CG) method was used for cell optimization. The limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) method was used for geometry optimization. The convergence criterion for a maximum geometry change was less than 3.0 × 10⁻³ Bohr, maximum force component less than 4.5 × 10⁻⁴ Hartree Bohr⁻¹, root-mean-square geometry change less than 1.5 × 10⁻³ Bohr, and root-mean-square force less than 3.0 × 10⁻⁴ Hartree Bohr⁻¹. The cell structure of the LLM-105 crystal was determined through the experiments with lattice parameters of a = 5.709 Å, b = 15.844 Å, and c = 8.416 Å.³¹ The lattice parameters calculated by the DFTB method were a = 5.58 Å, b = 16.422 Å, and c = 8.409 Å, indicating that the calculation deviations of the lattice parameters are less than 3.6%.

The density of states of systems can be obtained using eqn (1):

where E_n(k) describes the dispersion of the given band and the integral is determined over the Brillouin zone. An alternative representation of the density of states is based on the fact that N_n(E)dE is proportional to the number of allowed wave vectors in the nth band in the energy E to (E + dE).

To ensure the accuracy of the calculation results, we constructed a 4 × 1 × 3 supercell model. Fig. 1 shows the molecular structure, unit cell, and 4 × 1 × 3 supercell of LLM-105. The numbers of molecules in the unit cell and supercells of LLM-105 are 4 and 16, respectively. Among them, the LLM-105 supercell has three layers with 16 molecules per layer.

Fig. 1 Molecular structure (a), unit cell (b), and 4 × 1 × 3 supercell of LLM-105. Gray, blue, red, and white spheres stand for C, N, O, and H atoms, respectively.

Considering the concentration of vacancies, we constructed different vacancy models by deleting different numbers of LLM-105 molecules in different positions on the same molecular layer in the ideal supercell model. Table 1 presents the schematic plots, total energy, formation energy, and binding energy of the LLM-105 crystals with different vacancy concentrations. In the LLM-105 supercell model, one, two, and three LLM-105 molecules were deleted to obtain vacancy models with vacancy concentrations of 2.08%, 4.17%, and 6.25%, respectively. Here, when the vacancy concentration is the same, the different positions of the vacancy atoms will also affect the calculation results. Therefore, we designed the models with different vacancy positions based on the same vacancy concentration. For the model containing one vacancy molecule, we labeled it as V1. For the models containing two vacancy molecules, we labeled it as V2-1, V2-2, V2-3, and V2-4, respectively. For the models containing three vacancy molecules, we labeled it as V3-1, V3-2, V3-3, and V3-4, respectively. The Cartesian coordinates of all the optimized geometries are listed in the ESI.†

Table 1 Schematic plots, total energy, formation energy, and binding energy of the vacancy-containing LLM-105 crystals

Model	Projection along the b-axis^a	E (eV)	E f (eV)	E b (eV)
a Each square represents one LLM-105 molecule, and green represents a vacancy.
V0		−51 068.72	0
V1		−50 003.55	0.26
V2-1		−48 938.63	1.14	0.62
V2-2		−48 938.40	1.37	0.84
V2-3		−48 938.37	1.40	0.88
V2-4		−48 938.59	1.19	0.66
V3-1		−47 873.53	2.20	1.41
V3-2		−47 873.77	1.96	1.17
V3-3		−47 873.60	2.13	1.34
V3-4		−47 873.35	2.37	1.59

---

DFTB-MD simulations were performed with the NVT ensemble to study the thermal decomposition process of LLM-105 at 2500 K. Firstly, ideal supercell, single-vacancy, double-vacancy, and triple-vacancy models were relaxed at 298.15 K for 10 ps to obtain their equilibrium structure. Then, their relaxed structures were used as the initial structures for DFTB-MD simulations. The total simulated time is 15 ps with a time step of 0.5 fs. Based on DFTB-MD simulation, it is possible to obtain the change of products with time in order to study the decomposition mechanisms.

3. Results and discussion

The existence of defects may affect the stability of the crystal.^5,20 The vacancy formation energy (E_f) is often used to evaluate the stability of these systems. In addition, the vacancy binding energy (E_b)²¹ can be introduced to describe the change of the system energy of the formation of double-vacancy and triple-vacancy systems by the single-vacancy system, in order to explore the differences in vacancy formation energy between these systems. E_f and E_b can be obtained using formulas (2) and (3):

E_f = E + nE_LLM-105 − E₀ (2)

E_b = E_f + nE_f(V1) (3)

where E_f is the vacancy formation energy, E_b is the vacancy binding energy, E is the optimized energy of the system, n is the numbers of LLM-105 molecules in the system, E_LLM-105 is the energy of a single LLM-105 molecule, E₀ is the energy of the ideal system, and E_f(V1) is the single vacancy formation energy. The vacancy formation energy and vacancy binding energy for these different vacancy systems are presented in Fig. 2.

Fig. 2 Formation and binding energies of the vacancies in the LLM-105 crystals with different vacancy concentrations.

It can be seen in Fig. 2 that the vacancy formation energies for all the double-vacancy systems are larger than that of the single-vacancy system, and the vacancy formation energies of all the triple-vacancy systems are larger than those of the single vacancy system and double-vacancy systems, which indicates that as the concentration of the vacancy increases, the vacancy formation energies of these systems increase gradually and the formation of the vacancy systems becomes more difficult. When the concentration of the vacancy is the same, the position of the vacancy will also affect the energy of the system. In the multi-vacancy systems, the values of E_f and E_b of V2-1 and V3-2 systems are the lowest ones among the double-vacancy and the triple-vacancy systems, respectively, indicating that these two systems are most favorable in thermodynamics and the easiest to form. Compared with different systems with the same concentration of the vacancy among the double-vacancy and triple-vacancy systems, it can be seen that the structures of the vacancies in V2-1 and V3-2 are similar, whose vacancies are aggregated. Therefore, the aggregated arrangement of the vacancies in the LLM-105 crystal is the most favorable in thermodynamics.

To disclose the origin of the electronic structure caused by the defects, we performed the electronic structure analysis of the most stable configurations for corresponding vacancy-containing prototypes (denoted as V0, V1, V2-1, and V3-2). Table 2 summarizes the values of band gaps of these models. It is found that the band gaps of the vacancy systems are significantly decreased, causing the electron to transfer more rapidly between the conduction band minima (CBM) and valence band maxima (VBM) and thus promoting the decomposition of the explosive molecules, compared with the ideal system V0.

Table 2 Band gaps of theV0, V1, V2-1, and V3-2 systems

Models	V0	V1	V2-1	V3-2
Band gap (eV)	1.89	1.40	1.34	1.30

---

The total densities of states (TDOS) of the vacancy-containing LLM-105 crystals are displayed in Fig. 3. As the concentration of the vacancy increases, the TDOS distribution of the LLM-105 crystal not only moves to the CBM, but also their corresponding peak values have a declining trend. For compensation, Fig. 4 also presents the frontier molecular orbital distribution of the conduction band bottom (CBM) and the valence band top (VBM) in these systems. There is an obvious net charge aggregation occurring at diverse defects in the LLM-105 crystal compared with the V0 system, further suggesting that the frontier energy levels in the vacancy systems are all localized in the local region near the vacancies, and these lead to the decrease in the band gap of the system.

Fig. 3 Total density of states (TDOS) of LLM-105 in the V0, V1, V2-1, and V3-2 systems.

Fig. 4 Distributions of the frontier orbitals in the ideal and vacancy-containing systems. The red orbitals stand for the CBM, while the blue ones for the VBM.

To explore the effects of the vacancy on the decomposition mechanisms of LLM-105 at high temperatures, we selected V0, V1, V2-1, and V3-2 as the object systems to simulate their decomposition process at 2500 K in combination with the DFTB-MD method. The results from the dynamics simulations show that the presence of a vacancy does not change the decomposition mechanisms of LLM-105, but it changes the time of the complete decomposition of LLM-105 and the position of the bond-breaking in the LLM-105 crystal. Among them, the initial decomposition snapshots of these vacancy systems are given in Fig. 5. It can be seen that the positions where LLM-105 begins to decompose in the ideal system are relatively random, but the LLM-105 molecules that first decompose in the vacancy-containing systems are mainly distributed near the vacancies, which is consistent with the hot spot hypothesis.³² This is consistent with the conclusion based on the electronic structure analysis above that the VBM and CBM in the vacancy-containing systems are mainly distributed near the vacancies, which leads to a decrease in the band gaps of the systems and thus results in the instability of the LLM-105 crystals.

Fig. 5 Snapshots of the LLM-105 molecules that firstly decomposed in the V0, V1, V2-1, and V3-2 systems. The green area indicates the location of the vacancies.

The initial decomposition mechanisms of the LLM-105 molecule are illustrated in Fig. 6. Obviously, the initial decomposition mechanisms of the LLM-105 molecule mainly include four paths: C–NO₂ bond cleavage, nitro–nitrite rearrangement to eliminate NO, hydrogen transfer to eliminate OH, and hydrogen transfer to eliminate HONO. It is generally believed that the decomposition of LLM-105 is initiated by the cleavage of C–NO₂ to release NO₂.³³ In terms of the reaction energy barriers, the C–NO₂ bond cleavage in LLM-105 has an energy barrier of 63.8 kcal mol⁻¹. The reaction energy barrier of the hydrogen transfer is 12.7 kcal mol⁻¹. The nitro–nitrite rearrangement, on the other hand, has a reaction energy barrier of 57.3 kcal mol⁻¹.³⁴ Although the hydrogen transfer exhibits a lower energy barrier compared to the NO₂ removal, it typically occurs below 580 K, whereas the NO₂ removal tends to take place above this temperature.³⁵ Furthermore, the rates of the nitro–nitrite rearrangement and subsequent NO elimination are consistently lower than those of other reactions, indicating that the C–NO₂ bond dissociation is the predominant reaction path during the initial thermal decomposition of LLM-105.

Fig. 6 Thermal decomposition channels of the LLM-105 molecule.

Fig. 7 shows the main products during the decomposition of LLM-105 with different concentrations of vacancies at 2500 K. It is found that the main products are distributed throughout these systems. Additionally, the evolution trend of the decomposition products over time remains unchanged in different systems. This finding indicates that the presence of the vacancies does not alter the decomposition mechanisms of LLM-105. Among the main products, NO₂ is generated first, which indicates that the C–NO₂ bond cleavage occurs initially. In addition, the number of NO₂ initially increases and then decreases during the dynamic simulation process, suggesting that NO₂ participated in the generation of other small molecules in the later stage of the reaction. Interestingly, the formation of NO is relatively late compared to other products, while the formation times of HONO, H₂O, and OH are slightly later than that of NO. Despite LLM-105 having two nitro groups, which would imply a higher production of NO₂, the amount of NO generated is actually the highest one throughout the entire simulation process. This could be attributed to two possible reasons. Firstly, during the dynamic simulation process, the initially generated NO₂ undergoes cleavage of the N–O bond to form NO. Additionally, the cleavage of the HONO group also contributes to the formation of NO.

Fig. 7 Time evolution of the numbers of the products during the decomposition of the V0, V1, V2-1, and V3-2 systems.

Fig. 8 illustrates the trend of the number of the LLM-105 molecules over time in various vacancy-containing systems within 8 ps. It is worth noting that the numbers of LLM-105 in the vacancy systems decrease more sharply compared to those in the ideal system. Furthermore, the complete decomposition of LLM-105 occurs earlier in the vacancy systems. This indicates that the presence of the vacancies accelerates the complete decomposition of LLM-105.

Fig. 8 Time evolution of the numbers of the LLM-105 molecules during the decomposition of the V0, V1, V2-1, and V3-2 systems.

Since effective collision between two atoms from free radicals and small molecules produced during the decomposition of LLM-105 is the premise of their further reaction between atoms, it is necessary to explore the diffusion of different atoms in these systems. Here, we take C, H, O, and N atoms as the research object and explore the influence of vacancies on their diffusion degree in these systems. The diffusion coefficient is obtained in two steps. The first step is to find the root mean square displacement (MSD) by molecular dynamics simulation. The second step is to find the diffusion system according to the relationship between the diffusion coefficient given in the Einstein relation and the MSD.³⁶ The diffusion coefficient (D) can be calculated from the Einstein's eqn (4):

where is the position of atom n at time t and D is the diffusion coefficient in the unit of m² s⁻¹. The MSD curve can give information about the diffusion of atoms in the systems. Fig. 9 shows the MSD curves for the four atoms, C, H, O, and N during the thermal decomposition of different systems. The MSD curves for each atom show a constant upward trend with time. This is the fact that LLM-105 has a low degree of decomposition and the motion of the atom is bound at the beginning of the dynamic simulation. However, as the simulation progresses, the molecules will decompose into atoms and free radicals, the movement of the atoms accelerates and the MSD values of the atoms become correspondingly larger. In addition, the rising trend of the MSD curves of the atoms in the vacancy systems is faster than that of the ideal system, regardless of any atoms. To obtain the diffusion coefficients of the atoms, we fit the MSD curves by the least-squares method, and the results are listed in Table 3. It is found that in the same system, the diffusion rates of the four atoms can be ranked from the largest to the smallest as H > C > N > O. Moreover, the diffusion rate of each atom in the vacancy systems are larger than that in the ideal system, which indicates that the existence of the vacancy changes the diffusion rates of the atoms in the systems.

Fig. 9 MSD curves of the C, H, N, and O atoms in the V0, V1, V2-1, and V3-2 systems.

Table 3 Diffusion coefficients (D/10⁻⁸ m² s⁻¹) of the H, C, N, and O atoms during the decomposition of the V0, V1, V2-1, and V3-2 systems

Models	V0	V1	V2-1	V3-2
H	2.04	2.17	2.23	2.33
C	1.06	1.25	1.27	1.27
N	1.20	1.30	1.51	1.33
O	1.63	1.82	1.71	1.72

---

The thermal decomposition process of LLM-105 at high temperature involves the breaking and forming of chemical bonds. In order to explore the effects of the vacancies on the chemical bonds of the LLM-105 molecules during thermal decomposition, the numbers of specific chemical bonds over time for four systems were counted as shown in Fig. 10. It can be seen that the LLM-105 molecules undergo complete decomposition within 8 ps as shown in Fig. 8. After 8 ps, the influence of the vacancies on the system can be neglected. This is because a large number of molecular products lead to an increasing disorder of the system during the dynamic simulation process. Therefore, the statistical analysis of the numbers of chemical bonds over time for these four systems was focused on the initial 8 ps. The numbers of C–C, C–N, N–H, and N–O bonds over time in the different systems are shown in Fig. 10. It is evident that the numbers of the C–C, C–N, N–H, and N–O bonds decrease during the dynamic simulation. The decrease in these bond numbers are the consequence of the cleavage of the six-membered ring, removal of the NO₂ group, hydrogen transfer, and cleavage of NO₂ to form NO or nitro–nitrite rearrangement, respectively. The numbers of these bonds decrease in the order of C–N, N–H, C–C, and N–O, suggesting that NO₂ elimination occurred first, followed by hydrogen transfer, and subsequently, six-membered ring cleavage with nitro–nitrite rearrangement and NO₂ cleavage to form NO occurring last. In addition, the decreasing trends of the C–C, C–N, N–H, and N–O bonds in the V0 system are gentler than those in the V1 system, indicating that the presence of the vacancies promotes the decomposition of LLM-105.

Fig. 10 Time evolution of the numbers of the C–C, C–N, N–H, and N–O bonds during the decomposition of theV0, V1, V2-1, and V3-2 systems.

4. Conclusions

In this study, we investigated the impacts of vacancies on the formation energies, binding energies, electronic properties, and initial thermal decomposition paths of the LLM-105 crystal using DFTB and DFTB-MD methods. The presence of the vacancies can influence the energy of the LLM-105 crystal, thereby affecting its stability. The vacancies make the VBM and CBM in the vacancy systems become localized near the vacancy site, resulting in a reduction in the band gap of the systems and facilitating its electron transitions. The dynamic simulations reveal that the LLM-105 molecules near the vacancies tend to decompose first and the presence of the vacancies does not alter the decomposition mechanisms and types of products of LLM-105. However, the vacancies can promote its complete decomposition. In addition, the presence of the vacancies alters the diffusion rate of atoms in the LLM-105 crystals during the thermal decomposition process. This work may present new ideas for understanding the behaviors of explosive crystals with defects.

Conflicts of interest

There are no conflicts to declare.

References

1. J. Zhou, L. Hao, H. Hao, X. Ji, J. Li and L. Zhou, CrystEngComm, 2021, 23, 4214–4228.
2. W. Ding, X. Zhao, Y. Li, S. Qiu and J. Huang, CrystEngComm, 2023, 25, 5641–5649.
3. Z. Liu, Q. Wu, W. Zhu and H. Xiao, Phys. Chem. Chem. Phys., 2015, 17, 10568–10578.
4. Z. H. He, J. Chen, G. F. Ji, L. M. Liu, W. J. Zhu and Q. Wu, J. Phys. Chem. B, 2015, 119, 10673–10681.
5. P. Yan, X. Zhao, J. Rui, J. Zhao, M. Xu and L. Zhai, Crystals, 2021, 11, 329.
6. C. Li, W. Yang, Q. Gan, Y. Wang, L. Liang, W. Zhang, S. Zhu and C. Feng, Def. Technol., 2024, 31, 364–374.
7. T. T. Zhou and F. L. Huang, J. Phys. Chem. B, 2011, 115, 278–287.
8. J. Ji and W. Zhu, CrystEngComm, 2022, 24, 4030–4040.
9. J. Ji, K. Wang, S. Zhu and W. Zhu, CrystEngComm, 2021, 23, 2455–2468.
10. M. Li, C. Zhang, M. Li, F. Liu, L. Zhou, Z. Gao, J. Sun, D. Han and J. Gong, Chem. Eng. J., 2022, 429, 132450.
11. R. A. Austin, N. R. Barton, W. M. Howard and L. E. Fried, J. Phys.: Conf. Ser., 2014, 500, 052002.
12. P. Ye, H. Zhai, Y. Li, W. Zhu, G. Yang and C. Guo, J. Anal. Appl. Pyrolysis, 2022, 163, 105487.
13. T. Luo, Y. Wang, H. Huang, F. Shang and X. Song, Nanomaterials, 2019, 9, 854.
14. X. Wang, W. Hu, Y. Wu and F. Huang, RSC Adv., 2019, 9, 16095–16105.
15. Q. Xiao, H. Sui, Y. Yin, L. Liao, H. Wang, X. Fu, J. Wang, J. Ai and X. Ju, Radiat. Phys. Chem., 2022, 199, 110352.
16. Y. Ma, A. Zhang, C. Zhang, D. Jiang, Y. Zhu and C. Zhang, Cryst. Growth Des., 2014, 14, 4703–4713.
17. C. Zhang, X. Wang and H. Huang, J. Am. Chem. Soc., 2008, 130, 8359–8365.
18. M. R. Manaa, I. F. W. Kuo and L. E. Fried, J. Chem. Phys., 2014, 141, 064702.
19. H. X. Ma, J. R. Song, F. Q. Zhao, H. X. Gao and R. Z. Hu, Chin. J. Chem., 2008, 26, 1997–2002.
20. G. Hang, T. Wang, W. Yu, J. Wang and H. Shen, Theor. Chem. Acc., 2022, 141, 39.
21. M. M. Kuklja and A. B. Kunz, J. Phys. Chem. B, 1999, 103, 8427–8431.
22. W. Zhu, T. Wei, W. Zhu and H. Xiao, J. Phys. Chem. A, 2008, 112, 4688–4693.
23. W. Zhu, X. Zhang, W. Zhu and H. Xiao, Phys. Chem. Chem. Phys., 2008, 10, 7318–7323.
24. B. Aradi, B. Hourahine and Th. Frauenheim, J. Phys. Chem. A, 2007, 111, 5678–5684.
25. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7260–7268.
26. M. Elstner, Theor. Chem. Acc., 2006, 116, 316–325.
27. J. G. Brandenburg and S. Grimme, J. Phys. Chem. Lett., 2014, 5, 1785–1789.
28. Z. H. He, J. Chen, G. F. Ji, L. M. Liu, W. J. Zhu and Q. Wu, J. Phys. Chem. B, 2015, 119, 10673–10681.
29. T. Qi, C. W. Bauschlicher, J. W. Lawson, T. G. Desai and E. J. Reed, J. Phys. Chem. A, 2013, 117, 11115–11125.
30. M. R. Manaa, E. J. Reed, L. E. Fried and N. Goldman, J. Am. Chem. Soc., 2009, 131, 5483–5487.
31. B. B. Averkiev, M. Yu. Antipin, I. L. Yudin and A. B. Sheremetev, J. Mol. Struct., 2002, 606, 139–146.
32. F. P. Bowden, A. D. Yoffe and G. E. Hudson, Am. J. Phys., 1952, 20, 250–251.
33. Q. Lan, H. Zhang, Y. Ni, J. Chen and H. Wang, J. Energ. Mater., 2023, 41, 269–290.
34. N. Cheng, Q. Gan, Q. Yu, X. Zhang, R. Li, S. Qian and C. Feng, Molecules, 2019, 24, 125.
35. J. Wang, Y. Xiong, H. Li and C. Zhang, J. Phys. Chem. C, 2018, 122, 1109–1118.
36. J. Topping, Wuli Tongbao, 1956, 7, 281.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ce00371c

---