Face-Dependent Solvent Adsorption: A Comparative Study on the

Jul 11, 2017 - Figure 1. Schematic diagram on the construction of crystal-solvent model. ... density functional with DNP basis set in Dmol3 code of Ma...
0 downloads 0 Views 665KB Size
Subscriber access provided by AUSTRALIAN NATIONAL UNIV

Article

Face-Dependent Solvent Adsorption: A Comparative Study on the Interfaces of HMX Crystal with Three Solvents Yingzhe Liu, Weipeng Lai, Yiding Ma, Tao Yu, Ying Kang, and Zhongxue Ge J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b04470 • Publication Date (Web): 11 Jul 2017 Downloaded from http://pubs.acs.org on July 12, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Face-Dependent Solvent Adsorption: A Comparative Study on the Interfaces of HMX Crystal with Three Solvents Yingzhe Liu,* Weipeng Lai, Yiding Ma, Tao Yu, Ying Kang, Zhongxue Ge State Key Laboratory of Fluorine & Nitrogen Chemicals, Xi’an Modern Chemistry Research Institute, Xi’an 710065, P. R. China

* To whom correspondence should be addressed. E-mail: [email protected]. ACS Paragon Plus Environment

1

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 19

ABSTRACT: To understand the crystal-solvent interfacial interactions at the molecular scale, the interfaces between three solvents, i.e. acetone (AC), γ-butyrolactone (BL), and cyclohexanone (CH), and three growth faces of 1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX) crystal have been investigated with the aid of theoretical chemistry. The results show that the structural features of crystal faces play a critical role in the energetic, structural, and dynamic properties at the interfaces. For each solvent, the same change trend of some properties among the three faces of HMX crystal is observed, including adsorption affinity, local mass density, and solvent diffusion. For example, the rate of solvent diffusion at the three faces ranks as (011) > (110) > (020) regardless of solvent species. This can be attributed to the similar adsorption sites for solvent incorporation at the same face, which are concentrated at the cavities formed by surficial HMX molecules.

ACS Paragon Plus Environment

2

Page 3 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

INTRODUCTION Crystallization from solutions is a crucial avenue of separation and purification for the synthesis of energetic materials (EMs). Especially, the majority of EMs generally decomposes through melt crystallization. The crystal qualities of EMs, including polymorph, morphology, size, purity, defect, etc., can greatly influence their power and safety properties. Taking polymorph as an example, the most stable crystal form at ambient conditions can possess both the largest density and the lowest impact sensitivity for the same EM, such as HMX (1,3,5,7-tetranitro-1,3,5,7-tetrazocane).1 HMX is one of the most powerful EMs, which is widely utilized in civil as well as military fields, e.g. airbag inflators, explosives, rocket propellants. Three polymorphs (α, β, and δ) and one hemihydrate (γ) have been identified for the HMX crystal at ambient temperature and pressure.2-4 The stable crystal form, β-HMX, is the most valuable in practical usage. A tendency to form twins for the β-HMX crystal is observed experimentally, which depends on the solvent employed in the crystallization process. For instance, twins are practically formed in every crystal when β-HMX is grown from cyclohexanone, while less twining can be found in acetone.5 Recent molecular simulation showed that the twined HMX crystal can enhance the shock sensitivity and hotspot-forming ability due to large free volume near the twinning plane, leading to the faster decomposition in contrast to the perfect HMX crystal.6 Additionally, diverse crystal morphologies of HMX can be obtained by using different solvents in the course of crystallization.7 The sphere-shaped EM crystals can eliminate sharp edges that may act as detonation hot spots,8-10 and thus generally have lower mechanical and thermal sensitivities than the needle-like ones at the same level of size. In view of the time and money cost on crystallization experiments, the solvent-dependent morphologies of EM crystals could be theoretically predicted in advance. Previous computational simulations estimated the morphologies of HMX crystals grown from solutions, and probed several factors ranging from solvent, temperature, to supersaturation.7,11-13 The results revealed that the theoretical prediction of HMX crystal morphology is in reasonable agreement with the crystallization experiment. Despite extensive explorations of solvent-induced crystallization, the molecular-scale understanding ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 19

