Molecular Modeling on the Role of Local ... - ACS Publications

Article pubs.acs.org/crystal

Molecular Modeling on the Role of Local Concentration in the Crystallization of L‑Methionine from Aqueous Solution Hong-Min Shim,† Allan S. Myerson,‡ and Kee-Kahb Koo*,† †

Department of Chemical and Biomolecular Engineering, Sogang University, Seoul 04107, Korea Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States

‡

S Supporting Information *

ABSTRACT: In the present work, the importance of local concentration of the growth of L-methionine (L-Met) was explored by the spiral growth model and interfacial structure (IS) analysis. The spiral growth model shows that the decrease in the local concentration of growth units at crystal faces plays a key role in an appearance of thin hexagonal plate-like morphology of L-Met crystals from water. The (001)_2 layer with hydrogen bonds on its under surface was found to require a relatively large amount of energy for detachment of growth units from kink sites. This situation results in lower local concentration and leads to the decrease in net flux of solute molecules at kink sites compared with the (001)_1 layer interacting with its under surface by van der Waals force. Furthermore, the interfacial structure analysis provides important information on the growth of L-Met. Hydrophilic NH3+ and COO− groups of L-Met on the (001)_1 surface were found to have a tendency to form hydrogen bonding with water, and thus L-Met molecules should overcome the energy barrier required to detach water molecules adsorbed on the crystal surface for continuous growth of the (001) face. From those results, it can be concluded that the limitation of approaching growth units for the formation of the (001)_2 layer onto the (001)_1 surface is responsible for very slow growth of the [001] direction of L-Met crystals.

1. INTRODUCTION

emphasized enough that a mechanistic understanding of crystal growth should be considered. L-Met is comprised of non-centrosymmetric growth units that do not have a center of inversion. It accordingly brings on different kink rates with respect to the edges on the same crystal face.2 Moreover, the complex periodic bond chains (PBCs) exacerbate the approximation of kink density. The mechanistic spiral growth model able to deal with non-centrosymmetric growth units was reported by Kuvadia and Doherty2 and Kim et al.11 They highlighted the kink rate model for the step velocity in solution growth along with revealing that the kink rate could be determined from the required work for kink detachment. However, the importance of anisotropic local concentration on the growth shape of crystals has not been discussed in depth. In effect, the local concentration plays a key role in the crystal growth because the step velocity is affected by not only kink density but also local concentration considerably.12,13 In this paper, we aim at revealing the reason behind a hexagonal plate shape of L-Met crystals in water at the atomic level. In particular, we focus on the local concentration term contained in the kink rate expression. The (001) face, morphologically the important face of L-Met, is characterized by the existence of double layers in an interplanar distance. It

The morphology of crystals is mainly determined by the relative growth rate of crystal faces that is governed by the kinetics of surface integration processes. In general, when a large amount of supersaturation is consumed by nucleation, crystals can grow still at very low supersaturation where an infinite source of steps is favorably provided by screw dislocations on crystal face.1 In this respect, the spiral growth model is regarded to be useful for dealing with crystal growth at the atomic level, and the prediction of growth shape for numerous organic crystals has been successfully achieved.2−6 L-Methionine (L-Met), which is an initiating amino acid in the synthesis of most eukaryotic proteins, is a principal sulfurcontaining amino acid.7 L-Met crystals are usually produced as its salt form with HCl by reaction crystallization because of its lower solubility, and they grow as thin hexagonal plates and have a tendency to form agglomerates with ease.8 The morphology of LMet crystals influences the bulk density,9 which is one of the most important physical properties of amino acid crystals because products with higher bulk density reduce transportation load. However, there has been no report on the growth habit of LMet until a recent date. The change in crystal shape occurs unequivocally owing to the altered values of anisotropic terms such as kink density and kink rate which strongly affect the rate of movement of layers on the crystal face.10 Therefore, for the rational control of crystal shape of L-Met, it cannot be © XXXX American Chemical Society

Received: March 15, 2016 Revised: April 12, 2016

A

DOI: 10.1021/acs.cgd.6b00418 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

gives rise to the emergence of a rate-determining layer with much lower local concentration that reduces the flux of solute molecules into kink sites. In the present work, the local concentration determined by the spiral growth model is compared with the result obtained by the interfacial structure (IS) analysis,14 which may correlate the relationship between the local concentration and the growth shape of L-Met. Furthermore, a refined model dealing with the complex PBC networks is proposed.

Table 1. Bond Energies in the Crystal Graph of L-Met bond direction a b c d e f g h i

2. COMPUTATIONAL METHODS 2.1. Periodic Bond Chain (PBC) Analysis. The molecular structure of L-Met is given in Figure 1. The geometry of an L-Met

no.

dispersive energy, Ed (kcal/mol)

Coulombic energy, Ec (kcal/mol)

total energy, Etot (kcal/mol)

1 1 1 1 1 1 1 1 1

5.62 5.26 6.99 8.35 5.20 8.69 5.20 8.01 9.06

0.07 −0.04 0.56 −0.56 −0.92 −0.11 −1.00 −0.58 −0.23

−13.23 −12.03 −11.85 −7.81 −2.85 −3.61 −1.83 −0.37 −3.38

−13.16 −12.07 −11.29 −8.37 −3.77 −3.72 −2.83 −0.95 −3.61

