Prediction of Crystal Morphology of Cyclotrimethylene Trinitramine in

Nov 17, 2014 - Crystal morphology prediction of 1,3,3-trinitroazetidine in ethanol solvent by molecular dynamics simulation. Wenyan Shi , Yuting Chu ,...
0 downloads 5 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCA

Prediction of Crystal Morphology of Cyclotrimethylene Trinitramine in the Solvent Medium by Computer Simulation: A Case of Cyclohexanone Solvent Gang Chen,† Mingzhu Xia,*,† Wu Lei,† Fengyun Wang,† and Xuedong Gong‡ †

Institute of Industrial Chemistry and ‡Chemistry Department, Nanjing University of Science & Technology, Nanjing 210094, China S Supporting Information *

ABSTRACT: The crystal morphology of the energetic material cyclotrimethylene trinitramine (also known as RDX) influenced by the solvent effect was investigated via molecular dynamics simulation. The modified attachment energy (MAE) model was established by incorporating the growth parameter−solvent term. The adsorption interface models were used to study the adsorption interactions between solvent and RDX surfaces. The RDX crystal morphology grown from the cyclohexanone (CYC) solvent as a case investigation was calculated by the MAE model. The calculation results indicated that, due to the effect of CYC solvent, (210) and (111) faces had the greatest morphological importance on the final RDX crystal, while the morphological importance of (020), (002), and (200) faces were reduced. The predicted RDX morphology was in reasonable agreement with the observed experiment result.

1. INTRODUCTION Cyclotrimethylene trinitramine (common name RDX), one of the most famous single-compound explosives, has been widely applied in the military weapons and civilian fields since it was synthesized.1,2 It has found that the RDX from industrial preparation has high sensitivity and low load density.3,4 This seriously impacts the product performance of RDX, because the high sensitivity can cause unintended explosion and the low packing density results in poor detonation property.5,6 The research has shown that the sensitivity and detonation properties of RDX are largely influenced by the diverse crystal morphologies.7−9 The spherical RDX crystal can improve loading density and impact sensitivity to a great extent, compared with other crystal shapes such as needles and plates.1 Therefore, to reduce sensitivity and improve performance, the research on the crystal morphology of RDX has an important practical significance. The control of crystal morphology is a complex and difficult process, which depends on both the crystal internal structures and external growth conditions such as crystallization technologies,4,10 the solvents used,11,12 and the additives.13,14 Among them, it is believed that the solvent plays an important role in the crystal morphology. Compared to the traditional timeconsuming experimental methods, computational simulation has been a powerful and important tool for crystal workers to explore the crystallization mechanism at the atomistic/molecular level. Clydesdale et al.15,16 improved the attachment energy model considering the growth environment, and designed a computer program (HABIT95) for the prediction of crystal habit. Doherty et al.17,18 proposed a method taking into account the crystallization process to model the shape of solution-grown © 2014 American Chemical Society

organic crystals, which reasonably predicted the shape of adipic acid grown from water, etc. Lin et al.19 attempted to predict the crystal morphology in the aqueous solution. They successfully calculated the habits of glycine (α-form) crystal in water with the help of GenMol software. Stoica et al.20 investigated the solvent effect on the polar crystal habits by a molecular dynamics simulation method. They found that the habit modification was caused by special surface−solvent interactions, which affected the growth rate of the polar faces differently. Gale et al.21,22 studied the dissolution and growth of the molecular crystal of urea using the dynamic atomistic simulation approach, which offers insights into the role of the solvent, the supersaturation, and the contribution that extended defects (such as screw dislocations) make to crystal growth. In recent years, energetic materials scholars have focused on the investigation of crystalline morphologies of explosives by means of molecular simulation tools. Duan et al.23 discussed the solvent effect on the change of HMX crystal habits via molecular dynamics simulations. They modified the attachment energy model by introducing the energy correct term, and successfully predicted the morphology of HMX crystal in acetone solvent. Lee et al.24 attempted to explain the cosolvent effect on the shape evolution of ADNBF crystal by a molecular modeling technique. The crystallization morphology of ADNBF from cosolvents was predicted by the modified attachment energy model introducing the binding site densities of crystal surfaces. Zhang et al.25 proposed a so-called occupancy model for modeling the crystal morphologies of energetic Received: August 29, 2014 Revised: November 15, 2014 Published: November 17, 2014 11471

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

requires the removal of the adsorbed solvent molecules on the surface, that is, the desolvation of the crystal surface; thereby solute molecules can reach an incorporation (kink) site by surface diffusion. Such a desolvation process apparently costs energies, and the reduced energies of the system are the solvent interaction energies, which are determined by the adsorption strength between the solvent and the surface. From the analogy of the definition of the attachment energy, we introduce the sol‑att variable named the solvent attachment energy (Ehkl ) to represent the solvent effect, defined as the energy released by attaching a solvent slice of thickness dhkl on the (hkl) surface, schematically depicted in Figure 1. Due to the solvent effect, the

