Crystal morphology prediction of CL-20 and 1,4-DNI co-crystal at different temperatures

The morphologies of hexanitrohexaazaisowurtzitane (CL-20) and 1,4-dinitroimidazole (1,4-DNI) co-crystal under vacuum or solvent at different temperatures were predicted. The CL-20/1,4-DNI co-crystal has six important growth crystal planes: (002), (011), (101), (11‒1), (110), (111). The areas of (002), (101), and (011) planes account for a relatively large proportion, which are important crystal planes that affect the crystal morphology. The crystal habits at different temperatures were simulated. The simulation results showed that the crystal plane attachment energy of CL-20 and 1,4-DNI co-crystal increases with the increase of temperature, indicating that the increase of temperature is conducive to the growth of crystal planes. The aspect ratio decreases with the increase of temperature and the morphology of co-crystal becomes more spherical at a higher temperature. The theoretical predictions are in good agreement with the experiment. The simulation results can provide guidance for the crystallization of CL-20/1,4-DNI to obtain a nearly spherical crystal morphology. The CL-20/1,4-DNI unit cell structure was geometrically optimized by the COMPASS force field. The AE model was used to predict the morphology of CL-20/1,4-DNI under vacuum, resulting in the most morphologically important growth planes. Ethyl acetate was selected as the solvent. The interaction energy between the solvent and the crystal plane, and the attachment energies in solvent at 298 K, 320 K, 340 K, 360 K, and 380 K were predicted. The NVT ensemble is used in the molecular dynamics calculation process. The simulation step is 1 fs and the total simulation time is 500 ps. The Andersen thermostat is selected as the temperature control method. In the potential energy calculation, the atom-based and Ewald methods were selected to calculate the van der Waals force and the electrostatic interaction force, respectively.


Introduction
Energetic materials are a class of substances that can rapidly release a large amount of energy under external stimuli, but the high energy of energetic materials is often accompanied by high sensitivity, which will seriously hinder the development and application of energetic materials [1]. Hexanitrohexaazaisowurtzitane (CL-20) is a widely used high-energy-density energetic material [2]. However, it has a high sensitivity. In order to overcome this defect, CL-20 is co-crystalized with some low sensitive energetic compounds such as 1,4-Dinitroimidazole(1,4-DNI) [3,4]. High-energy co-crystal techniques typically combine two or more explosives into a single lattice through non-covalent interactions such as hydrogen bonding and van der Waals forces [5]. Therefore, the obtained high-energy co-crystals have different structures and properties compared to pure components and physical mixtures [6]. Co-crystallization of CL-20 and 1,4-DNI can yield high-energy energetic materials while maintaining low sensitivity and increasing safety.
The morphology of the crystal is the combined result of internal structural factors and external conditions, and the crystal presents different morphologies because of the different growth rates of each crystal plane. The growth rate of the crystal plane is related to the type of solvent, temperature and other factors [7]. Crystal morphology will affect the energy and safety of energetic materials [8]. Crystals of the same size, with spherical morphology, have lower impact strength and friction sensitivity than needle-like or plate-like morphologies. The security is better [9].
On the basis of these studies, this paper uses the attachment energy (AE) model to predict the morphology of CL-20/1,4-DNI co-crystal in vacuum. By obtaining the attachment energies and main growth crystal planes, the modified attachment energy (MAE) model and the molecular dynamics method were used to predict the crystal morphology of CL-20/1,4-DNI co-crystal at different temperatures. The aspect ratio of the crystal was obtained. The effect of different temperatures on the crystal plane was analyzed and compared with the experimental crystal morphology. The crystal morphology predicted by theoretical calculation provides guidance for obtaining the ideal crystal morphology through experiments.

