Noncovalent Interactions and Crystal Structure Prediction of Energetic Materials

The crystal and molecular structures, intermolecular interactions, and energy of CL-20, HATO, and FOX-7 were comparatively predicted based on molecular dynamic (MD) simulations. By comparison, the 2D fingerprint plot, Hirshfeld surface, reduced density gradient isosurface, and electrostatic potential surface were studied to detect the intermolecular interactions. Meanwhile, the effects of vacuum and different solvents on the crystal habit of CL-20, HATO, and FOX-7 were studied by AE and MAE model, respectively. The energy calculation was also analysed based on the equilibrium structures of these crystal models by MD simulations. Our results would provide fundamental insights for the crystal engineering of energetic materials.


Introduction
Energetic materials (EMs), mainly including explosives, propellants, and pyrotechnics are a kind of important compounds that can store and release considerable energy, and have occupied an important place in mining, military equipment, space exploration, and fireworks [1][2][3][4][5][6].
Dihydroxylammonium 5,5 -bistetrazole-1,1 -diolate (HATO or TKX-50) is a newly synthesized explosive with excellent comprehensive properties, which possesses lower impact sensitivity and higher detonation performance than the common explosives. Apart from the higher energy level and lower sensitivity, low toxicity is another advantage of HATO. However, poor crystal habit, irregular crystal morphology, and smaller crystal particle make an extremely constraint in its application [9][10][11].
The 1,1-diamino-2,2-dinitroethene (FOX-7) is synthesized in 1998 and is a push-pull symmetrical structure. There is a nitro-enamine structure with two nitro group at one carbon and two amino group at other in its molecular structure. Due to the special chemistry structure, FOX-7 possesses an excellent detonation performances and safety property, and is honoured as outstanding low sensitivity and highly energetic materials. However, poor crystal quality, monotonous crystal particle, more crystal defect, and poor solubility make its widespread application enormous challenging [12][13][14][15].
To meet the high requirement, the study on Ems with high crystal quality has always been an ultimate goal that researchers strive for. The morphology of energetic materials is a vital component and plays an important role in the crystal quality. It can be To meet the high requirement, the study on Ems with high crystal quality has always been an ultimate goal that researchers strive for. The morphology of energetic materials is a vital component and plays an important role in the crystal quality. It can be affected by intermolecular interaction, solvent, additive, and temperature in the process of crystallization [16,17]. Nevertheless, the crystallization mechanism of CL-20, HATO, and FOX-7 remain unclear, which significantly motivates us to probe the crystal habit in three explosives.
Herein, CL-20, HATO, and FOX-7 were chosen for comparatively simulating their structures and predicting their properties by quantum chemistry calculation. By comparison, the intermolecular interactions of CL-20, HATO, and FOX-7 were studied by twodimensional (2D) fingerprint plots, Hirshfeld surfaces, reduced density gradient (RDG) isosurface, and electrostatic potential surface analysis. Then, the effects of different solvents and additives on the crystal structures of CL-20, HATO, and FOX-7 were evaluated. This work aims to deepen the understanding of crystal growth and provide a better insight for the crystallization process of energetic material.
The chemical diagrams of CL-20, HATO, and FOX-7 are shown in Figure 1.

Hirshfeld Surface Analysis
The surface and 2D fingerprint plot, which was formed by the location of (di de) point and their relative frequency, can confirm distance and intensity of intermolecular interaction [18,19]. The 2D fingerprints and the associated Hirshfeld surfaces of CL-20, HATO and FOX-7 were shown in Figure 2, in which the intermolecular interaction was visualized.