materials. This model considering the solvent and temperature factors was a help in understanding the effect of the external growth conditions on the crystal growth. Shi et al.26 explored the influence of trifluoroacetic acid on the growth habit of ANPyO crystal by molecular simulation. They found that the change of ANPyO crystal morphology was related to the diffusion of solvent molecules onto specific surfaces. For the crystal morphology of RDX, the work also has made some progress.27−30 The investigation was mainly focused on the influence of solvents on the crystal habits of RDX via Monte Carlo or molecular dynamics methods. They attempted to explain the solvent effect at the molecular level, and indicated the solvent-specific RDX surface adsorption interactions were an important factor for the solvent-induced change of RDX growth habits. In addition, Shim et al.31 developed the Burton− Cabrera−Frank (BCF) model to predict the RDX shapes grown from solvents. They thought that the kink rate of a growth unit was a critical factor for the growth of RDX surface. This model was a help in understanding the RDX growth mechanism in solution. Although previous work has achieved many important and meaningful results, the ultimate relationships between the solvent interactions and crystal morphology are still unclear. Therefore, it is very necessary to further study the solvent effect on the crystal morphology of RDX. The subject of our work is to find the quantitative relation of solvent effect and the crystal morphology of RDX at the molecular level and to attempt to establish a reliable prediction model for RDX crystal grown from the solvent. In this paper, cyclohexanone (CYC) solvent is chosen as the objective solvent, which is usually used as the crystallization solvent for the preparation of RDX.32 The RDX surfaces−CYC solvent adsorption interface models with molecular dynamics simulations are employed to study the adsorption behaviors of CYC molecules on the RDX habit faces. A modified attachment energy model is developed to discuss the solvent effect and predict the RDX crystal morphology in the CYC solvent. We hope this work can provide some theoretical supports for high-quality RDX preparation.

Figure 1. Principle diagram for the calculation of the solvent attachment energy.

attachment energies of a crystal face decrease in solution. The solvent attachment energy therefore may be regarded as the energy correction term for the vacuum attachment energy, and it is calculated by the following equation: sol‐att ads Ehkl = ΔEhkl

sol Nhkl Nhkl

(1)

Nsol hkl

where Nhkl and are the molecule numbers in the oncoming crystal and solvent (hkl) slices, respectively; ΔEads hkl is the specific adsorption energy between solvent molecules and the (hkl) crystal face (the calculation is seen in Computation Simulations). The modified attachment energy (Emod‑att ) of a crystal face in hkl solution then can be calculated by the vacuum attachment energy minus the solvent attachment energy, as shown in eq 2.

2. MODEL AND COMPUTATIONAL SIMULATIONS 2.1. Model. The attachment energy (AE) model is a powerful tool to predict the crystal morphology, which takes into account the anisotropic energies in the crystal unit cell.33 The attachment energy (Eatt hkl) is defined as the energy released by adding a growth unit of one interplanar thickness dhkl onto the crystal face. It assumes that the growth rate of a crystal face is proportional to its attachment energy; that is, the face with the lower attachment energy is the face with slower growing rates and hence has the higher morphological importance. The AE model can give a reasonable prediction for crystal growth in the vapor phase.34,35 However, most crystal habits grown from solution are not in agreement with those predicted by the AE model. The difference may be attributed to the fact that the model mostly takes into account the internal structures of a given crystal but disregards the external crystallization factors such as the solvent. In solution, the adsorption of solvent and solute on the surface is determined by competition. It supposes that solvents can delay the crystal growth.11,28,36 If the adsorption interactions between the crystal face and solvent are strong, the adsorbed solvent molecules can occupy the bonding sites on the surface and block the incorporation of solute molecules; thus the solvation surface layer exists at the crystal−liquid interface and inhibits the growth of the corresponding crystal face. Obviously, crystal growth

mod‐att att sol‐att Ehkl = Ehkl − Ehkl

(2)

The relative growth rate of an (hkl) face is proportional to its attachment energy;33 hence correspondingly in the solvent we have sol mod‐att R hkl ∝ |Ehkl |

(3)

where Rsol hkl is the relative growth rate of an (hkl) crystal face in the solvent. Finally, the modified attachment energy (MAE) model incorporating the solvent effect term is established. Next, the MAE model with molecular dynamics (MD) simulations is applied to predict the crystal morphology of RDX in CYC solvent. 2.2. Computational Simulations. All computation simulations were completed in Materials Studio 4.2 software.37 The initial structure of the RDX unit cell is displayed in Figure 2a, which has eight irreducible RDX molecules with the lattice parameters a = 13.182 Å, b = 11.574 Å, c = 10.709 Å, and α = β = γ = 90°.38 An RDX molecule consists of three CH2−N−NO2 units arranged in a six-member ring as shown in Figure 2b, in which 11472

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

Figure 2. RDX (a) unit cell structure and (b) molecular configuration.