where ΔHsub is the enthalpy of sublimation, R is the gas constant, T is the temperature, and ΔEpt is the proton transfer energy. The experimental ΔHsub at 298.15 K is 29.89 kcal/mol.20 However, the ΔEpt cannot be determined experimentally and has not yet been revealed theoretically for L-Met. It has been reported that the ΔEpt determined by the density functional theory (DFT) for some amino acids lies in the range of −34.32 and −21.32 kcal/mol.19 The term of 2RT indicating the translational and rotational contributions yields 1.18 kcal/mol at 298.15 K. In this respect, the Elatt for L-Met could be cautiously presumed to be between −65.39 and −52.39 kcal/mol, and thus it seems that the calculated value of −59.77 kcal/mol in the present work is quite reliable. Figure 3 shows the molecular packing geometry for the (001), (100), (110), and (11̅0) faces of L-Met that are confirmed as F (flat) faces containing two or more PBCs in one layer. The (001) face is shown to have double layers within an interplanar distance of d(001). The (001)_1 and (001)_2 layers have the same bonding structures with each other as shown in Figure 4a, and they are characterized by the existence of hydrogen bonding with the under layer. The (001)_1 layer is found to interact with its under layer extensively by van der Waals force stemming from methyl (CH3) group of L-Met, whereas the (001)_2 layer interacts strongly with its under layer by hydrogen bonding between ammonium (NH3+) and carboxylate (COO−) groups. Likewise, the (100) face has double layers in an interplanar distance with providing the same bonding structure as shown in Figure 4b. The (110) and (11̅0) faces are symmetric with respect to the x-axis, and four different growth units exist in one layer. 2.2. Molecular Dynamics (MD) Simulation. The zwitterionic model for L-Met was chosen in the present work. The general AMBER force field was used to describe the intra- and intermolecular interactions of L-Met and water molecules. For a water molecule, SPC/E model was employed.21 The MD simulations were carried out using Gromacs (version 4.6.5)22−25 to investigate molecular behavior of L-Met near the crystal faces. The cutoff distances were 9 Å for Coulombic interactions and 14 Å for Lennard−Jones interactions. The time step was set at 2 fs, and periodic boundary conditions were applied in three directions. To improve the performance of calculation on long-range electrostatic interactions, the particle-mesh Ewald (PME) method was employed.26 Temperature and pressure coupling were done by using V-rescale27 and Berendsen,28 respectively. To build a crystal surface of L-Met exposed to an aqueous solution, the unit cell was cleaved corresponding to the Miller indices, and then the bulk crystal was created by producing a supercell. Every supercell was reorientated as a xy plane of the simulation box where the z direction is normal to the surface. The bulk solution was initially prepared by random insertion of L-Met and water molecules at XA = 0.0244. The constructed simulation system was geometrically optimized, and an isothermal−isobaric (NPT) simulation was carried out during 1 ns at the pressure of 1.013 bar and the temperature of 363.15 K. In these conditions, L-Met molecules are hydrated in aqueous solution and phase transition does not occur. The supersaturated bulk solution was added to both upper and lower regions of the simulation box of which a bulk crystal is in the middle. To induce crystal growth, an equilibrated solution at 363.15 K was quenched to the temperature of 303.15 K. According to the solubility data of L-Met in aqueous solution,29 the

Figure 1. Molecular structure of L-Met. molecule was optimized by the MNDO semiempirical method,15 and the atomic partial charges were obtained from the Mulliken population by using the INDO1 wave functions16 (Table S1). The crystal structure of L-Met is a monoclinic with space group P21 and Z = 4 in the unit cell (a = 9.493 Å, b = 5.201 Å, c = 14.831 Å, and β = 99.84°).17 The bond energies in the crystal graph of L-Met (Figure 2) are summarized in

Figure 2. Crystal graph of L-Met. Different colored lines represent different bonds (red, a; orange, b; cyan, c; green, d; blue, e; black, f; pink, g; yellow-green, h; brown, i).

Table 1. In this table, the bond energy between growth units was calculated by the MORPHOLOGY (Material Studio, version 7.0) where the dispersive energy between molecules was obtained by employing the general AMBER force field.18 The Elatt for a crystal comprising of zwitterionic molecules is given by19

E latt = −ΔHsub − 2RT + ΔEpt

[1̅00] [01̅0] [000] [01̅0] [01̅0] [020̅ ] [01̅0] [001̅] [010]

bond distance (Å)

(1) B



Article

Figure 3. Crystal packing structure for the (a) (001), (b) (100), (c) (110), and (d) (11̅0) faces of L-Met.

Figure 4. Bonding structure of the F-faces of L-Met: (a) (001)_1 (up) and (001)_2 (down) layers, (b) (100)_1 (up) and (100)_2 (down) layers, and (c) (110) (up) and (11̅0) (down) layers. relative supersaturation σ is 2.25, where σ is the supersaturation defined by (XA/X0A) − 1, XA is the concentration of solute in solution, and X0A is the equilibrium concentration (0.0075 at 303.15 K). The production run of MD simulations for the growth of L-Met on the crystal surface was carried out by an NPT ensemble at 303.15 K and 1.013 bar for 10 ns with an anisotropic pressure coupling in the z direction (Table 2). The

nL‑Met, cryst

nL‑Met, sol

nwat

(001)_1, (001)_2 (100)_1, (100)_2 (110)

600 640 640

160 192 160

6400 7680 6400

(2)

where h is the step height and τ is the characteristic rotation time required for one full turn of a spiral. For a convex N-sided spiral, the rotation time is given by30 N

Table 2. Details on the Simulation Conditions of L-Met layer

h τ

G=

τ=

∑

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

i=1

(3)

where lc,i+1 is the critical length of an edge i+1, αi,i+1 is the angle between edges i and i+1, and vi is the step velocity of an edge i. The critical length was calculated by the thermodynamic condition where the net free energy of the system becomes zero by adding a new row of growth units with length lc,i.

pressure in the x and y directions was set to 0 bar to keep dimension of a bulk crystal; thus, the fluctuation is allowed to the vertical direction. Initial and final configurations of L-Met molecules on each surface are shown in Figures S1−S3.

lc , i =

2ae , i⟨ϕi kink ⟩ Δμ

(4)

where ae,i is the distance between molecules along the edge i, ⟨ϕkink i ⟩ is the average value of kink energies, and Δμ is the chemical potential between the liquid and solid state defined as kBTln(σ+1). The step velocity is given by2