of solvent effect on the HMX crystal growth to an extent is still limited, such as the ordered arrangement of solvent molecules at the crystal face. The fundamental insights into the interactions and properties at the crystal-solvent interface are essential so that the crystal growth can be controlled to acquire the desired quality. In the present contribution, the interface between an HMX crystal and a solvent was theoretically investigated from a microscopic point of view. As a comparative study, three solvents, namely, acetone (AC), γ-butyrolactone (BL), and cyclohexanone (CH) were chosen because of their common use in the HMX crystallization. The growth faces of HMX crystal with morphological importance, including (011), (110), and (020) faces, were taken into account. The interactions as well as the energetic, structural, and dynamic properties at the interfaces were examined.

METHODS Models. The initial structure of β-HMX crystal was taken from the available experimental data (CCDC 79230), 14 which has two irreducible HMX molecules in one unit cell and belongs to a monoclinic crystal system with P21/n symmetry. The lattice parameters of β-HMX crystal are a = 6.526 Å, b=11.037 Å, c=7.364 Å, and β = 102.67º. After the optimization with the COMPASS force field,15 these lattice parameters became a = 6.468 Å, b= 10.336 Å, c= 7.629 Å, and β = 101.50º, indicative of an acceptable accuracy for molecular simulations. To further explore the crystal-solvent interactions, two models, i.e. adsorption model and interfacial model, were built and the corresponding construction processes are schemed in Figure 1. To determine the growth faces with morphological importance, the growth morphology of HMX crystal in vacuum was predicted on the basis of the attachment energy method,16,17 which is widely used in EMs crystals.18-26 The prediction result showed that the growth morphology of HMX crystal under vacuum condition is composed of five faces, i.e. (011), (110), (020), (101), and (101) faces. The (011), (110), and (020) faces are the most morphologically important because they constitute more than 97 percent of the entire crystal surface area. In this study, therefore, only these three faces were ACS Paragon Plus Environment

4

Page 5 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

investigated for the crystal-solvent interactions. The structure of growth face was constructed by cleaving the HMX crystal along the (h k l) index, and the normal to the growth face was aligned along the z direction. To avoid the free diffusion of surficial HMX molecules, the HMX face was constrained along a, b, and c directions in all the calculations and simulations. Crystal Surface

Adsorption Model

MC Anneal

AE Model

Cleave

Solvent Packing HMX Crystal

Growth Morphology

2D Crystal Facet

Crystal Layer

Interfacial Model

Figure 1. Schematic diagram on the construction of crystal-solvent model.

To verify the accuracy of COMPASS force field for three solvents, the solvent densities were estimated by molecular dynamics (MD) simulations. The simulation box was consisted of 400 randomly distributed molecules for each solvent. After geometry optimization, the MD simulation was carried out for 500 ps in the NPT ensemble at 298 K and 1 bar. The average densities of AC, BL, and CH during the trajectories of last 100 ps were calculated as 0.793, 1.082, and 0.931 g/cm3, which are in good accordance with the corresponding experimental ones of 0.785, 1.128, and 0.948 g/cm3 (ref.27). In the adsorption model, the growth face has three unit cell of HMX crystal in depth. The areas of (011), (110), and (020) faces are 660.03, 554.95, and 580.20 Å2, respectively, and are large enough to accommodate a solvent molecule. A thickness of 20 Å vacuum was added above the growth face to eliminate possible boundary effect. The conformations of a solvent molecule adsorbed at the crystal face were sampled by means of Monte Carlo (MC) method combined with the annealing technique. To increase the sampling of HMX-solvent configuration and reach the global minimum of potential surface, ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 19

the temperature range was set from 102 to 105 K in the annealing process and total ten cycles with 105 steps per cycle were performed. The conformation sampled with the lowest energy was regarded as the stable adsorption structure. To build the interfacial model, the growth face has six unit cell of HMX crystal in depth and is equally separated by a vacuum layer with a thickness of 40 Å. The vacuum layer was then packed with solvent molecules according to the predicted densities. Finally, the structural relaxation of solvent molecules was conducted with a 10-cycle simulated annealing for total 500 ps, where the temperature range is from 300 to 500 K. The detail of each interfacial model is summarized in Table 1.

Table 1. Detail of interfacial models investigated in this work* Interface

a/Å

