First-Principles-Based Force Field for 2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105).

2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105) is a highly promising energetic material (EM) with high safety. Understanding its microscopic response mechanisms within the external stimulus is meaningful for the design of EMs. In order to comprehend the complicated phenomena, it is necessary to employ molecular simulation methods to investigate the response mechanisms with the force field (FF) at an atomic level. In this work, we developed a tailored FF for LLM-105 based on first-principles calculations. The validity of the FF was evaluated by molecular dynamics simulations. The structural parameters of LLM-105 predicted by FF are in good agreement with the experimental values, such as lattice constant, bond length, bond angle, dihedral angle and center of mass, and so forth. Moreover, the FF possesses good performance to describe the structural response on pressure accurately. In general, our work not only builds a balanced FF in gas and condensed phases, but also provides a useful tool to study the properties about LLM-105 at a large scale, which is helpful to improve the understanding about the balance between energy and safety in EMs.

Introduction

Energetic materials (EMs), including pyrotechnics, propellants, and explosives, are widely applied for civil, industrial, and military purposes.[1,2] With the progress and development of the society, the requirements for EMs increase, such as high energy density, thermal stability, thermal and n class="Disease">shock, storability, and low handling hazards.[3,4] Because of a delicate balance between high energy and high safety, a great challenge is to search a new balance point for EMs using chemical synthesis and crystal engineering technology. Among EMs, 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) and 1,1-diamino-2,2-dinitroethylene (FOX-7) are always considered to be representatives at the balance point.[5,6] Besides TATB and FOX-7, 2,6-diamino-3,5-dinitropyrazine-1-oxide (LLM-105) are a highly promising candidate. Compared with main traditional EMs like 2,4,6-trinitrotoluene, 1,3,5-trinitro-1,3,5-triazine (RDX), and octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX), LLM-105 possesses outstanding thermal stability and insensitivity to friction, spark, impact, and shock.[7]

It is interesting to reveal how to keep the balance in high energy materials with great safety. In addition, under the Hugoniot relation, which governs the detonation process, the energy, pressure, and volume of an explosive before and after a shock are involved. Prediction of the lattice structure is therefore, important to locate the detonation performances of EMs. As we know, energy release and stability maintenance of EMs should be both related to the response under external stimuli. Since the discovery of n class="Chemical">LLM-105 in the mid-1990, numerous studies have been attracting attention to thermal studies, equation of state (EOS) measurements, and sensitivity of LLM-105.[8−11] Xu et al.[8] measured the mechanical and thermal sensitivity of LLM-105. Zhang et al.[9] investigated the morphology, particle size distribution, internal structure, sensitivity, and thermal decomposition properties of submicrometer-sized LLM-105 using scanning electron microscopy and X-ray diffraction. Gump et al.[10] employed synchrotron angle-dispersive X-ray diffraction and diamond anvil cells to ascertain the isothermal equations of the state of LLM-105 at static high pressure and temperature. Williamson et al.[11] presented experiments combined with a range of diagnostic instrumentation to describe the impact response of LLM-105, demonstrating the low-level reaction in it, consisting of hotspot ignition, comparatively little light output, and substantial levels of gas production. Recently, Xu et al.[12] have studied the structural stability of LLM-105 under different pressures and temperatures by X-ray diffraction, pressure-dependent Raman and infrared spectra, indicating that a structure phase transition occurs at approximately 30 GPa. Although many experiments have been performed, it is vague to comprehend the microscopic response mechanisms within the external stimulus.

It is necessary to apply molecule simulation methods to investigate the response mechanisms at the microscopic level. Recently, multiple studies have employed density functional theory (DFT) methods to simulate the structure and EOS under high pressure.[13−17] At the B3LYP/6-31G** level, He et al.[13] demonstrated systematically the existence of the intramolecular hydrogen bonds (HBs) in n class="Chemical">LLM-105 and revealed the influence of H-bonds on its properties. Wu et al.[14] reported the effect of high-pressure on the geometric, electronic, and absorption properties of LLM-105 at the PW91 level, proposing that there exists multiple phase transitions at different pressures. Manaa et al.[15] performed dispersion-corrected DFT calculations on the unreacted EOS of crystal LLM-105 under a hydrostatic compression of up to 45 GPa. Stavrou et al.[16] implemented dispersion-corrected first-principles molecular dynamics (MD) simulations to explore the EOS and the behavior of LLM-105 at various pressures, finding that the ambient pressure phase remains stable up to 20 GPa. Moreover, Zong et al.[17] investigated the structure, mechanical properties, vibrational spectra, and EOS of LLM-105 under a hydrostatic pressure of up to 100 GPa with the PBE level.

