Prediction of the Crystal Morphology of Î²-HMX using a Generalized

Subscriber access provided by Kaohsiung Medical University

Article

Prediction of the Crystal Morphology of #-HMX using a Generalized Interfacial Structure Analysis Model Bumjoon Seo, Seulwoo Kim, Minhwan Lee, Taewan Kim, Hyoun Soo Kim, Won Bo Lee, and Youn-Woo Lee Cryst. Growth Des., Just Accepted Manuscript • DOI: 10.1021/acs.cgd.7b01764 • Publication Date (Web): 08 Feb 2018 Downloaded from http://pubs.acs.org on February 28, 2018

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 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 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.

Crystal Growth & Design 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 36 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

Crystal Growth & Design

Prediction of the Crystal Morphology of β-HMX using a Generalized Interfacial Structure Analysis Model Bumjoon Seo†,#, Seulwoo Kim†,#, Minhwan Lee†, Taewan Kim†, Hyoun-Soo Kim‡, Won Bo Lee*,†, and Youn-Woo Lee*,† †

School of Chemical and Biological Engineering, and Institute of Chemical Processes, Seoul

National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826, Republic of Korea ‡

Agency for Defense Development, Yuseong, Daejeon, 34186, Republic of Korea

ABSTRACT

At sufficiently low supersaturations such that the spiral growth mechanism dominates, βcyclotetramethylenetetranitramine (HMX) grows from acetone into a polyhedron surrounded mainly by the (020) and (011) faces. In order to elucidate the morphology, a generalized form of the interfacial structure analysis model is suggested. In this method, the molecular order parameters of crystals are defined to identify the orientation and conformation of the adsorbed growth unit at the interface. This presents a robust method to calculate the orientational and

ACS Paragon Plus Environment

1

Crystal Growth & Design 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 36

conformational free energy surfaces that are utilized for the spiral growth model of centrosymmetric growth units with polygonal spiral edges. From the metadynamics simulation using these order parameters as collective variables, the free energy surfaces with respect to the collective variables revealed that high conformational free energy of the chair conformation discouraged pre-ordering of the growth units into crystal-like orientation and conformation. The resulting morphology was consistent with the previous experimental and theoretical results, indicating that the anisotropic local concentrations of the growth units at the interface play a critical role in the different relative growth rates of the slow-growing faces.


2

Page 3 of 36 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


1. INTRODUCTION Cyclotetramethylenetetranitramine (C4H8N8O8), usually called

HMX (High

Melting

eXplosive), belongs to a group of high explosives called nitramines which include RDX (C3H6N6O6) and HNIW (C6H6N12O12). Compared to the other widely used explosives such as TNT (C7H5N3O6), these explosives have higher densities which provide greater detonation properties but they also have higher sensitivities toward external stimuli resulting in higher likelihood of accidental ignitions. These properties vary greatly with the HMX polymorphs, and hence, the β form is highly favored owing to its high density and low friction/impact/shock sensitivities2. As for the β-HMX crystal, the sensitivities also depend largely on the crystallographic directions3 and consequently the morphology4, which encourage a thorough understanding of its growth morphology. As early attempts, modified attachment models5,6 were utilized to account for the interactions between the crystal face and solvents for predicting the morphology. The drawbacks of these models are their lack of considerations of the growth kinetics such as the growth mechanism and driving force. More recently, remarkable agreements with the experimental morphology were achieved by applying mechanistic growth models.7,8 These models explored the polyhedron morphology at low supersaturations and elongation in the [100] direction at large supersaturations by including the transition of the dominant growth mechanism. However, there are anticipated issues yet to be addressed, such as the anisotropic local concentration of the growth units at the interface9,10 and the orientational/conformational free energy barrier for the ordering (or reorientation) of the growth units.11,12 HMX provides a significant challenge in this respect owing to its molecular structure. Apart from the complexity of defining the relative orientations for different conformations, the sampling of the


3


Page 4 of 36

conformations poses difficulties even without considering the external solvents. Therefore, in previous works, the assumption of rigid molecule was implied or conformations were not considered at all for the prediction of morphology. Smith and Bharadwaj13 addressed the issue of different conformations by considering the flexibility of the HMX molecule for their proposed forcefield, but this was only utilized in solid or melt state calculations. As the HMX polymorphs are also comprised of different conformations (crown conformations for α, γ, and δ form and chair conformation for β form2), the issue of conformation is closely related to the overall growth of the HMX crystals. Metadynamics can be used to efficiently calculate the free energy surface with regard to the orientations and conformations of the adsorbed solute molecule at the interface.14-16 By adding a history-dependent bias to the space of chosen degrees of freedom, or the collective variables (CVs), metadynamics discourage the system from visiting the states that have been already visited. This property enables the algorithm to offer efficient exploration of the CV space and the calculation of the free energy surface. The key issue in metadynamics is that a small number, namely less than 3, of appropriate CVs that can distinguish all different states of interest should be defined. Based on the order parameter method for molecular crystals17-21, we suggest appropriate CVs for this purpose. In the present work, we develop a generalized interfacial structure analysis model to address the issues of anisotropic local concentrations at the interface and free energy barriers for the reorientation of the growth units for crystals comprised of centrosymmetric growth units. By using metadynamics with molecular order parameters as the CVs, it is possible to obtain the free energy surfaces with regard to the orientations and conformations of the adsorbed growth unit at


4

Page 5 of 36 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 interface. Moreover, the comparison of the morphology from seemingly different approaches shows the significance of the anisotropic local concentration at the interface.