b/Å

c/Å

α/(deg)

β(deg)

γ(deg)

NHMX

NAC

NBL

NCH

(011)

38.81

38.54

75.96

90.00

90.00

83.20

216

442

407

307

(110)

38.14

36.58

73.08

90.00

90.00

96.07

180

419

386

291

(020) 38.14 38.81 71.47 90.00 90.00 101.50 180 438 403 NHMX, NAC, NBL, NCH represents the number of HMX, AC, BL, and CH molecules, respectively

304

*

Electrostatic Potentials Calculation. To understand the interactions of HMX crystal face with solvent molecules, the electrostatic potentials on the 0.001 au (electrons/bohr3) molecular surfaces were computed by the PBE28,29 of generalized gradient approximation (GGA) density functional with DNP basis set in Dmol3 code of Materials Studio 8.0 software.30 A global orbital cutoff of 3.7 Å for plane-wave basis set and a convergence tolerance of 10-6 Hartree for electron self-consistency were set. The k-points were sampled on a grid of spacing 0.07 Å-1 in Brillouin zone for each crystal face. Adsorption Affinity Estimation. To quantify the adsorption affinity of the solvent molecule at the crystal face, the binding energy of adsorption motif was calculated by Gaussian09 program.31 The adsorption motif was composed of a solvent molecule and several HMX molecules within 5 Å of it, which was extracted from the stable adsorption model. The geometry of each adsorption motif was optimized and characterized to the relative energy minimum of the potential surface at the B3LYP-D3/6-31G* level of theory.32-34 ACS Paragon Plus Environment

6

Page 7 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The binding energy (Ebind) of adsorption motif was determined by the following formula: Ebind = − (Emotif – (Esol + EHMX))

(1)

where Emotif is the electronic energy of the adsorption motif, Esol and EHMX are the electronic energies of solvent and HMX molecules taken from the adsorption motif without optimization, respectively. The basis set superposition error was also considered using the standard counterpoise correction.35 The bigger the binding energy is, the stronger the adsorption affinity is. Interfacial System Simulation. The MD simulations on the interfacial models were performed with the Forcite module of Materials Studio 8.0 software in the NVT ensemble. The periodic boundary conditions were applied in three directions of Cartesian space. The time scale of each MD simulation is 10 ns, which is adequate for system equilibrium reflected by steady fluctuation of energy and temperature profiles. The temperature was kept at 300 K with the help of Andersen thermostat. The electrostatic forces were evaluated by Ewald summation approach and the van der Waals interactions were truncated by a smoothed of 12.5 Å spherical cutoff. The equation of motion was integrated by the time step of 1 fs. The trajectory data were collected every 2500 time steps. Partial analysis and visualization of simulation trajectories were conducted by means of VMD package.36

RESULTS AND DISCUSSION Structural features. As shown in Figure 2, distinct electrostatic potentials of three HMX crystal faces are observed as a result of different molecular packing. For the (110) face, the majority of negative potentials are located at peaks while the positive potentials are mostly concentrated at basins. It is evident from color depth that the range of potentials at the (110) face is larger than the other two faces, indicative of stronger electrostatic interactions with solvent molecules.

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 19

Figure 2. Electrostatic potential (unit in kcal/mol) on the 0.001 au molecular surface

The pattern of potential distribution is remarkably different between HMX and solvent molecules. The strongly positive potentials are focused above the central part of HMX molecular surface and the negative potentials are located around the periphery of HMX molecule, which is the representative characteristic of CHNO EMs37. By contrast, there are negative potentials around the electron-rich oxygen atom in solvent molecules as expected, exhibiting stronger than the positive ones, especially for CH. From the viewpoint of electrostatic attractions, the solvent molecules would adopt the optimal orientations to match the potential distribution of HMX crystal faces, e.g. the carbonyl groups of the solvent molecules contact with the positive-potential basin area of crystal face. Previous studies reported that the structural feature of crystal face, such as flatness and roughness, plays a significant role in the solvent adsorption since the adsorption sites mainly originate from the cavities formed at the basin area of crystal face.38,39 Accordingly, the crystal face with relative flatness is unfavorable for solvent adsorption. To quantitatively describe the structural feature of HMX crystal face, the parameter S is defined as S = Aacc/Ahkl

