The Crystal Structure and Morphology of 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20) p-Xylene Solvate: A Joint Experimental and Simulation Study

The crystal structure of 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaiso-wurtzitane (CL-20) p-xylene solvate, and the solvent effects on the crystal faces of CL-20 were studied through a combined experimental and theoretical method. The properties were analyzed by thermogravimetry-differential scanning calorimetry (TG-DSC), Fourier transform infrared spectroscopy (FTIR) and X-ray diffraction (XRD).The growth morphology of CL-20p-xylene solvate crystal was predicted with a modified attachment energy model. The crystal structure of CL-20p-xylene solvate belonged to the Pbca space group with the unit cell parameters, a = 8.0704(12) Å, b=13.4095(20) Å, c = 33.0817(49) Å, and Z = 4, which indicated that the p-xylene solvent molecules could enter the crystal lattice of CL-20 and thus the CL-20 p-xylene solvate is formed. According to the solvent-effected attachment energy calculations, (002) and (11−1) faces should not be visible at all, while the percentage area of the (011) face could be increased from 7.81% in vacuum to 12.51% in p-xylene solution. The predicted results from the modified attachment energy model agreed very well with the observed morphology of crystals grown from p-xylene solution.

the growth model of explosive CL-20 crystals in p-xylene solution at a given temperature. These results provided the fundamental information to predict the potential morphology changes during the formation of solvated CL-20 or other explosives.

Computational Theory
The growth morphology algorithm is based on the attachment energy (AE) method, which is proposed by Hartman and Bennema based upon period bond chain (PBC) theory [26]. The attachment energy, Eatt, is defined as the released energy of a growth slice to a growing crystal surface. Eatt is computed as [27]: (1) where Elatt is the lattice energy of the crystal, Eslice is the energy of a growth slice of thickness dhkl, the relative growth rate in vacuum of the crystal surface, Rhkl, is assumed to the proportional to the attachment energy: (2) The modified morphology is generated based on the relation stated by Hartman and Bennema, where R'hkl represents the growth rate in a particular direction, which is directly proportional to the modified attachment energy [28]: The solvent can be assumed to reduce the growth rate. First, the solvent needs to be removed from the surface before the crystal face can grow, which costs energy and herewith the apparent attachment energy decrease. So modified AE can be introduced with the following formula: (4) where correction factor S reflecting the surface features and the effect of solvent volume, defined as the accessible solvent surface per unit crystal area [29]: (5) where Aacc is the area of accessible solvent surface of an (h k l) slice in the unit cell, Ahkl is the surface area of an (h k l) slice in the unit cell.
Since the effect of solvent depends on both the different surface chemistry and the associated topography (step structures) of the crystal face, the correction term Es describing the energy of solvent binding on the crystal surface (h k l) can be calculated using the following formula [30]: (6) where Eint is the interaction energy between the solvent layer and the surface, Aacc is the area of accessible solvent surface of an (h k l) slice in the unit cell, Amodelis the total surface area of modeling box in the (h k l) direction. N is the number of CL-20 in the crystal layer. Eintis the interaction energy between the solvent layer and the surface, which is evaluated using the relationship: (7) where Etotal is the total energy of the surface and the solvent layer, Esurf is the energy of the surface without the solvent layer, and Esolv is the energy of solvent layer without the surface.