two nitro groups occupy axial positions (A) and the remaining nitro group is in the pseudoequatorial position (E). The morphology of RDX crystal in a vacuum was predicted by the AE model, which determined the morphologically important (hkl) faces. Subsequently, RDX crystal was cleaved parallel to the predicted (hkl) faces and constructed as a periodic superstructure. The partial charges of RDX and CYC molecules were determined by the built-in ab initio COMPASS force field. The COMPASS force field is a powerful ab initio force field,39 which has been confirmed to be able to give accurate prediction of structures and properties for condensed-phase materials including energetic compounds.29,40−42 In spite of this, since the partial charges were important sensitive parameters for the molecular crystals, the COMPASS force field was still tested via calculating the lattice energy of RDX crystal and the liquid properties of CYC solvent such as density and solubility parameter. The former was calculated by the geometry optimization method, and the latter were obtained by simulating the CYC solvent phase. The three-dimensional periodic solvent box with 100 random distributed CYC molecules was constructed by the Amorphous Cell tool and refined by MD techniques. The dimensions of the solvent box were consistent with the lattice parameters of the selected crystal face. For CYC solvent density, MD simulations with NPT ensemble were carried out at 298 K and 0.1 MPa, controlled by an Andersen thermostat43 and a Berendsen barostat.44 In the case of the solubility parameter, MD simulations with NVT ensemble were performed at the same temperature by the Andersen thermostat. Geometry optimizations were used before respective MD simulations were run with the time of 500 ps. The crystal surface−solvent interface model was used to study the influence of CYC solvent on RDX crystal morphology, which is displayed in Figure 3. To prevent the surface RDX molecules from diffusing into the solvent layer, RDX surface was constrained along a, b, and c directions, and the CYC solvent layer was placed onto the crystal surface along the c axis, in which all CYC molecules were able to move freely in the simulation. An 80 Å thick vacuum layer was set above the solvent phase to eliminate the effect of additional free boundaries. All the geometry optimizations and MD simulations were performed in the COMPASS force field. The initial interface models were first carried out as geometry optimizations. Then, MD simulations with NVT ensemble were carried out at 298 K, which is the usual crystallization temperature in CYC solvent32 controlled by the Andersen thermostat. For the equilibration stage, a period of 1 ns dynamics simulations with a time step of

Figure 3. Schematic diagram of the RDX surface−CYC solvent adsorption interface model.

0.1 fs was run. After the system equilibrium, the production stage was performed with the time of 500 ps during which the dynamics trajectories were collected every 100 time steps. For the potential energy calculation, the atom based method with a cutoff distance of 12.5 Å was set to calculate van der Waals (vdW) interactions and the Ewald summation method was applied to calculate Coulomb interactions. It can be considered that the total crystal surface−solvent adsorption interactions are mainly attributed to the specific interactions within a distance of several angstroms close to the interplanar distance near the surface.17 In this work, the CYC molecules within the solvent slices of thickness dhkl near the RDX surface were regarded as the CYC surface adsorption layer, and they were used to calculate the CYC solvent−RDX surface specific adsorption interactions. The specific ΔEads hkl values were averaged on the final 500 frames of the production trajectory (last 50 000 steps), computed by the following equation: ads ΔEhkl =

1 500

500

∑ (Etot, i − E sur, i−Esol,i)

(4)

i=1

where Etot,i is the total energy of each frame; Esur,i and Esol,i are the energies of the isolated crystal surface and solvent layer of each frame, respectively. Obviously, ΔEads hkl is a statistical result derived from the ensemble average of MD simulations, and is more accurate than the calculation of the random one from the static energy minimization. For the calculation of Nsol hkl, a CYC molecule is counted if the distance of any CYC molecule to RDX surface were within a thickness dhkl. Such a distance (r) is defined as the distance from the centroid of a CYC molecule to the centroid of the nearest RDX surface molecule. Consequently, the CYC adsorption numbers Nsol hkl were given by eq 5. sol Nhkl =

1 500

500,frame 100,solvent

∑

∑

i=1

j=1

⎧1 ⎨ ⎩0 ⎪

⎪

r ≤ dhkl ⎫ ⎬ r > dhkl ⎭ ⎪

⎪

(5)

Nsol hkl,

For it is also the ensemble average result of MD simulations. At last, the modified attachment energies of RDX 11473

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

that the calculated RDX lattice energy is close to the experimental value. Besides, the computed density and solubility parameter of CYC solvent are also in agreement with the corresponding experimental values. As a result, the test results indicate that COMPASS force field parameters are a validation for RDX crystal and CYC solvent. 3.2. RDX Crystal Morphology in a Vacuum. The crystal morphology of RDX in a vacuum predicted by the AE model is shown in Figure 4, which comprises the (111), (020), (002), (200), and (210) faces. The calculated crystal habit parameters of RDX are presented in Table 2. It can be found that, among these

