The effect of crystal-solvent interaction on crystal growth and morphology

The effect of crystal-solvent interaction on crystal growth and morphology

Accepted Manuscript The Effect of Crystal-Solvent Interaction on Crystal Growth and Morphology Jing-Wen Li, Shu-Hai Zhang, Rui-Jun Gou, Gang Han, Ming...

2MB Sizes 1 Downloads 51 Views

Accepted Manuscript The Effect of Crystal-Solvent Interaction on Crystal Growth and Morphology Jing-Wen Li, Shu-Hai Zhang, Rui-Jun Gou, Gang Han, Ming-Hua Chen PII: DOI: Reference:

S0022-0248(18)30580-3 https://doi.org/10.1016/j.jcrysgro.2018.11.007 CRYS 24841

To appear in:

Journal of Crystal Growth

Received Date: Revised Date: Accepted Date:

26 July 2018 13 October 2018 10 November 2018

Please cite this article as: J-W. Li, S-H. Zhang, R-J. Gou, G. Han, M-H. Chen, The Effect of Crystal-Solvent Interaction on Crystal Growth and Morphology, Journal of Crystal Growth (2018), doi: https://doi.org/10.1016/ j.jcrysgro.2018.11.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

J.-W.Li et al / Manuscript

The Effect of Crystal-Solvent Interaction on Crystal Growth and Morphology Jing-Wen Li,1 Shu-Hai Zhang,1 Rui-Jun Gou,1 Gang Han,1 Ming-Hua Chen,2 1

School of Environment and Safety Engineering, North University of China, Taiyuan 030051, Shanxi China

2

New Technology Applications Institute of Shijiazhuang, Shijiazhuang 050000, Shijiazhuang Hebei

Abstract Diacetone diperoxide (DADP)/1,3,5-tribromo-2,4,6-trinitrobenzen (TBTNB) and diacetone

diperoxide

(DADP)/1,3,5-trichloro-2,4,6-trinitrobenzen

(TCTNB)

cocrystals were selected to reveal the effect of interaction at solvent-crystal interface on crystal growth and morphology, and analyze the interaction at the interface by molecular electrostatic potentials (ESP), mass density distribution and radial distribution function (RDF). It can be found that in DADP/TBTNB-acetonitrile interface model, the interaction at (0 2 0) face is the greatest and presents as repulsive interaction, so the growth rate of (0 2 0) face is the largest and the (0 2 0) face is the smallest. Conversely, the least and attractive interaction at interface tends to make (1 1 -1) to be the largest. These deductions are in accordance with the simulation of DADP/TBTNB cocrystal morphology in acetonitrile, and similar conclusion also can be obtained from the growth of DADP/TBTNB cocrystal in cyclohexanone. These are related to cavities of cocrystal face, which affect the distribution and orientation of solvent molecules. It can be reasoned out that the quantity of cavities on crystal face and the electrophilicity or nucleophilicity of solvent molecules may inhibit or facilitate the growth of crystal face.

Keywords: A1. Crystal morphology; A2. Crystal growth; A1. Interfacial interaction; A1. MD simulation 

Corresponding author. Fax: +86 351 3920588 Tel.: +86 13994272868 E-mail address: [email protected] (S.-H.Zhang). 1

J.-W.Li et al / Manuscript

1. Introduction Nowadays, energetic materials are widely used in filling explosives, propellants and pyrotechnics [1]. Not only do energetic materials possess the higher energy than traditional explosives, but they have better properties and performances as well. For the study of energetic materials, from excessively pursuing the high energy to hunting novel materials with high density of energy and favorable comprehensive performances, the performance and safety of energetic materials are paid more and more attention [2]. The energetic materials also need to meet a lot of requirements, such as storage, powder charge and usage, etc. These requirements are dependent on the properties of materials. Looking for new energetic materials with favorable properties or finding methods to improve the comprehensive performances are of high interest in the energetic materials community [3]. There are roughly three ways to upgrade the power and promote the performances now: synthesizing novel energetic entities, seeking out the materials of novel conception and changing the energetic materials from micro [4]. Considering the facts that synthesizing and finding the new materials are challenging and time-consuming, improving the material’s properties is a high-efficiency and clearly improved way among these three methods [5]. The performances of energetic materials are influenced by the qualities of crystal, such as crystal morphology, purity, and granularity, etc. If the crystal morphology of energetic materials is sharp, the ‘hot-spot’ is possible to appear by friction and collision among the sharp crystals, which leads to insecure explosive. Therefore the crystal morphology is crucial for the performances and security of energetic materials. The crystal morphology can be affected by the external factors (solvent, additives and supersaturation) and internal factors (molecular structure and functional group), and the solvent plays a more important role on crystal morphology [6-9]. For the influence of solvent on the crystal morphology, pioneers had proposed some theories about how growth of crystal could be enhanced or inhibited by the solvent. Bennema et al. indicated that favorable interactions between solvent and solute on specific faces reduced interfacial tension 2

J.-W.Li et al / Manuscript