However, the theoretical research studies mentioned are performed on the system with a small quantity of atoms. Virtually, many phenomena are always largely complicated and consist of thousands of atoms, such as shock response, thermal diffusion, hot points, and defects. Thus, the methods based on the empirical force field (FF) are usually employed to investigate the mechanisms. Generally, some general FFs were chosen to simulate the condensed phase of EMs, such as Amber[18] and Compass.[19] However, to derive more dependable results for a particular material property, the tailor-made FFs are always applied to the simulations about the specific systems. Smith and Bharadwaj constructed a FF for n class="Chemical">HMX according to the DFT energies.[20] Mathew and Sewell used MD with a tailor-made FF to investigate the surface-initiated melting of TATB.[21] With the combination of symmetry-adapted perturbation theory and DFT [SAPT(DFT)],[22−26] a series of intermolecular interaction potential have been parameterized for RDX,[27] TATB,[28] and FOX-7.[29] Furthermore, Song et al.[30] developed all-atom, nonempirical, and tailor-made FFs for α-RDX, which can be reliably applied into crystal and density predictions and the sensitivity and detonation estimations. Because of its good accuracy, it is successful in describing the phenomenon for the respective molecules. However, it is not found that molecule-tailored FFs for LLM-105 have been established.

While the prediction of the lattice structures of EMs remains challenging, the new FF is helpful for studying the behaviors of neat LLM-105 under high pressures, which is crucial for the prediction of its detonation performance. In this work, we developed a tailored FF for n class="Chemical">LLM-105 based on first-principles calculations. After parameterization has been carried out step by step, the validity of the FF has been evaluated by the MD simulations. The FF possesses good performance to describe the LLM-105 molecular and crystal structures and present structural response on pressure. Our work provides not only a tailored FF for LLM-105, but also a powerful tool to explore the response mechanisms of LLM-105 at the atomic scale, which is helpful to improve the understanding about the balance between energy and safety in EMs.

Computational Methods

Force Field

LLM-105, with four formula units per unit cell, is a highly dense molecular solid that crystallizes in a monoclinic, P21/n space symmetry group. It exhibits an extensive network of HBs and π–π stacking in its crystal structure as shown in Figure a, which was drawn by n class="Chemical">mercury 3.9.[31] The topology of the LLM-105 molecule has been displayed in Figure b, which possesses only one stable conformation in the gas phase. The corresponding atom types are presented in Figure c. It is noted that a few atom types are applied for the FF construction.

Figure 1

(a) LLM-105 crystal structure and intermolar hydrogen bond in the crystal; the atom indexes (b) and atomic type (c) in the force field parameters.

(a) LLM-105 crystal structure and intermolar n class="Chemical">hydrogen bond in the crystal; the atom indexes (b) and atomic type (c) in the force field parameters.