habit faces were obtained according to eq 1, and the corresponding crystal morphology of RDX in CYC solvent was predicted in the Morphology module that allows researchers both to study particle shape and to consider the effects of altering the growth rate of particular faces on crystal morphology.

3. RESULTS AND DISCUSSION 3.1. Validation of COMPASS Force Field. The validation of the COMPASS force field can be evaluated by calculating the lattice energy of RDX crystal and the solvent properties of CYC such as the density and solubility parameter. The lattice energy (Elatt) of the crystal is defined as the total internal (intermolecular plus intramolecular) energy of the molecule in the crystal minus the corresponding energy of the molecule in the gas phase.45 It is also determined from the experimental sublimation enthalpy (ΔHsub) as shown in eq 6.46 E latt = −ΔHsub − 2RT (6)

Table 2. Crystal Habit Parameters of RDX in a Vacuum Calculated by AE Model

where R and T are the gas constant and temperature, respectively; the 2RT factor is a commonly accepted correction to take into account the zero point energy and thermal corrections at 298 K. Therefore, the reliability of the force field is assessed by the comparison of the calculated lattice energy with the “experimental” lattice energy. The Hildebrand solubility parameter (δ) is calculated as follows:47 δ=

ΔE = V

(7)

in which ΔE is the cohesive energy, ΔHvap is the vaporization enthalpy, and V is the molar volume of the solvent. The density and solubility parameter of CYC solvent as well as the RDX lattice energy calculated by the COMPASS force field at 298 K and 0.1 MPa are shown in Table 1. From Table 1, it is seen Table 1. Density and Solubility Parameter of CYC Solvent As Well As the RDX Lattice Energy Calculated by COMPASS Force Field at 298 K and 0.1 MPa Elatt/kcal·mol−1 expt COMPASS

−30.10 −27.90

48

ρ/g·cm−3 0.948 0.952

49

dhkl/Å

Eatt/kcal·mol−1

(210) (200) (002) (020) (111)

5.73 6.59 5.35 5.79 6.75

−133.26 −131.47 −108.90 −108.10 −103.32

important growth faces, the (111) face has the minimum attachment energy with 103.32 kcal·mol−1 and the maximum attachment energy of 133.26 kcal·mol−1 belongs to the (210) face. In a vacuum, the morphological importance of the (111) face is the strongest, and the (210) face has the weakest morphological importance. The RDX crystal surfaces have different molecular packing orientations since the molecular structure of RDX is noncentrosymmetric.31 The molecular packing diagrams of RDX habit faces are displayed in Figure 4. It is observed that the molecular arrangements of the (111) face are relatively flat at the molecular level, while the (020), (002), (210), and (200) faces show open and rough surface topographies with large voids. To describe the crystal surface features quantitatively, we introduce a parameter (S) defined as the area ratio of the solvent-accessible surface to the corresponding (hkl) face.23 The calculated S values of different RDX habit faces are listed in Table 3. It is found that the (111) face has the minimum S value of 1.02, which indicates morphological smoothness. The (200), (210), (002), and (020) faces have larger S values, of which the S value of the (200) face is the maximum with 1.44, indicating morphological relative roughness. The larger S value means the more complex surface

ΔH vap − RT V

(hkl)

δ/(J/cm3)1/2 20.2150 18.43

Figure 4. Crystal morphology and surface structures of RDX in a vacuum predicted by AE model. 11474

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

Table 3. S Values of RDX Habit Facesa (hkl) Aacc hkl Ahkl S

(111) 273.14 267.80 1.02

(020) 162.89 141.17 1.15

(002) 181.55 152.57 1.19

(200) 178.82 123.95 1.44

RDX surface, which implies the existence of strong adsorption interactions at the interface. From Figure 5, the adsorption orientations between CYC molecules and RDX surface can be observed: the ketone O atoms of CYC molecules point to the H atoms of RDX surfaces, while the RDX nitro O atoms attract the H atoms of CYC molecules. The solvent−surface adsorption interactions are essential intermolecular interactions. In general, intermolecular interactions include short-range interactions and long-range interactions. The short-range interactions are mainly hydrogen bonding (H-bonding) interactions and vdW interactions, whose effective distances are 2.6−3.1 Å and 3.1−5.0 Å, respectively. The long-range interactions are electrostatic interactions (also Coulomb interactions) whose distance range is usually above 5.0 Å. The composition of adsorption interactions can be analyzed by means of the radial distribution function (RDF). Taking the adsorption systems of the (210) and (020) faces as examples, the RDFs between CYC and RDX atoms are plotted in Figure 6. Figure 6a depicts the RDFs of CYC_O1 atoms∼RDX_H2 atoms. It is seen that the highest peak exists in the range 2.6−3.1 Å, implying the existence of hydrogen bonds. Some higher peaks appear in the range 3.1−5.0 Å, indicating strong vdW interactions; besides, it also has small peaks beyond the 5.0 Å region, suggesting weak electrostatic interactions. Figure 6b displays the RDFs of RDX_O2∼CYC_H1 atoms. It can be found that the higher peaks appear beyond the 5.0 Å region. Representing the system are main long-range interactions; the weak shoulder peaks are also seen in both the ranges 2.6−3.1 and 3.1−5.0 Å, indicating the existence of Hbonding interactions and vdW interactions among O2 and H1 atoms. It should be pointed out that the formation of hydrogen bonds among CYC and RDX molecules are the manner of C− H···O, which belongs to the weak hydrogen bond category.53