and caused a transition from smooth to rough interface, hence the growth of surface would be fast; on the other side, the growth of specific crystal face would be inhibited by the preferential adsorption [10]. Leiserowitz et al. pointed out that the growth of crystal face would be inhibited by the exposed functional groups which have a strong affinity with solvent, and the surface tends to be large. Otherwise the growth of crystal face would be enhanced by the exposed functional groups which have a weak affinity or strong repulsion with solvent, and the surface tended to be small [11]. Recently, Gou et al. reported that the density of solvent at the ‘basin’ of the rough crystal face is larger than the smooth one [12]. Accordingly, it can be speculated that the strength and affinity of the interaction between the solvent and crystal would be judged by the distribution of solvent molecules. By this token, it is important to explore the effect of the strength and affinity of interactions between solvent and crystal molecules on crystal growth and morphology. Diacetone (TBTNB,

diperoxide

(DADP,

C6N3Br3O6)

C6H12O4)/1,3,5-tribromo-2,4,6-trinitrobenzen and

diacetone

diperoxide

(DADP)/1,3,5-trichloro-2,4,6-trinitrobenzen (TCTNB, C6N3Cl3O6) are two cocrystals of energetic materials which were reported by Kira B. Landenberger in 2013. DADP is an acetone peroxide of energetic compound which is richly in electrons with accessible peroxide oxygen. TCTNB, 1,3,5-triamino-2,4,6-trinitrobenzen (TATB) precursor, is an energetic compound with electron-deficient aromatic ring which makes it ideally for cocrystallizing with DADP. For TBTNB, the brominated analogue of TCTNB, also could be designed as cocrystal with DADP. Both DADP/TBTNB and DADP/TCTNB are synthesized as 1:1. The molecular structures of DADP, TBTNB and TCTNB are shown in Fig. 1. The fundamental information of DADP/TBTNB and DADP/TCTNB could be found in Supporting Information (SI). [13]. They are two isostructural explosive cocrystals and their function groups are similar. Both DADP/TBTNB and DADP/TCTNB cocrystal are soluble in solvents of acetonitrile, cyclohexanone and heptane [14], and will be chosen to explore the solvent-cocrystal interaction effects on crystal growth and morphology. In general, experiment is a kind of complicated and expensive way. Nevertheless, computational chemistry is a sort of 3

J.-W.Li et al / Manuscript

research method for calculating the structures and performances of molecules, supramolecules, crystals and composites of various energetic materials [15]. The modified attachment energy (MAE) has been used to simulate the crystal morphologies which were grown from variety solvents, and the simulation results are in agreement with the experimental results [9, 14, 16-17]. Besides, Gao et al. have figured out that the cocrystal morphologies of DADP/TBTNB and DADP/TCTNB in acetonitrile, cyclohexanone, and heptane could be simulated by MAE model, and they managed to prove the feasibility of simulation [14]. In this study, the MAE model was selected to further explore the interaction effects on cocrystal growth and morphology of DADP/TBTNB and DADP/TCTNB.

Fig. 1 The molecular structures of DADP, TCTNB and TBTNB

2. Calculation Details 2.1 Molecular Dynamics Simulation 2.1.1 Theory The attachment energy (AE) model which is founded on the theory of periodic bond chain (PBC) was established by Hartman and Bennma. This model defines the energy released by a slice with thickness of d hkl when it was adsorbed on the (h k l) att growing surface as attachment energy Ehkl , it is used to simulate the crystal

morphology in vacuum [18]. For the growth of crystal in solvent, it needs more energy to remove the solvent molecules which absorbed on the crystal surface. Based on the AE model and in consideration of the effect of solvent, the interaction energy of solvent and morphologically important surface is taken to modify the attachment energy. So modified attachment energy (MAE) model is used to predict the crystal morphology in solvent condition [16, 19-20]. For the study, MAE model is chosen to 4

J.-W.Li et al / Manuscript

simulate the growth morphology of the DADP/TBTNB and DADP/TCTNB cocrystal in solvents. The modified attachment energy can be calculated by Eq. (2.1) [17]: mod  att att int Ehkl  Ehkl  Ehkl

S acce . S surf

(2.1)

mod  att att is the modified attachment energy of (h k l) face, kJ/mol. Ehkl is the Ehkl

attachment energy of (h k l) face in vacuum, kcal/mol. S acce is the surface area of a unit cell which can be touched by a solvent, nm2. S surf is the supercell surface area, int nm2. Ehkl is the interaction energy between solvent and (h k l) face molecule,

kcal/mol, it can be expressed by Eq. (2.2) [17]: int tot surf solv Ehkl  Ehkl  Ehkl  Ehkl

(2.2)

tot Among the formula of (2.2), Ehkl is the total energy of solvent-interface model surf which is the total system of (h k l) face layer and solvent layer, Ehkl is the energy of solv (h k l) face layer, Ehkl is the energy of solvent layer on the (h k l) face. All of the

units are kcal/mol. According to the definition of AE model, with the increase of the absolute value att of attachment energy Ehkl , time for crystal plate attachment is decreasing. So the

growth rate of (h k l) face in vacuum,

, is directly proportional to the absolute

att value of attachment energy Ehkl [18]. Thus the relatively growth rate of (h k l) face sol in solvent Rhkl is directly proportional to the absolute value of modified attachment mod  att energy Ehkl [21]. That is :

sol mod  att Rhkl  Ehkl .

(2.3)

2.1.2 MD Calculation In this research, Materials Studio (version 7.0) [22] software was used for Molecular Dynamics (MD) calculations. According to the cocrystal cell parameters of DADP/TBTNB and DADP/TCTNB [13], unit cell models were constructed for 5

J.-W.Li et al / Manuscript