(2)

where Aacc is the solvent accessible surface area of crystal face and Ahkl is the corresponding cross-sectional area. The larger S value means the more roughness of crystal face. The S values of the (011), (110), and (020) faces are 1.21, 1.30, and 1.31, respectively. It may be inferred that the solvent adsorption is more facilitative at the (110) and (020) faces than the (011) face. ACS Paragon Plus Environment

8

Page 9 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Adsorption

affinity.

The

optimized

structures

of

different

adsorption

motifs

at

the

B3LYP-D3/6-31G* level of theory are displayed in Figure 3. Expectedly, solvent molecules were adsorbed at the cavities of HMX crystal face, which can maximize their contact with HMX molecules. The capacity of surface cavity is limited so that it cannot accommodate the entire solvent molecule, particularly for BL and CH. The carbonyl groups of solvent molecules are primarily buried at the surface cavities while the alkyl groups are exposed outside. The solvent molecules adopt similar arrangements adsorbed at the (110) and (020) faces. Taking the (020) face as an example, the solvent molecules tilt towards the same direction with a nearly equal angle. Conversely, the adsorption arrangements of solvent molecules are different at the (011) face due to the relatively smooth surface structure.

Figure 3 The geometries of adsorption motifs optimized at the B3LYP-D3/6-31G* level of theory: (a) AC-(011), (b) AC-(110), (c) AC-(020), (d) BL-(011), (e) BL-(110), (f) BL-(020), (g) CH-(011), (h) CH-(110), and (i) CH-(020)

The adsorption affinity of each motif was obtained from different calculation levels, including PBE/DNP, B3LYP/6-311+G**, and M06-2X/6-311+G** (ref.40). To probe the importance of dispersion interaction, the dispersion correction was also considered in each theoretical level. As shown in Figure ACS Paragon Plus Environment

9

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 19

4, the binding energies were greatly underestimated by PBE and B3LYP methods when the dispersion correction was ignored, suggestive of the leading contribution of dispersion to adsorption interaction. In contrast, M06-2X level of theory is less sensitive to the dispersion correction and can give relatively reasonable results even without dispersion correction. The binding energies obtained from all the theoretical methods with dispersion correction show a similar change trend, namely, the adsorption affinity in the three crystal faces identically ranks as (110) > (020) > (011) for each solvent.

Figure 4 The adsorption affinities of HMX-solvent motifs calculated at different theoretical levels: (011) face in black, (110) face in red, and (020) face in blue.

Interfacial properties. The structural arrangement of solvent molecules in the interfacial region was firstly examined via the mass density profile along the normal to the crystal face. As exhibited in Figure 5, layered distribution of solvent molecules is observed in all the interfacial systems. The solvent density at the first peak, the nearest peak to the crystal face, is the highest as expected because of the strongest crystal-solvent interfacial interactions. With the solvent molecules gradually moving away from the crystal face, the density profile fluctuates weakly and the local density approaches the bulk density.

ACS Paragon Plus Environment

10

Page 11 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Mass density profiles of the solvent along the normal to the crystal face

At the same crystal face, the solvent density of the first peak ranks as BL > CH > AC, which is in consistent with the order of experimental density. Besides, the same change trend of first-peak density can be found as (011) > (110) > (020) for each solvent, revealing that the crystal face plays a critical role in the solvent density at the first peak regardless of the solvent species. To further understand the solvent structural arrangement, the molecular orientation of solvent in the interfacial region is quantified by the orientation angle θ, which is defined as the angle formed between the dipole vector of solvent molecule and the normal to HMX crystal face, i.e. z axis. Figure 6 depicts the solvent-molecule populations as a function of the cosine of orientation angle θ and the centroid position of solvent molecule projected to z axis. Evidently, preferred orientations of solvent molecule are found in all the interfaces and are mainly located at the first solvent peak, demonstrating that the solvent molecules are arranged in an ordered fashion at the crystal face.

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 19

Figure 6. Contour plots for the populations of solvent molecules in the interfacial region as a function of cosθ and z coordinate of the centroid of solvent molecule