(210) 341.22 285.27 1.20

a acc Ahkl

is the solvent-accessible area of an (hkl) face in the unit cell,51 calculated by the Connolly surface model;52 Ahkl is the surface area of the corresponding (hkl) face. Area unit is Å2.

topography, which is thus more convenient for solvent molecules incorporation. It is also observed from Figure 4 that the strong polar nitro groups (NO2) are exposed to (210), (200), (002), and (111) surfaces, which may contribute to attracting any oncoming solvent molecules, whereas the (020) face only exposes the nonpolar methylene groups (CH2). From surface chemistry, the (210), (002), and (200) faces as well as the (111) face may be classified as polar surfaces due to the exposed polar groups (e.g., NO2 groups), and the (020) face is considered as a nonpolar surface because of the appearance of nonpolar CH2 groups. 3.3. Effect of CYC Solvent. 3.3.1. Interface Structure. Due to the different surface structures and surface chemistry of each crystal face of RDX, the adsorption behaviors of CYC solvent on the different RDX interfaces may be complicated. The snapshots of RDX surface−CYC solvent interface structures after MD equilibrium are shown in Figure 5, which displays the different adsorption configurations between CYC solvent and RDX planes. As seen from Figure 5, CYC molecules have attached on RDX surfaces and formed CYC surface adsorption layers at the interface zones. Also, the distributions of CYC molecules are mainly concentrated on the solvent-accessible surfaces of RDX crystal faces. This indicates that the CYC surface adsorption layers within thickness dhkl have already closely attached on the

Figure 5. Snapshots of different RDX surface−CYC solvent interface models after MD equilibrium. Colors: the light blue represents the solventaccessible surface of the RDX habit face; the blue and purple arrows represent the adsorption orientations between the CYC molecule and the RDX surface molecule. 11475

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

Figure 6. RDFs between CYC atoms and RDX atoms for the adsorption systems of (210) and (111) faces. The O and H atoms in the CYC molecule are named O1 and H1, while the H and O atoms of the RDX molecule are named O2 and H2, respectively.

Therefore, it can be concluded that the adsorption interactions between RDX surface and CYC solvent are primarily the weak C−H···O interactions, strong vdW interactions, and electrostatic interactions. 3.3.2. Solvent Attachment Energy. The CYC solvent attachment energies of RDX habit faces calculated by eq 1 are shown in Table 4 (the calculation details are listed in the

of RDX crystal is stronger than on the nonpolar (020) face, in which that on the (210) face should be the strongest. 3.4. Crystal Morphology of RDX in CYC Solvent. Due to the CYC solvent interactions, the vacuum attachment energies of RDX crystal faces will be reduced to a different extent. The modified attachment energies of RDX habit faces in CYC solvent calculated by eq 2 are summarized in Table 5. As seen from Table

Table 4. Calculated CYC Solvent Attachment Energies of RDX Habit Facesa

Table 5. Modified Attachment Energies and the Relative Growth Rates of RDX Habit Faces in CYC Solvent Calculated by the MAE Modela

(hkl) ΔEads hkl Nsol hkl/Nhkl Esol‑att hkl a

(111) −57.16 0.88 −50.30

(020) −29.40 0.95 −27.93

(002) −38.38 0.98 −37.62

(200) −37.59 1.10 −41.35

(210) −74.39 1.04 −77.37

(hkl) Eatt hkl Emod‑att hkl Rvac hkl Rsol hkl

Energy unit is kcal·mol−1. a

Supporting Information). It is seen that the (210) face has the maximum Esol‑att of −77.37 kcal·mol−1, while the (020) face has hkl the minimum Esol‑att of −27.93 kcal·mol−1 (the minus sign hkl represents an exothermic process). In addition, the Esol‑att order hkl of RDX crystal faces is written as (210) > (111) > (200) > (002) > (020). That is, the CYC solvent attachment energies of RDX polar faces are much larger than that of the nonpolar (020) face. This is in accordance with the empirical fact that, during crystal growth, polar solvents (such as CYC) can preferentially adsorb on the crystal polar surfaces, whereas nonpolar solvents tend to bind on the nonpolar faces.54,55 The solvent attachment energy quantitatively reflects the solvent effect. If the solvent attachment energy is larger, the effect of CYC solvent on the RDX surfaces is stronger. Therefore, the effect of CYC solvent on the polar faces

(111) −103.32 −53.02 1 1

(020) −108.10 −80.17 1.05 1.51

(002) −108.90 −71.28 1.05 1.34