calculation. For the force field, the condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) would be selected which has been verified by our group. [14, 23]. They validated it among three force fields, condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS), the consistent-valence forcefield (CVFF), and the polymer consistent forcefield (PCFF), by calculating the solvents and cocrystals’ densities and found that the simulation densities in COMPASS forcefield are the closest to the real densities among three forcefields. Besides, Quality was chosen Ultra-Fine (the cutoff distance is 18.5 Å). The initial unit cells of DADP/TBTNB and DADP/TCTNB cocrystal were first minimized by Geometry Optimization of Forcite Module, and then their cocrystal morphologies in vacuum were predicted by Growth morphology of Morphology Module. In accordance with the law of Bravais, the importance of crystal face is usually judged by the lattice networking density because the crystal faces are often parallel to the lattice networking density. The more lattice networking density a crystal face has, the more important the face is. For the crystal faces which are remained in solvents, they all have big lattice networking densities. They could be selected by the percent of total facet area, the percentage of each facetted face’s area in total crystal area, which is shown in the result of the growth morphology simulation. After that, the morphologically important faces of the crystal were cleaved, and the thickness of each slice was set as three because it’s the most stable setting with the least number of molecules [17, 24-25], then expanding each slice in different multiples. Meanwhile, resetting the orientation and origin of slice to make the V axis along the Y axis, and the U axis was in the XY plane. The multiples of DADP/TBTNB and DADP/TCTNB’s morphologically important faces are shown in Table 1. Table 1 The expanded multiples for cleaved morphologically important faces of DADP/TBTNB and DADP/TCTNB DADP/TBTNB

DADP/TCTNB

Morphologically important face

Multiple

Morphologically important face

Multiple

(1 1 0)

6×3

(1 1 −1)

3×4

(1 1 −1)

3×3

(1 1 0)

4×4 6

J.-W.Li et al / Manuscript

(0 2 0)

6×2

(0 2 0)

4×3

(2 0 0)

3×6

(1 1 1)

4×3

Owing to the lengths of a cocrystal slice sides on U and V axis, the solvent boxes with a certain number of randomly distributed solvent molecules were built. Both DADP/TBTNB and DADP/TCTNB are soluble in acetonitrile, cyclohexanone and heptane solvents. The number of randomly distributed solvent molecules of acetonitrile, cyclohexanone and heptane were set as 250, 100 and 100, respectively. Later, the solvent boxes were optimized under the fixed a, b side. For building models, the solvent layer was placed on the cocrystal layer along the Z axis and kept a gap of 5 Å vacuum above the cocrystal layer at the same time. Then put a 100 Å vacuum over the solvent layer to avoid the border effects [15]. The lattice parameters of new model were in line with the lattice parameters of cocrystal face layer. The first layer of the cocrystal face was restricted along the X, Y and Z axis. The ensemble with constant temperature and constant volume (NVT) was employed for MD. The conditions of NVT ensemble were set as 298 K and 500.0 ps (500,000 fs) for the total simulation time. The Anderson was used for thermostat [26]. Set 1.0 ps as a time step, and output a frame every 50 time steps during the running time. When the operation ended and the energy of the entire system reached equilibrium, the most stable structure was sought out from the last 3000 frames to calculate the modified attachment energy by using the equations (2.1) and (2.2) in Part 2.1.1. The values of

S acce

and

S surf

are shown in Table

2.

Choose (0 2 0) surface of

DADP/TBTNB-acetonitrile interface models as an example, the schematic diagram of interface model establishing is shown in Fig. 2. Finally, the cocrystal morphologies in three solvents were predicted by Morphology Module. Table 2 The

and

of DADP/TBTNB and DADP/TCTNB

DADP/TBTNB

DADP/TCTNB

Cocrystal face

Cocrystal face

(0 2 0)

132.74

1329.44

(0 2 0)

158.75

1566.72

(1 1 −1)

332.18

1393.41

(1 1 −1)

185.44

1728.60

(1 1 0)

114.23

1332.98

(1 1 0)

154.24

1521.29

(2 0 0)

117.96

1769.37

(1 1 1)

276.97

1759.90 7

J.-W.Li et al / Manuscript

Remark: Both

and

are in Å2.

Fig. 2 The schematic diagram of building DADP/TBTNB-acetonitrile interface model. Red, white, blue, green, rusty red and grey small ball represent oxygen, hydrogen, nitrogen, chlorine, bromine and carbon atom, respectively

2.2 Quantum-chemical Calculation Quantum-chemical calculations were operated on Gaussian 09 [27] software. Using the density functional theory, B3LYP [28, 29], with 6-311++g (d,p) basis setting to optimize the acetonitrile, cyclohexanone and heptane molecules. The density functional theory at the B3LYP/6-311+g (d) level was utilized to optimize the morphologically important faces. At last, the coloring isosurface graph of the electrostatic potentials of all the optimized morphologically important faces were drawn to show the distribution of ESP at different areas.

3. Results and Discussion 3.1 Crystal Habits The simulation cocrystal morphologies of DADP/TBTNB and DADP/TCTNB in vacuum and in solvents are shown in Fig. 3. The attachment energies and modified 8

J.-W.Li et al / Manuscript

attachment energies on morphologically important faces of DADP/TBTNB and DADP/TCTNB are shown in Tables 3-5, respectively.

Fig. 3 The simulation cocrystal morphologies of DADP/TBTNB and DADP/TCTNB in solvents and in vacuum Table 3 The attachment energies for morphologically important faces of DADP/TBTNB and DADP/TCTNB DADP/TBTNB (h k l)

DADP/TCTNB % Total facet area

(h k l)

% Total facet area