Theory
The attachment energy (AE) model is a model proposed by Hartman and Bennema based on the Periodic Bond Chain theory [19,20]. The attachment energy (E att ) is defined as the energy released when a crystal face with a thickness of dhkl attaches to the crystal surface, and can be expressed as formula (1) [21].
where E latt is the lattice energy of the crystal; E slice is the energy released by growing a wafer with a thickness of d hkl .
The AE model considers that the relative growth rate (R hkl ) of the crystal face in the crystal is proportional to the absolute value of the corresponding attachment energy (E att ) [22], see Eq. (2). The low crystal face attachment energy corresponds to a slow growth rate, which is important in morphology.
Based on the given crystal structure, the growth morphology under vacuum conditions can be predicted by calculating the attachment energies of different growth crystal planes. In the solvent condition, the solvent will be attached to the growing crystal plane, and the solute must overcome the resistance of the solvent layer when growing on the crystal plane, resulting in the inhibition of the crystal plane growth rate. To this end, the attachment energy needs to be corrected [23]: where E s is the energy correction term introduced by the solvent adsorption, which represents the inhibitory effect of the solvent on the growth of the crystal plane. S is the correction factor that reflects the structural characteristics of the growing crystal plane, which can be understood as the weight of different crystal plane structures on E s and is defined as: where A acc is the solvent accessible surface area of the crystal plane unit (hkl). A hkl is the cross-sectional area of the crystal plane unit. Obviously, the larger the S value, the rougher the crystal plane. The smaller the S value, the flatter the crystal plane. E s can be calculated from the interaction energy (E int ) of the solvent with the crystal plane: where A box is the cross-sectional area of the growing crystal plane (hkl) in the simulation box. E int is the difference between the total energy (E tot ) of the crystal plane-solvent system, the crystal plane energy (E cry ) and the solvent energy (E sol ), namely: In the revised AE model, the crystal plane with a lower absolute value of the corrected attachment energy ( E ′ att ) corresponds to a slower relative growth rate ( R ′ hkl ), and the two are proportional, see Eq. (7)
The COMPASS force field is the most widely used in the field of energetic materials and is suitable for the simulation of energetic compound molecules [25]. The CL-20/1,4-DNI unit cell structure was geometrically optimized, and the experimental lattice parameters were compared with the optimized values ( Table 1). The relative error between the optimized lattice parameters and the experimental parameters is within 5%, indicating that the COMPASS force field is suitable for the simulation of CL-20/1,4-DNI.
Then the AE model was used to predict the morphology of CL-20/1,4-DNI under vacuum, resulting in the most morphologically important growth planes. Then, the important growth crystal planes are cut, and the supercell is used to construct a periodic structure with a length and width of about 40 Å. Ethyl acetate was selected as the solvent system, and a solvent box composed of randomly distributed solvent molecules was constructed.
The crystal plane and solvent bilayer structure model of CL-20/1,4-DNI was constructed by the Build Layer module. A vacuum thickness above the solvent layer is 50 Å to eliminate the influence of other free boundaries, and the vacuum thickness of the solvent layer and the crystal layer is 3 Å. The NVT ensemble is used in the molecular dynamics calculation process. The simulation step is 1 fs and the total simulation time is 500 ps. The data were collected every 5000 steps. The Andersen thermostat [26] is selected as the temperature control method. The temperature is set to 298 K, 320 K, 340 K, 360 K and 380 K. In the potential energy calculation, the atom based and Ewald methods were selected to calculate the van der  Waals force [27] and the electrostatic interaction force [28], respectively. Molecular dynamics simulation was carried out on the geometrically optimized bilayer model, so that the solvent molecules were distributed uniformly and the configuration reached an equilibrium. Finally, the stable configurations of E sol , E cry and E tot were calculated respectively, and obtained E ′ att was used to predict the possible crystal habit of CL-20/1,4-DNI under the influence of solvent. The simulation procedure was shown in Fig. 2.

Crystal morphology in vacuum
The morphology of CL-20/1,4-DNI in vacuum was simulated by the AE model, and the main crystal plane parameters were shown in Table 2. It can be seen in Fig. 3 that the morphology of CL-20/1,4-DNI in vacuum is a block-like structure with edges and corners, with 6 main growth crystal planes: (002), (011), (101), (11-1), (110), (111). The morphology of CL-20/1,4-DNI was predicted under vacuum, and the obtained aspect ratio is 2.359, in which the crystal plane area of (002) accounts for the largest 39.75%. The surface with fast growth rate would disappear, while a surface with a slow growth rate appears on the final morphology. The attachment energy affects the growth rate of the crystal plane. The larger the attachment energy, the faster the growth rate of the crystal plane, and the smaller the corresponding crystal plane. On the contrary, the smaller the crystal face attachment energy, the larger the crystal face, the (002) face has the smallest attachment energy, and the crystal face grows at the slowest speed, so the crystal face is the largest, which is consistent with the results of literature [29,30]. However, (111), (110) and (11-1) crystal planes have larger attachment energy and the smallest crystal plane area, which have little effect on the growth of CL-20/1,4-DNI co-crystal.
To understand the interaction of solvent molecules with different CL-20/1,4-DNI crystal planes, the Connolly surface of CL-20/1,4-DNI co-crystal is shown in Fig. 4. The blue mesh on top represents the accessible solvent surface. The surfaces of (002), (011), and (101) in Fig. 4 have relatively large concavities and convexities. They are relatively rough and have many voids, which are conducive to the adsorption of solvent molecules. The (111) crystal plane is relatively flat.
The calculated solvent accessible area (A acc ), surface area (A hkl ) and corresponding S (A acc /A hkl ) values of different CL-20/1,4-DNI surfaces were listed in Table 3. The larger the S value, the rougher the crystal plane [31], which makes the molecules in the solution more easily adsorbed. The S values of the (002), (011) and (101) crystal planes are larger ( Table 3), so that the growth of these crystal planes is inhibited, and the area of the crystal planes are larger. The smaller S values of (11-1) and (111) indicate that the surface is relatively smooth, the growth Table 2 The main growth planes and related parameters of CL-20/1,4-DNI in vacuum a a d hkl is interplanar distance in Å, A is surface area in Å 2 , n is multiplicity, E att is attachment energy in kcal·mol -1 , A/% is the percentage of a corresponding surface area to the total area  of the crystal plane is less inhibited by the solvent, and the area of the crystal plane is smaller, which is consistent with Fig. 4.

