Crystal Morphology Prediction of Hexahydro-1,3,5 ... - ACS Publications

Article pubs.acs.org/crystal

Crystal Morphology Prediction of Hexahydro-1,3,5-trinitro-1,3,5triazine by the Spiral Growth Model Hong-Min Shim and Kee-Kahb Koo* Department of Chemical and Biomolecular Engineering, Sogang University, Seoul 121-742, Korea S Supporting Information *

ABSTRACT: The crystal morphology of hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX) was predicted by the advanced Burton− Cabrera−Frank (BCF) model with consideration of non-centrosymmetric growth units. The present modeling showed that the advanced BCF model provides reliable results and understanding of the growth habit of RDX crystals grown from acetone and γbutyrolactone, in which the {210} and {111} faces are dominantly developed. In the present work, the kink rate, which is the net flux input and output from kink sites along the edge, is found to be a critical factor for the crystal growth of RDX, where the energy required for detaching the growth unit from a kink site determines the existence probability of a growth unit. The {210} and {111} faces were found to show slow kink rates originated from the kink sites with a low existence probability and consequently became morphologically important faces. Moreover, the effect of solvent on the crystal growth habit of RDX, the difference in morphology of RDX grown from acetone and γ-butyrolactone, was confirmed by utilizing local concentration contained in the kink rate term. dominantly generated from acetone and γ-butyrolactone. The disagreement between experimental and predicted morphology comes from the selected models including the Bravais− Friedel−Donnay−Harker (BFDH) model11 and AE model, in which the crystal growth kinetics is not considered.12−14 Even though there is a technique available for treating the growth kinetics, such as spiral growth model of BCF,15 the morphology prediction of RDX has been referred to as an extraordinarily difficult problem due to the non-Kossel growth unit and complex bonding network. Many efforts have been made in recent years to deal with non-Kossel growth unit.16−18 The crystal containing Kossel growth unit exhibits the same energy change when growth units are detached from kink sites regardless of crystal faces, and it leads to the same kink rate. Therefore, the kink density dominantly affects the step velocity. However, the crystal containing non-Kossel growth units has energetically different growth units for each kink site. For this reason, the step velocity is markedly affected by not only kink density but also kink rate. Fortunately, if the non-Kossel growth unit is centrosymmetric, it behaves like a Kossel growth unit. However, if the non-Kossel

1. INTRODUCTION Hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX) is one of the high energetic materials widely used in military and industry. In general, the crystal size and morphology of RDX are controlled by various recrystallization techniques such as cooling crystallization, drowning-out, and evaporation crystallization because they considerably affect the final performance and shock sensitivity.1−6 To control the crystal shape efficiently, the crystal growth mechanism at the molecular level should be understood because crystals are grown by integration events of solute molecules onto crystal faces. A few studies on crystal growth of RDX at the molecular level have been reported so far. Ter Host et al. tried to explain the morphology change of RDX grown in different solvents using Monte Carlo and molecular dynamics (MD) simulations.7,8 They suggested that the adsorption of solvent molecules on a crystal face influences the morphology of RDX. Zhang et al.9 performed morphology prediction of RDX using a modified attachment energy (AE) model10 with relative occupancy and studied the interactions between a solute or solvent molecule and a crystal face. Those approaches provided a good explanation for the relationship between the solvent/ crystal face interaction and crystal morphology. However, the predicted shape of RDX does not match the experimental results where the {210} and {111} faces are © 2014 American Chemical Society

Received: December 27, 2013 Revised: February 18, 2014 Published: March 12, 2014 1802

dx.doi.org/10.1021/cg401928m | Cryst. Growth Des. 2014, 14, 1802−1810

Crystal Growth & Design

Article

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

growth unit is non-centrosymmetric, various kink rates should be considered for each edge direction. Snyder and Doherty19 predicted the morphology of a crystal constructed by centrosymmetric growth units such as naphthalene and anthracene. The morphology of a crystal consisting of noncentrosymmetric growth units was also predicted by assuming that the growth units form a dimer with each other and behave as centrosymmetric growth units. However, this approach cannot be directly applied to the morphology prediction of RDX because the RDX molecules are non-centrosymmetric growth units and cannot form a dimer. Recently, a very remarkable theoretical technique available for treating the non-centrosymmetric growth units was developed by Kuvadia and Doherty.20 They provided an important theoretical approach to calculate kink density exposing different bonds along the edge and various kink rates for paracetamol and lovastatin crystals. In the present work, we address the detailed procedure for crystal morphology prediction of RDX grown from acetone and γ-butyrolactone. First, the spiral growth model, which is amenable to dealing with non-centrosymmetric growth units, is introduced. Second, for morphologically important faces of RDX, the kink density and kink rate are calculated, and consequently the relative growth rates of faces are estimated. Lastly, the predicted morphologies by the spiral growth model including the BFDH and AE models are compared with the experimental data.

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