(200) −131.47 −90.12 1.22 1.70

(210) −133.26 −55.89 1.22 1.05

Energy unit is kcal·mol−1.

5, the (111) face still has the minimum Emod‑att of −53.02 kcal· hkl mol−1, and the Emod‑att value of the (210) face decreases to hkl −55.89 kcal·mol−1 close to the (111) face, while the (200) face has the maximum Emod‑att of −90.01 kcal·mol−1 (the minus sign hkl represents an exothermic process). The rank of modified attachment energies of RDX habit faces in CYC solvent becomes (111) < (210) < (002) < (020) < (200). According to eq 3, the growth rate of an RDX crystal face relative to the (111) face is also presented in Table 5. It is found that, compared with that in a vacuum, the relative growth rate of the (210) face is slower and the relative growth rates of (020), (002), and (200) faces become faster, in which the growth of the (200) face is the fastest. If RDX crystallizes from CYC solvent, the morphological importance of 11476

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

work can provide some theoretical support for explosive morphology control technology. The next step is to test the applicability of the MAE model in other RDX crystallization systems.

both (210) and (111) faces will be the largest, while the morphological importance of (002), (020), and (200) faces will reduce, in which the decrease of the (020) face is the most significant due to its weakest CYC solvent effect. The RDX crystal morphology grown from CYC solvent predicted by the MAE model is shown in Figure 7a, and the

■

ASSOCIATED CONTENT

S Supporting Information *

System details of RDX surface−CYC solvent adsorption interface models calculated for the modified attachment energy. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: 0086-25-84315190. Fax: 0086-25-84315190. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Teipel, U. Energetic Materials: Particle Processing and Characterization; Wiley-VCH: Weinheim, Germany, 2005. (2) Lochert, I. J.; Dexter, R. M.; Hamshere, B. L. Evaluation of Australian RDX in PBXN-109; Systems Sciences Laboratory: Edinburgh, South Australia, 2002. http://hdl.handle.net/1947/3586. (3) Halfpenny, P.; Roberts, K.; Sherwood, J. Dislocations in energetic materials: IV. The crystal growth and perfection of cyclotrimethylene trinitramine (RDX). J. Cryst. Growth 1984, 69 (1), 73−81. (4) Tillotson, T.; Hrubesh, L.; Simpson, R.; Lee, R.; Swansiger, R.; Simpson, L. Sol−gel processing of energetic materials. J. Non-Cryst. Solids 1998, 225, 358−363. (5) Doherty, R. M.; Watt, D. S. Relationship between RDX properties and sensitivity. Propellants, Explos., Pyrotech. 2008, 33 (1), 4−13. (6) van der Heijden, A. E.; Creyghton, Y. L.; Marino, E.; Bouma, R. H.; Scholtes, G. J.; Duvalois, W.; Roelands, M. C. Energetic materials: crystallization, characterization and insensitive plastic bonded explosives. Propellants, Explos., Pyrotech. 2008, 33 (1), 25−32. (7) Roberts, C. W.; Hira, S. M.; Mason, B. P.; Strouse, G. F.; Stoltz, C. A. Controlling RDX explosive crystallite morphology and inclusion content via simple ultrasonic agitation and solvent evaporation. CrystEngComm 2011, 13 (4), 1074−1076. (8) van der Heijden, A. E.; Bouma, R. H. Crystallization and Characterization of RDX, HMX, and CL-20. Cryst. Growth Des. 2004, 4 (5), 999−1007. (9) Spyckerelle, C.; Eck, G.; Sjöberg, P.; Amnéus, A. M. Reduced sensitivity RDX obtained from Bachmann RDX. Propellants, Explos., Pyrotech. 2008, 33 (1), 14−19. (10) Kim, J.-W.; Park, D. B.; Shim, H.-M.; Kim, H.-S.; Koo, K.-K. Crystallization of RDX by drowning-out combined with fines dissolution and cooling process. Ind. Eng. Chem. Res. 2012, 51 (9), 3758−3765. (11) Lahav, M.; Leiserowitz, L. The effect of solvent on crystal growth and morphology. Chem. Eng. Sci. 2001, 56 (7), 2245−2253. (12) Bhat, M. N.; Dharmaprakash, S. Effect of solvents on the growth morphology and physical characteristics of nonlinear optical γ-glycine crystals. J. Cryst. Growth 2002, 242 (1), 245−252. (13) Thompson, C.; Davies, M. C.; Roberts, C. J.; Tendler, S. J.; Wilkinson, M. J. The effects of additives on the growth and morphology of paracetamol (acetaminophen) crystals. Int. J. Pharm. (Amsterdam, Neth.) 2004, 280 (1), 137−150. (14) Siegfried, M. J.; Choi, K.-S. Elucidating the effect of additives on the growth and stability of Cu2O surfaces via shape transformation of pre-grown crystals. J. Am. Chem. Soc. 2006, 128 (32), 10356−10357. (15) Clydesdale, G.; Roberts, K. J.; Docherty, R. HABIT95a program for predicting the morphology of molecular crystals as a function of the growth environment. J. Cryst. Growth 1996, 166 (1−4), 78−83.

