Prediction of Growth Habit of Î²-Cyclotetramethylene-tetranitramine

Article pubs.acs.org/crystal

Prediction of Growth Habit of β‑Cyclotetramethylene-tetranitramine Crystals by the First-Principles Models Hong-Min Shim and Kee-Kahb Koo* Department of Chemical and Biomolecular Engineering, Sogang University, Seoul 121-742, Korea S Supporting Information *

ABSTRACT: Experimentally, β-cyclotetramethylene-tetranitramine (βHMX) crystals were found to dramatically elongate to the [100] direction when a relatively high supersaturation was imposed. A sudden growth of β-HMX to the [100] direction is closely associated with a mechanistic transition from spiral growth to two-dimensional (2D) nucleation for the (110) face. The onset supersaturation for the growth by 2D nucleation, σ2D, was found to play a key role in the growth of βHMX. The present simulation results based on first-principles models such as the spiral growth model and the 2D nucleation model show that the values of σ2D on the (101) and (101̅) faces are smaller than those on the (020), (110), and (011) faces. This leads to the prediction of rapid growth rates for the (101) and (101̅) faces by 2D nucleation at low supersaturation and the appearance of a typical shape of βHMX. On the other hand, the needle-like shape of β-HMX begins to prevail when the supersaturation exceeds the σ2D for the (110) face because its growth mechanism is transformed from the spiral growth mechanism to the 2D nucleation mechanism which accompanies rapid growth of the (110) face. As a result, the present predictions are in remarkable agreement with the experiments. Furthermore, the kinetic Monte Carlo (KMC) simulation also shows that the σ2D for the (110) face is lower than that for the (011) face because the (011) face provides the surface topology on which growth units are unfavorably incorporated into the lattice sites. It evidently shows that the relative positions of σ2D bring on the advent of needle-like growth of β-HMX. level, but the effect of supersaturation on growth habit of β-HMX has not been tackled yet. In the present work, the needle-like crystals of β-HMX were obtained at high supersaturation induced by fast evaporation of solvent. They grew dramatically into the [100] direction where the (011) faces dominated over all crystal faces. As aforementioned, the growth of β-HMX to the [100] direction is significantly important with respect to its shock sensitivity. Therefore, much work should be devoted to the study on the growth habit of β-HMX. Currently, molecular modeling techniques on the prediction of crystal shape have been extensively developed. In particular, the predictive modeling of crystal shape based on a firstprinciples mechanistic model, the spiral growth mechanism or the 2D nucleation mechanism, is regarded as a phenomenally successful approach with providing a great coincidence with experimental observations.7−10 However, the usage of it was limited to the case where all F-faces of crystals, defined as faces with two or more periodic bond chains (PBCs), could be assumed to grow by the same growth mechanism. Therefore, it was not allowed to predict the crystal shape in moderate supersaturation where some F-faces were dominated by the spirals growth mechanism and others proceeded by the 2D nucleation mechanism. Recently, Lovette and Doherty demon-

1. INTRODUCTION Cyclotetramethylene-tetranitramine (HMX) has been known to have four polymorphs to date,1 and among them, β-HMX has been extensively investigated as a high energetic material in the formulation of polymer bonded explosives (PBXs), solid propellants, and so on because of its higher density and relatively good thermal stability.2 The importance on crystal shapes of βHMX has been emphasized until today because the shock sensitivity is strongly influenced by the anisotropy of crystal faces. The multiscale modeling carried out by Zamiri and De indicates that the shock response differs with crystallographic directions.3 Loading along the [100] direction leads to plastic deformation and results in high local temperatures whereas loading along the [010] direction shows brittleness that does not contribute to a temperature increase. Particularly, the substantial plastic deformation is observed along the [110] direction. In this regard, β-HMX can show different mechanical and thermal sensitivities under various morphologies.4 Recently, Jiang et al. obtained the 2D plate-like β-HMX crystals on the hydrophilic substrate and revealed the growth mechanism of the recess on the (011) face.5 They stressed that the emergence of the (011) face is caused by the lowest attachment energy. Zhang et al. carried out the molecular dynamics (MD) simulations for the β-HMX crystals grown from various solvents, where a high relative occupancy was ascribed to the generation of the (020) and (011) faces.6 They gave a good explanation about the crystal growth of β-HMX at the atomic © XXXX American Chemical Society

Received: May 3, 2015 Revised: July 2, 2015

A

DOI: 10.1021/acs.cgd.5b00605 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

Figure 1. (a) Unit cell of β-HMX and (b) its crystal graph.

Table 1. Crystal Graph of the β-HMX

a

bond

direction

no.

bond distance (Å)

Ed (kcal/mol)a

Ec (kcal/mol)b

Ep (kcal/mol)c

Etot (kcal/mol)d

a b c d

[1̅00] [1̅01] [001]̅ [1̅00]

2 4 2 4

6.52 7.74 7.36 7.02