Hirshfeld Surface Analysis
The surface and 2D fingerprint plot, which was formed by the location of (d i d e ) point and their relative frequency, can confirm distance and intensity of intermolecular interaction [18,19]. The 2D fingerprints and the associated Hirshfeld surfaces of CL-20, HATO and FOX-7 were shown in Figure 2, in which the intermolecular interaction was visualized.
As shown in Figure 2a-c, there are the Hirshfeld surfaces of CL-20, HATO and FOX-7, respectively. The shape of surface of HATO and FOX-7 is plate, which indicates an insensitivity energetic material than CL-20. Furthermore, there are red, white, and blue dots on the surface, which represents the close contact around molecules. The red dot represents the distance d less the Van der Waals, white dot represents equal and blue dot represents exceeds. Compared with HATO and FOX-7, more red dots, which represent the predominant intermolecular interaction, are located on the surface of CL-20, particularly at edges, indicating that high close contact around CL-20 molecule.
In Figure 2d-f, we can see the intermolecular interaction in 2D fingerprints plots of crystals of CL-20, HATO, and FOX-7, respectively. We can see that there are O···H, O···N, and O···O intermolecular interactions in both of CL-20 and FOX-7 crystal. Compared with CL-20 and FOX-7, intermolecular interactions are O···H, O···N, and O···O in HATO. In addition, we can know that there are O···H intermolecular interactions in all of CL-20, HATO, and FOX-7, and it is a pair of remarkable spikes on the bottom left in 2D fingerprints.

Hirshfeld Surface Analysis
The surface and 2D fingerprint plot, which was formed by the location of (di de) point and their relative frequency, can confirm distance and intensity of intermolecular interaction [18,19]. The 2D fingerprints and the associated Hirshfeld surfaces of CL-20, HATO and FOX-7 were shown in Figure 2, in which the intermolecular interaction was visualized.  As shown in Figure 2a-c, there are the Hirshfeld surfaces of CL-20, HATO and FOX-7, respectively. The shape of surface of HATO and FOX-7 is plate, which indicates an insensitivity energetic material than CL-20. Furthermore, there are red, white, and blue dots on the surface, which represents the close contact around molecules. The red dot represents the distance d less the Van der Waals, white dot represents equal and blue dot represents exceeds. Compared with HATO and FOX-7, more red dots, which represent the predominant intermolecular interaction, are located on the surface of CL-20, particularly at edges, indicating that high close contact around CL-20 molecule.
In Figure 2d-f, we can see the intermolecular interaction in 2D fingerprints plots of crystals of CL-20, HATO, and FOX-7, respectively. We can see that there are O···H, O···N, and O···O intermolecular interactions in both of CL-20 and FOX-7 crystal. Compared with CL-20 and FOX-7, intermolecular interactions are O···H, O···N, and O···O in HATO. In addition, we can know that there are O···H intermolecular interactions in all of CL-20, HATO, and FOX-7, and it is a pair of remarkable spikes on the bottom left in 2D fingerprints.

Reduced Density Gradient Analysis
To further reveal the information on inter-and intra-molecular effects and comprehensively study their influence on crystal structure, on the basis of the electron density and its RDG, the NCI plots of CL-20, HATO, and FOX-7 were determined in real space. In this method, the differences between hydrogen bonds, Van der Waals interactions, and repulsive steric clashes were observed by the study of the relationship between quantummechanical electron density and the reduced density gradient [20][21][22][23]. As shown in Figure  3, it is clear that there are abundant elliptical green slabs around the nitro group of each molecular edge, which denotes that weak hydrogen bonds and strong Van der Waals interactions exist in CL-20, HATO, and FOX-7. In addition to this, we can see that there are

Reduced Density Gradient Analysis
To further reveal the information on inter-and intra-molecular effects and comprehensively study their influence on crystal structure, on the basis of the electron density and its RDG, the NCI plots of CL-20, HATO, and FOX-7 were determined in real space. In this method, the differences between hydrogen bonds, Van der Waals interactions, and repulsive steric clashes were observed by the study of the relationship between quantum-mechanical electron density and the reduced density gradient [20][21][22][23]. As shown in Figure 3, it is clear that there are abundant elliptical green slabs around the nitro group of each molecular edge, which denotes that weak hydrogen bonds and strong Van der Waals interactions exist in CL-20, HATO, and FOX-7. In addition to this, we can see that there are red slabs among five-membered ring of CL-20 and HATO, which indicates the steric effect in CL-20 and HATO. Due to more red slabs, CL-20 possess more greater steric effect than HATO. Due to red region only on the edge of elliptical slabs, FOX-7 possess the weakest steric effect in three molecule structure.

Electrostatic Potential Surface Analysis
Electrostatic potential surfaces (ESP) illustrate the necessary energy that a unit positive charge moves from infinity to a certain position in a chemical system [24]. It clearly shows the charge density distribution and the charge variation region of molecule, which can be used to determine the electrophilicity and nucleophilicity at different positions in the system, and also to reveal the properties of non-covalent molecular interaction [25,26]. Hence, the ESP was carried out in order to understand the electronic property and the bond strength variation in CL-20, HATO, and FOX-7. The ESP mapped surfaces of CL-20, HATO, and FOX-7 were presented in Figure 4.