Figure 7. (a) RDX crystal morphology grown from CYC solvent predicted by the MAE model; (b) corresponding experimental shape.47

corresponding experimental shape obtained by the cooling crystallization56 is shown in Figure 7b. It is seen that the predicted RDX morphology is in reasonable agreement with the experimental result. This indicates that the proposed MAE model is valid for the prediction of RDX crystal morphology in solvents.

4. CONCLUSION In this paper, molecular dynamics simulations have been performed to investigate the effect of CYC solvent on the growth of RDX crystal. The RDX surface−CYC solvent adsorption interface models are built to study the adsorption interactions between RDX surfaces and CYC molecules. A modified attachment energy model incorporating the solvent effect term is established and applied to predict the RDX crystal morphology in CYC solvent, and it is compared with the experiment. The major results presented in this work are summarized as follows: The calculation of the AE model indicates that the crystal habits of RDX in a vacuum are dominated by the (111), (020), (200), (002), and (210) faces, in which the morphological importance of (111) is the strongest and that of the (210) face is the weakest. The analysis of radial distribution functions shows that the adsorption interactions between CYC solvent and RDX surface comprise weak C−H···O interactions, vdW interactions, and electrostatic interactions. The calculation of the solvent attachment energy suggests that the effect of CYC solvent on the polar faces of RDX crystal is stronger than that on the nonpolar (020) face, among which that on the (210) face is the strongest. The calculation of the MAE model indicates that, due to the CYC solvent effect, (210) and (111) faces have the largest morphological importance, while the morphological importance of (002), (020), and (200) faces reduces. The predicted RDX crystal morphology is in reasonable agreement with the experimental result. In a word, this paper explores a reliable theoretical model for the systematic examination of the solvent effect on the RDX crystal habits with the aid of computer simulation. We hope this 11477

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478

The Journal of Physical Chemistry A

Article

(16) Clydesdale, G.; Roberts, K. J.; Saunders, V. R.; Pugh, D.; Jackson, R. A.; Meenan, P. Prediction of the polar morphology of sodium chlorate using a surface-specific attachment energy model. J. Phys. Chem. B 1998, 102 (36), 7044−7049. (17) Winn, D.; Doherty, M. F. A new technique for predicting the shape of solution-grown organic crystals. AIChE J. 1998, 44 (11), 2501− 2514. (18) Winn, D.; Doherty, M. F. Modeling crystal shapes of organic materials grown from solution. AIChE J. 2000, 46 (7), 1347−1368. (19) Lin, C. H.; Gabas, N.; Canselier, J. P.; Pèpe, G. Prediction of the growth morphology of aminoacid crystals in solution: I. α-Glycine. J. Cryst. Growth 1998, 191 (4), 791−802. (20) Stoica, C.; Meekes, H.; van Hoof, P. J. C. M.; Kaspersen, F. M.; Vlieg, E. Understanding the Effect of a Solvent on the Crystal Habit. Cryst. Growth Des. 2004, 4 (4), 765−768. (21) Piana, S.; Gale, J. D. Understanding the barriers to crystal growth: Dynamical simulation of the dissolution and growth of Urea from aqueous solution. J. Am. Chem. Soc. 2005, 127 (6), 1975−1982. (22) Piana, S.; Reyhani, M.; Gale, J. D. Simulating micrometre-scale crystal growth from solution. Nature 2005, 438, 70−73. (23) 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 (1), 175−180. (24) Lee, H.-E.; Lee, T. B.; Kim, H.-S.; Koo, K.-K. Prediction of the growth habit of 7-amino-4, 6-dinitrobenzofuroxan mediated by cosolvents. Cryst. Growth Des. 2009, 10 (2), 618−625. (25) 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. 2012, 13 (1), 282−290. (26) Shi, W.; Xia, M.; Lei, W.; Wang, F. Solvent effect on the crystal morphology of 2, 6-diamino-3, 5-dinitropyridine-1-oxide: A molecular dynamics simulation study. J. Mol. Graphics Modell. 2014, 50, 71−77. (27) Ter Horst, J.; Geertman, R.; Van der Heijden, A.; Van Rosmalen, G. The influence of a solvent on the crystal morphology of RDX. J. Cryst. Growth 1999, 198, 773−779. (28) Ter Horst, J.; Geertman, R.; Van Rosmalen, G. The effect of solvent on crystal morphology. J. Cryst. Growth 2001, 230 (1), 277−284. (29) Chen, G.; Xia, M.; Lei, W.; Wang, F.; Gong, X. A study of the solvent effect on the morphology of RDX crystal by molecular modeling method. J. Mol. Model. 2013, 19 (12), 5397−5406. (30) Chen, G.; Xia, M.; Lei, W.; Wang, F.; Gong, X. D. Molecular dynamics investigation of the effect of solvent adsorption on crystal habits of Hexogen. Can. J. Chem. 2014, 92 (9), 849−854. (31) Shim, H.-M.; Koo, K.-K. Crystal Morphology Prediction of Hexahydro-1, 3, 5-trinitro-1, 3, 5-triazine by the Spiral Growth Model. Cryst. Growth Des. 2014, 14 (4), 1802−1810. (32) Kim, D.-Y.; Kim, K.-J. Solubility of cyclotrimethylenetrinitramine (RDX) in binary solvent mixtures. J. Chem. Eng. Data 2007, 52 (5), 1946−1949. (33) Hartman, P.; Bennema, P. The attachment energy as a habit controlling factor: I. Theoretical considerations. J. Cryst. Growth 1980, 49 (1), 145−156. (34) Hartman, P. The attachment energy as a habit controlling factor II. Application to anthracene, tin tetraiodide and orthorhombic sulphur. J. Cryst. Growth 1980, 49 (1), 157−165. (35) Hartman, P. The attachment energy as a habit controlling factor: III. application to corundum. J. Cryst. Growth 1980, 49 (1), 166−170. (36) Chen, J.; Trout, B. L. Computer-aided solvent selection for improving the morphology of needle-like crystals: A case study of 2, 6dihydroxybenzoic acid. Cryst. Growth Des. 2010, 10 (10), 4379−4388. (37) Materials Studio 4.0; Accelrys Software Inc.: San Diego, CA, USA, 2007. (38) Choi, C. S.; Prince, E. The crystal structure of cyclotrimethylenetrinitramine. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1972, 28 (9), 2857−2862. (39) Sun, H. COMPASS: an ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds. J. Phys. Chem. B 1998, 102 (38), 7338−7364.