Computational Details
All of the MD and morphology predictions were done with the commercial molecular modeling software package Materials Studio 5.5 [31] with the Compass force-field. The crystal structure of γ-CL-20 (CCDC number 117,778) [32] was optimized with the forcite algorithm. Subsequently, the crystal morphology of CL-20 in vacuum was predicted with AE model to acquire the main stable crystal faces withdifferent miller indices (h k l).Then, the CL-20 crystal was cleaved along the main stable faces (10−1), (101), (011), (002), (110) and (11−1) with a depth of 2dhkl. The corresponding crystal layer was built as a periodic superstructure of about 30 × 30 Å 2 unit cell. A solvent layer filled with 100 random distributed p-xylene molecules was constructed with the amorphous cell tool with a layer density of 0.8611g/cm 3 and whose size was set appropriate to the crystal layer. The two-layer interfacial model was employed to investigate solvent effect on the crystal morphology ( Figure 1).  In this model, the two layers were a crystal slice and asolvent layer, i.e., the solvent layer docked CL-20 crystal surface along c axis and moved freely in the simulation. And a vacuum slab thickness of 30 Å was built above the solvent layer to eliminate the effect of additional free boundaries on the structure. Energy minimization for the interfacial model was carried out before the dynamics simulation. MD simulation was performed for 500 ps with a time step of 1 fs (T = 298 K and P = 1 atm) with the NVT ensemble [33].Coupling to the heating bath was carried out with the Andersen thermostat [34]. After the system was equilibrated at the target temperature, the interfacial configurations were sampled for 100 ps with a sampling interval of 1 ps. The Compass force-field, a powerful ab initio force-field parameterized and validated using condensed phase properties, enabled the accurate and simultaneous prediction of structural, conformational, vibrational, and other physical properties for a broad range of covalent molecules including most common organics [35]. For potential-energy calculations, the Coulomb interactions and Van der Waals forces were calculated by employing the standard Ewald method with a calculation accuracy of 0.001kcal/mol and cutoff distance was set as 12.5 Å [36,37].
In the CL-20 p-xylene solvate crystal packing, the γ-conformation of the CL-20 molecules was observed in the crystal structure of the solvate, in which there was one p-xylene molecule for every CL-20 molecule. Alternation of two CL-20 layers and one p-xylene layer could be found, as shown in Figure 2. In the CL-20 layers, individual molecules with both the molecular conformation and the relative spatial distribution observed in the γ-polymorph of CL-20 crystals (as shown in Figure 2). The C-H…O intermolecular interactions resulted in a network of weakly linked CL-20 molecules along two dimensions. The "face-to-face" arrangement (all six-membered rings of the cage) was oriented to c-axis, which leaded to the significative CL-20/CL-20 intermolecular interactions in the structure.
In the crystal structure, the most important intermolecular interactions between the hydrogen atoms (H) of the p-xylene molecules and the NO2 groups of CL-20 molecules were the primary driving force for the formation of CL-20 p-xylene solvate, which also made contribution to the stability of the crystal structure ( Figure 3). The possible hydrogen bond distances and angles are listed in Table 2 with the labeled atoms in Figure 3. In addition to the intermolecular hydrogen bond interactions, the O…C, and O…N stacking interactions between CL-20 and p-xylene played an important role in the stabilization of the solvated structure, which brought the two molecules together to form the CL-20 p-xylene solvate (Table 1). For example, the very strong C10-H10A…O5 hydrogen bond between the hydrogen atom of p-xylene and the nitro group of CL-20 was one of the main forces to form CL-20 p-xylene solvate, which had a bond distance of 2.6995 Å and bond angle of 166.6°. Table 1. Crystallographic data and structure refinement parametersfor CL-20 p-xylene solvate.
Limiting indices

Solvate Properties
The TG-DSC curve of CL-20-p-xylene solvate (dried at room temperature) is presented in Figure 4. There is a gradual weight loss around 105 °C due to the desolvation and a steep weight loss from 225 to 240.8 °C owing to decomposition. Two peaks were found on the DSC curve. The first peak at 105 °C was attributed to desolvation, in conformity with the slight weight loss. The second sharp peak at 240.8 °C could be assigned to the CL-20 decomposition [38]. FTIR spectroscopy is normally used for the phase determination of CL-20 and works very well in determination of the predominant phase and the changes caused by intermolecular hydrogen bonds are relatively apparent in the infrared spectrum [39]. The FTIR spectrum showed that several regions were highly sensitive to the structure changes of CL-20 p-xylene solvate ( Figure 5). The 813 cm −1 band could be assigned to C-H flexural out-plane vibration of p-xylene. Bands at 1600 cm −1 , 1580 cm −1 and 1500 cm −1 due to C=C stretching of benzene ring appear in the CL-20 p-xylene solvate structure. Among the three nitramine group stretching bands of CL-20, including as (NO2) (antisymmetric stretching vibration), s (NO2) (symmetric stretching vibration) and (N-N), only as (NO2) could be considered as a characteristic one to study [40]. In CL-20 p-xylene solvate, the 1628-1553 cm −1 bands could be assigned to asymmetrical NO2 stretching of CL-20, which was shifted to 1632-1564 cm −1 due to the interactions between CL-20 and p-xylene. The symmetrical NO2 stretching modes were mixed with other vibrations in the same region, especially with C-H bonds. Based on TG-DSC and FT-IR spectra, it could be confirmed that theCL-20p-xylene solvate was formed.