The initial molecular structure of LLM-105 (Cambridge Structural Database: refcode YEKQAG03) was obtained from the work by Liu et al.[32] After optimization with the PBE0/aug-cc-pVDZ[33,34] level, more than 1000 dimer structures were generated randomly referring to the method introduced by refs (27) and (28). To obtain the accurate intermolecular energy, the combination with PBE0/aug-cc-pVDZ and exchange-hole dipole moment dispersion correction (n class="Chemical">XDM)[35−39] were employed. The XDM dispersion model has demonstrated good performance in sublimation enthalpies and absolute lattice energies of molecular crystals.[40] Moreover, some studies[41−43] obtained accurate energy ranking in the crystal-structure prediction of molecular crystals with the DFT–XDM method. As a whole, the FF form is similar with Amber FF.[44] The nonbonding interaction is expressed with the electrostatic and van der Waals (vdW) terms. As shown in eq , the point charge model and the vdW potential were chosen to describe the two terms, respectively. Instead of the Lennard-Jones (LJ) potential in Amber, the vdW interaction is expressed using the Buckingham (BKH) potential. This FF form has been applied to describe HMX by Smith and Bharadwaj.[20]where, q is the charge on atom i, and r is the distance between atom i and j, A, B, and C are the parameters in the BKH potential. A revised electrostatic potential fitting charge (RESP) scheme was used to obtain the point charge.[45] RESP charges always show good agreement with the electrostatic component of the ab initio interaction energies. BKH has more reasonable short-range interactions and is more suitable to describe the EMs in high pressure. In order to obtain reasonable BKH parameters, the FIT potential parameters[46,47] were chosen for the initial values directly, and the fitting bounds were confined to be 0.9 and 1.1 times of the parameter. In addition, the least-square method was used to optimize the parameter referring to the dimer energies in the condition of constrained charge.

For bonding interactions, a series of perturbed structures were generated based on one of the equilibrium conformations according to the redundant internal coordinates, and optimized with the constrained structural parameters at the PBE0/aug-cc-pVDZ level, and the harmonic style functions were chosen for bond and angle terms (Ebond and Eangle), as followwhere kb and kθ are the force constants of bond and angle terms, b and b0 are the actual and balanced bond lengths, while θ and θ0 are the actual and balanced bond angles. A fourier style potential was chosen to describe dihedral terms (Edihedral)where kϕ, n, ϕ, and d are the force constants, scale factor, and actual and balanced dihedral angels. To fix the nitro group planer, a harmonic improper term was chosen for n class="Chemical">nitro groups (Eimproper)where kψ and ψ are the force constants and improper angle, and the initial bonding parameters were obtained from Amber FF parameters. In order to get reliable results, the bounds of parameters were set to be 0.8–1.2 of the initial parameters, and the evolutionary algorithm was employed to optimize the parameters with minimization of the following score function (Escore)where K is the slope determined by the linear fitting between the referenced DFT energy Eref and FF calculated energy EFF with the optimizations at several structural parameter constraints. The Escore can be applied to balance between the slope and root mean square error. It is worth noting that all improper force constants kψ were fixed at 1.12 kcal mol–1 rad–2 during the parameterization. In all DFT/MM optimization and energy calculation were performed with Gaussian 09[48] and LAMMPS[49] packages.

MD Simulation

In our MD simulations, the calculations were carried out with the NPT ensemble. The time step length and cutoff distance of vdW interaction were set to 0.5 fs and 12 Å, respectively. In order to ensure that all the edge lengths of the cell were at the least twice cutoff distance, a bulk supercell of a 5 × 2 × 3 crystallographic unit cell was used. All structure relaxations were performed for 100 ps from 0 to 20 GPa at 300 K. After equilibration at each pressure, atomic configurations were recorded every 1000th step, and the last 20 000 step trajectories were gathered. By statistical average, we obtained the mean lattice parameters, atom coordinates, volume, center of mass, etc. In addition, a gaseous molecular structure was also optimized with the FF. All calculations in this part were implemented in the LAMMPS program.

Results

Geometry in the Gaseous Phase

The LLM-105 molecule is a planar structure with C point group symmetry in the gaseous phase. Thus, atom types were defined referring to the chemical groups and their connections. Figure a shows the structural overlap between FF and DFT. Overall, the structure from the FF is very similar with that from DFT. In detail, Table compares the structural parameters optimized at the PBE0/aug-cc-pVDZ and fitted FF levels. Because of the planar structure, dihedral and improper angles are not included. Obviously, the parameters from the FF are in good agreement with those from DFT. The mean error (ME) of bond lengths is about 0.006 Å, and the ME of bond angles is about 0.1°. The root mean squared error (RMSE) of bond lengths and angels are 0.015 Å and 1.0°. The maximum error (MAXE) about bond lengths existing in the bond between C1 and N12 is about 0.20 Å. Among bond angles, the MAXE is about 1.9° in the angle C2–N3–C4.