where ap,i is the distance which is propagated by adding a monolayer to the ith edge, v0 is the frequency of attachment and detachment trials, and δ is the average spacing between kinks. The value (ae/δ)i is the kink density, ρi, which is the probability of finding a kink along the ith edge, ΔE is the energy barrier for the incorporation of a molecule into a kink site, kb is the Boltzmann constant, and Vm is the molecular volume. The kink energy and surface energy for edge i have a relationship as below

h τ

Table 1. Solubility Parameters for Acetone and γButyrolactone24

(2)

∑ i=1

solvent

δd (cal/cm3)1/2

δp(cal/cm3)1/2

Vm (cm3/mol)

acetone γ-butyrolactone

8.23 9.42

6.98 10.36

74.0 76.8

(3)

2.2. Kink Density. In the crystal consisting of centrosymmetric growth units, the kink rates have the same value independent of kink sites because the bonds exposed at any kink sites are the same. This term drops out in the calculation of the relative growth rates, and the step velocity is dependent only on ap,i and ρi. The probability of finding a kink site in a step follows the Boltzmann distribution.15

where ae,i is the distance between molecules along edge i, ⟨ϕkink,i⟩ is the average kink energy, Δμ is the chemical potential between the liquid and solid state, which is equal to RT ln(S + 1) where R is the ideal gas constant, T is absolute temperature, and S is the supersaturation, defined by S = (XA/X0A) − 1 where XA is the concentration of solute molecules in solution and X0A is the equilibrium concentration. For a convex N-sided spiral, the rotation time is given by τ=

(9)

where δz is the solubility parameter related to the component z, where z represents d and p for the dispersive and polar bonding components, respectively, and Vm is the molar volume (Table 1). In the present work, these energies were calculated by using eq 9 assuming f = 0.14.

2ae, i⟨ϕkink, i⟩

N

(8)

γz = fδz2Vm1/3

The critical lengths for all the spiral edges are determined by thermodynamic approximation.21

Δμ

(7)

where the γcryst,d and γcryst,c are interfacial energies for dispersive and Coulombic bonding components of a crystal, and the γsolv,d and γsolv,p are the interfacial energies for dispersive and polar bonding components. The cohesion energies for each solvent were estimated by empirical correlations.

⎪

lc, i =

Akink, i = a p, ih

− 2 γcryst,cγsolv,p

where h is step height, which is usually regarded as an interplanar distance of a face, and τ is characteristic rotation time of the spiral, which is the time required for one full turn of a spiral. To start to grow for an ith edge, the length of an edge has to reach a critical length lc.

⎪

(6)

γcryst,solv = γcryst,d + γsolv,d + γcryst,c + γsolv,p − 2 γcryst,dγsolv,d

(1)

v = 0, l ≤ lc ⎫ ⎬ v = v∞ , l > lc ⎭

ϕkink, i = γkink, iAkink, i

Since a crystal is growing from solution contacting a liquid phase, the interfacial energy (γcryst,solv) of two immiscible phases is approximated by the cohesion energies for each phase and their adhesion energy which is approximated by the geometric mean rule.23

2. SPIRAL GROWTH MECHANISM 2.1. Theoretical Backgrounds. In the spiral growth model, the growth rates of crystal faces depend on the step height and rotation time of a spiral. The growth rate G can be expressed as G=

(5)

ρ=

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

2 exp( −ϕkink /k bT ) 1 + 2 exp( −ϕkink /k bT )

(10)

where ϕkink is the work or energy change to create a kink along the edge.

(4) 1803



Article

In the present work, the approach proposed by Kuvadia and Doherty.20 was employed. If two different bonds, ϕ1 and ϕ2, are connected alternatively along the edge, the kink sites possible to be formed from the straight edge are classified into five cases as follows: (1) Four kink sites exposing ϕ1,up, ϕ1,down, ϕ1,up, and ϕ1,down (2) Four kink sites exposing ϕ1,up, ϕ1,down, ϕ2,up, and ϕ2,down (3) Four kink sites exposing ϕ2,up, ϕ2,down, ϕ1,up, and ϕ1,down (4) Four kink sites exposing ϕ2,up, ϕ2,down, ϕ2,up, and ϕ2,down (5) No kink sites Each probability was weighted by Boltzmann distribution. The average expense energy per kink are written as ε1 = (ϕ1 + ϕ1)/4, ε2 = (ϕ1 + ϕ2)/4, ε3 = (ϕ2 + ϕ1)/4, ε4 = (ϕ2 + ϕ2)/4, respectively. The density is generally written as the summation of all possible probability to form a different kink site. The denominator of the equation is summation over all probability to form kink sites, and the numerator is summation of probability to form kink sites except for the zero kink probability.