Electrostatic Potential Surface Analysis
Electrostatic potential surfaces (ESP) illustrate the necessary energy that a unit positive charge moves from infinity to a certain position in a chemical system [24]. It clearly shows the charge density distribution and the charge variation region of molecule, which can be used to determine the electrophilicity and nucleophilicity at different positions in the system, and also to reveal the properties of non-covalent molecular interaction [25,26]. Hence, the ESP was carried out in order to understand the electronic property and the bond strength variation in CL-20, HATO, and FOX-7. The ESP mapped surfaces of CL-20, HATO, and FOX-7 were presented in Figure 4.

Electrostatic Potential Surface Analysis
Electrostatic potential surfaces (ESP) illustrate the necessary energy that a unit positive charge moves from infinity to a certain position in a chemical system [24]. It clearly shows the charge density distribution and the charge variation region of molecule, which can be used to determine the electrophilicity and nucleophilicity at different positions in the system, and also to reveal the properties of non-covalent molecular interaction [25,26]. Hence, the ESP was carried out in order to understand the electronic property and the bond strength variation in CL-20, HATO, and FOX-7. The ESP mapped surfaces of CL-20, HATO, and FOX-7 were presented in Figure 4.

CL-20
HATO FOX-7  In Figure 4, it can be seen that the surface regions of ESP are coloured by different colours and the sites of maxima and minima are depicted by orange and cyan spheres, respectively. The blue region indicates negative ESP, which is the attraction of the proton by the concentrated electron density [27,28]. The blue region unexceptionally distributes on the edges of the molecules. Particularly, it situates on -NO 2 group (in CL-20 and FOX-7) and -NO group (in HATO). The red region indicates positive ESP, which is the repulsion of the proton by the atomic nuclei in regions where low electron density exists and the nuclear charge is incompletely shielded [27,28]. The red region distributes on the edge of the -NH group (in HATO and FOX-7) and the central regions of the CL-20 molecule. When the polar of ESP is the same on the solvent and crystal face, they are manifested as a repulsive interaction, while the attractive interaction exists between the solvent and crystal face, when they are the opposite polar. The ESP qualitatively expresses the charged areas of molecules and the site of molecular reactively [29].

Crystal Morphology in Vacuum and Dominant Crystal Facets Analysis
The morphologically important facets of CL-20 crystal were simulated by the growth morphology algorithm (see Figure 5 and Table 1). This process was calculated based on the AE model, which was proposed by Hartman and Bennema according to the PBC theory in 1980. The attachment energy (E att ) describes an amount of energy release when a crystal slice of thickness d hkl attached to the (hkl) facet, and this model does not consider the solvent media, temperature, pressure, and other factors during the simulation process. It can be seen from Figure 6   In Figure 4, it can be seen that the surface regions of ESP are coloured by different colours and the sites of maxima and minima are depicted by orange and cyan spheres, respectively. The blue region indicates negative ESP, which is the attraction of the proton by the concentrated electron density [27,28]. The blue region unexceptionally distributes on the edges of the molecules. Particularly, it situates on -NO2 group (in CL-20 and FOX-7) and -NO group (in HATO). The red region indicates positive ESP, which is the repulsion of the proton by the atomic nuclei in regions where low electron density exists and the nuclear charge is incompletely shielded [27,28]. The red region distributes on the edge of the -NH group (in HATO and FOX-7) and the central regions of the CL-20 molecule. When the polar of ESP is the same on the solvent and crystal face, they are manifested as a repulsive interaction, while the attractive interaction exists between the solvent and crystal face, when they are the opposite polar. The ESP qualitatively expresses the charged areas of molecules and the site of molecular reactively [29].