(1 1 0)

−53.52

59.64

(1 1 −1)

−58.08

36.25

(1 1 −1)

−120.31

20.40

(1 1 0)

−61.50

35.91

(0 2 0)

−56.18

13.03

(0 2 0)

−53.79

16.82

(2 0 0)

−71.90

6.92

(1 1 1)

−87.19

11.02

Remark:

is in kcal/mol.

Table 4 The modified attachment energies between morphologically important faces of DADP/TBTNB and solvents Acetonitrile

Cyclohexanone %Total facet area

(h k l)

Heptane

%Total facet area

%Total facet area

(1 1 0)

−24.58

−340.76

18.65

−25.39

−331.27

0.01

−31.29

−261.80

65.32

(1 1 −1)

11.73

−554.95

69.60

3.03

−518.38

95.03

−49.62

−297.13

27.21

(0 2 0)

−30.34

−262.22

0.31

−29.05

−275.33

0

−35.47

−210.24

7.48

(2 0 0)

−26.80

−676.47

11.44

−28.97

−643.83

4.96

−51.70

−302.91

0

Remark: Both

and

are in kcal/mol.

9

J.-W.Li et al / Manuscript

Table 5 The modified attachment energies between morphologically important faces of DADP/TCTNB and solvents Acetonitrile

Cyclohexanone

%Total facet area

(h k l)

Heptane

%Total facet area

%Total facet area

(1 1 −1)

−20.95

−346.62

24.40

−19.33

−361.74

17.69

−30.24

−259.85

31.93

(1 1 0)

−22.05

−389.16

15.70

−25.73

−352.93

0

−35.19

−259.61

15.71

(0 2 0)

−22.64

−307.52

5.64

−24.47

−289.48

0

−30.75

−227.47

13.11

(1 1 1)

13.69

−641.34

54.26

6.81

−597.63

82.31

−26.71

−384.54

39.25

Remark: Both

and

are in kcal/mol.

Fig. 4 The simulation and experiment of DADP/TBTNB cocrystal morphology in cyclohexanone

(A) DADP/TBTNB-cyclohexanone system From Table 3 and Table 4, the facet area order of DADP/TBTNB in vacuum is (1 1 0) (59.64%) > (1 1 –1) (20.40%) > (0 2 0) (13.03%) > (2 0 0) (6.92%), while the order of DADP/TBTNB faces in cyclohexanone is (1 1 –1) (95.03%) > (2 0 0) (4.96%) > (1 1 0) (0.01%) > (0 2 0) (about 0). It can be found that the facet area of (1 1 –1) tends to be larger in cyclohexanone, while the facet areas of (0 2 0), (2 0 0) and (1 1 0) tend to sol be smaller. Because Rhkl and

mod  att are positive relevant, the smaller Ehkl

mod  att Ehkl

may result in slower surface growth and larger surface areas. From Table 4, the order mod  att of Ehkl in DADP/TBTNB-cyclohexanone is (1 1 –1) (3.03 kcal/mol) < (1 1 0)

(–25.39 kcal/mol) < (2 0 0) (–28.97 kcal/mol) < (0 2 0) (–29.05 kcal/mol). Therefore, the facet area of (1 1 –1) tends to be the largest and (0 2 0) tends to be the smallest, which can be verified by the facet area order of DADP/TBTNB cocrystal important faces in cyclohexanone and vacuum. Because the difference of facet area between (1 10

J.-W.Li et al / Manuscript

1 –1) and (2 0 0) is too large and (0 2 0) and (1 1 0) almost disappear (faces areas close to 0), it formed a needle-like cocrystal morphology in simulation, in agreement with the experimental morphology of DADP/TBTNB in cyclohexanone by Landenberger et al [13]. The above discussions also prove that MAE model is a proper way to simulate the cocrystal morphology. (B) DADP/TBTNB- and DADP/TCTNB-other solvents system From Table 3 and Table 4, the facet area order of DADP/TBTNB in acetonitrile is (1 1 –1) (69.60%) > (1 1 0) (18.65%) > (2 0 0) (11.44%) > (0 2 0) (0.31%). Comparing with the facet area in vacuum, (1 1 –1) and (2 0 0) are larger while (1 1 0) mod  att Ehkl

and (0 2 0) become smaller. For the

in DADP/TBTNB-acetonitrile, (0 2 0)

is the largest while (1 1 –1) is the smallest, resulting in the change tendency of facet area. The similar conclusion can be found in other solvents, such as (1 1 0) and (2 0 0) of DADP/TBTNB in heptane, whose

mod  att are the smallest and largest, tending Ehkl

to be larger and smaller, respectively. It’s also workable for the changes of DADP/TCTNB morphologies in solvents. For the (1 1 1) and (0 2 0) in acetonitrile, they have the smallest and largest

mod  att and will be large and small, respectively. Ehkl

For (1 1 1) and (1 1 0) in cyclohexanone, they will be large and small with the smallest and largest

mod  att , respectively. The (1 1 1) and (1 1 0) in heptane head to Ehkl

be large and small, because they have the smallest and largest

mod  att . It should be Ehkl

noted that the above discussions are based on MD simulation and may guide for the future experiments. According to the solvent effect, the cocrystal morphologies grew from different solvents are different. Moreover, the cocrystal morphologies of DADP/TBTNB and DADP/TCTNB in the same solvent are similar. Theoretically, crystal growth is that crystal molecular stacked layer by layer [30]. If crystal face is favorable for solvent adsorption, the process of stacking crystal molecular layers will be difficult because it needs more energy to get rid of the adsorption of solvent molecules, and the growth of crystal will be inhibited along certain direction and form a large crystal face and thus change the morphology. From Table 4, the absolute value 11