VmX A0 =

⎛ ΔWkink ⎞ VmX A0 exp⎜ − ⎟= ⎝ k bT ⎠ 1 − VmX A0

The VmXA in a supersaturated solution is The step velocity can be expressed as

⎛ V X 0 ⎞⎞ m A ⎟⎟ ⎜ 0 ⎟ ⎝ 1 − VmX A ⎠⎠

For a dilute solution, ≪ 1, eq 19 becomes eq 5. For centrosymmetric growth units, the ΔWkink is half the lattice energy. However, the ΔWkink is anisotropic for noncentrosymmetric growth units that gives rise to the local concentration dependent on crystal face, XA(hkl), as suggested by Schwinn et al.27 and Curry and Cushman.28 There exists the local concentration of the effective growth units near the crystal surface participating efficiently in dynamic equilibrium with surface-located molecules. For a vacuum, the ΔWkink could be directly obtained by summation over all bonds connected to the growth unit at a kink site. However, the crystal grown from solution is surrounded by solvent molecules. Therefore, the bond energy at the crystal−solvent interface was estimated using eqs 6 and 8 assuming that the area per growth unit exposed to a solvent is the same regardless of direction. In this work, the average area of a RDX molecule, A, can be written by A=

(12)

j− = (1 − VmXA )k−

(14)

⎛ ΔE ⎞ k+ = v0 exp⎜ − ⎟ ⎝ k bT ⎠

(15)

⎛ ΔE + ΔWkink ⎞ k− = v0 exp⎜ − ⎟ k bT ⎝ ⎠

(16)

⎛ V ⎞2/3 ⎜ ⎟ ⎝Z⎠

(20)

where V is the volume of a RDX unit cell and Z is the number of molecules in the unit cell. The bond energy at the crystal− solvent interface was used in the calculation of ΔWkink. On the other hand, the activation energy cannot be directly obtained from bond chains without MD simulation. In this work, it was assumed that the activation energy has little effect on the energy barrier. The assumption is not valid for all crystals. However, it has been reported that the anisotropic physical properties on crystal faces can be compensated by the anisotropic ΔWkink, which dominantly affects the crystal morphology.20 Hence, in the present work, the kink rate was reduced as u′ = u/[v0 exp(−ΔE/kbT)] leaving the only anisotropic terms. For the kink rate calculation of non-centrosymmetric growth units, we employed the deposition mechanism of the AB crystal in which each kink site has a different probability of a growth unit.29 For example, two kink sites exist along the edge at the equilibrium state; the kink rate can be written by

with (13)

(19)

VmX0A

(11)

j+ = VmXA k+

(S + 1).

⎛ ΔE ⎞⎛ = a p, iρi v0 exp⎜ − ⎟⎜⎜VmXA − (1 − VmXA ) ⎝ k bT ⎠⎝

2.3. Kink Rates. In the solution crystallization where diffusion is very fast and the diffusion boundary layer is diminishing, the rate of incorporation of a solute molecule at kinks becomes a limiting factor. In a statistical mechanical formalism, the net flux of entering growth units can be expressed as25,26

u = j + − j−

(18)

VmX0A

vi = a p, iρi ui

⎡ ⎛ |ε | ⎞ ⎛ |ε | ⎞ ⎢1 + exp⎜ − 1 ⎟ + exp⎜ − 2 ⎟ ⎢⎣ ⎝ k bT ⎠ ⎝ k bT ⎠

⎛ |ε | ⎞ ⎛ | ε | ⎞⎤ + exp⎜ − 3 ⎟ + exp⎜ − 4 ⎟⎥ ⎝ k bT ⎠ ⎝ k bT ⎠⎥⎦

(17)

thus

⎡ ⎛ |ε | ⎞ ⎛ |ε | ⎞ ⎛ |ε | ⎞ ρ = ⎢exp⎜ − 1 ⎟ + exp⎜ − 2 ⎟ + exp⎜ − 3 ⎟ ⎢⎣ ⎝ k bT ⎠ ⎝ k bT ⎠ ⎝ k bT ⎠ ⎛ | ε | ⎞⎤ + exp⎜ − 4 ⎟⎥ ⎝ k bT ⎠⎥⎦