3. THEORETICAL BACKGROUNDS In the spiral growth model, the growth rate of a F-face, G, is expressed as C



Article

vi = ap , iuiρi

The kth kink rate, uk, is defined as2

(5)

uk = j+ Pk − jk−+ 1 Pk + 1

where ap,i is the distance the edge propagates by adding a monolayer into the edge i, ui is the kink rate, and ρi is the kink density that is the probability of finding a kink along the edge i. For centrosymmetric growth units, the kink density can be determined by using Boltzmann distribution. ρi =

2 exp( −ϕi

kink

1 + 2exp( −ϕi

Pk + 2

/kBT )

(6)

where is the required work for generation of a kink from a straight edge i, and kB is Boltzmann constant, respectively. ϕkink is i equal to − Etot,i/2 for a vacuum.3 For the non-centrosymmetric growth units, the kink density can be determined by counting all the different microstates residing along the edge.2 When two kinks exposing E1 and E1 or E2 and E2 are generated by removing the growth units from the straight edge, two kinks exposing E1 and E1 or E2 and E2 can be formed by deposition of themselves on that edge. In that case, all the possible ways of rearrangement are type I (1,1 → 1,1), II (1,1 → 2,2), III (2,2 → 1,1), and IV (2,2 → 2,2) where values of the kink kink energy expense per kink are denoted as ϕkink 1,1 , ϕ1,2 , ϕ2,1 , and kink ϕ2,2 , respectively. The kink density is given by

u1 = u 2 = ··· = un = u/n

⎛ ϕ kink ⎞ ⎛ ϕ kink ⎞ ⎛ ϕ kink ⎞⎫ ⎪ 1,2 ⎟ ⎜ − 2,1 ⎟ + exp⎜ − 2,2 ⎟⎬ + exp⎜⎜ − + exp ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ kBT ⎠ ⎝ kBT ⎠ ⎝ kBT ⎠⎪ ⎭ (7) kink kink where ϕkink 1,1 = −(E1 + E1)/4, ϕ1,2 = −(E1 + E2)/4, ϕ2,1 = −(E2 + kink E1)/4, and ϕ2,2 = −(E2 + E2)/4, respectively. The surface energy of two immiscible phases, γi, is given by3

Wad, i = 2

Wad, i =

|γd,cryst | i

⎛ γ crystγ solv ⎞ d, i d ⎟ 4⎜⎜ cryst solv ⎟ | γ | + γ ⎝ d, i ⎠ d

+

+2

|γc,cryst | i

|γc,cryst |γasolv i

⎛ γ crystγ solv ⎞ c, i a ⎟ 4⎜⎜ cryst solv ⎟ | γ | + γ ⎠ ⎝ c, i a

(15)

= (1 − xeq)v0 exp(− (ΔU + ΔWk)/kBT )

ϕssi

1 ss (ϕ + ϕi ff ) 2 i

(16)

ϕffi

where and are interaction energies between solid−solid and fluid−fluid phases, respectively; ϕssi = −γssi si and ϕffi = −γffi si, and the ϕsfi is a gain in surface energy by forming one solid−fluid contact; ϕsfi = −Wad,isi/2 where Wad,i can be estimated by the geometric mean or the harmonic mean between γssi and γffi , and si = (Vunit/Z)2/3 where Vunit is the volume of a L-Met unit cell and Z is the number of molecules in the unit cell. All the values of ϕi are summarized in Table 3. The quantity xeq can be obtained using the mass balance of fluxes at equilibrium (j+ = j−),

solv where γcryst tot,i is the total surface energy at the kink site i, γtot is the surface energy of a solvent, and Wad,i is the work of adhesion at the crystal−solvent interface. Practically, polarization and hydrogen bonding interactions are inseparable in common force fields where the hydrogen bonding is assumed to be contained in the Coulombic interactions. Therefore, Wad can be determined by mixing the Coulombic interactions of the crystal and the associative interactions of the solvent. Wad can be estimated by the geometric mean equation (eq 9)3,31 or the harmonic mean equation (eq 10).32,33

γc,cryst i

(14)

jk−,eq

ϕi = ϕisf −

(8)

|γd,cryst |γdsolv i

j+ = xeqv0 exp( −ΔU /kBT )

where xeq is defined as VmCeq, where Vm, and Ceq are the molar volume and molar concentration of solute molecules in the solution, respectively, v0 is the vibrational frequency of attachment and detachment trials, and ΔWk is the required energy for removing a solute molecule from the kth kink site that can be obtained by summing over all bonds connected to the m growth unit at the kth kink site, ΔWk = ∑ j (ϕk , i)j , where ϕi is the energy of formation of the ith solid−fluid bond, and ΔU the energy barrier for transition state of a solute molecule into a kink site that can be assumed to be constant regardless of crystal faces because the attachment barrier is equal to that of breaking of the solvation shell near the growth unit.11 Therefore, the kink rate can be reduced as u′ = u/[v0 exp(− ΔU/kBT)] leaving only anisotropic terms.5 In the present work, the ϕi can be estimated by the relation35

⎧ ⎛ ϕ kink ⎞ ⎪ 1,1 ⎟ ⎨1 + exp⎜⎜ − ⎟ ⎪ k BT ⎠ ⎝ ⎩

cryst solv γi = γtot, + γtot − Wad, i i

(13)

The attachment and detachment rates are given as11

⎧ ⎛ ϕ kink ⎞ ⎛ ϕ kink ⎞ ⎛ ϕ kink ⎞ ⎪ 1,1 ⎟ ⎜ − 1,2 ⎟ + exp⎜ − 2,1 ⎟ ρ = ⎨exp⎜⎜ − + exp ⎜ k BT ⎟ ⎜ kT ⎟ ⎪ ⎝ kBT ⎟⎠ B ⎠ ⎝ ⎠ ⎝ ⎩

γd,cryst i

(12)

where j is the attachment rate of growth units into kink sites that is independent of the type of kink, jk−+ 1 is the detachment rate of growth units from the (k+1)st kink site, and Pk and Pk+1 are the probabilities of kth and (k+1)st kink sites terminated along the n edge i, respectively, where ∑k = 1 Pk = 1. To maintain the crystal stoichiometry, the overall kink rate, u, defined as number of molecules added per unit time, has the relationship34

ϕkink i

⎛ ϕ kink ⎞⎫ 2,2 ⎟⎪ + exp⎜⎜ − ⎟⎬ ⎝ kBT ⎠⎪ ⎭

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

+

/kBT )