2. THEORETICAL BACKGROUND 2.1 Spiral Growth Model In contrast to non-mechanistic models such as the Bravais-Friedel-Donnay-Harker model22 and the attachment energy model by Hartman and Perdok23, mechanistic models consider the external environments including supersaturation, growth mechanism, and solvents. Based on the periodic bond chain (PBC) theory by Hartman-Perdok approach23, crystal faces are classified into three categories: flat (F)-face with more than two PBCs, stepped (S)-face with one PBCs, and kinked (K)-face with no PBCs parallel to the faces. Among them, the F-faces grow relatively slowly by either the spiral growth model by Burton-Cabrera-Frank (BCF)24 or the twodimensional nucleation model25,26 resulting in crystals surrounded by F-faces. The dominant mechanism between these two layered growth models is dependent on the supersaturation. At low supersaturations, crystals are grown in the spiral growth dominant regime whereas as the supersaturation increases, a crossover to the two-dimensional nucleation growth dominant regime occurs. In both of these layered growth mechanisms, the kink sites are of utmost importance among all the incorporation sites. The reason is that kinks are the only sites for the attachment and detachment of the growth units, that (1) maintain the surface free energy of the edge as constant and (2) create a new kink site upon incorporation12. For solution growth of small organic molecules with layered growth mechanisms, ignoring the surface diffusion and assuming the


5


Page 6 of 36

kink incorporation as the rate determining step27 whose barrier corresponds to that of the desolvation of the solute molecule28 have become a consensus in mechanistic models.29-40. This desolvation barrier is usually assumed isotropic for all the faces. In the approaches referred to as the interfacial structure analysis model29-30,38-40 two additional concepts were addressed: (1) the local concentration of the growth units at the interface and (2) the barrier related to the reorientation of the growth units at the interface, introduced as “the orientational and conformational entropy barrier”.41 Before desolvation for incorporation into a kink, the adsorbed growth units are assumed to go through an reorientation process with a free energy barrier between its current orientational/conformational state (called an F2-unit state) and the crystal-like orientational/conformational state (called an F1-unit state).29-30 A generalized interfacial structure analysis model for the Kossel-like crystal systems of centrosymmetric growth units,33,42 based on the kink rate expressions by Shim and Koo37 is presented here. In the spiral growth model, the growth rate is expressed as

=

ℎ

(1)

where the step height ℎ is typically equal to the interplanar distance (or a fractional

multiple of depending on the identity of the growth unit) and the rotation time is the time required for one full turn of the spiral. The rotation time for N-sided spiral is11,31

=

, sin ,

(2)


6

Page 7 of 36 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


where , is the critical length of edge i, , is the angle between edge i and i+1, and the step velocity is the rate of lateral advancement of steps. The critical length , is obtained from

thermodynamic approximations as32

, =

! 2,

% − 1 "# $

(3)

where the supersaturation ratio is equal to &' /&',) , , is the distance between growth units along edge i, and the kink energy

!

is the work or free energy change required to form a

̅ kink on edge i. The average critical length on the face , is obtained by using the average

! values over the face for the edge length +, and kink energy + .

The average step velocity is proportional to the kink density and the concentration of the effective growth units at the interface as suggested in the Supporting Information, ! ̅ ∝ +-, &'() 1 00

(4)

where +-, is the average step propagation length on the face, &'() is the concentration of the 00

! effective growth units at the interface, and 1 is the kink density. As expected, the

dependence of step velocity on the thermodynamic driving force ( − 1) is also isotropic for the

̅ spiral growth model.11,12,35 Parabolic law24 is obtained from the scaling of 1/ , and ̅ with ( − 1).

The concentration of the effective growth units is related to the equilibrium value of the

concentration at the interface and the average effectiveness factor 2 , which is determined from

∗ two free energy barriers: the free energy barriers for desolvation Δ4567 and reorientation ∗ Δ6:4: .


7


&'() = 2 &'() 00

Page 8 of 36

∗ Δ4567 @ "# $ ≈ &'() ∗ ∗ Δ6:4:() Δ4567 exp ? @ + exp % "# $ "# $

exp ?

(5)