k− k+ + k−

where j+ and j− are the fluxes into and out of a kink site, respectively, k+ and k− are the rate constants for the attachment and detachment of growth units into and out of a kink site, respectively, v0 is the vibrational frequency of attachment and detachment attempts, ΔE is the energy barrier for transition state of a molecule into a kink site, and ΔWkink is the energy required for removing a solute molecule from a kink site. The VmXA is the value involved in the probability of a solute molecule available for incorporating into a kink site. Therefore, (1 − VmXA) is the value representing the solvation ability of a solvent molecule from a kink site. The quantity VmX0A at equilibrium can be obtained using the mass balance of fluxes (j+ = j−) 1804

u1 = j+ p1 − j2− p2

(21)

u 2 = j+ p2 − j1− p1

(22)

u1 = u 2

(23)

p1 + p2 = 1

(24)



Article

Figure 1. (a) Molecular structure of RDX. (b) The unit cell and crystal graph of RDX. (Different colored lines represent different bond chains.)

Table 2. Crystal Graph of the RDX bond

direction [uvw]

multiplicity

bond distance (Å)

Coulombic energy (kcal/mol)

dispersive energy (kcal/mol)

total energy (kcal/mol)

a b c d e f

[1̅00] [001]̅ [01̅0] [1̅00] [11̅ 0̅ ] [1̅01̅]

1 2 2 2 2 1

4.32 7.18 7.46 6.77 6.28 7.20

−0.09 −0.09 −0.49 −0.35 −0.50 −0.86

−4.61 −1.58 −1.31 −1.62 −2.81 −1.64

−4.70 −1.67 −1.80 −1.97 −3.31 −2.50

Figure 2. The bonding networks for (a) the {200}_1 and (b) {200}_2 faces of RDX. (Circle and cubic represent different growth units.)

⎛ j+ + j − ⎞ k+1 ⎟ Pk + 2 = ⎜⎜ ⎟Pk + 1 − − j ⎝ k+2 ⎠

where u1 and u2 are the kink rates at kink site 1 and 2, j+ is the input flux into the kink site, j−1 and j−2 are the output fluxes from the kink site 1 and 2, p1 and p2 are the probabilities that the chain is terminated by growth unit 1 and 2. In a similar way, the governing equation dealing with the edge consisted of various kink sites could be expanded by imposing the probability to each kink site. If n different types of kink sites are incorporated in series in a single edge direction i, the governing equation relating the probability and rate of the kth kink is written by20 ui = (j+ pk − jk−+ 1 pk + 1 )i

⎛ j+ ⎞ ⎜⎜ − ⎟⎟Pk ⎝ jk + 2 ⎠

(26)

n

∑ Pk = 1 k=1

(27)

3. COMPUTATIONAL METHODS The molecular structure of RDX is shown in Figure 1a. The crystal structure of RDX is an orthorhombic space group Pbca with Z = 8 in the unit cell, and cell dimensions are a = 13.18 Å, b = 11.57 Å, and c = 10.71 Å.30 In the present work, atomic point charges were obtained from electrostatic potential (ESP) using the PM6 semiempirical method,31 and the dispersive

(25) 1805



Article

Figure 3. The bonding networks for (a) the {020} and (b) {002} faces of RDX.

Figure 4. The bonding networks of the {210} face; (a) top view of the face, (b) side view of the face, (c) schematic image of incorporation of a growth unit into a kink site.

energy was calculated by AMBER force field. From the experimental sublimation enthalpy of 30.1 kcal/mol,32 the lattice energy was estimated as about −26.1 kcal/mol.33 In the present work, the calculated lattice energy, which was obtained by energy minimization, was −24.7 kcal/mol. The crystal graph of RDX in the solid state was constructed by using the program Morphology as shown in Figure 1b. The total energy including the multiplicity, bond distance, Coulombic energy, and dispersive energy are written in Table 2. The F-faces containing two or more PBCs in one layer were chosen as the {200}, {020}, {002}, {210}, and {111} faces.