kink

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

(11)

xeq =

exp( −ΔWk /kBT ) 1 + exp( −ΔWk /kBT )

(17)

In the present work, xeq is obtained by substituting ΔWk for the arithmetic mean of them, ΔWk,avg ,because ΔWk varies with respect to kink sites for non-centrosymmetric growth units. In a supersaturated solution, xeq becomes xeq(σ+1).

(9)

(10)

4. CALCULATION OF GROWTH RATES 4.1. Kink Density. As shown in Figure 3a, the (001) faces have double layers within an interplanar distance. The bonding structures of the (001)_1 and (001)_2 layers in Figure 4a show

where the subscripts d, c, and a denote the dispersive, Coulombic, and associative components of the surface energies. 2 For water, γsolv = 7.5, γsolv = 63.5, and γsolv d a tot = 71.1 erg/cm , respectively.3 D



Article

Table 3. Interaction Energy between Solid−Fluida bond

ϕgi (kcal/mol)

ϕhi (kcal/mol)

a b c d e f g h i

2.15 1.62 1.83 0.60 0.07 0.05 0.15 0.62 0.01

3.55 2.90 3.02 1.14 0.13 0.10 0.29 1.02 0.03

kink energy of ϕkink was considered to be reasonable for g,g estimating the kink density on the [01̅0] edge. In Figure 3b, there are double layers in an interplanar distance of the (100) face, and the bonding structures of the (100)_1 and (100)_2 layers are the same as shown in Figure 4b. Two bonds of Ec and Eh reside along the [001] edge where the kink density can be determined by eq 7. Meanwhile, it should be noted that the advance of growth units on the [010] edge changes the bonds perpendicular to the chain. In that case, the kink density cannot be directly obtained by the classical Boltzmann formula. In the same manner of the [010̅ ] edge on the (001) face, four different types of rearrangement of growth units can be classified as shown in Figure 6; type I (e,e → g,g), II (g,g → e,e), III (e,e → e,e), and IV (g,g → g,g), respectively, in which the energy expense per kink, kink = −[(Ee + Eg) + nedge{(Ec + Ed + Ei) − Eh}]/4, ϕkink ϕe,g g,e = −[(Eg kink + Ee) + nedge{Eh − (Ec + Ed + Ei)}]/4, ϕe,e = −(2Ee + 2Eg)/4, and kink = −(2Ee + 2Eg + 2Ed + 2Ei)/4, respectively. The ϕg,g rearrangement by type I or II is found to change the bonds exposed to the [001] direction by thermal rearrangement. For kink is found to become larger as nedge increased, the type I, ϕe,g which means that the rearrangement by type I is energetically unfavorable and rarely occurs. However, for the type II, ϕkink g,e tends to decrease as nedge increases and it results in a negative value when even nedge = 1. It implies that the rearrangement by type II is energetically more favorable and it exceedingly takes place. Meanwhile, type III and IV have constant values of kink kink = 3.30 and ϕkink energy, ϕe,e g,g = 9.29 kcal/mol, respectively, among which a rather smaller kink energy was used for the calculation of kink density. In Figure 3c,d, the (110) and (11̅0) faces are symmetry with respect to the x-axis, and they have the same bonding structure (Figure 4c). For the (110) face, the kink density on the [001] edge where two bonds of Ec and Eh alternate can be determined by eq 7. However, the only two growth units are found to be connected with the chain consisting of Ea and Eb along the [11̅0] direction, but the others remained without any bonds to that direction. Possible rearrangements of growth units on the [11̅0] edge of the (110) face are classified as given in Figure 7. The rearrangement by type I (0,0 → b,b) results in the kink energy ϕkink 0,b = −{Eb + nedge(Ec − Eh)}/4, while the rearrangement by type II (b,b → 0,0) consumes the kink energy ϕkink b,0 = −{Eb + nedge(Eh − Ec)}/4. The rearrangement by type I rarely occurs because Ec is stronger than Eh; hence, the kink energy becomes larger as nedge increases. However, the rearrangement by type II vigorously takes place because the kink energy is a negative value at nedge > 2. Type III (0,0 → 0,0) and IV (b,b → b,b) are found to kink consume the same energy where ϕkink 0,0 = ϕb,b = −(2Eb)/4. In the present work, the kink density on the [1̅10] edge of the (110) layer can be determined by eq 7. The kink energies in vacuum, , and water, ϕkink ϕkink,cryst i i , were summarized in Table 4, where the is the kink area used for the estimation of γkink,cryst and γkink Akink i i i .

a

Superscript g and h denote the geometric mean and harmonic mean, respectively.