For the same solvent, the preferred orientations of solvent molecules are dependent on the crystal faces, which can be reflected by the shape and position of red areas in Figure 6. Taking the AC solvent as an example, the distribution of preferred orientations at the (011) face is relatively spread due to the weak adsorption interactions. At the (110) face, in contrast, the preferred orientations of AC solvent are highly concentrated. Moreover, only one preferred orientation of AC solvent is found at the (110) and (020) faces. For the same crystal face, the preferred orientations of the three solvents have some similarities. Especially at the (020) face, the shape and position of preferred-orientation distributions are nearly the same for AC, BL, and CH solvents. The surface diffusion of solvent in the interfacial region is characterized by the residence probability function (P) that is defined as ACS Paragon Plus Environment

12

Page 13 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

P

N sol (t ) N sol (0)

(3)

where Nsol(t) and Nsol(0) are the number of solvent atoms present in the interfacial layer (i.e. within 5 Å of the crystal face) at time t and time zero, respectively. The longer the solvent stay in the interfacial layer, the more slowly P decays. It is notable that the re-entering of solvent into the interfacial layer is not considered in the P calculation. As shown in Figure 7, P gradually decays to 0 for all the systems, demonstrating that no solvent can continuingly stay in the interfacial layer during the time scale of 0.4 ns. At the (011) and (020) faces, the AC solvent has a faster diffusion than the other two solvents, and P rapidly reaches 0 after about 0.1 ns. Conversely, the differences of P profiles in the three solvents at the (020) face are marginal. For each solvent, interestingly, the same order of diffusion rate among the three faces is discovered, namely, (011) > (110) > (020).

Figure 7. Residence probability function (P) of solvent molecules at the HMX crystal face

To further probe the solvent incorporation at the crystal face in the course of adsorption, the solvent occupancy was calculated from the simulation trajectory. Primarily, the three-dimensional space of simulation system was fully divided into cubic grids with the size of 0.5 Å along each axis direction. For ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 19

each grid, the occupancy was then set to either 1 or 0, depending on whether the grid contained solvent atoms or not. Finally, the fractional occupancy was obtained by averaging over the simulation trajectory for the last 5 ns. A high occupancy corresponds to a region that could be continuously occupied by solvent molecules for a long simulation time, i.e. the adsorption sites. For the sake of clarity, the adsorption sites occupied by solvent molecules are defined as the regions encompassed by the isosurface with an occupancy of 0.8, which are shown as the areas colored in pink in Figure 8.

Figure 8. Snapshots for the adsorption sites (pink color) of solvent molecules at the HMX crystal faces

As clearly seen from Figure 8, the adsorption sites for solvent incorporation are concentrated at the cavities formed by the surficial HMX molecule, indicative of a dominant dependency on the structural feature of the crystal face. At the (020) face, the adsorption sites are arranged in an ordered and uniform pattern regardless of the solvent species, which is in accordance with the unique preferred orientation of solvent molecules (see Figure 6). In contrast, there are multiple distributions of the adsorption sites at the (011) and (020) faces. Furthermore, the regions of adsorption sites are larger at the (020) face than those at the other faces, suggesting that the (020) face is more favorable for the solvent residence. ACS Paragon Plus Environment

14

Page 15 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

CONCLUSIONS The molecular-scale understanding of the HMX crystal-solvent interfaces were comparatively studied via computational chemistry. The adsorption sites for solvent incorporation are concentrated at the surface cavities formed by HMX molecules, leading to some similar adsorption behaviors at the same face, such as the preferred orientations of solvent molecules at the (020) face. The adsorption affinity obtained from quantum chemistry calculations ranks as (110) > (020) > (011) for each solvent. The layered distribution of solvent molecules along the normal to the HMX crystal face was found in all the systems, and the local solvent density at the first peak was higher than the bulk density. For each solvent, the order of first-peak density is (020) < (110) < (011). The preferred dipole orientation of solvent molecules was observed for each interface, indicative of the ordered arrangements of solvent molecules. Furthermore, the residence probability function showed that the solvent diffusion is the fastest at the (011) face and is the slowest at the (020) face, which is independent on solvent species. The results reported here would provide a fundamental comprehending of HMX crystal-solvent interactions, especially the dominant effect of crystal face on the interfacial properties.