determining a stable face that was generated in the process of spiral growth, the additional consideration of a lateral energy was employed. It was assumed that the layer with a strong lateral energy was maintained to be stable during spiral growth. In a vacuum, the calculated lateral energies per growth unit for the {200}_1 and {200}_2 faces are −6.94 and −13.82 kcal/ mol, respectively. However, the lateral energies of both faces are −2.16 and −2.03 kcal/mol from acetone, and −7.13 and −5.42 kcal/mol from γ-butyrolactone, respectively. Therefore, the {200}_1 face with larger lateral energies from both solvents was selected as a stable F-face. The kink density for each edge direction was calculated using eq 10. However, when estimating the kink rate, the additional efforts should be made to consider various kink sites. As shown in Figure 2, the growth units do not exist as a flat surface.

4. GROWTH RATE CALCULATION 4.1. The {200} Faces. Two different types of F-faces for the {200} faces exist as shown in Figure 2. As a means of 1806



Article

even in the same direction, the averaged kink rate of the ⟨001⟩ edge was used in the calculation of step velocity. 4.4. The {111} Faces. As shown in Figure 5, three different types of PBCs are located on the {111} faces. The strongest

Hence, it leads to the energy difference between the upper and lower levels even in one layer. Along the ⟨001⟩ or ⟨010⟩ direction, the energies required for removing a growth unit from each kink site in a vacuum were calculated as −5.44 kcal/ mol and −19.26 kcal/mol, respectively. In that case, the net kink rate was calculated by using the eqs 21−24. However, when trying to use these equations, the big problem immediately arises because the VmXA varies with the kink sites related to the rate constant k−, and it brings about a different local concentration even in the same edge direction. In the present work, the average local concentration XA,avg was used instead of XA. The local concentration is originated from the energy required to detach a growth unit from a specific kink site. Therefore, the arithmetic mean value of the consumed energies for each kink site, ΔWkink,avg, was put into eq 17. 4.2. The {020} and {002} Faces. In the {020} and {002} faces, one type of F-face with four PBCs (Figure 3) exists. For the {020} faces, to determine the edge of a spiral motion, the average energy of bond chains for each PBC was calculated. First, the ⟨001⟩ PBCs with bond chain a and f (average energy per bond, −3.60 kcal/mol) was selected as the strongest PBC, and it can make a stable edge relative to the ⟨001⟩ PBCs consisting of bond chain b (energy, −1.67 kcal/mol). On the basis of the Hartman theory, two different bond chains cannot coexist to form an edge. Therefore, the ⟨102⟩ PBCs (average energy per bond, −2.64 kcal/mol) were not considered in the spiral growth due to the coexistance of bond chain a and f. For this reason, it is reasonable to assume that four-sided spiral motion consisting of the ⟨100⟩ (bond chain d) and ⟨001⟩ (bond chains a and f) PBCs occurs. The kink density was calculated by eqs 10 and 11, and the kink rate was calculated by eqs 21−24. For the {002} faces, the strongest PBC was initially selected as the ⟨010⟩ PBCs with bond chain e. In a same manner, the ⟨120⟩ PBCs along which bond chain c coexists are excluded. 4.3. The {210} Faces. Along the ⟨12̅0⟩ edge even in the same bonding network, there are two different types of PBCs, stable and unstable PBCs, as shown in Figure 4a. The energy changes at the kink sites along the unstable [120̅ ] edge (ΔW1 = 11.78 kcal/mol, ΔW3 = 8.75 kcal/mol, ΔW5 = 9.58 kcal/mol, ΔW7 = 8.75 kcal/mol) are relatively smaller than those at kink sites along the stable [12̅0] edge (ΔW2 = 12.92 kcal/mol, ΔW4 = 15.95 kcal/mol, ΔW6 = 15.12 kcal/mol, ΔW8 = 15.95 kcal/ mol). It stems from the bonding structure, in which the growth units at the lower level are connected more strongly to the bottom layer with many bonds (Figure 4b). Hence, the consumed energies for the growth units located along the stable edge are much bigger than those along the unstable edge. In this situation, the growth units at the unstable edge with low existence probabilities (P1 = 0.08790, P3 = 0.00058, P5 = 0.00222, P7 = 0.00058) are easily dissolved into solution than at the stable edge with high existence probabilities (P2 = 0.23587, P4 = 0.22228, P6 = 0.23223, P8 = 0.21834). Consequently, in the process of spiral growth, the only stable edge can be maintained. To deal with the unstable edge, formation of double kinks originally proposed by Kuvadia and Doherty20 was employed. In Figure 4c, once a growth unit is incorporated into the kink site 1, two different types of the kink sites 2 and 3 are consequently generated. Since the existence probability at the kink site 2 is larger than at the kink site 3, the growth units are incorporated in the kink site 2 maintaining a stable edge. The kink rate with eight different types of kink sites was calculated by eqs 25−27. Since the ⟨001⟩ edge has two different kink rates