that the interaction in the lateral direction is identical for both layers where the kink density is equal to each other. For the (001) _1 layer, two different bonds, Ea and Eb, are found to exist along the [11̅0] edge, and the kink density can be determined by eq 7. However, the estimation of the kink density on the [010̅ ] edge is complicated because the energy perpendicular to the chain is not constant when growth units are rearranged on that edge. Figure 5a illustrates one of the possible rearrangements of growth units on the [010̅ ] edge of the (001)_1 layer, labeled as type I (g,g → e,e) that generates two kinks exposing Eg and Eg with accompanying two kinks exposing Ee and Ee. The energy expense for this thermal rearrangement per kink, ϕkink g,e , is given by − [(Eg + Ee) + nedge{(Eb + Ef) − Ea}]/4, where nedge is the number of growth units removed from the straight edge. In the same manner, the Figures 5b−d show the rearrangements of growth units, type II (e,e → g,g), III (g,g → g,g), and IV(e,e → e,e); the kink energy expense per kink, ϕe,g = −[(Ee + Eg) + nedge{Ea − (Eb + Ef)}]/4, ϕkink = −(2E + 2E )/4, and ϕkink g,g e g e,e = −(2Ee + 2Eg + 2Ef)/ 4, respectively. Types I and II are found to change the bonds perpendicular to the edge by rearrangement. The rearrangement by type I removes Ea from a straight edge and makes Eb and Ef, whereas the rearrangement by type II eliminates Eb and Ef from a straight edge and forms Ea. However, the rearrangement by type III or IV does not show any change in the bonds perpendicular to the edge direction. The rearrangement type I is energetically unfavorable because the kink energy becomes larger as nedge increases. However, the kink energy associated with the rearrangement by type II becomes smaller as nedge increases and even shows a negative value at nedge > 3, which means that the formation of kink sites is energetically more kink favorable. Meanwhile, the values of ϕkink g,g and ϕe,e are constant as 3.30 and 5.16 kcal/mol, respectively. Fundamentally, thermal rearrangement is achieved by minimizing energy consumption associated with forming a new surface of a kink site. Therefore, the rearrangement by type I hardly occurs while the rearrangement by type II is aggressively accomplished; hence, they rarely affect the distribution of kink sites. For this reason, a rather lower

Figure 5. Rearrangements of growth units to form different types of kink sites on the [01̅0] edge of the (001)_1 layer of L-Met; (a) type I (g,g → e,e), (b) type II (e,e → g,g), (c) type III (g,g → g,g), and (d) type IV(e,e → e,e). E



Article

Figure 6. Rearrangements of growth units to form different types of kink sites on the [010] edge of the (100)_1 layer of L-Met; (a) type I (e,e → g,g), (b) type II (g,g → e,e), (c) type III (e,e → e,e), and (d) type IV(g,g → g,g).

Figure 7. Rearrangements of growth units to form different types of kink sites on the [1̅10] edge of the (110) layer of L-Met; (a) type I (0,0 → b,b), (b) type II (b,b → 0,0), (c) type III (0,0 → 0,0), and (d) type IV (b,b → b,b).

Table 4. Kink and Surface Energies for the (001)_1, (100)_1, and (110) layersa layer

edge

Akink (Å2) i

ϕkink,cryst (kcal/mol) i

γkink,cryst (erg/cm2) i

γkink,g (erg/cm2) i

ϕkink,g (kcal/mol) i

γkink,h (erg/cm2) i

ϕkink,h (kcal/mol) i

(001)_1

[01̅0] [1̅10] [010] [001]̅ [001] [1̅10]

79.06 37.99 69.33 24.31 24.57 67.33

3.30 6.58, 6.03 3.30 5.64, 0.47 5.64, 0.47 6.58, 6.03

28.97 120.25, 110.29 33.04 161.20, 13.56 159.52, 18.65 67.85, 62.23

11.93 20.50, 10.98 10.01 40.48, 32.27 39.81, 22.54 10.61, 5.33

1.36 1.12, 0.60 1.00 1.42, 1.13 1.41, 0.80 1.03, 0.52

22.10 27.18, 19.00 18.83 63.10, 49.41 61.87, 38.18 8.79, 7.02

2.59 1.49, 1.04 1.88 2.21, 1.73 2.19, 1.35 0.85, 0.68

(100)_1 (110) a

All values were obtained at T = 300 K. Superscript g and h denote the geometric mean and harmonic mean, respectively.

Figure 8. Progress of kink sites on the [01̅0] edge of the (001)_1 layer (a) and those on the [010] edge of the (100)_1 layer (b) when the thermal rearrangement occurs by type III, respectively.

4.2. Kink Rate. The progress of kink sites on the [01̅0] edge of the (001)_1 layer is illustrated in Figure 8a. When double layered kink sites are generated by thermal rearrangement, growth units of type 2 are located outermost on the edge. First, the growth unit migrates to a crystal face and is incorporated into kink site 1. Thereafter, two kink sites 1 and 2 are possible for the

incorporation of growth units. In that case, the event of attachment of a growth unit occurs mainly at kink site 2 (Pg1 = 0.42 and Pg2 = 0.58; Table S2) because a growth unit is strongly bonded at that site. However, a different progress of kink sites is found to take place along the [010] edge of the (100)_1 layer (Figure 8b). The F



Article

Figure 9. Thermal rearrangement and progress of kink sites for the [010] edge of the (100)_1 layer of L-Met, respectively.

Figure 10. Progress of kink sites for the [1̅10] edge of the (110) layer.

Figure 11. L-Met crystals grown from water by cooling (a) with agitation (cooling rate: 0.25 °C/min) and (b) without agitation.

because of their same kink energy. Once the double layered kink sites are generated by the rearrangement of type IV, the incorporation of a growth unit into the kink site 1 will be followed by the formation of two different kink sites 2 and 3. As given in Table S2, the probability at the kink site 2 is found to be higher than at the kink site 3; hence, a growth unit tends to be incorporated into the kink site 2 with retaining a stable edge. In the same manner, since the probability of a growth unit at the kink site 4 is larger than at the kink site 1, the edge is able to proceed with filling kink sites 2 and 4.