J.-W.Li et al / Manuscript

of interaction energy on (0 2 0) of DADP/TBTNB is the smallest in three solvents, so it may be inferred that the adsorption affinity on (0 2 0) face is the weakest, and the (0 2 0) has a trend to be smaller. On the contrary, the absolute value of interaction energy on (1 1 −1) face is the largest, so the adsorption affinity on (1 1 −1) face is the strongest, and the (1 1 −1) tends to be larger. The growth of DADP/TCTNB cocrystal face is affected by the interaction and the adsorption affinity as well. 3.2 Molecular Electrostatic Potentials Analysis The molecular electrostatic potential (ESP) refers to the energy of moving a unit positive charge from infinity to a certain position in a chemical system. It shows the electrophilicity and nucleophilicity at different positions in the system, and also has contributed to reveal the properties of non-covalent molecular interaction [31, 32]. The ESP in different region of a molecule structure is shown as negative or positive, which means the region is expressed as electrophilicity or nucleophilicity. If the polar of ESP on solvent and crystal face are the same, there will be repulsive interaction between that solvent and crystal face. Otherwise, there will be attractive interaction between solvent and crystal face. It would be qualitatively expressed by the coloring isosurface graphs of electrostatic potentials, which visualized the charged areas of molecules. The coloring isosurface graphs of the molecular electrostatic potentials on cocrystal morphologically important faces of DADP/TBTNB and DADP/TCTNB are illustrated in Fig. 5, and the coloring isosurface graphs of solvents molecules are shown in Fig. 6. In these coloring isosurface graphs, the positive and negative potential regions are depicted by the different colors. Dark blue is symbolized of the most positive potential and dark red is symbolized of the most negative potential. For the morphologically important faces, the positive potentials are mostly located at the basins, it can be considered as nucleophilic in the basins. For the solvent molecules, they selected the optimal orientations to match the cocrystal faces’ potential distribution in accordance with the theory of electrostatic attractions, which inhibit or facilitate the growth of cocrystal faces.

12

J.-W.Li et al / Manuscript

Fig. 5 The coloring isosurface graphs of the molecular electrostatic potentials for the morphologically important faces of DADP/TBTNB and DADP/ TCTNB (Red, white, blue, green, rusty red and grey small ball represent oxygen, hydrogen, nitrogen, chlorine, bromine and carbon atom, respectively)

Fig. 6 The coloring isosurface graphs of the molecular electrostatic potentials for the acetontrile, cyclohexanone and heptane molecule

As shown in Fig. 7, for instance, it can be found that a certain amount of acetonitrile molecules oriented with the positions of negative potential to the basins of DADP/TBTNB’ s cocrystal faces. The orientations of acetonitrile molecules are different between different cocrystal faces. By roughly statistic, the number of acetonitrile molecules oriented with negative potential to the basin of (1 1 −1) face is 13

J.-W.Li et al / Manuscript

the most (more than half), but on (0 2 0) face is the least (less than half). It means that there are more attractive interactions between acetonitrile molecules and (1 1 −1) face than repulsive interactions due to the orientation and other factors. It will inhibit the growth of (1 1 −1) face and so acquires a larger face. For (0 2 0) face, repulsive interactions between acetonitrile molecules and (0 2 0) face are more than that of the attractive interactions due to the orientation of acetonitrile, so there is repulsive interaction between acetonitrile and (0 2 0) face. It will facilitate the growth of (0 2 0) face and gets a small face. The orientations of cyclohexanone on (1 1 −1) face and (0 2 0) face are similar to the orientations of acetonitrile. But the condition of the heptane orientations is different because almost all the adsorption sites on heptane molecule are the positive potentials that show as nucleophilicity, which affects the electrostatic attraction between heptane and basins to a large extent.

Fig. 7 The molecular arrangement at the interfaces between morphological important faces of DADP/TBTNB and acetonitrile

The reason for why the arrangement of solvent molecules orientations is different on different faces can be explained by structural features of cocrystal faces. Liu et al. hold that the flatness and roughness of crystal face are of vital importance to the solvent adsorption because the adsorption sites would form the cavities which born 14

J.-W.Li et al / Manuscript

from the ‘basin’ [33]. It can be characterized by the parameter S, the definition of S is: (3.2) Among the definition, the value of S represents the roughness of crystal face, the solvent accessible area of crystal face, and

is

is the corresponding

cross-sectional area of crystal face. The values of S of DADP/TBTNB and DADP/TCTNB are showed in Table 6. Table 6 The roughness S of DADP/TBTNB and DADP/TCTNB DADP/TBTNB

DADP/TCTNB

(h k l)

(0 2 0)

(1 1 −1)

(1 1 0)

(2 0 0)

(0 2 0)

(1 1 −1)

(1 1 0)

(1 1 1)

S

1.18

1.98

1.52

1.19

1.21

1.28

1.61

1.81

From Table 6, it may be inferred that in DADP/TBTNB the solvent adsorption is more facilitative at the (1 1 −1) face than the (0 2 0) face, and in DADP/TCTNB the solvent adsorption is more facilitative at the (1 1 1) face than the (0 2 0) face. The solvent accessible surface of DADP/TBTNB and DADP/TCTNB are shown in Fig. 8.