Crystal Morphology in Vacuum and Dominant Crystal Facets Analysis
The morphologically important facets of CL-20 crystal were simulated by the growth morphology algorithm (see Figure 5 and Table 1). This process was calculated based on the AE model, which was proposed by Hartman and Bennema according to the PBC theory in 1980. The attachment energy (Eatt) describes an amount of energy release when a crystal slice of thickness dhkl attached to the (hkl) facet, and this model does not consider the solvent media, temperature, pressure, and other factors during the simulation process. It can be seen from Figure 6       The morphologically important facets of HATO (see Figure 5, Table 2) and FOX-7 (see Figure 5, Table 3) crystals can be seen from the below figures and tables. In Figure 6, we can know that with an aspect ratio of 1.88, the predicted morphology of HATO crystal is mainly dominated by the {1 0 0}, {1 1 0}, {0 2 0}, {0 1 1}, and {1 1 −1} crystal facets. In Table 2 Figure 6). In Table 3, the largest area percentage of FOX-7 is {0 1 1} facet, and this facet is the most morphological important in FOX-7 crystal morphology. With a larger interplanar spacing {d hkl }, {1 1 0} facet is the slow-growing facet whereas other facets occupy smaller {d hkl } are the fast-growing facets. Apart from, there is E att mainly composed with electrostatic interaction, Van der Waals force and hydrogen bonding interaction in three crystals. In CL-20 and FOX-7, the E att is mainly involved Van der Waals force and electrostatic interaction, whereas, in HATO, the E att is mainly involved in Van der Waals force and hydrogen bonding interaction.

Effect of Ethyl Acetate, Acetonitrile, Acetone, DMF, DMSO, and NMP
The interaction of solvents in the morphologically important crystal facets was simulated by the molecular dynamic method (shown in Tables 4-6). The interaction energy (E int ) is the energy between each crystal face and solvent molecule which is negative value in three crystals, indicating that the adsorption of solvent in the crystal face was exothermic, during which the solvent will spontaneously adhere to the crystal interface [30]. In Tables 4-6, we can see that different solvents have different effects on each crystal facet, which leads to the change of the importance of each crystal facet during the process of crystal growth. The crystal morphologies of CL-20, HATO, and FOX-7 in different solvents predicted and experimented were shown in Figures 6-8, respectively. Compared with the morphology in vacuum, the crystal morphology predicted by MAE model in solvent is very obviously different.        In Table 4  In Figure 6, we can see that the morphology of CL-20 in ethyl acetate and DMF is a column, whereas, in acetonitrile and acetone the morphology tends to a rhombus. Meanwhile, in Figure 7, the morphology of HATO is rodlike in DMSO and NMP. In Figure 8, the morphology of FOX-7 is irregular polyhedron in NMP, DMSO, and NMP. Moreover, we can know that the morphology predicted of three crystals in different solvents is quite consistent with the corresponding experimental results.

Hirshfeld Surface
The Hirshfeld surface (HS) and corresponding 2D fingerprints were calculated by CrystalExplorer [31]. In CrystalExplorer, the electron density, which comes from tabulations of atomic wavefunctions, was used to compute the HS. The pro-molecule electron density ρ promol (r) and the pro-crystal electron density ρ procryst (r) were constructed based on the spherical atomic electron densities ρ A (r). Then, the weight function was defined based on ρ promol (r) and ρ procryst (r). The weight function is as follows: In all space, this continuous scalar function has 0 < W(r) < 1, and the 0.5 isosurface of this function is the HS. The 2D fingerprint was displayed to its HS by CrystalExplorer.

Electrostatic Potential Surface
The calculation function of the electrostatic potential surface (ESP) is as follows: The electrostatic interaction between a unit point charge placed at r and the system of interest can be measured by this function. ESP was calculated by Multiwfn [35] and visualized using the VMD program [36]. The cartesian atomic coordinates can be found in Supplementary Materials: Tables S1-S3.

Construction of Models
All computation simulations were performed using the Materials Studio software. We constructed the initial unit cell structure of CL-20, HATO, and FOX-7 based on the Cambridge Structural Database available at Cambridge Crystallographic Data Centre (CCDC) using Material Studio software. The crystal cell model of three crystals is shown in Figure 9.