probability of a growth unit at the kink site 1 is much higher than that at the kink site 2 (Pg1 = 0.81 and Pg2 = 0.19; Table S2). As shown in Figure 6, the growth unit of type 1 is strongly attached on the edge with the bonds of Ec, Ed, and Ei that reside within the slice, while the growth unit of type 2 is only combined with Eh. Therefore, the growth unit incorporated into the kink site 2 is readily annihilated, and it gives rise to the stable edge comprised of the growth unit of type 1. For paracetamol, the stable edge is found to be coincident with that formed by thermal rearrangement.2 However, the edges for the (100) layers of L-Met generated by rearrangement are unstable, and the outermost growth units are readily annihilated. A stable edge that is bound to the crystal with more bonds has a tendency to maintain its structure of an edge, whereas an unstable edge that is weakly interacted with solute molecules within the crystal readily disintegrates. A new model for dealing with the unstable edge for the L-Met is illustrated in Figure 9. The double layered kink sites can be generated by the thermal rearrangement of type III, where the growth unit of type 2 is weakly bonded along the edge compared with the growth unit of type 1. Previously, it was confirmed that the growth units could be mainly rearranged by the type II along the edge because of its negative value of kink energy. Therefore, the initial structure of an edge generated by the rearrangement of type III has a tendency to be transformed to that by the rearrangement of type II, so that the stable edge consisting of growth unit 1 becomes prevalent. As a result of additional rearrangement by type II, the kink rate can be determined by eqs 11−15. In Figure 10, four different kink sites was found along the [1̅10] edge of the (110) layer. The thermal rearrangements of type III and IV bring on the same distribution of kink sites

5. RESULTS AND DISCUSSION 5.1. Growth Shape of L-Met Predicted by the Spiral Growth Model. In general, L-Met is difficult to grow as a single crystal from an aqueous solution and prone to form flakes (Figure 11a). In the present experiments, L-Met crystals were obtained from an aqueous solution without agitation to prevent the formation of agglomerates. A saturated L-Met/water solution of 2 mL at temperature of 60 °C was prepared in a 5 mL vial and kept at 25 °C until L-Met crystals grow to an appreciable extent. Figure 11b shows that a faceted L-Met crystal captured by a scanning electron microscope (SEM, JEOL/JSM-6010LA) is a hexagonal plate with the {001} faces that are morphologically important. Figure 12 presents the growth shape of L-Met predicted by the spiral growth model where the (001)_2, (100)_1, and (110) layers were considered, assuming that the net growth rate is governed by the layer of which growth rate is much slower at each face. In the present work, the growth rate of the (001)_2 layer is much smaller than that of the (001)_1 layer. Thus, the growth of the (001) face is controlled by the step of incorporation of G



Article L-Met crystals becomes thinner owing to a larger ΔWk at the (001)_2 layer when the work of adhesion is determined by using the harmonic mean. The work required for detaching a growth unit from kink sites determines the local concentration which considerably influences the kink rate. In particular, the solute molecules contained in the (001)_2 layer are strongly bound to the surface by forming hydrogen bonds, which brings about the lowest value of local concentration (Table S2). The kink rates of the (001)_2 layer are much slower than those of the (100) and (110) layers. As a result, the lowest local concentration of the (001)_2 layer reduces the flux into the kink site; thus, the kink rate becomes much slower and an extensive appearance of the (001) face follows. 5.2. Interfacial Structure Analysis. To understand the relationship between the surface structures and their local concentration in the crystallization of L-Met, the IS analysis14 was carried out. The effective growth units (named F1 units) are defined as those whose orientation is similar to that of solid molecules, SI, at the surface. The F1 units are very important because they can be directly incorporated into the kink sites at solid phase. Among F2 units with an orientation different from that of SI, there exist F1-like F2 units whose ΔG* is negligibly small so that they can be regarded as effective growth units. Therefore, the effective growth units are considered as both F1 and F1-like F2 units. In the present work, the orientation of the C5 → C2 dipole vector of the L-Met molecule was considered to identify the effective growth unit with providing the altitude angle (θC−C) and the azimuthal angle (ϕC−C). The interfacial concentration of effective growth units Xeff A(hkl) is given by14

Figure 12. Growth shapes of L-Met in water predicted by the spiral growth model: (a) the geometric mean and (b) the harmonic mean used in the calculation of the work of adhesion.

growth units into the (001)_2 layer. The relative growth rates of faces are summarized in Table S3. When the work of adhesion is estimated by the geometric mean, the spiral growth model produces a regular hexagonal plate much thicker than that from experiments. However, when the harmonic mean is used, the predicted morphology is in a remarkable agreement with the experiments. The harmonic mean gives generally lower values for the work of adhesion than the geometric mean,36 and therefore it gives rise to relatively large values of energies as given in Tables 3 and 4. The hydrogen bonding interaction typically belongs to the Coulombic interactions for common force field. Therefore, when Coulombic interactions within in a crystal and associative interactions within a solvent are mixed, the work of adhesion associated with the polar interaction becomes larger than it actually is. Accordingly, it may give rise to an overestimate of the work of adhesion. The intermolecular interactions of L-Met are greatly governed by strong Coulombic interactions. Therefore, the predicted shape of

X Aeff(hkl) = δXA(hkl)

(18)

where the factor δ is defined as

Figure 13. 3D plots of Gibbs free energy distribution for the orientation of C−C for the (a) (001)_1, (b) (001)_2, (c) (100)_1, (d) (100)_2, (e) (110), and (f) (110̅ ) layers, respectively. Colored circles indicate each position of SI unit. H



Article