−4.77 −3.52 −3.07 −1.80

−1.91 −0.57 −0.36 −0.43

−0.16 −0.11 −0.11 −0.10

−6.68 −4.09 −3.43 −2.23

Ed = dispersive energy. bEC = Coulombic energy. cEp = polarization energy. dEtot = total energy.

Figure 2. β-HMX crystals grown by evaporative crystallization with (a) slow evaporation (method 1) and (b) fast evaporation (method 2).

2. MATERIALS AND METHODS

strated a new model able to foresee the possible shapes of critically sized 2D nuclei,11 by which the onset supersaturation for 2D nucleation can be approximated. Hence, it is possible to predict the crystal shape over the whole range of supersaturations. Our aim with this work is at revealing the reason why the βHMX crystals grow with dominance of the {011} faces and they remarkably elongate into the [100] direction at higher supersaturation. For this, the prediction of a β-HMX crystal shape is carried out by the Lovette−Doherty (L−D) method and then the effect of supersaturation on the anisotropic growth behavior of β-HMX is considered. This provides a stunning qualitative agreement with experimental outcome and supports a better understanding of the growth mechanism. At last, the predicted growth shapes of β-HMX by the L−D method were compared with those by the kinetic Monte Carlo (KMC) simulation.

2.1. Experimental Details. β-HMX (99.9%, Hanwha Co., Daejeon, Korea) was used as received. Acetone (≥99.9%, Sigma-Aldrich, St. Louis, MO, USA) was used as a solvent. To obtain β-HMX crystals from different supersaturation conditions, an acetone was evaporated by two different ways. Method 1: Crystallization at Low Supersaturation Induced by Slow Evaporation. A saturated HMX/acetone solution of 10 mL at temperature of 20 °C was prepared in a 50 mL beaker, and the beaker was sealed by an aluminum foil with punched holes with diameter of about 0.8 mm. The beaker was kept at an ambient temperature of 20 °C for 2 days. In this condition, it was found that at least 24 h were required to detect HMX crystals. Method 2: Crystallization at High Supersaturation Induced by Fast Evaporation. A 4 mL saturated HMX/acetone solution was prepared in a 20 mL vial without a cap at temperature of 30 °C. The vial was placed on the rotator (IKA Vibrax VXR) for 30 min to promote evaporation at an ambient temperature of 20 °C. 2.2. Computer Modeling. In the present work, the crystal structure of β-HMX determined recently by Deschamps et al.12 was used instead of the previous one reported by Choi and Boutin.13 It should be noted B



Article

that the (111̅), (102̅), and (100) faces are newly indexed as the (110), (101), and (101̅) faces, respectively. The unit cell of β-HMX is illustrated in Figure 1a. The crystal structure of β-HMX is a monoclinic with space group P21/n and Z = 2 in the unit cell (a = 6.5255 Å, b = 11.0369 Å, c = 7.3640 Å, and β = 102.67°). In the present work, the intermolecular interactions were calculated by the GAFF force field14 and the atomic point charges were obtained by the MNDO/C semiempirical calculation.15 The polarization energy was calculated by using the AA-CLP program.16 Table 1 summarizes the bond energy between growth units contained in the β-HMX unit cell. The calculated lattice energy is −45.50 kcal/mol that gives a good agreement with the experimental value of −45.34 kcal/mol.17 The crystal graph of β-HMX in Figure 1b is displayed by the MORPHOLOGY (Material studio, version 7.0), where each colored line indicates different bonds. Through analyzing the structure of PBCs, F-faces were found to be the (020), (011), (110), (101), and (101̅) faces.

The attachment energy (AE) model considers the intermolecular interactions between growth units contained in a crystal.20 The attachment energy, Eatt hkl, is defined as the energy released when a new layer is attached on the hkl face. Rhkl is proportional to the magnitude of Eatt hkl: att R hkl ∝ |Ehkl |

(2)

From Table 2, it can be anticipated that the (020), (011), (110), and (101̅) faces with larger dhkl dominantly appear by the Table 2. Interplanar Distance and Attachment Energy of βHMX

3. RESULTS AND DISCUSSION 3.1. Experimental Morphology. Figure 2a shows the typical shape of β-HMX grown in acetone where the {011} faces dominate and they are elongated to the [100] direction.18 Furthermore, morphologically different growth shapes were observed by the fact of existence or nonexistence of the {020} face. On the other hand, the β-HMX with a high aspect ratio was obtained at higher supersaturation as shown in Figure 2b. Figure 3 displays the powder X-ray diffraction (PXRD) patterns of β-

face

dhkl (Å)

Eatt hkl (kcal/mol)

{020} {011} {110} {101} {101̅} {111}̅

5.51 6.02 5.51 4.31 5.38 4.84

−25.28 −19.50 −26.00 −29.14 −36.58 −36.95

Figure 4. Growth shapes of β-HMX predicted by the (a) BFDH and (b) AE models, respectively.