Crystal Morphology Prediction
The external morphology of crystals was predicted with the AE model. The predicted crystal morphology in vacuum was much different from the observed morphology in solvent as shown in   The crystal packing diagrams (Figure 7) revealed the surface chemistry and topography of CL-20 crystal morphology faces. The (101) and (011) faces were rather open and rough on the molecular level with nitro groups, hydrogen and oxygen atoms of the CL-20 exposed on the surface. Oxygen atoms were only observed exposed on the (10−1) face. The surface chemistry of (110) and (11−1) faces were mainly dominated by the exposed oxygen atoms and nitro groups. The (002) face was found to be smooth on the molecular level. Surface areas of six crystal faces were listed in Table 3. In the Connolly surface calculation, the grid interval was 0.4000 Å at a probe radius of 1.0 Å. The accessible solvent surfaces area of unit crystal denoted as S in Table 3. The Aacc and S were calculated quantitatively to compare the step structures of different crystal faces. As listed in Table 3, it was found that the (101) face had the highest S value, and the (002) face had the least S value on the contrary among all the exhibiting surfaces.   The solvent effect Es and modified attachment energy mod Eatt were calculated with formulas (1-6). The solvent effects depended on the different surface chemistry and the step structure of the crystal faces, which were quantified as the energy correction term Es in this study. With a stronger interaction, it was found that the absolute value of the interaction energy was larger. For example, the low growth rate of (011) face in the p-xylene solution could be attributed to the polar nitro groups, hydrogen atoms and oxygen atoms exposed on the surface, and the large accessible solvent surface, which would be conducive to the easy solvent binding. Therefore, the desolvation process became very difficult with the high energy barrier, and the solute molecules were hampered to reach the surface. In the (002) face, its absolute value of the Es was the lowest, resulting in the fastest growth in this direction. Inspite of the nearly equal interaction energies of (101) and (110) faces with the p-xylene molecules, the difference in the accessible solvent surface of unit area (see the S value in Table 2) had brought faster growth of (110) face than (101) face.
It is well known that the morphological importance of crystal faces is inversely proportional to their growth rates, which is affected significantly by the solvent. In a particular direction, the more negative the attachment energy is, the faster the growth rate will be and the less morphological importance the crystal face will have. As shown in Table 4, the attachment energies of all the habit faces were decreased with the solvent effect taken into consideration. This decrease indicated that interactions of the solvent molecules with the crystal layers were more favored compared to the CL-20 molecule itself. Therefore, the solvent effect restrained the growth rate. The comparison between Figure 6a,b revealed that the faster relative growth rates of (002) and (11−1) faces led to their disappearance in the final crystal habit. Conversely, the area percentage of (011) was increased up to 12.51% in p-xylene solution from 7.81% in vacuum due to strong solvent interactions. The changes of relative growth rates leaded to crystal shape which was close to a hexagon with the aspect ratio of 1.95, which was consistent with the observed morphology in p-xylene solvent. These results of the morphological predictions ( Figure 6b) were in good agreement with the experimental morphology (Figure 6c), which indicated the successful prediction of the modified AE model. Table 4. Calculated attachment energies for dominant crystal morphology faces together with corrected attachment energies and relative growth rates for main faces.