ACKNOWLEDGEMENTS The authors greatly appreciate the financial support from the National Natural Science Foundation of China (Nos. 21403162 and 21503160).

REFERENCES (1) Palmer, S. J. P.; Field, J. E. The Deformation and Fracture of β-HMX. Proc. R. Soc. London, Ser. A 1982, 383, 399−407. (2) Cady, H. H.; Larson, A. C.; Cromer, D. T. The Crystal Structure of α-HMX and a Refinement of the Structure of β-HMX. Acta Cryst. 1963, 16, 617−623. (3) Cobbledick, R. E.; Small, R. W. H. The Crystal Structure of the δ-Form of 1, 3, 5, 7-Tetranitro-1, 3, 5, 7-Tetraazacyclooctane (δ-HMX). Acta Cryst. B 1974, 30, 1918−1922. ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 19

(4) Main, P.; Cobbledick, R. E.; Small, R. W. H. Structure of the Fourth Form of 1, 3, 5, 7-Tetranitro-1, 3, 5, 7-Tetraazacyclooctane (γ-HMX), 2C4H8N8O8∙0.5H2O. Acta Cryst. C 1985, 41, 1351−1354. (5) van der Heijden, A. E. D. M.; Bouma, R. H. B. Crystallization and Characterization of RDX, HMX, and CL-20. Cryst. Growth Des. 2004, 4, 999−1007. (6) Wen, Y. S.; Xue, X. G.; Zhou, X. Q.; Guo, F.; Long, X. P.; Zhou, Y.; Li, H. Z.; Zhang, C. Y. Twin Induced Sensitivity Enhancement of HMX versus Shock: A Molecular Reactive Force Field Simulation. J. Phys. Chem. C 2013, 117, 24368−24374. (7) Zhang, C. Y.; Ji, C. L.; Li, H. Z.; Zhou, Y.; Xu, J. J.; Xu, R. J.; Luo, Y. J. Occupancy Model for Predicting the Crystal Morphologies Influenced by Solvents and Temperature, and Its Application to Nitroamine Explosives. Cryst. Growth Des. 2013, 13, 282−290. ( 8 ) Armstrong, R. W.; Ammon, H. L.; Elban, W. L.; Tsai, D. H. Investigation of Hot Spot Characteristics in Energetic Crystals. Thermochim Acta 2002, 384, 303−313. (9) Tarver, C. M.; Chidester, S. K.; Nichols, A. L. Critical Conditions for Impact-and Shock-Induced Hot Spots in Solid Explosives. J. Phys. Chem. 1996, 100, 5794−5799. ( 10 ) Kim, J. H.; Yim, Y. J. Effect of Particle Size on the Thermal Decomposition of ε-Hexanitrohexaazaisowurtzitane. J. Chem. Eng. Jpn. 1999, 32, 237−241. (11) Duan, X. H.; Wei, C. X.; Liu, Y. G.; Pei C. H. A Molecular Dynamics Simulation of Solvent Effects on the Crystal Morphology of HMX. J. Hazard. Mater. 2010, 174, 175–180. (12) Yan, T.; Wang, J. H.; Liu, Y. C.; Zhao, J.; Yuan, J.M.; Guo, J. H. Growth and Morphology of 1, 3, 5, 7-Tetranitro-1, 3, 5, 7-Tetraazacy-Clooct Ane (HMX) Crystal. J. Cryst. Growth 2015, 430, 7−13. (13) Shim, H. M.; Koo, K. K. Prediction of Growth Habit of β-Cyclotetramethylene-Tetranitramine Crystals by the First-Principles Models. Cryst. Growth Des. 2015, 15, 3983−3991. (14) Deschamps, J. R.; Frisch, M.; Parrish, D. Thermal Expansion of HMX. J. Chem. Cryst. 2011, 41, 966−970. (15) Sun, H. COMPASS: An ab initio Force-Field Optimized for Condensed-Phase Applications Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338−7364. (16) Hartman, P.; Bennema, P. The Attachment Energy as a Habit Controlling Factor: I. Theoretical Considerations. J. Cryst. Growth 1980, 49, 145–156. (17) Hartman, P. The Attachment Energy as A Habit Controlling Factor: III. Application to Corundum. J. Cryst. Growth 1980, 49, 166–170. (18) Ter Horst, J. H.; Geertman, R. M.; van der Heijden, A. E.; van Rosmalen, G. M. The Influence of a Solvent on the Crystal Morphology of RDX. J. Cryst. Growth 1999, 198, 773−779. (19) Lee, H. E.; Lee, T. B.; Kim, H. S.; Koo, K. K. Prediction of the Growth Habit of 7-Amino-4, ACS Paragon Plus Environment