Figure 2

Overlap for molecules from DFT and FF in the gaseous phase (a); a front (b) and a lateral view (c) from the experiment and FF in the crystal form.

Table 1

Structural Parameters Optimized at the DFT and FF Levels

structural parameters	DFT	FF	Δ(FF–DFT)
Bond Lengths (Å)
B(C1–C2)	1.419	1.402	–0.018
B(C2–N3)	1.308	1.323	0.016
B(C1–N12)	1.371	1.391	0.020
B(C1–N14)	1.325	1.345	0.020
B(C2–N6)	1.459	1.464	0.005
B(N14–H19)	1.015	1.025	0.010
B(N14–H18)	1.010	1.022	0.013
B(N6–O8)	1.233	1.224	–0.009
B(N6–O7)	1.206	1.221	0.015
B(N12–O13)	1.290	1.275	–0.015
ME			0.006
RMSE			0.015
Angles (deg)
A(H18–N14–H19)	126.0	127.8	1.8
A(C2–C1–N12)	116.2	116.1	–0.1
A(C1–N12–C5)	122.0	122.5	0.5
A(N3–C2–N6)	116.3	116.8	0.5
A(C2–C1–N14)	129.5	130.2	0.7
A(C1–N14–H19)	115.4	113.9	–1.5
A(C1–N14–H18)	118.6	118.3	–0.3
A(C2–N3–C4)	119.5	121.3	1.9
A(N12–C1–N14)	114.3	113.7	–0.6
A(C1–C2–N3)	123.1	122.0	–1.1
A(O7–N6–O8)	125.0	123.3	–1.6
A(C2–N6–O7)	118.5	119.3	0.8
A(C2–N6–O8)	116.5	117.3	0.8
A(C1–C2–N6)	120.6	121.2	0.6
A(C1–N12–O13)	119.0	118.8	–0.3
ME			0.1
RMSE			1.0

Lattice Structure

After the 200 000 step NPT simulation, the equilibrium crystal structure was obtained with averaging lattice parameters and atom coordinates in last 20 000 step trajectories. Figure shows the overlap between simulated and experimental structures. In the 15/15 similarity standard, the displacement between two structures is calculated be 0.16 Å by pymatgen,[50] which reflects that the structure from the FF is very similar with that from the experiment. Table compares the lattice parameters based on DFT,[14,16,17] FF simulations, and experimental measure.[10,15,16,32,51,52] Compared to the experimental and DFT calculated values, our FF gives a slightly smaller lattice cell. The lattice parameters were predicted to be 5.64, 15.55, and 8.24 Å for a, b, and c axis, respectively. The lengths of a, b, and c decreases by about 0.08, 0.30, and 0.18 Å, or 1.4, 1.9, and 2.1%, respectively, and the lattice angle β varies little, which is predicted to be 100.96°. Because of the little shorter axis, the lattice volume is reduced from 748 Å3 to 711 Å3, by about 5.0%. Although our FF gives a smaller volume, the molecular position and orientation are in good agreement with the experimental measurements.[32] Because of the symmetry, the position is described with the fraction coordinate of the centroid of a chosen molecule, and the orientation is expressed with three angles between x, y, and z axis and the vector from n class="Chemical">N3 to N4 on the molecule. In the result from the FF, the position is predicted to be (0.560, 0.869, and 0.549) for the molecule, which are very similar with (0.559, 0.867, and 0.553) from the experiment.[32] Also the orientation angles are simulated to be (145.65, 67.38, and 114.38°), which are consistent with those from experimental measures.[32]

Figure 3

Overlap between simulated (normal color) and experimental structures (green).

Table 2

Structural Parameters, Bulk Modulus and Its Derivative of the LLM-105 Crystal from FF, DFT, and Experiments

	FF	LDA[14]	PW91[14]	PBE[17]	PBE-D2[16]	EXP[10,15,16,51,52]