∫ ζ(τ)P(τ) dτ * − ΔG **(τ )]/2kBT } = 1/2 ∫ P(τ ) × exp{[ΔG kink

δ=

* − ΔG **(τ )]/2kBT } dτ × sech{[ΔG kink

(19)

where P(τ) is the probability for an adsorbed growth unit at the state τ, and ζ is an effective factor to characterize the effectivity of F2 units. ΔGkink * is the potential barrier for transition from an F1 units to an SI units. ΔG**(τ) = 0 if ΔG*(τ) ≤ 0, and ΔG**(τ) = ΔG*(τ) if ΔG*(τ) > 0. The Gibbs free energy of the solute molecule at state τ can be obtained from the probability distribution function by the statistical thermodynamics. G*(τ ) = −kBT ln[P(τ )] + const

(20) Figure 14. Interfacial structure of L-Met molecules on the (a) (001)_2 and (b) (001)_1 surfaces, respectively. Water is indicated by oxygen (red) and hydrogen (white) atoms. For the sake of clarity, an oxygen atom of L-Met is only colored as red and the others as yellow.

The corresponding Gibbs free energy distributions were obtained as shown in Figure 13 (Configurations of the L-Met molecules on each crystal surface, relative concentration distribution profiles of L-Met, and their probability density distributions are provided in Figures S1−S3, Figures S4−S6, and Figure S7, respectively.). For the (001)_1 and (001)_2 layers (Figure 13a,b), the Gibbs free energy surface is shown to be very rugged, and the free energy barrier ΔG* is smaller than kBT. When ΔG* ≤ ΔGkink * , where ΔGkink * can be assumed to be comparable to the energy of thermal vibration, it leads to ΔG*kink ≈ kBT, ζ ≈ 1, and δ ≈ 1. In that case, Xeff A(hkl) ≈ XA(hkl), which implies that F2 unit as well as F1 unit can be regarded as the effective growth units.14 However, for the (100)_1 and (100)_2 layers, the lowest free energy local minimums are found to be around regions of SI units as shown in Figure 13c,d. In that case, the conformations of F2 units can favorably roll down into the corresponding traps of F1 unit where ΔG* < 0. Therefore, L-Met molecules orientated in SI unit become prevalent. Likewise, the (110) and (110̅ ) layers show that the lowest free energy local minimum is located around regions of SI units (Figure 13e,f) and furthermore ΔG* is negative. Accordingly, F2 unit can be spontaneously transformed into F1 unit with ease by diverse pathways. As a result of the MD simulations, the interfacial concentration of effective growth units, Xeff A(hkl), was determined as follows: eff eff eff Xeff A(001)_1 = 0.0703, XA(001)_2 = 0.0371, XA(100)_1 = 0.0465, XA(100)_2 eff = 0.0650, Xeff = 0.0476, and X = 0.0479, respectively. A(110) A(11̅0) Those results underpin the fact that the (001)_2 layer has the lowest value of xeq in the spiral growth model. As discussed earlier, the solute molecules within the (001)_2 layer are strongly bound to the under surface by forming hydrogen bonds. It brings about the lowest value of local concentration followed by reducing the kink rate. Therefore, the growth of the (001)_2 layer becomes a rate-determining step. Furthermore, the MD simulations give a new understanding of how the interfacial concentration of the (001)_2 layer is remarkably reduced compared with that of the (001)_1 layer. The hydrophobic CH3 group of L-Met is exposed to the solution on the (001)_2 surface, which exerts little impact on the approach of growth units belonging to the (001)_1 layer as shown in Figure 14a. However, Figure 14b shows that the hydrophilic NH3+ and COO− groups of L-Met on the (001)_1 surface are contacted with water molecules by the formation of hydrogen bonding, which hampers the approach of growth units appertaining to the (001)_2 layer. Therefore, L-Met molecules should surmount the energy barrier required for detaching water molecules adsorbed on the surface to continue to grow.

The decrease in interfacial concentration mainly occurs when F2 units showing a complete misorientation have larger ΔG*.14 However, in the present work, the ΔG* is very small or negative, and the lowest free energies are located near the regions of SI unit. It means the orientation of F2 unit is readily transformed into that of F1 unit. In that case, the local concentration determined by the spiral growth model, although which does not consider the entropy of a solute molecule on the crystal face, becomes comparable with the interfacial concentration obtained by the MD simulations. 5.3. Morphology of L-Met Predicted by the AE Model. In the present work, the morphology prediction of L-Met was also performed by the attachment energy (AE) model37 with two different set of faces: set I: (001)_1, (100)_1, (110) and (11̅0) layers, and set II: (001)_2, (100)_2, (110) and (110̅ ) layers. Attachment energies for L-Met were given in Table 5. The Table 5. Attachment Energy for L-Met layer

Eatt hkl (kcal/mol)

(001)_1 (001)_2 (100)_1 (100)_2 (110), (110̅ )

−0.95 −23.27 −13.16 −15.79 −33.98

morphology predicted by the AE model with set I seems to give a good agreement with the experimental results as shown in Figure 15a. On the other hand, the AE model with set II predicts the morphology of L-Met as a much thicker hexagonal plate, which is attributed to the relatively larger attachment energy of the (001) _2 layer (Figure 15b). In that case, it could be concluded that the (001)_1 layer plays an important role in determining the growth rate of the (001) face, which goes against the simulation results. In the AE model, the growth rate of a crystal face is simply assumed to be proportional to the magnitude of the attachment energy so that it might be misunderstood that the growth step of the (001)_1 layer is rate-determining rather than that of the (001)_2 layer. However, it should be noted that the AE model does not take into account of water molecules intensely interacting with the NH3+ and COO− groups of L-Met on the (001) crystal face. In this respect, to predict the crystal shape I



Article

hindrance of L-Met molecules approaching to the crystal face caused by hydration shell was not observed.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.cgd.6b00418. Table of atomic charge of zwitterionic L-Met. Table of work, existence probability, and local concentration for LMet in water. Table of kink density, critical length, kink rate, and relative growth rate for L-Met. Figures of initial and final configuration of the L-Met molecules on various surface; figures of relative concentration distribution profiles of L-Met as a function of distance; 3D plots of probability density distributions for the orientations of C− C (PDF)

Figure 15. Crystal shapes of L-Met predicted by AE model with (a) set I and (b) set II, respectively.

using attachment energies, the important anisotropic factors such as XA(hkl) and ΔWk encountered in the spiral growth model should be considered. Furthermore, in the AE model, the kinetics of growth units for att L-Met cannot be incorporated. Usually, the face with larger |Ehkl| produces a relatively weaker tangential force named the slice energy, Eslhkl, which yields a relatively lower kink energy followed by an appearance of an edge with a high kink density. Therefore, the faces with larger |Eatt hkl| tend to grow rapidly. However, this approach is not applicable for the prediction of crystal shape of LMet owning to the existence of two different layers in an interplanar distance. In that case, a smaller |Eatt hkl| did not guarantee the face with a larger |Eslhkl| and a rather layer-dependent local concentration determined by ΔWk plays an important role in the growth of L-Met as a hexagonal plate shape.

■

AUTHOR INFORMATION

Corresponding Author

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

The authors declare no competing financial interest.

■

REFERENCES

(1) Winn, D.; Doherty, M. F. AIChE J. 1998, 44, 2501−2514. (2) Kuvadia, Z. B.; Doherty, M. F. Cryst. Growth Des. 2011, 11, 2780− 2802. (3) Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2013, 13, 3341− 3352. (4) Winn, D.; Doherty, M. F. Chem. Eng. Sci. 2002, 57, 1805−1813. (5) Shim, H.-M.; Koo, K.-K. Cryst. Growth Des. 2014, 14, 1802−1810. (6) Shim, H.-M.; Koo, K.-K. Cryst. Growth Des. 2015, 15, 3983−3991. (7) Brosnan, J. T.; Brosnan, M. E. J. Nutr. 2006, 136, 1636S−1640S. (8) Matsuoka, M.; Yamanobe, M.; Tezuka, N.; Takiyama, H.; Ishii, H. J. Cryst. Growth 1999, 198-199, 1299−1306. (9) Timofeeva, E. V.; Routbort, J. L.; Singh, D. J. Appl. Phys. 2009, 106, 014304. (10) Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2012, 12, 656− 669. (11) Kim, S. H.; Dandekar, P.; Lovette, M. A.; Doherty, M. F. Cryst. Growth Des. 2014, 14, 2460−2467. (12) Gnanasambandam, S.; Rajagopalan, R. CrystEngComm 2010, 12, 1740−1749. (13) Liu, X. Y.; Bennema, P. J. Cryst. Growth 1996, 166, 117−123. (14) Liu, X. Y.; Boek, E. S.; Briels, W. J.; Bennema, P. J. Chem. Phys. 1995, 103, 3747−3754. (15) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899−4907. (16) Pople, J. A.; Beveridge, D. L. Approximate Molecular Orbital Theory; McGraw-Hill: New York, 1970. (17) Dalhus, B.; Gorbitz, C. H. Acta Chem. Scand. 1996, 50, 544−848. (18) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (19) No, K. T.; Cho, K. H.; Kwon, O. Y.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1994, 98, 10742−10749. (20) Svec, H. J.; Clyde, D. D. J. Chem. Eng. Data 1965, 10, 151−152. (21) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. (22) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43−56. (23) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Mod. 2001, 7, 306− 317. (24) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (25) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447.

6. CONCLUSIONS In the present work, we attempted to explain the reason why LMet crystals grow as a hexagonal plate from an aqueous solution, through the periodic bond chain and interfacial structure analysis. At the atomic level, the structure of L-Met is markedly characterized by the fact that two layers reside in an interplanar distance for the (001) and (100) faces, which leads to the emergence of different growth rates even in the same direction. The present simulation results clearly show that the (001)_2 and (100)_1 layers determine the growth shape of L-Met as a whole due to their very slow growth rates. In an aqueous solution, L-Met molecules in the (001)_2 layer are difficult to detach from their kink sites compared with those in the other layers, which leads to a dramatic decrease in the local concentration of the (001)_2 layer. In that case, a considerably reduced step velocity for the (001)_2 layer becomes a rate determining factor for the appearance of plate-like L-Met crystals. The existence of different local concentrations appearing in the spiral growth model was also underpinned by the present MD simulation results. On the (001)_1 surface, the NH3+ and COO− groups of L-Met are exposed to the solution, which enables water molecules to form hydrogen bonds with the crystal surface. However, on the (001)_2 surface, the CH3 group of L-Met is mainly exposed to the solution. The sulfur of L-Met connected to a methyl group is not so highly nucleophilic, and the side chain is most hydrophobic.7 Therefore, a van der Waals force dominates on the surface, and the energy involved in dehydration becomes quite smaller. For this reason, L-Met molecules can readily approach the (001)_2 surface by providing a relatively high interfacial concentration for the formation of the (001)_1 layer. On the other hand, on the {100}, {110}, and {110̅ } faces, the J



Article

(26) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (27) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101. (28) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (29) Polenske, D.; Lorenz, H. J. Chem. Eng. Data 2009, 54, 2277−2280. (30) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. London, Ser. A 2009, 465, 1145−1171. (31) Owens, D. K.; Wendt, R. C. J. Appl. Polym. Sci. 1969, 13, 1741− 1747. (32) Wu, S. J. Polym. Sci., Part C: Polym. Symp. 1971, 34, 19−30. (33) Wu, S. J. Adhes. 1973, 5, 39−55. (34) Zhang, J.; Nancollas, G. H. J. Colloid Interface Sci. 1998, 200, 131− 145. (35) Bennema, P. J. Cryst. Growth 1996, 166, 17−28. (36) Hansen, F. K. The Measurement of Surface Energy of Polymers by Means of Contact Angles of Liquids on Solid Surfaces. A Short Overview of Frequently Used Methods. University of Oslo, Oslo, 2004. (37) Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145−156.

K


Molecular Modeling on the Role of Local ... - ACS Publications

Recommend Documents