Effect of Temperature on Morphology of CL-20/1,4-DNI co-crystal
The CL-20/1,4-DNI co-crystal data at different temperatures calculated by molecular dynamics simulation using the modified attachment energy model were listed in Table 4. The interaction energies of E int between the solvent and the crystal face are all negative, indicating that the adsorption of the solvent is an exothermic process. Therefore, the larger the absolute value of the interaction energy, the stronger the adsorption effect of the solvent and the slower the growth rate of the crystal plane. Correcting the absolute value of the attachment energy determines the area of the corresponding crystal plane.
Through the value of E′ att , the growth rate of the crystal plane can be obtained, then the crystal morphology can be predicted. The predicted crystal morphology at different temperatures is shown in Fig. 5. The crystal morphology at 320 K is prismatic that is similar to the experiment [18].
By comparing the E′ att values of different crystal planes as shown in Fig. 6, it can be found that the E′ att value of the (002) crystal plane does not change much with temperature, indicating that the temperature has little effect on the (002) crystal plane. The E′ att value of other crystal planes changes greatly with the change of temperatures, and these crystal planes are greatly affected by temperature. With the increase of temperatures, the absolute value of E′ att increases, indicating that the increase of temperature will promote the growth of crystal planes. These three crystal planes have the largest area, so they are the most important crystal planes that affect the crystal morphology of CL-20/1,4-DNI.
By comparing the crystal morphology of the CL-20/1,4-DNI co-crystal at different temperatures and the calculated aspect ratio (Fig. 7), it is found that the aspect ratio decreases with the increase of temperature, and the crystal morphology is closer to a spherical shape at higher temperature. The minimum ratio is 1.88.

Non-covalent interactions analysis
Radial distribution function (RDF) is the interactions between molecules based on the measured probability g(r) of the particle appearing at a distance r from a reference atom. Intermolecular interactions mainly include hydrogen bonding and van der Waals interactions. Based on the Fig. 4 Connolly surface of each main growth crystal plane of CL-20/1,4-DNI Table 3 The parameter values for the crystal habit surfaces of CL-20/1,4-DNI a a A acc is is the solvent accessible surface area, A hkl is the cross-sectional area of the crystal plane unit, S is A acc /A hkl   Fig. 8. The first peak in Fig. 8(a) is around 2.7 Å, indicating that there is van der Waals force between the ethyl acetate molecules and the co-crystal plane. The first peaks in Fig. 8(b) and Fig. 8(c) are around 3 Å and 3.3 Å, indicating the existence of van der Waals force between solvent molecules and co-crystal planes. The first peak at 298 K in Fig. 8(d) is around 3.5 Å. There is a van der Waals force between the solvent molecule and the co-crystal plane, and the van der Waals force disappears as the temperature increases. The curve in Fig. 8 shows that the peak at 298 K is the largest and the intermolecular van der Waals forces are the strongest in comparison with higher temperatures.

Conclusions
The crystal morphology and main surface of CL-20/1,4-DNI in vacuum and solvent were simulated using the AE model. The effect of temperature on each crystal plane of CL-20/1,4-DNI was studied by simulating the adsorption of ethyl acetate on the habitual plane of CL-20/1,4-DNI. The crystal morphologies of CL-20/1,4-DNI at different temperatures were predicted. The main results can be summarized as follows: 1. The simulated morphology was a prismatic structure with six important growth crystal planes: (002), (011), (101), (11-1), (110) and (111). Among them, the areas of (002), (101) and (011) planes account for a relatively large proportion, which are the main crystal planes that affect the crystal morphology. 2. The radial distribution functions show that the interaction of solvent molecules with the CL-20/1,4-DNI co-crystal surface includes van der Waals forces, and the interaction strength of solvent molecules with the co-crystal surface varies with temperatures. The absolute value of the attachment energy of each crystal plane increases with the increase of temperature, indicating that the increase in temperature will promote the growth of crystal planes. Temperature has the greatest influence on (002), (011) and (101) crystal planes. The crystal planes with the largest areas control CL-20/1,4-DNI morphology. 3. The crystal morphology is in good agreement with the experiment at room temperature. The aspect ratio becomes smaller with the increase of temperature, which facilitates the crystallization of a spherical morphology at higher temperature.  (2). c N(1)‧‧‧O (2). d H(1)‧‧‧O (2)