a (Å)	5.64	5.84	6.01	5.64	5.81	5.71–5.75
b (Å)	15.55	15.58	18.28	15.96	16.05	15.60–15.87
c (Å)	8.24	8.22	8.71	8.46	8.40	8.41–8.48
β (deg)	100.96	99.51	100.75	100.93		101.01–101.22
V0 (Å3)	711.00	737.24	939.20	747.38	769.51	746.21–750.08
center of mass	0.560					0.559
	0.869					0.867
	0.549					0.553
orient (deg)	145.65					145.65
	67.38					67.38
	114.38					114.38
B0	13.8			16.5	12.7	11.19–16.2
B′	11.7			9.4	9.4	8.30–18.54

Further, the charges in the molecular structure between those predicted by the FF and experimental measure[32] are focused on. Because of crystal packing and HB formation, the symmetry on molecular structures has been broken, and the molecule is not kept with a planar structure longer. Figure b,c also gives a front view and a side lateral view of the overlap structure. The structure from the FF is consistent with that from the experiment.[32]Table lists the main deformation which is found around the interaction sites such as nitro and amino groups. Apparently, the FF can also be applied to give similar molecular structural parameters with the experimental values.[32] The dihedral angle of C1–C2–N6–O8 varied from 0.00° in gaseous to 18.82° in crystal phase, and the FF gave an accurate prediction about the dihedral angle with 18.82° in the crystal. The other dihedral angle (D(N3–C4–N9–O11)) about C–NO2 is estimated to be −3.48° by the FF and in good agreement with the experimental value.[32] After the crystal packing, all N–H bonds are shortened to 0.860 Å from 1.010 Å, and the angles of n class="Chemical">nitrogen atoms in amino groups increase to about 120°. Also the bond N12–O13 with the last hydrogen bond O13 bond is lengthened from 1.275 to 1.321 Å. The experimental values[32] of B(N–H), A(H–N–H), and B(N12–O13) are 0.860 Å, 120.01°, and 1.321 Å, respectively.

Table 3

Structural Parameters Optimized at the FF Levels and Experimental Values

structural parameter	EXP	FF	Δ(FF-EXP)
Bond Lengths (Å)
B(N15–H16)	0.860	0.860	0.000
B(N15–H17)	0.860	0.860	0.000
B(N14–H18)	0.860	0.860	0.000
B(N14–H19)	0.860	0.860	0.000
B(N12–O13)	1.321	1.321	0.000
Bond Angles (deg)
A(C1–N14–H18)	120.01	120.01	–0.01
A(C1–N14–H19)	119.98	119.98	0.00
A(H18–N14–H19)	120.01	120.01	0.01
A(C5–N15–H16)	120.00	120.01	0.00
A(C5–N15–H17)	120.02	120.01	0.00
A(H16–N15–H17)	119.98	119.98	0.00
Dihedral Angles (deg)
D(C1–C2–N6–O8)	18.82	18.82	0.00
D(N3–C2–N6–O8)	–159.43	–159.43	0.00
D(C5–C4–N9–O10)	–2.67	–2.67	0.00
D(N3–C4–N9–O11)	–3.48	–3.48	0.00

Because of less charge in the molecule, the cell difference between the simulation and experiment should be attributed to the description about the intermolecular interaction. Figure shows four main dimer structures with the great intermolecular interaction. Figure a,b is typical π–π stacking dimers, called π–π1 and π–π2, respectively. In addition, Figure c,d shows doublet HB systems, marked with HB1 and HB2. Three structural parameters are applied to compare the simulated and experimental π–π stacking dimers:[32] the distance between the two geometry centers of the six-member rings (Rgc), the distance between the two π-planes constructed by the six-member rings (Rpl), and the angle between the normal vector of a plane and the vector from the center to the other (Anc). Overall, few differences can be found on Anc, while both distances are shortened slightly. In π–π1, the experimental Rgc and Rpl are 4.113 and 3.212 Å, while the simulated ones are 3.960 and 3.096 Å. The Rgc and Rpl were decreased by 0.153 and 0.116 Å, respectively. For π–π2 in Figure b, the differences on Rgc and Rpl were less with 0.102 and 0.052 Å. Compared with the structures with π–π stacking, the HB dimers have slightly larger structural charges as shown in Figure c,d. Here, the length (RHB) and angle (AHB) of the HB were applied to describe the HB dimer. The RHB between H18(A) and O8(B) is shorten by 0.141 Å and that between O13(A) and H19(B) is decreased by 0.201 Å. However, the AHB of N14–H19–O13 and N14–H18–O8 varied little. Thus, one can conclude that shorter equilibrium distances between molecules should contribute to the smaller lattice crystal.