(40) Bunte, S. W.; Sun, H. Molecular modeling of energetic materials: the parameterization and validation of nitrate esters in the COMPASS force field. J. Phys. Chem. B 2000, 104 (11), 2477−2489. (41) Qiu, L.; Xiao, H.-M.; Zhu, W.-H.; Xiao, J.-J.; Zhu, W. Ab initio and molecular dynamics studies of crystalline TNAD (trans-1, 4, 5, 8tetranitro-1, 4, 5, 8-tetraazadecalin). J. Phys. Chem. B 2006, 110 (22), 10651−10661. (42) Xu, X.-J.; Xiao, H.-M.; Xiao, J.-J.; Zhu, W.; Huang, H.; Li, J.-S. Molecular dynamics simulations for pure ε-CL-20 and ε-CL-20-based PBXs. J. Phys. Chem. B 2006, 110 (14), 7203−7207. (43) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72 (4), 2384−2393. (44) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (45) Singh, M. Predicting lattice energy and structure of molecular crystals by first-principles method: Role of dispersive interactions. J. Cryst. Growth 2014, 396, 14−23. (46) Gavezzotti, A. Structure and intermolecular potentials in molecular crystals. Modell. Simul. Mater. Sci. Eng. 2002, 10 (3), R1. (47) Hildebrand, J.; Scott, R. The Solubility of Nonelectrolytes; Reinhold: New York, 1950; p 411. (48) Rosen, J. M.; Dickinson, C. Vapor pressures and heats of sublimation of some high-melting organic explosives. J. Chem. Eng. Data 1969, 14 (1), 120−124. (49) Radhamma, M.; Venkatesu, P.; Rao, M.; Lee, M.-J.; Lin, H.-m. Excess molar volumes and ultrasonic studies of dimethylsulfoxide with ketones at T = 303.15 K. J. Chem. Thermodyn. 2008, 40 (3), 492−497. (50) Liang, B. Polymer Physics; Chinese Textile Press: Shanghai, 2000. (51) Lee, B.; Richards, M. F. The interpretation of protein structures: Estimation of static accessibility. J. Mol. Biol. 1971, 55 (3), 379−400. (52) Connolly, M. L. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983, 221 (4612), 709−713. (53) Desiraju, G. R. Hydrogen bridges in crystal engineering: interactions without borders. Acc. Chem. Res. 2002, 35 (7), 565−573. (54) Berkovitch-Yellin, Z. Toward an ab initio derivation of crystal morphology. J. Am. Chem. Soc. 1985, 107 (26), 8239−8253. (55) Nokhodchi, A.; Bolourtchian, N.; Dinarvand, R. Crystal modification of phenytoin using different solvents and crystallization conditions. Int. J. Pharm. (Amsterdam, Neth.) 2003, 250 (1), 85−97. (56) Yu, X.-H.; Sima, T.-L.; Sun, K.-D. Analysis of Mechanical Properties of RDX Crystals Obtained from Different Solvents. Chin. J. Energ. Mater. 2004, 12 (2), 78−81.

11478

dx.doi.org/10.1021/jp508731q | J. Phys. Chem. A 2014, 118, 11471−11478