Fig. 8 The solvent accessible arrangement on morphologically important faces of DADP/TBTNB and DADP/TCTNB

In Fig. 8, it can be found that the (0 2 0) face of DADP/TBTNB is the most flat surface and has the least number of cavities, so it possesses the least number of available adsorption sites. This would give rise to the least solvent molecules fill in the cavities of (0 2 0) face and finally get the weakest adsorption affinity. For the (1 1 15

J.-W.Li et al / Manuscript

−1) face, it is the most rough surface and has the most number of cavities, so it holds the most number of available adsorption sites. This would cause the most solvent molecules set in the cavities and get the strongest adsorption affinity at last. It would explicates why the orientations of solvent molecules at (0 2 0) and (1 1 −1) face are completely different. From Table 6, it shows that the value S of (0 2 0) face from DADP/TCTNB, 1.21, is close to the value of (1 1 −1) face from DADP/TCTNB, 1.28. It means that they are in similar roughness. But from Table 5, it is found that their growth rates are quite different, the growth rate of (0 2 0) face is larger than (1 1 −1) face. In Fig. 8, it can be known that although the roughness of (0 2 0) face and (1 1 −1) face are very close, the number of cavities is different. Obviously, the number of cavities on (1 1 −1) face is more than (0 2 0) face, which increases the adsorption affinity of (1 1 −1) face to a certain extent and relatively inhibits the growth of (1 1 −1) face. 3.3 Mass Density Distribution After the qualitative analyses by ESP and the structural features of DADP/TBTNB and DADP/TCTNB, quantitative analyses of the interaction between the crystal face and the solvents are needed. The mass density distribution of solvent molecules, the solvent molecules which along the normal of crystal face, is influenced by the interaction between solvent and crystal face and the plot of solvent mass density distribution can be seen in Fig. 9. In Fig. 9, the abscissa z is the value of solvent-interface model, and the ordinate represents the mass density distribution at the corresponding z positions. Taking DADP/TBTNB as an example, after comparing the mass density distribution of the three solvent molecules, it is found that the values of the first prominent peak on the (1 1 0) face are always the lowest and most of them are between 0.8 g/cm−3and 0.9 g/cm−3. The value of the first prominent peak on the (0 2 0) face is extremely close to the (2 0 0) face, and (0 2 0) face is nearly always the highest value among these morphologically important faces in the solvents. According to Fig. 9, there are several small peaks before the first prominent peak on the (1 1 −1) face of DADP/TBTNB, it can be inferred that (1 1 −1) face is relatively rough.

16

J.-W.Li et al / Manuscript

Fig. 9 The functions of mass density distribution of solvent molecules on the morphologically important faces of DADP/TBTNB and DADP/TCTNB

From the mass density distribution of (0 2 0) face and (2 0 0) face in DADP/TBTNB, it may be supposed that the densities of solvent molecules around the (0 2 0) face and (2 0 0) face are larger than (1 1 −1) face and (1 1 0) face, and the interactions on (0 2 0) face and (2 0 0) face are also larger than that on (1 1 −1) face and (1 1 0) face. The solvent densities of (0 2 0), (2 0 0) and (1 1 −1) are more than 0.9 g·cm−3, it can be indicated that the interactions between solvent molecules and (0 17

J.-W.Li et al / Manuscript

2 0), (2 0 0) or (1 1 −1) are stronger than that of the (1 1 0) [34]. According to ESP analysis, there is repulsive interaction between the (0 2 0) interface and attractive interaction between the (1 1 −1) interface. So it can be inferred that the stronger interactions on (0 2 0) face and (1 1 −1) face would facilitate and inhibit the growth, and the size of (0 2 0) face would be the minimum and the size of (1 1 −1) face would be the maximum. It’s in conformity with the results of qualitative analyses. 3.4 Radial Distribution Function (RDF) Radial distribution function is the statistically average distribution of other atoms in the radial direction around a certain atom [35]. The abscissa r indicates the distance from a particle to the reference atom, and the ordinate g(r) indicates the probability of appearance for a particle which at a distance to reference atom. In this study, it would be used to explore the hydrogen bonding effect between the solvent and cocrystal molecules at the interface to explore the interactions between the solvent and cocrystal quantitatively. The oxygen atoms from the DADP/TBTNB cocrystal layer and the hydrogen atoms from the solvent molecules layer were picked out to analyze the H···O interaction, and the plots of radial distribution function about hydrogen and oxygen atoms are shown in Fig. 10. The abscissa of the first peak represents the strength of intermolecular interaction and the ordinate of the first peak stands for the probability of that intermolecular interaction. In general, if the first peak is high and pointed, the combination between the atoms is tight and the distribution of the atoms is regular. For the abscissa of first peak, the numerical range from 2.0 Å to 3.1 Å could be regarded as hydrogen bonding. If the values are greater than 3.1 Å, they are called van der Waals force [36]. In Fig. 10, all the values on abscissa at first peak of morphological important faces are less than 3.0 Å, it can be concluded that all of them are considered to be hydrogen bond. The values on abscissa are completely similar, so the locations of H···O bonds at the interface of these morphologically important faces are approximately similar. For the DADP/TBTNB-acetonitrile interface, the values of ordinate at first peak of (0 2 0), (1 1 −1), (1 1 0) and (2 0 0) are 0.7, 1.8, 1.0 and 1.9, respectively. It can be concerned that the probability of H···O bonding also represents 18

J.-W.Li et al / Manuscript