Figure 4

Four main dimer structures in intermolecular interaction energy: (a,b) dimers with π–π stacking; (c,d) dimers with hydrogen bonds. Estimated values are outside of parentheses, and experimental values in the parentheses.

Although the FF tends to give a smaller lattice volume, we test the FF estimation on pressure response of the crystal structure under pressure in the range from 0 to 20 GPa. Additionally a series of NPT simulations were performed to obtain the equilibrium crystal structure at the different pressures. Using the lattice parameters, the density and relative density of LLM-105 under various pressures are derived as shown in Figure S4. The density increases rapidly under low pressures (<3 GPa) and increases slowly but in a nearly linear way under high pressures. The predicted density and relative density are more or less greater than the measured ones but the discrepancies are less than 0.2 g/cm3 or 3.4%. The pressure-dependent relative cell volumes and lattice parameters are shown in Figure a,b, respectively, together with the results of the experimental measurement and theoretical calculations. It is evident that the results agree well with the experimental data and DFT calculations.[16] Similar with the experiment, the b-axis is more compressible than a and c axes when pressure is applied. The EOS is critical in describing the properties of materials. Thus, we conducted unweighted n class="Disease">fits of the pressure–volume data using a third-order Birch–Murnaghan (B–M) EOS as followwhere, B0 is the bulk modulus and B′ is the pressure derivative of B0, V is the volume corresponding the pressure, and V0 is the equilibrium volume at zero pressure. The fitted B0 and B′ from the FF are 13.8 GPa and 11.7, respectively, which are close to the experimental values of 15 GPa and 9 in ref (16) and 14.6 GPa and 10.6 in ref (15). The B0 and B′ are also summarized in Table and compared to the experimental and calculated values. Apparently, the calculated results are in agreement with other experimental[10,15,16] or other theoretical values.[16,17] Compared with DFT estimations, the FF seem give a better prediction on B0 and B′.

Figure 5

Pressure-dependent relative cell volumes (a) and lattice parameters (b) from the FF, experiment, and DFT.

Discussion

The packing features of a crystal can be expressed using a cluster. As a crystal with P21/n and Z′ = 1, the cluster contains a center LLM-105 molecule and its 14 closed neighbors, which form 14 dimers as shown in Figure . Further, 14 dimers can be divided into eight groups according to the structure similarity. In Figure , the molecules with same color means the same dimer can be constructed with one of them and the center molecule. Table lists the intermolecular interaction energies of the dimers and the corresponding energy decomposition with CE-B3LYP[53] based on the experimental structure. In the energy decomposition scheme, the interaction energy (Eint) is expressed in terms of electrostatic (Eele), polarization (Epol), dispersion (Edisp), and exchange-repulsion (Erep) as the following equationwhere, ki (i = es, pol, disp, and rep) is the fitted scale factor for the corresponding energy component. Because of the lack of explicit physical meanings, total energies (Etot) without any scale factors were also provided in Table . Overall, π–π stacking and the HB should dominate the crystal packing. Among all dimers, π–π stacking interactions are slightly stronger than those of the HB in either Etot or Eint. π–π1, π–π2, HB1, and HB2 in Figure possess the greatest Eint. Although the number of π–π stacking dimers is less than HBs, the overall Etot of π–π stacking is similar to that of HBs because of their greater strength. From the energy decomposition, it is found that the interaction is contributed mainly by the electrostatic and dispersion. The dispersion effect plays a main role in π–π stacking, while the electrostatic term governs HB interaction.

Figure 6

LLM-105 molecule and its close neighbors. The center molecule and the same color neighbor can construct the same dimer in geometry. 14 neighbors can be divided into eight groups.

Table 4