BFDH model (Figure 4a), whereas the (020), (011), (110), and (101) faces with smaller |Eatt hkl| become more prominent on a crystal predicted by the AE model (Figure 4b). The predicted growth shapes of β-HMX were found to be quite different with the experiments in which the (020), (101), and (101̅) faces almost disappeared (Figure 2a). Furthermore, the aspect ratio of crystals is much smaller than those obtained at high supersaturation (Figure 2b). A large discrepancy in prediction and experiments presumably stems from their unrealistic assumptions that the growth rates are proportional to the normal forces at crystal faces and the crystal growth is achieved by the attachment of layers on crystal faces excluding the process parameters such as solvent and supersaturation that strongly influence on the growth habit of crystals. 3.3. Lovette−Doherty Method. The L−D method allows for predicting crystal shapes by determining key anisotropic parameters for F-faces that grow by either a spiral or a 2D nucleation mechanism.11 The procedure of the L−D method is followed: (1) The growth rates of F-faces are calculated by the spiral growth model where the step velocity is mainly influenced by the kink density. (2) The growth rates of F-faces are calculated by the 2D nucleation model where the critical shape of a 2D nuclei is predicted and the onset supersaturation for 2D nucleation is estimated.

Figure 3. PXRD patterns for β-HMX grown by crystallization method: (a) slow evaporation (method 1) and (b) fast evaporation (method 2), respectively.

HMX where Miller indices were calculated by Reflex (Materials Studio, version 7.0). The β-HMX grown at higher supersaturation was noticeably characterized by the (011) and (022) peaks. Hence, the shape transition in the growth of β-HMX seems to be caused by different supersaturation levels. However, there is still an absence of the explanation why β-HMX grows into the [100] direction with the {011} basal faces at higher supersaturation. 3.2. BFDH and AE Morphology. According to the Bravais− Friedel−Donnay−Harker (BFDH) model,19 the relative growth rate of a face, Rhkl, is inversely proportional to the interplanar distance, dhkl: 1 R hkl ∝ dhkl (1) C



Article

(3) The absolute growth rates encompassing two mechanisms are determined over the specified range of supersaturations. For the spiral growth, the growth rate of a face, GS, is expressed as

GS =

h τS

H̅c, i =

(3)

ae, i ⎛ ϕi kink ⎞ ⎟ ⎜ σ ⎜⎝ kBT ⎟⎠

lc, i =

∑

ϕe =

F≡

cn ≡ (6)

1 ⎛ cnϕe ⎞ ⎜ ⎟ 12 ⎝ kBT ⎠

(13)

2

(14)

P A1/2

σ2D = −

(15)

6 F 7 W − 6 FM 6/7 −1 7

(

)

(16)

where W−1 is the −1 branch of the Lambert W function and M is the dimensionless quantity consisting of anisotropic factors which can be fixed as M ≈ 5 × 10−5 for molecular organic crystals.11 The absolute growth rates for different faces, Ghkl(σ), are determined by ⎧ ⎛ σ ⎞2 ⎪GS, hkl(σ *)⎜ ⎟ for σ ≤ σ2D, hkl ⎝ σ* ⎠ ⎪ Ghkl(σ ) = ⎨ ⎛ σ2D, hkl ⎞2 ⎪ G2D, hkl(σ ) ⎜ ⎟ G ( ) σ * S, hkl ⎪ G (σ ⎝ σ* ⎠ ⎩ S, hkl 2D, hkl)

2 exp( −ϕi kink /kBT ) /kBT )

(12)

where P is the perimeter of the nucleus and A is the area. According to the L−D method, the onset supersaturation, σ2D, is given by

where ρi can be obtained by assuming that the probability of forming a kink follows the Boltzmann distribution.25