In recent mechanistic models, the effect of reorientation is expected to be much smaller than the ∗ ∗ effect of desolvation such that exp(Δ4567 /"# $) ≫ expΔ6:4:() /"# $ , resulting in

&'() = &'() . 00

! The kink density 1 is the probability of encountering a kink site on an edge and can be

approximated for centrosymmetric growth units as43 ! 1

! + 1 = 1 + exp %% 2 "# $

C

(6)

Subsequently, the relative growth rate can be calculated by collecting the anisotropic factors of

.

D ∝

+-,

+, E sin +

! + 1 2 & 1 + exp %% ! '() + 2 "# $

C

(7)

where E is the number of nearest neighbors in the two-dimensional lattice and + is the

angle between the regular spiral edges. As for the values of +-, and +, , the previously

used

assumptions

were +, ≈ +-, ≈ = (F )/G where H is

the

molecular

volume.29,30,38-40 In the present work, 1 and √3/2 were assigned as the values of +-, /+, for the cases of square and regular hexagon, respectively to account for the different polygonal spirals. Cancelling of the terms results in the final expression of


8

Page 9 of 36 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


D

! + 1 ∝ 2 & 1 + exp %% ! '() 2 "# $ E +

C

(8)

The average kink energy is obtained from the experimental dissolution enthalpy and other habit∗ controlling factors through the surface scaling factor K() and molecular orientational factor

L ,30

∗ K()

=

ln &'() 00

ln &'

L = − ln

4 55 ΔN = ΔN 4 55

&'() 00

&'()

∗ ! + = O K() ΔN 4 55 /E

(9)

(10)

(11)

5 The surface anisotropy factor O = P /PQRR and E are calculated from the PBC

∗ analysis. K() is calculated from the values of &'() and 2 obtained from molecular

dynamics simulations of the interface.

2.2 Molecular Order Parameters In a metadynamics simulation, coarse-grained dynamics is driven by the forces from the free energy of the system and the sum of history-dependent bias potentials in the space of a set of CVs.14 The CVs (1) should be explicit functions of the atomic positions and should also have continuous derivatives with respect to them, (2) should be able to distinguish the different states of the system by representing every slow evolving degree of freedom of the system. For the


9


Page 10 of 36

purpose of describing a crystallization process, traditionally used CVs include coordination numbers and Steinhardt parameters17 which are applicable to crystal systems with highly symmetric molecules or spherical geometry. For non-spherical cases e.g., most organic molecular crystals, however, the order parameters suggested by Santiso and Trout18 have been considered appropriate which have been modified and extended for various nucleation simulations.19-21 These parameters are based on the concept of a generalization of the pair distribution function. In addition to the center of mass distances between two molecules, the relative orientations between them and internal degrees of freedom of the molecule are considered. For a crystal at 0 K, the peaks of the parameter are completely localized with Y

S(T, U, V) = 2(T − TW )2(U − UX )2(V − VX ) X

(12)

where T and U are the center of mass distances and relative orientation between a molecule and its neighbor and V denotes the internal degrees of freedom such as conformation. For a crystal at finite temperatures, the delta functions are replaced by the probability densities as Y

Y

X

X

S(T, U, V) = ZX (T, U, V) ≈ ZXT (T)ZX) (U)ZX (V) [

(13)

where Z denotes the associated probability density. In recent applications of the extended order

parameters19-21, the order parameters discriminate between the disordered state at S = 0 and the

ordered state at S = 1. In this work, two order parameters were defined to investigate the free energy surface for the orientation and conformation of an HMX molecule. The molecular structure of HMX consists of an eight-membered ring whose conformation can be described equivalently by the Cartesian coordinates of the ring atoms, 8 dihedral angles,


10

Page 11 of 36 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


Strauss-Pickett coordinates44, and Cremer–Pople puckering coordinates45. Using the Cremer– Pople puckering coordinates, five puckering coordinates (three amplitudes and two phase angles) can be defined from the ring atom positions for an eight-membered ring. Equivalently, five mutually orthogonal Cartesian forms can be defined.46 For a vector of Cartesian form Cremer– Pople puckering coordinates CP are given as,

\] = (K^ , K^_ , K^G , K^` , K^a ) = (U_ bcd

_ , U_ deE _ , UG bcd G , UG deE G , U` )

(14)

The distribution of the puckering coordinates is approximated as a multiplication of five independent Gaussian terms [ ZX (\])

a

= f

1

! g2hσ_ jkl

mno p−

(K^! − ++++ CP! )_ 2σ_jkl

s

(15)

+++! and σjk are the mean and standard deviation of the corresponding probability where +CP l distribution of the Cremer–Pople puckering coordinates. There are eight equivalent peaks for the chair conformation.47 For the relative orientation between two molecules, molecular frames are defined on the ring nitrogen atoms that are connected to the nitro groups. The relative orientation between two molecular frames qi,j and qk is calculated as18,48,49

∗ U: ( ,t,) = U ,t U = (U:,u ( ,t,) , U:,v ( ,t,) , U:,w ( ,t,) , U:,x ( ,t,) )

(16)

where i = 1, 2 is the index of asymmetric units, j = 1, 2 is the index of the ring nitrogen atom with axial NO2 groups of a crystal molecule, and k = 1, 2, 3, and 4 is the index of the ring


11


Page 12 of 36

nitrogen atoms with both axial and equatorial NO2 groups of a solute molecule. The distribution of the relative quaternions can be approximated as a bipolar Watson distribution of the form18,

ZX) (U) ≈

1 _ exp zOX UX ∙ U: ( ,t,) | y (1/2, 2, OX ) 1 _ = exp}OX U:,u ( ,t,) ~ y (1/2, 2, OX )

(17)

where y (1/2, 2, OX ) is the hypergeometric function, qα is the mean relative orientation of the same asymmetric units of a crystal at 0 K and is equal to the quaternion multiplication identity (1, 0, 0, 0), OX is the concentration parameter of a crystal at the temperature of interest, and the

symbol ∙ denotes the standard 4D dot product. From the above two probability densities, the two order parameters are defined as

++++!,X K^! − CP = f exp − 2σ_jkl,

a

_

(18)

X !

_

_

exp OX U:,u ( t ,) exp OX U:,u ( t ,) _ = exp(OX ) exp(OX ) _

( , )

(19)

where i is the index of different asymmetric units, j is the index of the axial NO2 groups of a

crystal molecule, and (" , "_ ) is the pair of indices of NO2 groups at the opposite side of the ring

in the solute molecule.


12

Page 13 of 36 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


Table 1. Parameters used for the calculation of the molecular order parameters. Parameters Values

++++ ,W , ++++ ,W , ++++ ,W [Å]a

0.00000

W [Å]b

(0.078518, 0.074183, 0.040945, 0.050535, 0.046278)