Intermolecular Interaction Energies of the Neighbor Dimers and the Corresponding Energy Decomposion with CE-B3LYP (kcal/mol)

To analyze the discerption on intermolecular interaction further, the potential energy curves about the four main dimers are represented in Figure . As a whole, the FF can be used to give a similar description on no-bonding interaction with that obtained at the PBE0-XDM/aug-cc-pVDZ level. Except HB2, the equilibrium distances and well depths are in great agreement with those by DFT. The curves of π–π2 from FF and DFT are nearly overlapping, and the equilibrium distance of π–π2 should be shorter by less than 0.1 Å, and the well depth of HB1 are increased by less than 0.6 kcal/mol. Although our FF can give three reasonable curves, the greater error can be found on HB2. Obviously, the equilibrium distance is decreased by about 0.2 Å, which is very similar with that in the crystal. Despite this, the curve from the FF possesses a good shape like that from DFT as shown in Figure d. This might be why our FF can give a good bulk modulus at the same time of obtaining a smaller volume.

Figure 7

Potential energy curves from the FF and DFT for the four main dimers: (a,b) the dimers with π–π stacking; (c,d) dimers with hydrogen bonds.

Moreover, the interaction with π–π stacking is stronger than that with HBs, which agrees with results from CE-B3LYP. In π–π1 and 2, the interaction energies with PBE0-XDM are estimated to be −11.53 and −9.28 kcal/mol, respectively, greater than the energy with the doublet HBs of −8.05 kcal/mol. The greater interaction energy prefers a greater weight during parameterization, which might make the no-bonding parameters more suitable to describe the π–π interaction. Because of the dominance by the electrostatic effect, the description accuracy of HB interaction should be dependent on the atomic charges. The atomic charges in the FF were obtained referring to the gaseous molecular structure rather than the structure in the crystal, and the effect of the crystal environment on the electrostatic potential were ignored. Moreover, few atomic types make a charge average. All of these might be responsible for the poorer description on HB interaction. Nevertheless, the parameters of the FF are obtained mainly from the unimolecular properties of n class="Chemical">LLM-105 calculated by DFT, which should describe LLM-105 systems without any preference. Overall, the FF with a general FF form can be applied to investigate the properties of LLM-105.

Conclusions

We have developed an all-atom FF for LLM-105 compounds. As a remarkable energetic crystal, its fundamental properties have been focused on using several ab initio methods but the FF. Considering the compatibility of the force filed, a modified Amber function form was chosen for describing the interaction in n class="Chemical">LLM-105, in which the BHK vdW potential is used instead of the LJ potential. The electrostatic energy was represented with the interaction among the RESP charges, and vdW parameters were fitted using the potential energies of more than 1000 dimers estimated at the PBE0/aug-cc-pVDZ level with the XDM dispersion correction. Then, the bonding parameters were fitted from a series of perturbation structures based on the equilibrium structure in the gas phase.

By the MD simulations, the molecular and crystal structures were relaxed, and the pressure response of the lattice was predicted. All results are in good agreement with the experimental measures. The molecular structure from the FF is consistent with that from DFT. Although the FF tends to give smaller lattice parameters, the structure, orient, and relative position of the molecules in the crystal concur with the experiment. Moreover, the bulk modulus and its pressure derivative are predicted to be 13.8 GPa and 11.7, which are closer to the experimental measure than several DFT results. Moreover, the packing features of the LLM-105 crystal were analyzed with the intermolecular interaction using an energy decomposition scheme. In n class="Chemical">LLM-105, the intermolecular interaction is contributed mainly by the electrostatic and dispersion effects. The π–π stacking and HB dominate the crystal packing, and π–π stacking dimers have greater strength than HBs. We found that the FF will give a better equilibrium distance of the π–π system, but underestimate slightly on the bond length of the HB. Because of the underestimation, the FF tends to give a decreased volume by about 5%. The developed FF is used to simulate the behaviors of neat LLM-105 under a shock. To predict its detonation parameters, other FFs have to be developed for the reactions and resulting products. In general, we provide a reliable FF tool for the complicated phenomenon of LLM-105 at a large scale.