16

Page 17 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

6-Dinitrobenzofuroxan Mediated by Cosolvents. Cryst. Growth Des. 2009, 10, 618−625. (20) Chen, G.; Xia, M.; Lei, W.; Wang, F. Y.; Gong, X. D. A Study of the Solvent Effect on the Morphology of Rdx Crystal by Molecular Modeling Method. J. Mol. Model. 2013, 19, 5397–5406. (21) Chen, G.; Xia, M. Z.; Lei, W.; Wang, F. Y.; Gong, X. D. Prediction of Crystal Morphology of Cyclotrimethylene Trinitramine in the Solvent Medium by Computer Simulation: A Case of Cyclohexanone Solvent. J. Phys. Chem. A 2014, 118, 11471−11478. (22) Chen, G.; Chen, C. Y.; Xia, M. Z.; Lei W.; Wang, F Y.; Gong, X. D. A Study of the Solvent Effect on the Crystal Morphology of Hexogen by Means of Molecular Dynamics Simulations. RSC Adv. 2015, 5, 25581−25589. (23) Shim, H. M.; Koo, K. K. Crystal Morphology Prediction of Hexahydro-1, 3, 5-Trinitro-1, 3, 5-Triazine by the Spiral Growth Model. Cryst. Growth Des. 2014, 14, 1802−1810. (24) Shim, H. M.; Kim, H. S.; Koo, K. K. Molecular Modeling on Supersaturation-Dependent Growth Habit of 1, 1-Diamino-2, 2-Dinitroethylene. Cryst. Growth Des. 2015, 15, 1833−1842. (25) Shi, W. Y.; Chu, Y. T.; Xia, M. Z.; Lei, W.; Wang, F. Y. Crystal Morphology Prediction of 1, 3, 3-Trinitroazetidine in Ethanol Solvent by Molecular Dynamics Simulation. J. Mol. Graph. Model. 2016, 64, 94–100. (26) Liu, N.; Li, Y. N.; Zeman, S.; Shu, Y. J.; Wang, B. Z.; Zhao, Q. L.; Wang, W. L. Crystal Morphology of 3, 4-Bis (3-Nitrofurazan-4-Yl) Furoxan (DNTF) in a Solvent System: Molecular Dynamics Simulation and Sensitivity Study. CrystEngComm 2016, 18, 2843−2851. (27) Lide, D. R. CRC Handbook of Chemistry and Physics (84th edn.); CRC; Boca Raton, 2004 (28) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. ( 29 ) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. (30) Material Studio 8.0; Acceryls Inc.: San Diego, 2014. (31) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Gaussian 09, Revision E.1; Gaussian, Inc.: Wallingford, CT, 2009. (32) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into A Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. (33) Becke, A. D. Density‐Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (34) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 19

Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. (35) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. (36) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics, 1996, 14, 33−38. (37) Politzer, P.; Murray, J. S. High Performance, Low Sensitivity: Conflicting or Compatible? Propell. Explos. Pyr. 2016, 41, 414–425. (38) Liu, Y. Z.; Lai, W. P.; Yu, T.; Ma, Y. D.; Kang, Y.; Ge, Z. X. Understanding the Growth Morphology of Explosive Crystals in Solution: Insights from Solvent Behavior at the Crystal Surface. RSC Adv. 2017, 7, 1305–1312. (39) Liu, Y. Z.; Yu, T.; Lai, W. P.; Ma, Y. D.; Kang, Y.; Ge, Z. X. Adsorption Behavior of Acetone Solvent at the HMX Crystal Faces: A Molecular Dynamics Study. J. Mol. Graph. Model. 2017, 74, 38–43. (40) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functional. Theor. Chem. Acc. 2008, 120, 215–241.

ACS Paragon Plus Environment

18

Page 19 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC Graphic

ACS Paragon Plus Environment

19