1 + 2 exp( −ϕi

N

∑i = 1 ni

where cn is the shape factor that is defined as

where ap,i is the distance edge i that is propagated by adding a new monolayer into the edge, v0 is the frequency factor of attachment and detachment trials, and δ is the average spacing between kinks. The value (ae,i/δ)i becomes the kink density ρi that is the probability of finding a kink along the ith edge, Vm is the molecular volume, and ΔE is the barrier height for the incorporation of a molecule into a kink site, which is taken to be an isotropic because the attachment barrier is profoundly related to the breaking of the solvation shell around the growth unit.24 The growth units of β-HMX are centrosymmetric, and hence they exhibit Kossel-like behavior where XA is independent of surface. In that case, the step velocity was reduced as vi = a p, iρi (7)

kink

∑i = 1 niϕiedge

where F is defined as

where αi,i+1 is the angle between edges i and i+1 and vi is the step velocity as follows23

ρi =

sin(αi − 1, i)

G2D ∝ σ 5/6 exp( −F /σ )

(5)

⎛ ΔE ⎞ ⎛ ae, i ⎞ vi = a p, iv0⎜ ⎟ exp⎜ − ⎟Vm(XA − X A0 ) ⎝ δ ⎠i ⎝ kBT ⎠

H̅c, i − 1 − H̅c, i cos(αi − 1, i)

where ni is the number of molecules existing along the edge i and N is the number of edges on the nucleus. The growth rate of a face by 2D nucleation is given by

lc, i + 1 sin(αi , i + 1)

i=1

+

N

(4)

vi

sin(αi , i + 1)

where αi,i+1 and αi‑1,i are the angle between edges i and i+1 and between i−1 and i, respectively.11 Once the shape of a critical nucleus is obtained, the effective edge energy ϕe can be calculated by

where ae,i is the distance between molecules along edge i, is the kink energy that is defined as the energy change or work required to create one kink site along an edge i, σ is the supersaturation defined by σ = (XA/X0A) − 1, where XA is the concentration of solute molecules in solution and X0A is the equilibrium concentration, kB is the Boltzmann’s constant, and T is the absolute temperature. For a convex N-sided spiral, the rotation time is given by22 τS =

H̅c, i + 1 − H̅c, i cos(αi , i + 1)

(11)

ϕkink i

N

(10)

where H̅ c,i is the perpendicular distance from a center of the crystal to face i and vm is the molecular volume. The critical length of the ith edge can be represented as

where h is the step height that is usually regarded as an interplanar distance of a face and τS is the characteristic rotation time required for one full turn of a spiral. The edge starts to grow when the length of an edge reaches a critical value lc. The critical length of an ith edge is determined by thermodynamic approximation.21 lc, i = 2

edge vm ⎛ γi ⎞ ⎜⎜ ⎟ σ ⎝ kBT ⎟⎠

for σ > σ2D, hkl

(8)

(17)

In general, = −Etot,i/2 for centrosymmetric growth units. For the 2D nucleation, the edge energy is a key factor because it affects the free energy barrier of 2D nucleation and the shape of a critical nucleus. The edge energy γedge,cryst can be estimated by i

where σ* is an arbitrarily low supersaturation at which all faces are grown by the spiral growth model. The interfacial energy of two immiscible phases, γi, can be approximated by using the macroscopic relationship where the adhesion energy is obtained by the geometric mean rule.11,24

ϕkink i

γiedge,cryst = 2

26

energy of bonds broken area created by the two edges

cryst solvent γi = γtot, + γtot − 2 γd,cryst γ solvent − 2 γp,cryst γ solvent i i d i a

(9)

For an N-sided nucleus, Lovette and Doherty obtained

(18) D



Article

Figure 5. Bonding structures for the (a) (020), (b) (011), (c) (110), (d) (101), and (e) (101̅) faces of β-HMX, respectively, and their possible shapes for spiral and 2D nucleus.

Table 3. Simulation Results of β-HMX by the Spiral Growth Modela face

edge

ϕkink,cryst (kcal/mol) i

γkink,cryst (erg/cm2) i

γkink (erg/cm2) i

ϕkink (kcal/mol) i

ρ

lc (Å)

dhkl (Å)

Rs,hkl

(020)

[001] [100] [100] [111]̅ [11̅1] [001] [11̅1] [11̅ 1] [111̅] [111]

1.71 3.33 3.33 2.04 1.11 1.71 1.11 2.04 2.04 1.11

33.15 57.17 49.76 33.59 18.37 30.78 18.20 35.01 42.56 20.55

5.86 21.16 17.81 7.38 5.10 5.38 5.10 7.74 9.99 5.10

0.303 1.236 1.195 0.449 0.310 0.299 0.313 0.452 0.480 0.277

0.542 0.192 0.204 0.480 0.539 0.544 0.538 0.479 0.467 0.553

767 1770 2678 1195 748 758 756 1203 1276 669

5.51

1.00

6.02

0.82

5.51

1.88

4.31 5.38

1.46 3.75

(011)

(110)

(101) (101̅)

Growth shape of β-HMX was simulated by the spiral growth model at T = 293.15 K and σ = 0.01. Rs,hkl is a relative growth rate by the spiral growth model. a

Table 4. Simulation Results of β-HMX by the 2D Nucleation Modela face

edge

ϕedge,cryst (kcal/mol) i

γedge,cryst (erg/cm2) i

γedge (erg/cm2) i

ϕedge (kcal/mol) i

lc (Å)

ϕe (kcal/mol)

cn

σ2D

R2D,hkl

(020)

[001] [100] [100] [111̅] [111̅ ] [001] [11̅1] [1̅11] [111]̅ [111]

6.67 3.43 6.32 8.91 10.76 6.32 7.52 5.66 4.08 2.23

115.08 66.54 112.28 133.81 178.02 108.60 135.37 92.73 85.45 41.23

52.77 16.82 42.62 62.09 90.65 40.50 55.62 30.18 28.20 8.35

3.062 0.867 2.401 4.136 5.483 2.359 3.090 1.844 1.349 0.453

21 68 149 153 277 91 87 96 35 10

1.351

4.82

1.06

∼1020

4.301

3.55

7.63

1

2.430

3.46

1.90

∼1017

1.349 0.453

4 4

0.69 0.06

∼1021 ∼1024

(011)

(110)

(101) (101̅)

Growth shape of β-HMX was simulated by the 2D nucleation model at T = 293.15K and σ = 1.00. R2D,hkl is a relative growth rate by the 2D nucleation model. a

The associate component is comprised of polar and hydrogen bonding components with

where the subscripts d, c, and p denote dispersive, Coulombic, and polarization, respectively. The cohesion energies for solvent can be estimated by empirical correlations.27 γzsolvent

=

2

fδz Vm1/3

δa ≡ (19)

δ h 2 + δp2

(20)

For acetone, γsolvent = 17.33, γsolvent = 11.28, and γsolvent = 24.47 d a tot 2 erg/cm , where δd = 7.6, δp = 5.1, and δh = 3.4 (cal/cm3)1/2, respectively, with Vm = 74.0 cm3/mol and f = 0.0715.27 In the present work, γsolvent is estimated by a correlation on a three tot

where δz is the solubility parameter related to component z, with z = d, p, h, and a for the dispersive, polar, hydrogen bonding, and associate components, respectively, and Vm is the molar volume. E



Article

component solubility parameter for organic liquids (γsolvent = tot f V1/3[δd2 + 0.632δa2]).28 The bonding structures for all of the F-faces of β-HMX are represented in Figure 5. For the spiral growth model, the ϕkink,cryst i was obtained by analyzing the PBCs, and the ϕikink was approximated by eqs 18, 19, and 20, after which the ρi was calculated by eq 8. At given ae,i, ϕkink i , T, and σ, the lc was calculated by eq 4 and then τS could be estimated by eqs 5 and 7. Consequently, the relative growth rate, Rs, was determined by eq 3. The simulation results of β-HMX by the spiral growth model are summarized in Table 3. In this way, the γedge,cryst for the 2D i nucleation model was obtained from the bonding structures by eq 9. Considering the effect of a solvent on the edge energy, the γedge was approximated by eqs 18, 19, and 20. To predict the i shape of a critical nucleus, the critical length of the ith edge was calculated by eqs 10 and 11, which in sequence leads to the acquisition of ϕe by eq 12. The relative growth rate, R2D, was determined by eq 13. The simulation results of β-HMX by the 2D nucleation model are summarized in Table 4 where the σ2D was obtained by solving eq 16. The possible shapes for the spiral and 2D nucleus in theory are also shown in Figure 5. Figure 6 presents the observation of spirals for the (020) and (011) faces that are quite similar to the predicted spirals resulting

the aspect ratio is excessively higher than those shown in Figure 2b. The reason why the 2D nucleation model predicts a needlelike shape stems from the highest value of edge energy on the (011) face. F-faces with large edge energy (ϕe > 3RT) were found to be the (011) and (110) faces among which the (011) face results in the highest value. Therefore, the growth rate of (011) face becomes considerably slower in comparison with the other faces. Figure 8 demonstrates the growth shapes of β-HMX predicted by the L−D method where the Ghkl(σ) was determined by eq 17 at σ* = 0.01. The predicted growth shapes in the range of σ ≤ 2.0 were found to be in good agreement with the experimental morphology of β-HMX grown at a moderate supersaturation (Figure 2a). It was noticed that the (101) and (101)̅ faces of βHMX disappear when σ < 1.0 because they grow rapidly by 2D nucleation. The (020) face began to disappear near σ = 1.0 above which it tends to grow by 2D nucleation. Interestingly, the aspect ratio of β-HMX was found to increase when σ > 1.9. The L−D method accurately predicted the growth shape of β-HMX with a large aspect ratio when σ = 2.7. It implies that the elongation of βHMX at high supersaturation is caused by the commencement of the 2D nucleation mechanism on the (110) face. Experimentally, the {110} surfaces of β-HMX were also found to be rough at higher supersaturation as shown in Figure 9. Unequivocally, the change in growth habit can be observed only when the growth mechanism of F-faces is transformed from a spiral growth to a 2D nucleation. When all of the F-faces are grown by the spiral growth mechanism, the crystal shape is constant because the supersaturation term, σ, containing the expression of a critical length, lc, in eq 4, drops out in the calculation of relative growth rates, Rs, where the kink density significantly affects the growth shape. For the noncentrosymmetric growth units, albeit supersaturation has an impact on the kink rate, the effect on the change in a growth shape was found to be negligibly small.9,29 However, when supersaturation is large enough to overcome the energy barrier of forming 2D nuclei, it brings on the mechanistic transition from a spiral growth to a 2D nucleation. Such progression of growth mechanisms with supersaturation has been reported by Land et al.30 and Land and De Yoreo.31 The σ2D for β-HMX was estimated by eq 16 on which the effective edge energy, ϕe, exerted a strong influence. In particular, for organic crystals, the bonding structures consisting of growth units were found to evidently differ with respect to crystal faces and it leads to the different edge energy. Therefore, the edge energy has been referred to as a key anisotropic factor to determine the relative growth rates of faces growing by the 2D nucleation mechanism. The most important thing that should be remembered in this work is the dependency of σ2D upon the bonding structures of crystal faces. Through analyzing the simulation results, the σ2D seems to become larger as the ϕe

Figure 6. Optical observation of spirals for the (a) (020) and (b) (011) faces of β-HMX crystals obtained by slow evaporation (method 1).

in a high accordance of spiral edges comprising PBCs. According to the spiral growth model, the kink density is a key anisotropic factor to have impact on the step velocity of a layer, and it determines the growth shape of a crystal. The (020) and (011) faces exhibit a rather larger kink energy originated from the PBC retaining bond a, which leads to a lower kink density according to eq 8. Therefore, the (020) and (011) faces grow very slowly compared to the other faces. In Figure 7a, the growth shape of βHMX predicted by the spiral growth model resembles the shape predicted by the AE model. However, it does not coincide with the experimental outcome where the (101) face disappears. When employing the 2D nucleation model, a needle-like βHMX was predicted as shown in Figure 7b. The disappearance of the (101) faces is consistent with the experimental results, but

Figure 7. Growth shapes of β-HMX predicted by the (a) spiral growth model and (b) 2D nucleation model, respectively. F



Article

Figure 8. Predicted growth shapes for β-HMX at (a) σ < 1.0, (b) 1.0 < σ < 1.9, (c) σ = 2.0, and (d) σ = 2.7, respectively.

The detachment rate of a growth unit, which is dependent on the specific lattice site i where a growth unit exists, can be represented as ⎛ ΔWi ⎞ K i− = K 0− exp⎜ − ⎟ ⎝ kBT ⎠

where ΔWi is the energy required to remove a molecule from a site i. At the equilibrium (Δμ = 0), the attachment rate is K+0 , which becomes K−0 exp(−ΔW/kT), where K−0 is the frequency of detachment trials. The average ΔW is half the lattice energy which can be substituted by the dissolution enthalpy.

Figure 9. Surface morphology of the {110} faces of β-HMX grown by (a) method 1 and (b) method 2, respectively.

⎛ E ⎞ K 0+ = K 0− exp⎜ − latt ⎟ ⎝ 2kBT ⎠

increases. Thermodynamically, it is reasonable to prefer a lower edge energy to create a 2D nucleus because the total free energy of formation should reach a maximum value where the free energy relevant to the increase in the area should be minimized. Accordingly, the σ2D should become larger enough to form 2D nuclei when the crystal face retains a large value of ϕe. In the cooling crystallization of β-HMX,32 the typical growth shapes resemble those predicted by the L−D method at σ ≤ 2.0. However, the needle-like β-HMX could be observed only by the crystallization with rapid evaporation of solvent (Figure 2b), which means that method 2 reaches a relatively high supersaturation able to trigger 2D nucleation on the (110) face of βHMX. 3.4. Monte Carlo Simulation. At first, the total energies, Einf tot,i, at the crystal−solvent interface were calculated by using eqs 18, 19, and 20 assuming that the area per growth unit exposed to a solvent is the same regardless of directions,29 and then they were rescaled as the value, Escaled tot,i , on the basis of the dissolution enthalpy in acetone by eq 21. The dissolution enthalpy was estimated as 2.9 kcal/mol from our solubility data (Supporting Information Figure S1) by using the ideal solution theory.33 scaled inf Etot, i = Etot, i

(23)

(24)

The KMC simulations considered the only scenarios of the attachment and detachment of growth units on a crystal face without surface diffusion. The growth rate could be correlated by the sticking fraction defined as Shkl = (Natt − Ndet)/Natt, where Natt is the number of attached growth units on the (hkl) face and Ndet is the number of detached growth units from the (hkl) face. The used algorithm is based on the n-fold way proposed by Bortz et al.,36 where the number n of events is classified with their probability. All possible events were screened in advance by a random insertion, and a majority of situations where insertion and annihilation of a growth unit competitively occur on a crystal face was considered.37 Figure 10 represents the sticking fraction for each crystal face with respect to driving forces. The (101̅) face immediately begins to grow in a rough mode at the lowest driving force whereas the

ΔHdiss inf |E latt |

(21)

inf where Einf latt is the lattice energy obtained by summing over all Etot,i.

In the present work, the KMC simulations were carried out on the basis of the solid-to-solid (STS) model.34 According to Gilmer’s approach,35 the attachment rate of a growth unit was defined as

⎛ Δμ ⎞ K + = K 0+ exp⎜ ⎟ ⎝ kBT ⎠

(22)

Figure 10. KMC simulation results for the F-faces of β-HMX. The sticking fraction is plotted as a function of Δμ/kT. The simulations were carried out on a crystal face of 20 × 20 unit cells during 1,000 steps.

where is the frequency of attachment trials and Δμ is the chemical potential. K+0

G



Article

Figure 11. Growth shape of β-HMX predicted by the KMC simulations at Δμ/kT = 2.9.

In the present work, the molecular modeling based on a firstprinciples model encompassing the spiral growth mechanism and 2D nucleation mechanism is shown to give remarkable agreement with experimental data. It enables us to analyze the bonding structure for each crystal face and to predict the growth habit of β-HMX over the specified range in supersaturations. However, it still remains the uncertainty about the prediction of crystal shapes by a first-principles model although the predicted shapes are phenomenally in good agreement with experiments. The shape of a nucleus becomes unclear as the supersaturation increases, which exacerbates the approximation of the effective edge energy, ϕe. Moreover, there is an absence in the information on the growth by bulk transport of molecules to the surface that depends on the concentration, diffusion coefficient, and thickness of the depletion boundary layer.42 The present work does not consider the onset of bulk-transport-limited growth, σBT, above which all F-faces are rapidly grown by roughening regardless of growth mechanisms. The needle-like β-HMX crystals were produced by fast evaporation of solvent with facilitating a rotator. Fortunately, the supersaturation seems to not reach the σBT. Nevertheless it should be remembered that the simulation results are solely applicable for the range of supersaturations where F-faces are grown by layered growth.

(011) face appears to grow at much higher driving force of Δμ/ kT = 2.9. Figure 11 displays the predicted growth shape of βHMX that resembles the β-HMX crystals grown at relatively high supersaturation (Figure 2b). However, the KMC simulation could not predict the typical shapes of β-HMX at lower supersaturation. The KMC simulations show that it is extraordinarily difficult for the growth units to attach on the (011) face and this leads to a low sticking fraction. However, the (110) face begins kinetic roughening prior to the (011) face, and hence the advent of needle-like β-HMX crystals is observed. As a result, the predicted growth shape by the KMC simulations seems to be in good agreement with that by the L−D method at high supersaturation. Theoretically, the sticking fraction is considered as an important kinetic parameter to describe the insertion and annihilation of growth units at a flat surface or kinks that are energetically different. Therefore, the sticking fraction as well as the edge energy has been used for estimating relative growth rates of F-faces.38,39 The net flux of insertion and annihilation is fundamentally associated with the local concentration of solute on the crystal face, which is not equal to the bulk concentration from solution.40 The local concentration for the (hkl) face stems from the dissolution enthalpy of the (hkl) face, ΔHhkl diss, which is determined by the required work for detaching growth units from a (hkl) face. The larger the edge energy, the local concentration becomes typically lower,41 in which the KMC simulation gives a low sticking fraction. Resultantly, the (011) face provides a much lower local concentration, and it increases the threshold driving force for growth. Therefore, the β-HMX grows into the [100] direction with the {011} basal faces.

■

ASSOCIATED CONTENT

S Supporting Information *

β-HMX solubility in acetone and theoretical growth rates of Ffaces by the L−D method. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.cgd.5b00605.

■

4. CONCLUSIONS The present study successfully demonstrates the effect of supersaturation on the growth habit of β-HMX by employing a first-principles model. For the β-HMX, the (101) and (101̅) faces were found to have much lower edge energy; thus, the onset supersaturation for 2D nucleation, σ2D, becomes smaller. Therefore, those faces tend to more rapidly grow by the 2D nucleation mechanism and finally disappear. The advent of needle-like β-HMX was found to be caused by the transformation of growth mechanisms of the (110) face. Theoretically, the typical growth shapes of β-HMX can be observed only when the (011) and (110) faces are equally grown by the spiral growth mechanism. However, if the supersaturation reaches the σ2D of the (110) face, the β-HMX will grow much more rapidly to the [100] direction. This leads to the emergence of the needlelike β-HMX. Furthermore, according to the KMC simulation results, it was revealed that the (011) face provided a surface topology on which growth units are hardly attached on the crystal face. Because the crystal face providing a lower sticking fraction offers a lower local concentration, it is difficult for the (011) face to grow at lower supersaturation. Consequently, the KMC simulation also underpins the results predicted by the L−D method.

AUTHOR INFORMATION

Corresponding Author

*Phone: +82-2-705-8680. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by Hanwha and the Agency for Defense Development (Grant UC120019GD; Republic of Korea) and by the Next-Generation Converged Energy Materials Research Center.

■

REFERENCES

(1) Zhu, W.; Xiao, J.; Ji, G.; Zhao, F.; Xiao, H. J. Phys. Chem. B 2007, 111, 12715−12722. (2) Gallagher, H. G.; Miller, J. C.; Sheen, D. B.; Sherwood, J. N.; Vrcelj, R. M. Chem. Cent. J. 2015, 9, 22. (3) Zamiri, A. R.; De, S. Interact. Multiscale Mech. 2011, 4, 139−153. (4) Zhang, C.; Peng, Q.; Wang, L.; Wang, X. Propellants, Explos., Pyrotech. 2010, 35, 561−566. (5) Jiang, Y.; Xu, J.; Zhang, H.; Liu, Y.; Pu, L.; Li, H.; Liu, X.; Sun, J. Cryst. Growth Des. 2014, 14, 2172−2178. (6) Zhang, C.; Ji, C.; Li, H.; Zhou, Y.; Xu, J.; Xu, R.; Li, J.; Luo, Y. Cryst. Growth Des. 2013, 13, 282−290. H



Article

(7) Bisker-Leib, V.; Doherty, M. F. Cryst. Growth Des. 2003, 3, 221− 237. (8) Grimbergen, R. F. P.; Reedijk, M. F.; Meekes, H.; Bennema, P. J. Phys. Chem. B 1998, 102, 2646−2653. (9) Kuvadia, Z. B.; Doherty, M. F. Cryst. Growth Des. 2011, 11, 2780− 2802. (10) Dandekar, P.; Doherty, M. F. AIChE J. 2014, 60, 3720−3731. (11) Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2012, 12, 656− 669. (12) Deschamps, J. R.; Frisch, M.; Parrish, D. CCDC 792930: Experimental Crystal Structure Determination, 2011. (13) Choi, C. S.; Boutin, H. P. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1970, 26, 1235−1240. (14) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (15) Thiel, W. J. Am. Chem. Soc. 1981, 103, 1413−1420. (16) Gavezzotti, A. New J. Chem. 2011, 35, 1360−1368. (17) Qian, W.; Zhang, C.; Xiong, Y.; Zong, H.; Zhang, W.; Shu, Y. Cent. Eur. J. Energ. Mater. 2014, 11, 59−81. (18) Gallagher, H. G.; Sherwood, J. N.; Vrcelj, R. M. Chem. Cent. J. 2014, 8, 75. (19) Donnay, J. D. H.; Harker, D. Am. Mineral. 1937, 22, 446−467. (20) Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145−156. (21) Lovette, M. A.; Doherty, M. F. J. Cryst. Growth 2011, 327, 117− 126. (22) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. London, Ser. A 2009, 465, 1145−1171. (23) Markov, I. V. Crystal Growth for Beginners; World Scientific: Singapore, 1995. (24) Kim, S. H.; Dandekar, P.; Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2014, 14, 2460−2467. (25) Burton, W. K.; Cabrera, N.; Frank, F. C. Philos. Trans. R. Soc., A 1951, 243, 299−358. (26) Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2013, 13, 3341− 3352. (27) Barton, A. F. M. Chem. Rev. 1975, 75, 731−753. (28) Beerbower, A. J. Colloid Interface Sci. 1971, 35, 126−132. (29) Shim, H.-M.; Koo, K.-K. Cryst. Growth Des. 2014, 14, 1802−1810. (30) Land, T. A.; Malkin, A. J.; Kuznetsov, Yu. G.; McPherson, A.; De Yoreo, J. J. Phys. Rev. Lett. 1995, 75, 2774−2777. (31) Land, T. A.; De Yoreo, J. J. J. Cryst. Growth 2000, 208, 623−637. (32) Cao, X.; Duan, X.; Pei, C. Cryst. Res. Technol. 2013, 48, 29−37. (33) Prausnitz, J. M.; Lichtenthaler, R. N.; Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria; Prentice Hall: Englewood Cliffs, NJ, USA, 1986; Chapter 9. (34) Boerrigter, S. X. M.; Josten, G. P. H.; van de Streek, J.; Hollander, F. F. A.; Los, J.; Cuppen, H. M.; Bennema, P.; Meekes, H. J. Phys. Chem. A 2004, 108, 5894−5902. (35) Gilmer, G. H. J. Cryst. Growth 1976, 36, 15−28. (36) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10−18. (37) Shim, H.-M.; Kim, H.-S.; Koo, K.-K. Cryst. Growth Des. 2015, 15, 1833−1842. (38) Cuppen, H. M.; van Eerd, A. R. T; Meekes, H. Cryst. Growth Des. 2004, 4, 989−997. (39) Deij, M. A.; Meekes, H.; Vlieg, E. Cryst. Growth Des. 2007, 7, 1949−1957. (40) Winn, D.; Doherty, M. F. AIChE J. 2000, 46, 1348−1367. (41) Liu, X. Y.; Boek, E. S.; Briels, W. J.; Bennema, P. J. Chem. Phys. 1995, 103, 3747−3754. (42) Lovette, M. A.; Browning, A. R.; Griffin, D. W.; Sizemore, J. P.; Snyder, R. C.; Doherty, M. F. Ind. Eng. Chem. Res. 2008, 47, 9812−9833.

I


Prediction of Growth Habit of Î²-Cyclotetramethylene-tetranitramine

Recommend Documents