,W ,++++ ,W ) [Å]a (++++ W b

(±0.47800, ±0.76736) or (±0.76736, ±0.47800)

211.58

a

from 0 K MD simulation of a crystal

b

from 293.15 K MD simulation of a crystal

3. MATERIALS AND METHODS 3.1 Experimental Methods β-HMX powder was provided by Hanwha Co. (Republic of Korea). Acetone (99.7%; SigmaAldrich, USA) was used as the solvent. Carbon dioxide (99.99%; Hyoup-sin Gas Co., Republic of Korea) was used as the antisolvent. All materials were used without further purification. In the present work, Gas Anti-Solvent (GAS) process is used to induce low supersaturations in the HMX-acetone solution. Teipel and coworkers50 used this method to recrystallize β-HMX from acetone solution into 65 µm-sized particles with a truncated octahedron morphology that is similar to the morphologies from crystallization by evaporation.5-7 In this work, smaller sized particles with similar morphology were reproduced using a different experimental condition. Before crystallization, 20 mL of 2.17 wt% (90% of the saturation concentration at the operating temperature) HMX solution is placed inside the crystallizer. The antisolvent passes through a precooler whose temperature is kept at 259 K by a refrigerated bath circulator (model: RBC-20, JEIO TECH, Republic of Korea) before going into the pump to avoid any cavitation. Then, the antisolvent is introduced from the bottom of the crystallizer to build up the pressure using a high-pressure pump (model: HKS-12000, Hanyang Accuracy, Republic of Korea). A


13


Page 14 of 36

heating circulator (model: RW-1040G, Lab Companion, USA) is used to control the temperature of the preheater and crystallizer during the crystallization. The temperature, pressure, antisolvent addition rate, and agitation rate of the crystallizer are kept at 303 K, 100 bar, 200 ml/min, and 900 rpm for 30 minutes after the desired pressure is met by two back-pressure regulators (model: 26-1700, TESCOM, USA). Fresh antisolvent is supplied throughout the process. Finally, the pressure is lowered and the crystals are collected from the membrane filter (FLUOROPORETM membrane filters 0.45 µm FH, Millipore, Ireland). A schematic diagram of the apparatus is provided in the Supporting Information.

3.2 Computational Methods For an accurate description of the solution interactions, the generalized AMBER force field51,52 was used for the HMX and acetone molecules. It was reported that when the generalized AMBER force field was applied to the HMX crystal structure that is used in the present work, the simulation box symmetry and size (and other properties such as the density, isothermal compressibility, bulk modulus, and coefficients of thermal expansions) were preserved at ambient and high temperature-pressure conditions.53 The relaxed molecular structure of the crown and chair conformations of an HMX molecule were obtained using density functional theory (DFT) structural optimization performed at the B3LYP/6-31G(d,p) level starting from the conformations in the α-HMX and β-HMX crystals. The same procedure was performed on the acetone molecule. The partial atomic charges were computed using the restrained electrostatic potential (RESP) formalism54 on the basis of the electrostatic potential evaluated at the HF/631G* level of theory using R.E.D. tools.55 All DFT calculations were carried out using GAMESS-US.56,57


14

Page 15 of 36 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 initial surfaces of the (011), (020), (110), and (101) faces were cleaved from the unit cell using Materials Studio (version 7.0), in the appropriate orientations using the crystal structure of β-HMX reported by Deschamps et al.58 (CCDC 792930) instead of the structure reported earlier

by Choi and Boutin59 (CCDC 1225495). The faces (110), (101) , and (101+ ) of the crystal

structure by reported Deschamps et al. are equivalent to the faces (111+ ), (102+ ), and (100) in the crystal structure reported by Choi and Boutin, respectively.7 The unit cell parameters of the crystal structure are a = 6.5255 Å, b = 11.0369 Å, c = 7.364 Å and β = 102.67°. The β-HMX crystallizes in a monoclinic structure with the space group P21/n with Z = 2 per unit cell. The

values of O , E , and ΔN4 55 (= 2.9 kcal/mol) are adopted from the previous study by Shim

and Koo.7 The corresponding crystal graph and bonding structures are shown in Figure 1 and Figure 2, respectively.

Figure 1. Crystal graph of β-HMX.


15


Page 16 of 36

Figure 2. Bonding structures for the F-faces of β-HMX crystals: (a) (020), (b) (011), (c) (110), and (d) (101) faces.

All molecular dynamics simulations were performed using GROMACS 5.1.460-63 equipped with PLUMED 2.364 plugin. An acetone solution with the HMX mole fraction of 0.00724 at 313 K (56 HMX molecules and 7680 acetone molecules) was prepared using the genbox package of GROMACS. The system was cooled to 293.15 K, where the concentration corresponds to a supersaturation value of 0.613 defined as = &' /&' − 1.65 This is approximately 90% of the lowest onset supersaturation _ , the supersaturation at which the transition from the spiral

growth dominant regime to 2D nucleation growth dominant regime occurs34, among the reported _ values for the major faces.7

After insertion, energy minimization by steepest descent method and equilibration by NVT (isothermal-isochoric) simulation at 313 K and 1.013 bar for 1 ns were performed. The resulting solution was brought into contact with the specific faces, and subsequently, the entire system was


16

Page 17 of 36 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


equilibrated with a short 100 ps NVT simulation at 313 K followed by 10 ns NPT (isothermalisobaric) simulation at 313 K and 1 bar. After equilibration, production runs were carried out in NPT simulation for 200 ns at 293 K, 1 bar. For the temperature and pressure control, BussiDonadio-Parrinello thermostat66 and semiiisotropic Berendsen barostat67 followed by the semiisotriopic Parrinello-Rahman barostat68,69 were used. The coupling constants for the thermostat and barostat were 0.1 and 2 ps, respectively. All bonds were constrained by the LINCS algorithm with a timestep of 2 fs. The cutoff distances were 1.4 nm for both Coulombic interaction and Lennard-Jones interactions. The particle mesh Ewald (PME) method70 was used for the long-range electrostatic interactions. Three-dimensional periodic boundary conditions were applied. In order to investigate the free energy surfaces associated with the reorientation of an adsorbed HMX molecule in the acetone solution in contact with specific crystal faces, well-tempered metadynamics15 was used to achieve the convergence of the free energy estimates. NVT simulations were carried out with the reduced cut off distances of 1.2 nm (1.1 nm for 101 face) to match half the box size. Gaussian hills with initial height of 0.1 kJ/mol and width of 0.1 in each CV were deposited every 1 ps with the bias factor of 5. The resulting growth morphology was drawn using the WinXMorph software.71

4. RESULTS AND DISCUSSION From the molecular dynamics simulations, the concentration of the growth units at the interface is calculated from the mole fraction of the growth unit in the first fluid layer, which corresponds to the nearest layer parallel to each surface with its height equal to the interplanar


17


Page 18 of 36

distance . The density profiles along the perpendicular directions of the faces are shown in Figure 3.

Figure 3. Density profile along the z direction in the simulation for the (a) (020), (b) (011), (c) (110), and (d) (110) faces.

As the free energy surface along the two CVs is monotonic with a plateau, the corresponding barrier can be observed using the finite temperature string (FTS) method72,73. The contour of the two-dimensional free energy surface and the free energy along a possible minimum energy path S on the surface is shown in Figure 4. For all the faces, starting from a stable disordered state, an initial barrier was observed, followed by a plateau and further increase in the free energy at the most ordered state. The origin corresponds to the F2-unit whereas the plateau corresponds to the


18

Page 19 of 36 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


F1-unit. The deduced barrier was higher for the (011) face, which is known as the dominant face.74,75

Figure 4. Free energy surfaces from the two collective variables S1 and S2 and along a possible minimum energy path S according to the FTS method for the (a) (020), (b) (011), (c) (110), and (d) (101) faces.

Noticeable characteristics of the free energy surface are the high free energy at S1 = 1 and S2 = 1 and the plateau. These are related to the free energies of different conformations of the HMX molecule. As for the gas phase conformations from DFT calculations13,75, the chair conformation is considered as one of the low-energy conformers. The potential energies were determined to be in the order of boat-chair < chair < boat-boat < crown. Our structural optimization of the four low-energy conformers suggested by Molt et al.76 converged to crown and chair conformations. Since other initial structure could lead to different conformations, this suggests that chair is one


19


Page 20 of 36

of the low-energy conformers in gas phase. On the other hand, in the investigation of the conformational free energies for cyclooctane and its derivatives, the chair conformation is considered a transition state. Martin et al.77 calculated the MM3 energy of cyclooctane to determine the potential energy landscapes and transition pathways for the three conformation families: crown (CR: crown, CC: chair-chair, TCC: twist-chair-chair), boat (B: boat, BB: boatboat, TB: twist-boat), and boat-chair (BC: boat-chair, TBC: twist-boat-chair, C: chair, TC: twistchair). The BC conformation (and TBC) was observed as the global minimum whereas the CR conformation had slightly higher energy than BC. BB and TB in the boat family had relatively higher energy than the former ones. The B and C conformations were observed to be transition states. Spiwok and Kralova78 used metadynamics in the reduced dimension of path collective variables defined by Isomap projections to calculate the conformational free energy surface of trans,trans-1,2,4-trifluorocyclooctane. It was observed upon starting the simulation from BC conformation, that CR and B have higher free energies (>10 kJ/mol). Our simulation coincides with the latter findings in that the chair conformation is instantaneously converted to boat-chair and transition occurs between boat-chair/crown/boat-boat conformations. This is the reason for the high free energy values at S1 = 1 and S2 = 1 and a plateau for all the chair conformations at S1 > 0. Before predicting the morphology using the model, we first present two limiting cases with respect to the two free energy barriers. ∗ ∗ if Δ4567 ≫ Δ6:4: (δ = 1 and &'() = &'() ) 00

! ∗ + Δ4567 1 D ∝ exp ?− @ & 1 + exp %% ! '() "# $ 2 "# $ E +

C

(20)


20

Page 21 of 36 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


∗ ∗ ≪ Δ6:4: (δ ≪ 1 and &'() ≪ &'() ) if Δ4567 00

D ∝ exp −

∗ Δ6:4:()

"# $

! + 1 & + exp % 1 %% ! '() 2 "# $ Eℎ" +

C

(21)

where the desolvation free energy barrier term was the omitted isotropic part of the step velocity (See Supporting Information for the full equation). The first case assumes a relatively large desolvation barrier and ignores the effect of the reorientation barrier. As the desolvation corresponds to the breaking of the solvation shell on the solute, it is independent of the identity of the face, and therefore considered isotropic. The original interfacial structure analysis and further studies30,38-40,79 assumed that the value of this desolvation barrier is close to the thermal fluctuation, "# $, but this could be significantly different from the actual value. If the actual

value is much larger than the assumed value as suggested by Vekilov,12 where the expected values are approximately 28 kJ/mol for the enthalpy part and near zero for the entropy part80, the effect of reorientation barrier could be negligible. The second case is similar to the original

interfacial structure analysis except for the factor of E regarding the rotation time for different spiral geometries from Equation (2). For the case of β-HMX presented here, the Equation (20) is used for another reason. The free energy surfaces with respect to the order parameters suggest that the reorientation of the growth units at the interface is unfavorable because an ordered F1-unit is more unstable than an unordered F2-unit due to the significantly high conformational free energy of the chair conformation. Therefore, we expect that all the growth units at the interface will go through a similar process of kink incorporation starting from an F2-unit without any pre-ordering. This is similar to the previous study where 2 = 1 or, equivalently, L = 0 was assumed for


21


Page 22 of 36

paraffin.30 In this case, the concentration of the effective growth units &'() is the same as the 00

concentration of the growth units at the interface &'() . Therefore, only the effect of anisotropic local concentration of the growth units at the interface is evaluated whereas the effect of the reorientation barrier is ignored. It should be noted that since the desolvation barrier term is

considered isotropic, it is not considered for the calculation of the relative growth rates. ∗ From the above analysis, the internal ( , O , E ) and external L , K() habit-

controlling factors are obtained. The relative growth rates are calculated using Equation (19) with the habit-controlling factors. The calculated habit-controlling factors and relative growth rates are summarized in Table 2. The growth morphology of β-HMX from the Frank–Chernov condition28 is shown in Figure 5. The obtained morphology is in good agreement with the morphologies of crystals obtained from the GAS process in Figure 6. Compared to the results from other spiral growth models,7,8 the morphology is characterized by a relatively larger (020) face similar to the experimental morphology reported by Duan et al.5 However, it shares the same feature of definite (101) face which has low recurrences in the experimental morphologies.73,74 The effect of anisotropic concentration of the growth units at the interface, which is implicit in other methodologies, could have resulted in the similarity in the morphology. From the point of view of interfacial structure analysis models, more accurate determination and comparison of the kink incorporation barriers including those of the desolvation step and reorientation step are desired for further analysis of the growth process. In this work, only the free energy surface for reorientation related to the pre-ordering at the interface was considered which is different from that of reorientation during the actual incorporation into a kink.


22

Page 23 of 36 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


Table 2. Habit-controlling factors and relative growth rates for the F-faces of β-HMX crystals. Layer Miller index

a

(Å)

a

∗ \()

T

5.518

0.4438

4

0

0.6038

0.865

6.021

0.4940

6

0

0.5845

1.00

5.515

0.4048

6

0

0.5850

1.13

4.318

0.3682

4

0

0.5298

1.39

calculated from the slice and lattice energy values in Shim et al.7

Figure 5. Predicted β-HMX growth morphology by the spiral growth model.

Figure 6. β-HMX crystals grown by the GAS process.


23


Page 24 of 36

5. CONCLUSION In the present work, a generalized interfacial structure analysis model was developed to combine the reorientation effects of both orientation and conformations of the adsorbed growth units and adequate descriptions of the spiral rotation for the crystals of centrosymmetric growth units. The free energy surfaces with respect to the orientation and conformation of an adsorbed growth unit were constructed by defining the molecular order parameters to distinguish crystallike and solution-like states. The obtained free energy surfaces from metadynamics simulations revealed that the ordering of the growth units at the interfaces of the β-HMX crystal and acetone solution was unfavorable due to the high free energy of the chair conformations. The resulting morphology prediction from the computations of the anisotropic local concentration of the growth units at the interface exhibited consistency with the experimental morphology and with the morphologies from previous theoretical studies with different approaches.

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge on the ACS Publications website at Derivation of the step velocity and kink rates. Supplementary information on Cremer–Pople puckering coordinates and quaternions. Table of atomic charges of the HMX molecule. Figures of morphology prediction from BFDH model and attachment energy model. Figure of the concentration of the adsorbed HMX acid molecules at the crystal-solution interface. Figure of the GAS process apparatus.


24

Page 25 of 36 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


AUTHOR INFORMATION Corresponding Author * E-mail: [email protected]. * E-mail: [email protected]. Author Contributions #

These authors contributed equally

ORCID Youn-Woo Lee: 0000-0003-4793-6693 NOTES The authors declare no competing financial interest.

ACKNOWLEDGMENT The authors would like to thank Dr. Shim at the Agency for Defense Development (Republic of Korea) for the helpful discussions. This work was supported by Agency for Defense Development, National Institute of Supercomputing and Network, Korea Institute of Science and Technology Information (KSC-2015-C2-034), and the National Research Foundation Grant (2015R1A2A2A01007379).

REFERENCES (1) Bernstein, J.; Polymorphism in Molecular Crystals; Oxford University Press: Oxford, 2002.


25


Page 26 of 36

(2) Zhu, W.; Xiao, J.; Ji, G.; Zhao, F.; Xiao, H. First-Principles Study of the Four Polymorphs of Crystalline Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine. J. Phys. Chem. B 2007, 111, 12715-12722. (3) Zamiri, A. R.; De, S. Multiscale modeling of the anisotropic shock response of β-HMX molecular polycrystals. Interact. Multiscale Mech. 2011, 4, 139-153. (4) Song, X.; Wang, Y.; An, C.; Guo, X.; Li, F. Dependence of particle morphology and size on the mechanical sensitivity and thermal stability of octahydro-1,3,5,7-tetranitro-1,3,5,7tetrazocine. J. Hazard. Mater. 2008, 159, 222-229. (5) Duan, X; Wei, C.; Liu, Y.; Pei, C. A molecular dynamics simulation of solvent effects on the crystal morphology of HMX. J. Hazard. Mater. 2010, 174, 175-180. (6) Zhang, C.; Ji, C.; Li, H.; Zhou, Y.; Xu, J.; Xu, R.; Li, J.; Luo, Y. 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. (7) Shim, H.-M.; Koo, K.-K. Prediction of Growth Habit of β-Cyclotetramethylenetetranitramine Crystals by the First-Principles Models. Cryst. Growth Des. 2015, 15, 3983-3991. (8) Tilbury, C. J.; Doherty, M. F.; Modeling Layered Crystal Growth at Increasing Supersaturation by Connecting Growth Regimes. AIChE J. 2017, 63, 1338-1352. (9) Schwinn, T.; Gaub, H. E.; Rabe, J. P. Supramolecular structures and dynamics of organic adsorbate layers at the solid-liquid interface. Supramol. Sci. 1994, 1, 85-90.


26

Page 27 of 36 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


(10) Adobes-Vidal, M.; Shtukenberg, A. G.; Ward, M. D.; Unwin, P. R. Multiscale Visualization and Quantitative Analysis of L-Cystine Crystal Dissolution. Cryst. Growth Des. 2017, 17, 1766-1774. (11) Chernov, A. A. Notes on interface growth kinetics 50 years after Burton, Cabrera and Frank. J. Cryst. Growth 2004, 264, 499-518. (12) Vekilov, P. G. What determines the rate of growth of crystals from solution? Cryst. Growth Des. 2007, 7, 2796-2810. (13) Smith, G. D.; Bharadwaj, R. K. Quantum Chemistry Based Force Field for Simulations of HMX. J. Phys. Chem. B 1999, 103, 3570-3575. (14) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 49, 12562–12566. (15) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. (16) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the Equilibrium Boltzmann Distribution from Well-Tempered Metadynamics. J. Comput. Chem. 2009, 30, 1616-1621. (17) Steinhardt, P.; Nelson, D.; Ronchetti, M. Bond-orientational order in liquids and glasses. Phys. Rev. B 1983, 28, 784-805. (18) Santiso, E. E.; Trout, B. L. A general set of order parameters for molecular crystals. J. Chem. Phys. 2011, 134, 064109.


27


Page 28 of 36

(19) Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. Uncovering Molecular Details of Urea Crystal Growth in the Presence of Additives. J. Am. Chem. Soc. 2012, 134, 17221-17233. (20) Salvalaglio, M.; Vetter, T.; Mazzotti, M.; Parrinello, M. Controlling and Predicting Crystal Shapes: The Case of Urea, Angew. Chem. Int. Ed. 2013, 52, 13369-13372. (21) Giberti, F.; Salvalaglio, M.; Mazzotti, M.; Parrinello, M. Insight into the nucleation of urea crystals from the melt. Chemical Engineering Science 2015, 121, 51-59. (22) Donnay, J. D. H.; Harker, D. A new law of crystal morphology extending the law of Bravais. Am. Mineral. 1937, 22, 446-467 (23) Hartman, P.; Bennema, P. The attachment energy as a habit controlling factor: I theoretical considerations. J. Cryst. Growth 1980, 49, 145−156. (24) Burton, W. K.; Cabrera, N.; Frank, F. C. The growth of crystals and the equilibrium structure of their surfaces. Philos. Trans. R. Soc. A 1951, 243, 299−358. (25) Volmer, M.; Marder, M. Zur Theorie der linearen Kristallisationsgeschwindigkeit unterkühlte Schmelzen und unterkühlter fester Modifikationen. Z. Phys. Chem. A 1931, 154, 97112. (26) Kaischew, R.; Stranski, I. N. Zur Theorie der linearen Kristallisationsgeschwindigkeit, Z. Phys. Chem. A 1934, 170, 295-299. (27) Chernov, A. A.; Modern Crystallography III: Crystal Growth; Springer-Verlag: Berlin, 1984


28

Page 29 of 36 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


(28) Chernov, A. A. The Kinetics of the Growth Forms of Crystals. Sov. Phys. Cryst. 1963, 7, 728-730. (29) Liu, X. Y.; Boek, E. S.; Briels, W. J.; Bennema, P. Prediction of crystal growth morphology based on structural analysis of the solid-fluid interface. Nature 1995, 374, 342-345. (30) Liu, X. Y.; Bennema, P. Prediction of the growth morphology of crystals. J. Cryst. Growth 1996, 166, 117-123. (31) Snyder, R. C.; Doherty M. F. Predicting crystal growth by spiral motion. Proc. R. Soc. A 2009, 465, 1145-1171. (32) Lovette, M. A.; Doherty, M. F. Reinterpreting edge energies calculated from crystal growth experiments. J. Cryst. Growth 2011, 327, 117-126. (33) Kuvadia, Z. B.; Doherty, M. F. Spiral Growth Model for Faceted Crystals of NonCentrosymmetric Organic Molecules Grown from Solution. Cryst. Growth Des. 2011, 11, 27802802. (34) Lovette, M. A.; Doherty, M. F. Predictive Modeling of Supersaturation-Dependent Crystal Shapes. Cryst. Growth Des. 2012, 12, 656-669. (35) Kim, S. H.; Dandekar, P.; Lovette, M. A.; Doherty, M. F. Kink Rate Model for the General Case of Organic Molecular Crystals. Cryst. Growth Des. 2014, 14, 2460-2467. (36) Li, J.; Tilbury, C. J.; Joswiak, M. N.; Peters, B.; Doherty, M. F. Rate Expressions for kink Attachment and Detachment During Crystal Growth. Cryst. Growth Des. 2016, 16, 3313-3322.


29


Page 30 of 36

(37) Shim, H.-M.; Koo, K.-K. Molecular Approach to the Effect of Interfacial Energy on Growth Habit of ε‑HNIW. Cryst. Growth Des. 2016, 16, 6506-6513. (38) Gnanasambandam, S.; Rajagopalan, R. Growth morphology of α-glycine crystals in solution environments: an extended interface structure analysis. Cryst. Eng. Comm. 2010, 12, 1740-1749. (39) Jion, A. I.; Rajagopalan, R. On the study of crystal growth via interfacial analysis and string optimization. Cryst. Eng. Comm. 2014, 16, 6224-6233. (40) Seo,B.; Kim, T.; Kim, S.; Ryu, J. H.; Ryu, J.; Yoon, J.; Lee, W. B.; Lee, Y.-W. Interfacial Structure Analysis for the Morphology Prediction of Adipic Acid Crystals from Aqueous Solution. Cryst. Growth Des.2017, 17, 1088-1095. (41) Liu, X. Y. Interfacial Effect of Molecules on Nucleation Kinetics. J. Phys. Chem. B 2001, 105, 11550-11558. (42) Kossel, W. Zur Theorie des Kristallwachstums. Nachur. Ges. Gottingen 1927, 206, 135145. (43) Frenkel, J. On the surface motion of particles in crystals and the natural roughness of crystalline faces. J. Phys. U.S.S.R. 1945, 9, 392-398. (44) Strauss, H. L.; Pickett, H. M. Conformational structure, energy, and inversion rates of cyclohexane and some related oxanes. J. Am. Chem. Soc. 1970, 92, 7281-7290. (45) Cremer, D.; Pople, J. A. A General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97, 1354-1358.


30

Page 31 of 36 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


(46) Evans, D. G.; Boeyens, J. C. A. Mapping the Conformation of Eight-Membered Rings. Acta Cryst. 1988, B44, 663-671. (47) Kessler, M.; Perez, J. Equivalence properties of the Cremer & Pople puckering coordinates for N-membered rings. J. Math. Chem. 2012, 50, 187-209. (48) Karney, C. F. F. Quaternions in molecular modeling. J. Mol. Graph Model. 2007, 25, 595604. (49) Kuipers, J. B. Quaternions and Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality; Princeton University Press: Princeton, 1999. (50) Teipel, U.; Förter-Barth, U.; Krause, H. H. Crystallization of HMX-Particles by Using the Gas Anti-Solvent Process, Propellants Explos. Pyrotech. 1999, 24, 195-198. (51) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P.A. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. (52) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (53) Bergh, M.; Caleman, C. A Validation Study of the General Amber Force Field Applied to Energetic Molecular Crystals. J. Energ. Mater. 2016, 34, 62-75. (54) Bayly, C.; Cieplak, P.; Cornell, W.; Kollman, P. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269-10280.


31


Page 32 of 36

(55) Dupradeau F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. The R.E.D. tools: advances in RESP and ESP charge derivation and force field library building, Phys. Chem. Chem. Phys, 2010, 12, 7821-7839. (56) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery J. A. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 13471363. (57) Gordon, M. S.; Schmidt M.W. Advances in electronic structure theory: GAMESS a decade later. In Theory and Applications of Computational Chemistry: the first forty year; Dykstra C. E. ,Frenking, G., Kim, K. S., Scuseria, G. E., Eds; Elsevier: Amsterdam, 2005. (58) Deschamps, J. R.; Frisch, M.; Parrish, D. CCDC 792930: Experimental Crystal Structure Determination, 2014. (59) Choi, C. S.; Boutin, H. P. A study of the crystal structure of β-cyclotetramethylene tetranitramine by neutron diffraction. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1970, 26, 1235−1240. (60) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (61) Lindahl, E.; van der Spoel, D.; Hess, B. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7, 306−317. (62) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718.


32

Page 33 of 36 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


(63) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (64) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604−613. (65) Chen, L.-Z.; Zhang, J.; Wang, W.-Y.; Diao, Y. Solubility of β-HMX in Acetone + Water Mixed Solvent Systems at Temperatures from 293.15 K to 313.15 K. J. Solution Chem. 2012, 41, 1265-1270. (66) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (67) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (68) Parrinello, M.; Rahman, M. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (69) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (70) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (71) Kaminsky, W. WinXMorph: a computer program to draw crystal morphology, growth sectors and cross sections with export files in VRML V2.0 utf8-virtual reality format. J. Appl. Crystallogr. 2005, 38, 566−567.


33


Page 34 of 36

(72) Weinan, E.; Ren, W.; Vanden-Eijnden, E. Finite Temperature String Method for the Study of Rare Events. J. Phys.Chem. B 2005, 109, 6688-6693. (73) Vanden-Eijnden, E.; Venturoli, M. Revisiting the finite temperature string method for the calculation of reaction tubes and free energies. J. Chem. Phys. 2009, 130, 194103 (74) Gallagher, H. G.; Sherwood, J. N.; Vrcelj, R. M. Growth and Dislocation studies of βHMX. Chem. Cent. J. 2014, 8, 75. (75) Gallagher, H. G.; Sherwood, J. N.; Vrcelj, R. M. The growth and perfection of bcyclotetramethylene-tetranitramine (HMX) studied by laboratory and synchrotron X-ray topography. J. Cryst. Growth 2017, 475, 192-201. (76) Molt, R. W.; Watson, T.; Bazanté , A. P.; Bartlett, R. J. The great diversity of HMX conformers: probing the potential energy surface using CCSD(T). J. Phys. Chem. A 2013, 117, 3467−3474. (77) Martin, S.; Thompson, A.; Coutsias, E.; Watson, J.-P. Topology of cyclo-octane energy landscape. J. Chem. Phys. 2010, 132, 234115. (78) Spiwok, V.; Králová, B. Metadynamics on the conformational space nonlinearly dimensionally reduced by Isomap. J. Chem. Phys. 2011, 135, 224504. (79) Shim, H.-M.; Myerson, A. S.; Koo, K. K. Molecular Modeling on the Role of Local Concentration in the Crystallization of L‑Methionine from Aqueous Solution. Cryst. Growth Des. 2016, 16, 3454-3464.


34

Page 35 of 36 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


(80) De Yoreo, J. J. In 13th International Conference on Crystal Growth; Hibiya, T., Mullin, J. B., Uwaha, M., Eds.; Elsevier: Kyoto, Japan, 2001.


35


Page 36 of 36

For Table of Contents Use Only

Prediction of the Crystal Morphology of β-HMX using a Generalized Interfacial Structure Analysis Model Bumjoon Seo, Seulwoo Kim, Minhwan Lee, Taewan Kim, Hyoun-Soo Kim, Won Bo Lee*, and Youn-Woo Lee*

Table of Contents graphic and synopsis

The orientation and conformation of an HMX molecule at the interface can be identified using order parameters that distinguish crystal-like molecules. Using metadynamics with the order parameters revealed that all the adsorbed molecules act as equivalent growth units. The prediction of the morphology using the anisotropic concentration of these growth units with the spiral growth model exhibited consistency with the experimental morphology of β-HMX crystal.


36

Prediction of the Crystal Morphology of Î²-HMX using a Generalized

Recommend Documents