Figure 5. The bonding networks for the {111} face of RDX.

PBCs with a large bond energy averaged over the bond chains were determined as the ⟨213̅⟩, ⟨011̅⟩, and ⟨2̅31̅⟩ PBCs in a series. However, bond chain a and e of the ⟨011⟩̅ PBCs coexist in the ⟨213̅⟩ PBCs. Since two chains cannot coexist to form one spiral edge, the ⟨213̅⟩ and ⟨2̅31̅⟩ PBCs were taken into consideration. For the alternative bond chains, the kink density was calculated by eq 11 with averaging over bond chain a and e. There are four different kinds of kink sites along the edge directions. On the {111} faces, it was found that the kink rate varies dependent on not only edge direction but also perpendicular direction even in the same PBC. For instance, the [213̅] edge gives the different kink rates dependent on perpendicular directions, [2̅31̅] and [23̅1]. For this reason, due to four different step velocities, the spiral layer spread asymmetrically.

5. RESULTS AND DISCUSSION In the present work, cooling crystallization of RDX from acetone and γ-butyrolactone was carried out as follows: (1) A RDX (12.4 g)/acetone (80.0 g) solution at 50 °C was cooled down to 20 °C at a cooling rate of 0.4 °C/min with stirring. (2) A RDX (44.0 g)/γ-butyrolactone (80.0 g) solution at 100 °C was cooled down to 40 °C at a cooling rate of 0.1 °C/min with stirring. From the RDX/acetone solution, the RDX crystals have a slightly elongated prismatic morphology with the {210} and {111} faces dominantly appearing (Figure 6a). On the other hand, from the γ-butyrolactone solution, RDX crystals are

Figure 6. RDX crystals grown from (a) acetone and (b) γbutyrolactone. 1807



Article

0.04 and 0.10. Typically, the morphology change dependent on supersaturation can be observed when 2-D nucleation occurs on crystal faces. As supersaturation increases, the critical formation energy of nuclei, which is the required energy to form a stepped interface on crystal faces, decreases. Accordingly, a majority of nuclei are prone to be twodimensionally formed, where the growth rate of a face is proportionally correlated with supersaturation. However, at low supersaturation, in which 2-D nucleation is extremely slow, screw dislocations on crystal faces become another source of steps, where the solute molecules do not require the additional energy associated with creating a kink or step. Therefore, in contrast to the 2-D nucleation mechanism, the dramatic morphology change at low supersaturation did not occur, whereas the effect of temperature on the growth rate of RDX was found to be negligibly small. Theoretically, the spiral growth model produces reliable results only if crystals are grown at low supersaturation. Once nucleation occurs from a solution, supersaturation rapidly decreases and then crystals are grown at low supersaturation. Therefore, it is reasonable to assume that crystals are grown at low supersaturation when a solution is cooled down with a very slow cooling rate. Although the supersaturation varies with respect to time, the spiral growth model can be used in this work since the crystals are grown for a relatively long time at low supersaturation and under the high kink energy regime that the largest kink energy is larger than the value of kbT. Interestingly, when considering the dispersive and polar bonding terms of solvent, a distinctive change on crystal morphology of RDX was observed. The reason behind the morphological change comes from the local concentration. There exists the effective concentration of growth units efficiently participating in dynamic equilibrium with surface located molecules, which strongly affects the probability of growth units available for incorporation into the kink sites.34 Strictly speaking, the local concentration is a major anisotropic term determining the relative growth rates of crystal faces. However, in the present work, the anisotropic kink rate containing the probability of growth units at kink sites was employed, by which the effect of anisotropic local concentration on crystal growth was compensated for. For this reason, each average local concentration for both solvents leads to different kink rates and consequently influences the growth rate of crystal faces, where different morphologies of RDX are observed. This approach enables us to calculate the kink density, kink rate, and existence probability, which are very important terms of theoretical description of spiral motion. According to the spiral growth model, it was found that the {210} and {111} faces grow very slowly compared to the other faces. The edges with a noticeably low kink rate are found to have kink sites whose existence probability is almost zero. The growth unit that is loosely bound to the kink site is easily solvated into solution. Hence, the process of incorporation of a growth unit into a kink site is consequently delayed. This result is observed from the only crystals consisting of non-centrosymmetric growth units which are energetically different growth units. For the {210} and {111} faces, the kink sites at which existence probabilities are almost zero were found to be in common, which leads to a slower rate of spiral motion on crystal faces.