Choice of Force Field
The force field plays a very important role to obtain reliable molecular dynamics simulation result. The popular COMPASS (condensed-phase optimized molecular potentials for atomistic molecular dynamics studies) force field has proved to be a universal all-atom force field, which has been calibrated from empirical parametrization techniques and state-of-the-art ab initio methods [8]. The COMPASS is suitable for molecular dynamic simulations involving conformation, vibration, structure, energy, and physical properties simulation of nitro (−NO2) containing compounds, such as RDX, HMX, and CL-20 [37,38]. So, CL-20 and FOX-7 were calculated by the COMPASS force field. However, there was no reasonable parameter for the N-oxides on the azole ring of HATO based on COMPASS force field in the simulation process. This result can be calibrated when replace COMPASS force field with Dreiding force field. Dreiding force field is a purely diagonal forcefield with harmonic valence terms and a cosine Fourier expansion torsion term and is more suitable for the N-oxides on the azole ring [39]. So, HATO cell structure was calculated by the Dreiding force field.

Structure Optimization and Morphology in Vacuum
The atomic positions, lattice parameters and the geometry structure of CL-20, HATO, and FOX-7 were firstly optimized to minimize the total energy and obtain an optimized molecule geometry by the selected force field. Then, the crystal morphology of CL-20, HATO, and FOX-7 was simulated under vacuum with an attachment energy (AE) model, in which the growth morphology method was employed to identify the morphologically important facets of these crystals.

Construction of Interface Models
The morphologically important growth crystal facets of CL-20, HATO, and FOX-7 crystals were cut to obtain the 3D periodic structure. Next, the periodic supercell was constructed. In the Amorphous Cell module, a solvent box was built with 200 molecules randomly distributed solvent molecules and it was consistent with the 3D periodic structure sizes of the important facets. Due to a favourable solubleness of CL-20 in ethyl acetate, N,N-Dimethylformamide(DMF), acetonitrile and acetone, these solvents were choice during the computation [40]. In the same way, the solvent molecules of HATO were dimethyl sulfoxide (DMSO) and N-methylpyrrolidone (NMP), and the solvent molecules of FOX-7 were DMF, DMSO, and NMP. Then, optimizing the solvent layer by the selected force field. The two-layer structure model containing solvent layer and main facet was constructed, which has a 50 Å thickness of the vacuum layer to eliminate the effect of additional free boundaries with a periodic boundary.

Choice of Force Field
The force field plays a very important role to obtain reliable molecular dynamics simulation result. The popular COMPASS (condensed-phase optimized molecular potentials for atomistic molecular dynamics studies) force field has proved to be a universal all-atom force field, which has been calibrated from empirical parametrization techniques and state-of-the-art ab initio methods [8]. The COMPASS is suitable for molecular dynamic simulations involving conformation, vibration, structure, energy, and physical properties simulation of nitro (−NO 2 ) containing compounds, such as RDX, HMX, and CL-20 [37,38]. So, CL-20 and FOX-7 were calculated by the COMPASS force field. However, there was no reasonable parameter for the N-oxides on the azole ring of HATO based on COMPASS force field in the simulation process. This result can be calibrated when replace COMPASS force field with Dreiding force field. Dreiding force field is a purely diagonal forcefield with harmonic valence terms and a cosine Fourier expansion torsion term and is more suitable for the N-oxides on the azole ring [39]. So, HATO cell structure was calculated by the Dreiding force field.

Structure Optimization and Morphology in Vacuum
The atomic positions, lattice parameters and the geometry structure of CL-20, HATO, and FOX-7 were firstly optimized to minimize the total energy and obtain an optimized molecule geometry by the selected force field. Then, the crystal morphology of CL-20, HATO, and FOX-7 was simulated under vacuum with an attachment energy (AE) model, in which the growth morphology method was employed to identify the morphologically important facets of these crystals.

Construction of Interface Models
The morphologically important growth crystal facets of CL-20, HATO, and FOX-7 crystals were cut to obtain the 3D periodic structure. Next, the periodic supercell was constructed. In the Amorphous Cell module, a solvent box was built with 200 molecules randomly distributed solvent molecules and it was consistent with the 3D periodic structure sizes of the important facets. Due to a favourable solubleness of CL-20 in ethyl acetate, N,N-Dimethylformamide(DMF), acetonitrile and acetone, these solvents were choice during the computation [40]. In the same way, the solvent molecules of HATO were dimethyl sulfoxide (DMSO) and N-methylpyrrolidone (NMP), and the solvent molecules of FOX-7 were DMF, DMSO, and NMP. Then, optimizing the solvent layer by the selected force field. The two-layer structure model containing solvent layer and main facet was constructed, which has a 50 Å thickness of the vacuum layer to eliminate the effect of additional free boundaries with a periodic boundary.