Radial Distribution Functions of the Interfacial Model
The radial distribution function (RDF) gives a measure of the probability that, given the presence of an atom at the origin of an arbitrary reference frame, there would be an atom with its center in a spherical shell with infinite thickness and a distance, r, from the reference atom. In the radial function graph g(r)-r, the peaks were all within 3.0 Å, which indicated the hydrogen bonds between O atoms and H atoms [41]. Some high peaks existed in the region of 3.0-8.0 Å, implying the strong Van der Waals' forces and the electrostatic interactions between this atom pair. [42] In order to analyze the solvent-crystal interactions clearly, the interface structures between CL-20 (10−1) faces and p-xylene molecules were explored with RDF(s), as shown in Figure 8. The RDF(s) graphs of two pairs of atoms were plotted, which were hydrogen atoms in p-xylene molecules to oxygen atoms and nitrogen atoms of the nitro groups in CL-20 molecules, respectively denoted as Hp-xylene-OCL-20 and Hp-xylene-NCL-20 pairs. Similar RDF(s) could be obtained when analyzing other interfacial model systems comprising p-xylene molecules and the other five crystal faces.
The sharp peaks of O-H with r = 2.5 Å and 3 Å indicated the existing hydrogen bonds between the O and H atoms. With the interval of 3-8 Å, wide peaks could be observed in g(r)O-H-r and g(r) N-H-r, which indicated a strong Van der Waals' forces and the electrostatic interactions among CL-20 and p-xylene.

Crystallization of CL-20 p-Xylene Solvate by Evaporation
A glass vial was loaded with γ-CL-20 (2.19 mg, 5.0 µmol) and p-xylene (20 mL). Mild heating was utilized to completely dissolve all solids. About 20 mL of the solution was prepared and then carefully filtered using 0.45 µm PTFE filter into a clean glass vial. The solution was allowed to evaporate at room temperature over several days, single crystals were obtained and collected and then used for further characterization.

Characterization
Crystals suitable for X-ray crystallographic analysis were obtained as mentioned above. The CL-20 p-xylene solvate crystal structure was determined by a Rigaku RAXIS IP diffractometer and SHELXTL crystallographic software package of molecular structure. The single crystals were mounted on a Rigaku RAXIS RAPID IP diffractometer equipped with a graphite-monochromatized MoKα radiation (λ = 0.71073 Å). Data were collected by the ω scan technique. The structure was solved by direct method with SHELXS-97 [43] and expanded by using the Fourier technique. The non-hydrogen atoms were refined anisotropically. The hydrogen atom was determined with theoretical calculations and refined with an isotropic vibration factor.
Simultaneous TG-DSC of crystal was carried out on a thermal analysis system (Netzsch STA-449F3, Selb, Germany). The thermal behavior of sample was studied over a range at a heating rate 10 °C/min from 50 °C to 500 °C in nitrogen atmosphere.
FTIR spectra were recorded on a Bruker tensor 27 instrument with KBr pellets. The resolution was 4 cm −1 , and the scan range was 400-4000 cm −1 . The sample was ground with KBr powder and then pressed to obtain a suitably sized pellet for spectrum recording.

Conclusions
The CL-20 p-xylene solvate crystal was obtained through slow evaporation. The p-xylene molecules could enter the CL-20 crystal lattice by forming weak intermolecular interactions and stacking interactions, which affected not only the crystal structure but also the crystal morphology. The solvent effects on different crystal faces varied due to the diversity of functional groups, and then to change the crystal morphology significantly. With the modified AE, (002) and (11−1) faces were not visible at all due to the exposure of the less polar groups or atoms and hence a small solvent interaction, while the area percentage of (011) face was increased up to 12.51% in p-xylene solution from 7.81% in vacuum due to the strong solvent interaction. The predicted morphology was in good agreement with the observed experimental morphology, and this consistency validated the applicability of the modified AE model. Furthermore, the RDF analysis indicated that the p-xylene solvent molecules were adsorbed on the CL-20 faces mainly via interactions between hydrogen bonds, Van der Waals forces and electrostatic interactions. The molecular simulation could be an effective tool for the solvent selection to modify the morphology of the solvated CL-20 crystal and other solvated explosive crystals. In addition, molecular simulation could reduce the development cycle and the cost of explosive crystallization process.