shown to be hexagonal bipyramidal shape in which the {210} facial area decreases and the {111} faces are dominantly generated (Figure 6b). Figure 7a presents the morphology of RDX predicted by the BFDH model, where the area of the {210} faces was found to

Figure 7. Predicted crystal morphologies of RDX by the (a) BFDH and (b) AE model.

be much smaller than that of the other faces, and the {200} and {111} faces were dominantly developed. The reason behind the appearance of extensive {200} and {111} faces comes from a larger interplanar distance (Table 3). On the other hand, when Table 3. Interplanar Distance and Attachment Energy of RDX face

d (Å)

Eatt (kcal/mol)

{200} {020} {002} {210} {111}

6.59 5.78 5.35 5.73 6.75

−10.90 −10.23 −10.54 −13.30 −11.26

the AE model was employed, the {210} faces were found to grow much faster than the other faces, and their facial area consequently decreased due to the largest attachment energy of the {210} faces (Figure 7b). However, when the spiral growth model was employed, the predicted morphologies of RDX for each solvent were found to be in a good agreement with experimental data as shown in Figure 8. The {210} and {111} faces are dominantly generated from both solvents, whereas the {200} face slightly appears only from γ-butyrolactone. Table 4 summarizes the relative growth rates of RDX calculated by different conditions, solvent, supersaturation, and temperature. One can see that supersaturation has little effect on the morphology of RDX for S =

Figure 8. Predicted morphology of RDX grown from (a) acetone and (b) γ-butyrolactone using spiral growth model at S = 0.04 and T = 300 K. 1808



Article

Table 4. Relative Growth Rates of RDXa S = 0.04 face

Rvac

Ract

Rgbl

Rvac

Ract

Rgbl

300 K

{200} {020} {002} {210} {111} {200} {020} {002} {210} {111}

0.10 >100 35.39 1.00 0.01 0.14 >100 22.74 1.00 0.02

1.16 4.80 4.74 1.00 1.74 1.27 4.13 4.00 1.00 1.67

1.12 86.42 >100 1.00 0.97 1.15 42.78 72.68 1.00 1.05

0.10 >100 35.33 1.00 0.01 0.14 >100 22.68 1.00 0.02

1.16 4.76 4.72 1.00 1.74 1.27 4.11 4.22 1.00 1.67

1.12 85.48 >100 1.00 0.97 1.15 42.26 71.54 1.00 1.05

350 K

a

S = 0.10

temperature

The subscripts vac, act, and gbl denote vacuum, acetone, and γ-butyrolactone, respectively. (2) Doherty, R. M.; Duncan, S. W. Propellants, Explos., Pyrotech. 2008, 33, 4−13. (3) Song, X.; Li, F.; Wang, Y.; An, C.; Wang, J.; Zhang, J. J. Energ. Mater. 2012, 30, 1−29. (4) Van der Heijden, A. E. D. M.; Bouma, R. H. B. J. Cryst. Growth 2004, 4, 999−1007. (5) Kim, J.-W.; Park, D.-B.; Shim, H.-M.; Kim, H.-S.; Koo, K.-K. Ind. Eng. Chem. Res. 2012, 51, 3758−3765. (6) Kim, J.-W.; Shin, M.-S.; Kim, J.-K.; Kim, H.-S.; Koo, K.-K. Ind. Eng. Chem. Res. 2011, 50, 12186−12193. (7) ter Horst, J. H.; Geertman, R. M.; van der Heijden, A. E.; van Rosmalen, G. M. J. Cryst. Growth 1999, 198/199, 773−779. (8) ter Horst, J. H.; Geertman, R. M.; van Rosmalen, G. M. J. Cryst. Growth 2001, 230, 277−284. (9) Zhang, C.; Ji, C.; Li, H.; Zhou, Y.; Xu, J.; Xu, R.; Li, J.; Luo, Y. Cryst. Growth Des. 2013, 13, 282−290. (10) Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145−156. (11) Donnay, J. D. H.; Harker, D. Am. Mineral. 1937, 22, 446−467. (12) Pina, C. M.; Becker, U.; Risthaus, P.; Bosbach, D.; Putnis, A. Nature 1998, 395, 483−486. (13) Teng, H. H.; Dove, P. M.; Orme, C. A.; De Yoreo, J. J. Science 1998, 282, 724−727. (14) Rimer, J. D.; An, Z.; Zhu, Z.; Lee, M. H.; Goldfarb, D. S.; Wesson, J. A.; Ward, M. D. Science 2010, 330, 337−341. (15) Burton, W. K.; Cabrera, N.; Frank, F. C. Phil. Trans. R. Soc. A 1951, 243, 299−358. (16) Hollander, F. F. A.; Plomp, M.; van de Streek, J.; van Enckevort, W. J. P. Surf. Sci. 2001, 471, 101−113. (17) Chernov, A. A.; Petrova, E. V.; Rashkovich, L. N. J. Cryst. Growth 2006, 289, 245−254. (18) Cuppen, H. M.; Meekes, H.; van Enckevort, W. J. P.; Vlieg, E. Surf. Sci. 2004, 571, 41−62. (19) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. A 2009, 465, 1145− 1171. (20) Kuvadia, Z. B.; Doherty, M. F. Cryst. Growth Des. 2011, 11, 2780−2802. (21) Modern Crystallography III. Crystal Growth; Chernov, A. A.; Springer-Verlag: Berlin, Germany, 1984. (22) Markov, I. V. Crystal Growth for Beginners; World Scientific: Singapore, 1995. (23) Girifalco, L. A.; Good, R. J. J. Phys. Chem. 1957, 61, 904−909. (24) Belmares, M.; Blanco, M.; Goddard, W. A.; Ross, R. B.; Caldwell, G.; Chou, S.-H.; Pham, J.; Olofson, P. M.; Thomas, C. J. Comput. Chem. 2004, 25, 1814−1826. (25) Winn, D.; Doherty, M. F. AIChE J. 2000, 46, 1348−1367. (26) Lovette, M. A.; Doherty, M. F. Phys. Rev. E 2012, 85, 021604. (27) Schwinn, T.; Gaub, H. E.; Rabe, J. P. Supramol. Sci. 1994, 1, 85− 90. (28) Curry, J. E.; Cushman, J. H. J. Chem. Phys. 1995, 103, 2132− 2139.