for its strength. So the strength of H···O bonding at the interface of acetonitrile and DADP/TBTNB can be considered as (0 2 0) < (1 1 0) < (1 1 −1) < (2 0 0). The interaction energies between acetonitrile and (0 2 0), (1 1 −1), (1 1 0) and (2 0 0) are −262.22 kcal/mol, −554.95 kcal/mol, −340.76 kcal/mol and −676.47 kcal/mol, respectively. The order of their absolute values is consistent with the conclusion of RDF. According to the analysis of ESP and mass density distribution, it is inferred that the growth rate of (0 2 0) is the greatest and the (1 1 −1) is the least, so the size of (0 2 0) is the minimum and the size of (1 1 −1) is the maximum. Although the interaction on the (2 0 0) face is the strongest, the interaction does not express repulsion or attraction obviously so the effect on the face growth is not significant. For the DADP/TBTNB-cyclohexanone interface, the plot of RDF at the first peak on (1 1 −1) face is most sharp. It can be inferred that the H···O bonding on (1 1 −1) face is the strongest and verified by the growth of (1 1 −1) face. In general, the cocrystal growth could be explicated by the RDF analyses on the whole.

19

J.-W.Li et al / Manuscript

Fig. 10 The functions of radial distribution of H···O distance at the interface between cocrystal and solvent molecules

4. Conclusion In summary, the simulations and analyses confirmed that DADP/TBTNB and DADP/TCTNB cocrystal growth and morphology could be affected by the strength and affinity of interactions between cocrystal surface and solvent molecules. By analyzing the ESP and mass density distribution on different region of the solvent and cocrystal molecules at the interface, it can be found that the interactions between 20

J.-W.Li et al / Manuscript

solvent and (0 2 0) face and (1 1 −1) face show repulsion and attraction because of the electrophilicity and nucleophilicity at different molecular structure regions, and the growth of (0 2 0) face and (1 1 −1) face would be facilitated and inhibited and finally formed a small and a large face. The orientations of the solvent molecules, which could determine whether attraction or repulsion at the interface, are affected by the number and the size of the cavities on cocrystal face. The strength of interactions can be known by the mass density distribution of solvent molecules, and the strength of attractive interaction can be obtained from the radial distribution function. They can judge the strength and property of interactions and determine the size of grown face. The conclusions are helpful to select the solvents and additives for acquiring desired crystal morphologies, such as crystal morphologies spheroidization of energetic materials to reduce their sensitivity by using specific solvents or additives, or designing special crystal morphologies to meet particular needs for performances. It’s still not clear the accurate distribution of the cavities on cocrystal face. Perhaps it’s the direction of the next step. We will try to control the growth of crystals by controlling the properties and orientation of solvent molecules.

5. Acknowledgement We gratefully acknowledge the North University of China for financial support. Deep thanks to Prof. Zhang, Prof. Gou and Dr. Han’s help and support in this research. Sincerely thanks the correcting and helping by anonymous reviewers.

6. Statement This paper only represents the authors’ view. If there’s different opinion or a better sound to the content, hope to contact us for communication.

7. References [1] U. Teipel, Energetic materials: particle processing and characterization, WILEY-VCH, Weinheim, 2005. [2] J.P. Agrawal, High energy materials: propellants, explosives and pyrotechnics, WILEY-VCH, Weinheim, 2010. [3] S.R. Anderson, P. Dubé, M. Krawiec, J.S. Salan, D.J. Ende, P. Samuels, Promising 21

J.-W.Li et al / Manuscript

CL-20-based energetic material by cocrystallization, Propell. Explos. Pyrot. 41 (2016) 783-788. [4] L.J. Powell, Insensitive munitions-design principles and technology developments, Propell. Explos. Pyrot. 41 (2016) 409-413. [5] O. Bolton, A.J. Matzger, Improved stability and smart-material functionality realized in an energetic cocrystal, Angew. Chem. 123 (2011) 9122-9125. [6] M.J. Siegfried, K.S. Choi, Elucidating the effect of additives on the growth and stability of Cu2O surfaces via shape transformation of pre-grown crystals, J. Am. Chem. Soc. 128 (2006) 10356-10357. [7] J.X. Chen, J.K. Wang, J. Ulrich, Q.X. Yin, L.Z. Xue, Effect of solvent on the crystal structure and habit of hydrocortisone, Cryst. Growth Des. 8 (2008) 1490-1494. [8] M.N. Bhat, S.M. Dharmaprakash, Effect of solvents on the growth morphology and physical characteristics of nonlinear optical g-glycine crystals, J. Cryst. Growth 242 (2002) 245-252. [9] Q.L. Zhao, N. Liu, B.Z. Wang, W.L. Wang, A study of solvent selectivity on the crystal morphology of FOX-7 via a modified attachment energy model, Rsc Adv. 6 (2016) 59784-59793. [10] P. Bennema, Theory of growth and morphology applied to organic crystals; possible applications to protein crystals, J. Cryst. Growth 122 (1992) 110-119. [11] M. Lahav, L. Leiserowitz, The effect of solvent on crystal growth and morphology, Chem. Eng. Sci. 56 (2001) 2245-2253. [12] G. Han, Q.F. Li, R.J. Gou, S.H. Zhang, F.D. Ren, L. Wang, R. Guan, Growth morphology of CL-20/HMX cocrystal explosive: insights from solvent behavior under different temperatures, J. Mol. Model. 23 (2017) 360. [13] K.B. Landenberger, O. Bolton, A.J. Matzger, Two isostructural explosive cocrystals with significantly different thermodynamic stabilities, Angew. Chem. Int. Edit. 52 (2013) 1-5. [14] H.F. Gao, S.H. Zhang, F.D. Ren, R.J. Gou, C.L. Wu, Theoretical insight into the temperature-dependent acetonitrile (ACN) solvent effect on the diacetone diperoxide (DADP)/1,3,5-tribromo-2,4,6-trinitrobenzene (TBTNB)