6. CONCLUSIONS The morphology prediction of RDX crystal was performed employing spiral growth model describing spiral motions of a monolayer. The present work focused on understanding the reason why the {210} and {111} faces of RDX crystals are dominantly generated from acetone and γ-butyrolactone. The growth rates of crystal faces are largely influenced by both kink density and kink rate, but the kink rate was found to be a key factor in this work. The edges containing the kink sites with low existence probability appears to grow slowly because the loosely bound growth unit at kink sites hinders the integration of a step monolayer. This phenomenon is originated from the noncentrosymmetric growth unit of RDX, by which the energetically different kink sites are generated. Because of the kink sites with low existence probability, the {210} and {111} faces become morphologically important faces. Furthermore, the local concentration term contained in the kink rate was found to be an important factor describing the effect of solvent on crystal growth of RDX. The resulting difference in local concentrations for both solvents yields a different kink rate and consequently alters the crystal morphology of RDX.

■

ASSOCIATED CONTENT

S Supporting Information *

Detailed values used for estimating crystal morphology of RDX are tabulated. This material is available free of charge via the Internet at http://pubs.acs.org.

■

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 Agency for Defense Development (UC120019GD; Republic of Korea) and by Next-Generation Converged Energy Materials Research Center.

■

REFERENCES

(1) Van der Steen, A. C.; Verbeek, H. J.; Meulenbrugge, J. J. Proceedings of the 9th International Detonation Symposium on Detonation; August 28−September 1, 1989, Portland, Oregon, USA; Office of the Chief of Naval Research: Arlington, VA, USA; pp 83−88. 1809



Article

(29) Zhang, J.; Nancollas, G. H. J. Colloid Interface Sci. 1998, 200, 131−145. (30) Choi, C. S. Acta Crystallogr. 1972, B28, 2857−2862. (31) Stewart, J. P. J. Mol. Model. 2007, 13, 1173−1213. (32) Hazards of Chemical Rockets and Propellants; Chemical Propulsion Information Agency: Columbia, MD, 1984; Vol. 2, pp 11−13. (33) Zhu, W.; Xiao, J.; Zhu, W.; Xiao, H. J. Hazard. Mater 2009, 164, 1082−1088. (34) Gnanasambandam, S.; Rajagopalan, R. CrystEngComm 2010, 12, 1740−1749.

■

NOTE ADDED AFTER ASAP PUBLICATION This paper was published on March 25, 2014, with an error to equation 14. The corrected version was reposted with the issue on April, 2, 2014.

1810


Crystal Morphology Prediction of Hexahydro-1,3,5 ... - ACS Publications

Recommend Documents