cocrystallization,

Comp. 22

J.-W.Li et al / Manuscript

Mater. Sci. 121 (2016) 232-239. [15] J.J. Xiao, W.H. Zhu, W. Zhu, H.M. Xiao, Molecular dynamics of energetic materials (in chinese), Beijing: Science Press, Beijing, 2013. [16] N. Liu, Y.N. Li, S. Zeman, Y.J. Shu, B.Z. Wang, Y.S. Zhou, Q.L. Zhao, L.W. Wang, Crystal morphology of 3,4-bis (3-nitrofurazan-4-yl) furoxan (DNTF) in a solvent system: molecular dynamics simulation and sensitivity study, Crystengcomm 18 (2016) 2843-2851. [17] X.H. Duan, C.X. Wei, Y.G. Liu, C.H. Pei, A molecular dynamics simulation of solvent effects on the crystal morphology of HMX, J. Hazard. Mater. 174 (2010) 175-180. [18] P. Hartman, P. Bennema, The attachment energy as a habit controlling factor: I. Theoretical considerations, J. Cryst. Growth 49 (1980) 145-156. [19] M. Brunsteiner, S.L. Price, Morphologies of organic crystals:  sensitivity of attachment energy predictions to the model intermolecular potential, Cryst. Growth Des. 1 (2001) 447-453. [20] P. Hartman, W.G. Perdok, On the relations between structure and morphology of crystals. Ⅰ, Acta Crystallogr. 8 (1955) 49-52. [21] Z. Berkovitch-Yellin, Toward an ab initio derivation of crystal morphology, J. Am. Chem. Soc. 107 (1985) 8239-8253. [22] Accelrys Software Inc., Materials Studio Release Notes, Release 7.0, Accelrys Software Inc., San Diego, 2013. [23] H. Sun, COMPASS:  An ab Initio Force-Field Optimized for Condensed-Phase Applications

Overview with Details on Alkane and Benzene Compounds, J. Phys.

Chem. B 102 (1998) 7338-7364. [24] G. Chen, C.Y. Chen, M.Z. Xia, W. Lei, F.Y. Wang, X.D. Gong, A study of the solvent effect on the crystal morphology of hexogen by means of molecular dynamics simulations, Rsc Adv. 5 (2015) 25581-25589. [25] G. Chen, M.Z. Xia, W. Lei, F.Y. Wang, A study of the solvent effect on the morphology of RDX crystal by molecular modeling method, J. Mol. Model. 19 (2013) 5397–5406. 23

J.-W.Li et al / Manuscript

[26] H.C. Andersen, Molecular dynamics simulations at constant pressure and/or temperature, J. Chem. Phys. 72 (1980) 2384-2393. [27] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, O. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision A.01, Gaussian Inc, Wallingford, CT, 2009. [28] A.D. Becke, Density ‐ functional thermochemistry. III. The role of exact exchange, J. Chem. Phys. 98 (1993) 5648-5652. [29] C. Lee, W. Yang, R.G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 35 (1988) 785-789. [30] J.W. Mullin, Crystallization, Oxford, Butterworth-Heinemann, 2001 [31] T.P. Chiu, S. Rao, R.S. Mann, B. Honig, R. Rohs, Genome-wide prediction of minor-groove electrostatic potential enables biophysical modeling of protein-DNA binding, Nucleic Acids Res. 45 (2017) 12565-12576. [32] J.A. McCaughan, D.M. Turner, R.K. Battle, Electrostatic Potential Change in a Paired Epitope: A novel explanation for Bw4 antibodies in patients with B13 (Bw4) antigens, Transplantation 100 (2016) E32–E34. [33] Y.Z. Liu, W.P. Lai, Y.D. Ma, T. Yu, Y. Kang, Z.X. Ge, Face-Dependent Solvent Adsorption: A Comparative Study on the Interfaces of HMX Crystal with Three Solvents, J. Phys. Chem. B 121 (2017) 7140-7146. [34] Y.Z. Liu, W.P. Lai, T. Yu, Y.D. Ma, Y. Kang, Z.X. Ge, Understanding the growth 24

J.-W.Li et al / Manuscript

morphology of explosive crystals in solution: insights from solvent behavior at the crystal surface, Rsc. Adv. 7 (2017) 1305-1312. [35] R.V. Vaz, J.R.B. Gomes, C.M. Silva, Molecular dynamics simulation of diffusion coefficients and structural properties of ketones in supercritical CO2 at infinite dilution, J. Supercrit. Fluid. 107 (2016) 630–638.

25

J.-W.Li et al / Manuscript

Highlights:

1. Electrostatic potentials were used to analyze the cocrystal-solvent interaction. 2. Cocrystal surface growth is influenced by the strength & affinity of interaction. 3. Repulsive interaction will facilitate the surface growth and then form a small face. 4. Attractive interaction will inhibit the surface growth and then form a large face.

26

J.-W.Li et al / Manuscript

Graphical abstract

27