Nucleation Mechanisms of a Polymorphic Molecular Crystal: Solvent

Article pubs.acs.org/crystal

Nucleation Mechanisms of a Polymorphic Molecular Crystal: SolventDependent Structural Evolution of Benzamide Aggregates Philipp Ectors, Patrick Duchstein, and Dirk Zahn* Lehrstuhl für Theoretische Chemie/Computer-Chemie-Centrum, Friedrich-Alexander-Universität Erlangen-Nürnberg, Nägelsbach Str. 25, D-91052 Erlangen, Germany ABSTRACT: We report on molecular dynamics simulations dedicated to the aggregation and the early stage of nucleation of benzamide molecular crystals. As suggested by an earlier study of the molecular interactions in different benzamide polymorphs,1 we find the solvent to take a prominent role for polymorph control. The underlying driving forces are now rationalized by in-depth investigation of the formation of dimers, subsequent molecule association, and self-organization into aggregates of up to 100 benzamide molecules. Within the forming nuclei, we identify the development of ordered molecular motifs as characteristic for either the P1 or the P3 polymorph structure. From the different structural evolution related to the two nucleation scenarios, the preference of the stable P1 form is observed in aqueous solution, whereas precipitation from the vapor exhibits a tendency for the P3 polymorph of benzamide. Whereas benzamide aggregation from the vapor is mainly driven by hydrogen bonding between benzamide molecules, aggregation in aqueous solution is dominated by segregation of the hydrophobic moieties and π−π interactions of the benzene rings.

■

■

INTRODUCTION

The standard approach to tackling molecular crystal polymorphism is still trial and error. For industrial purposes, the testing of different solvents and thermodynamic conditions was even automated such that robots systematically probe suitable synthesis routes. Meanwhile, in-depth studies of the underlying molecular interactions, often based on computer simulations, are developing into an increasingly rational approach to control molecular crystal polymorphism. This applies to (i) the theoretical assessment of the manifold of different polymorphs2−4 and, more recently, (ii) molecular simulation as tools aiming at the rationalization of the nucleation mechanisms leading to the various crystal structures.5−11 By the example of an experimentally well-established benchmark system, benzamide, we recently demonstrated the analysis of molecular interactions in different polymorphic forms and its possible implications for polymorph control.1 From this, the favoring of either P1 or P3 forms of benzamide was attributed to different layer−layer interactions within periodic models of the bulk crystals. Guided by this postulated mechanism of P1 versus P3 polymorph selectivity, in the present study, we aim at the in-depth understanding of the whole morphogenesis of a forming aggregate, i.e., molecule-bymolecule association and subsequent self-organization in favor of crystalline motifs. To account for the role of the solvent, two different nucleation routes are investigated discriminating the amphiphilic nature of benzamide molecules during precipitation from polar solvent or nonpolar environment. © 2014 American Chemical Society

SIMULATION METHODS

Force Fields. The molecular interactions are mimicked by empirical potential energy functions that allow efficient sampling of aggregate structures at reasonable accuracy. The force field used for benzamide is based on fully flexible molecular models using the general AMBER force field (GAFF).12 We, however, refitted both the atomic charges and the amide−benzene torsion with respect to quantum calculations (see Scheme 1 and Table 1). For the partial

Scheme 1. Atomic Representation of Benzamide As Also Used in Table 1a

a

The torsion around the C(3)−C(5) axis has been reparameterized.

charge calculation, the benzamide dimer was used as reference, as all experimentally confirmed crystal structures13−15 show this very distinct motif. The dimer structure was optimized at the MP2 6-311+G(d,p) level and then used to calculate a restricted ESP-fit16 (Table 1). Using the same quantum approach, force-field parameters for the amide− benzene torsion were derived. The corresponding interactions are described by the below torsion potential plus the nonbonding Received: February 18, 2014 Revised: March 24, 2014 Published: April 16, 2014 2972

dx.doi.org/10.1021/cg500247c | Cryst. Growth Des. 2014, 14, 2972−2976

Crystal Growth & Design

Article

ambient conditions, aggregate relaxation is explored from 250 ps simulated-annealing runs at constant volume. In particular, during the relaxation runs for aggregates grown from the vapor, it is useful to not impose ambient conditions, but to start simulated annealing from a lower simulation temperature (200 K) to avoid excessive dissociation rates experienced for the newly added solute. This allows for more careful relaxation after the rather simplistic docking procedure and thus reduces failed attempts during the Monte Carlo-type steps mimicking putative aggregate growth steps. Several degrees of such “precooling” were tested for the annealing runs (260, 220, and 200 K) along with an increase in the sampling times up to 10 ns. We also draw the reader’s attention to the fact that a rigorous grandcanonical Monte Carlo procedure would imply accepting/rejecting growth attempts on the basis of the free energy of solute association to the forming aggregate. While, in principle, accessible to molecular simulation, here, we choose a (computationally much cheaper) acceptance criterion, namely, the direct inspection if the solute gets incorporated into the aggregate or is dissociated upon relaxation in solution. This intuitive criterion was found to perform excellently for a wide range of solutions11,19 with the exception of very high concentrations such as that of overconcentrated NaCl in water.19 Motif Recognition. We used a motif-recognition algorithm to identify crystalline features in the growing aggregates.20 The identification of crystalline motifs in only partially ordered nuclei is a nontrivial challenge, and particular care is required to avoid biasing the analysis by a too intuitive choice of the underlying order parameters. Indeed, the search for a specific molecular arrangement (such as the hydrogen-bonded dimers in the present work) is always subject to the chosen tolerance (root-mean-square deviation). Nonetheless, a much less biased analysis can be performed on a relative scale when comparing equivalent motifs of different polymorphs using identical tolerance. While the high-energy form P2 is of no relevance for the present study, the two most stable polymorphs of benzamide are illustrated in Figure 1. Both structures

Table 1. Atom Types and Force-Field Parameters of Benzamide atom

GAFF type

charge

ϵ [kJ/mol]

σ [Å]

N(1) H(2) C(3) O(4) C(5) C(6) C(7) C(8) H(9) H(10) H(11)

N HN C O CA CA CA CA HA HA HA

−0.716853 0.331928 0.588880 −0.542819 −0.038452 −0.089859 −0.168663 −0.118405 0.122741 0.145213 0.144929

0.71128 0.06569 0.35982 0.87864 0.35982 0.35982 0.35982 0.35982 0.06275 0.06275 0.06275

3.25 1.06908 3.39968 2.95993 3.39968 3.39968 3.39968 3.39968 2.59965 2.59965 2.59965

interactions as denoted in Table 1, using the usual 0.83 and 0.5 scaling factors for the Coulomb and the Lennard-Jones potentials, respectively.

Vtorsion = 10 kJ mol−1· [1 + cos(2φ − π )] + 1.4 kJ mol−1· [1 + cos(6φ)] Using this optimized force field, we were able to reproduce all crystal structures of benzamide polymorphs with an accuracy in potential energy of 1 kJ/mol per molecule (as compared to quantum calculations of the P1 and P3 structures).1,17 Note that the P1 and P3 polymorphs are energetically equivalent within these error margins. Moreover, the force field has been successfully used to predict metastable defects in form 1 of benzamide, in good agreement with NMR experiments.17 As the standard GAFF parameters were adopted for the nonbonded interaction models of the force field, mixing intermolecular interaction parametersas employed for the benzamide−water interactionsis possible. In what follows, water is described by the standard TIP3P model.18 Molecular Dynamics Simulations. The molecular dynamics simulations were carried out using a time step of 1 fs. Explicit Coulomb summation was used for the aggregates modeled in the gas phase, whereas a (shifted-force) cutoff of 1 nm was found appropriate for the aqueous solutions. Molecule-by-molecule association was investigated by means of a recently developed simulations scheme.11,19 The Kawska−Zahn approach reflects an iterative procedure that allows investigating aggregate evolution molecule-by-molecule starting with the formation of a dimer. Each growth step is devised into two processes, (i) the diffusion of a solute molecule to reach the aggregate and (ii) solute association and aggregate reorganization. The first step (i) reflects a random walk of initially dispersed solute molecules until close proximity to the forming aggregate. To avoid computationally demanding molecular dynamics simulation of this process, we mimic the outcome of solute diffusion to the aggregate by a simple docking procedure. Starting with the formation of a benzamide dimer, the Kawska−Zahn approach models each association step by placing a solute on a random position in proximity of the forming aggregate to account for unbiased solute diffusion leading to aggregate−solute contacts. Next, possible association or solute rejection is explored from simple potential energy minimization, which is implemented as rigid aggregate−solute docking in the absence of solvent molecules. Finally (ii), aggregate relaxation is explored from unconstrained molecular dynamics simulations allowing full flexibility of the aggregate and considering the solvent. Thus, an effective Monte Carlo-type procedure is used for docking new solutes to a forming aggregate, while, after each solute association, the relaxation of the aggregate is explored in detail from molecular dynamics simulations.11,19 When exploring aggregation from solution, the relaxation step (ii) is implemented by placing the aggregate in a bulk model of explicit water molecules treated as a cubic box of ∼5 nm cell edge length. After 50 ps of constant pressure, constant temperature simulations at

Figure 1. Illustration of the P1 and P3 structures of benzamide.13,14 While the units cells are shown as gray boxes, hydrogen-bonded dimers constituting the key building blocks are highlighted in red. To discriminate P1 and P3 motifs, the orientation of such dimers with adjacent benzamide molecules is analyzed (dashed green line). The tilt angles are 0° and 74° in P1 and P3, respectively. are composed of hydrogen-bonded dimers that form hydrogen-bonded sheets within the (001) plane. P1 and P3, however, differ in terms of the contact angle and distance of the neighboring molecules of adjacent sheets. Thus, trimers comprising a hydrogen-bonded dimer (highlighted in red in Figure 1) and an adjacent molecule “attached” via a benzene−benzene contact represent the smallest molecular motif that allows discrimination between P1- and P3-type ordering. Comparing a putative dimer/trimer motif in the aggregate to the corresponding dimer/trimer motifs in the bulk crystal structures as reference, a root-mean-square deviation of less than 1.3 Å was chosen as the threshold for crystalline motif recognition. This treshold was found to provide best discrimination between P1- and P3-type motifs.

■

RESULTS Guided by the suggested preference of P1- and P3-type polymorph growth in polar and nonpolar solutions, respectively, we explored two setups of benzamide precipitation. 2973



Article

Molecule-by-molecule association up to the aggregation of 125 benzamide molecules was explored (a) in the vapor as a proxy to nonpolar solution (several nonpolar solvents were tested for the relaxation of a small series of aggregates demonstrating vacuum as an excellent proxy for such solvation. An important exception is benzene, which does π−π stacking with the solute and thus interacts more selectively) and (b) in aqueous solution. For each setup, three independent growth runs were performed to provide at least qualitative insights into the statistical manifold of possible nucleation and growth routes. The mechanisms discussed in the following were found in all of these simulation runs. Benzamide aggregation from the vapor is dominated by hydrogen bonding, which is particularly evident for aggregates of less than about 10 molecules. A representative growth run is illustrated in Figure 2a. Upon further aggregate growth, Figure 3. (a) Snapshots of benzamide aggregate growth from aqueous solution. Solute−solute interactions are characterized by the interplay of hydrophobic segregation and hydrogen-bonded dimer formation. Colors: H (white), C (gray), N (blue), and O (red). The solvent molecules and the hydrogen atoms of the benzene ring are omitted for clarity. (b) Aggregate of 125 benzene molecules as precipitated from aqueous solution at different representations, but identical orientation. Left to right: all molecules, highlighting of dimer motifs, highlighting of parallel dimers (sheets) highlighting P1 (red) versus P3 (green) type motifs.

aggregate explored, comprising 125 benzamide molecules, is illustrated in more detail in Figure 3b. Whereas the onset of ordering is observed in the inner part of the aggregate, the outer boundaries are poorly ordered and exhibit amide−amide and amide−water hydrogen bonds at roughly equal rates. In contrast to the rigorous benzamide−vacuum interface, this leads to a weakly ordered and more extended benzamide− water interface region. Despite the still small number of ordered motifs observed during these early stages of aggregation, it is already possible to discriminate the different rates of hydrogen-bonded benzamide dimer formation as identified for aggregation from the vapor and in aqueous solution. Figure 4 shows the absolute number of such motifs as a function of aggregate size as sampled from three independent growth runs for precipitation both from the vapor and in water. Note that this ordering increases gradually and not by spontaneous disorder−order transitions. While the dimers reflect the smallest and most simple motif of the molecular crystal, more subtle discrimination between P1- and P3-type polymorph nucleation requires the analysis of trimer motifs, as outlined in Figure 1. Even for aggregates comprising more than 100 molecules, the overall number of such motifs is rather small (see Figures 2b and 3b for illustration). Nonetheless, for aggregate growth from the vapor, we can identify the biasing in favor of the metastable polymorph P3. The overall sampling from all aggregate runs up to 125 molecules shows 80% of all trimer motifs are of P3-type. On the other hand, the analogous analysis of aggregation from aqueous solution exhibits only a weak bias in favor of 60% P1-type motifs. Thus, even for the very early stages of aggregation explored in this study, we can relate the favoring of different polymorphic structures to the preferential formation of P1versus P3-type trimer motifs, in agreement to experimental and other theoretical evidence, both collected for macroscale crystallites.1,17 With growing aggregate size, the importance of

Figure 2. (a) Snapshots of aggregate growth from the vapor. Starting from the hydrogen-bonded dimer compact, almost spherical hydrogenbonded clusters are formed. Colors: H (white), C (gray), N (blue), and O (red). The hydrogen atoms of the benzene ring are omitted for clarity. (b) Aggregate of 125 benzamide molecules as precipitated from the vapor at different representations, but identical orientation. Left to right: all molecules, highlighting of dimer motifs, highlighting of parallel dimers (sheets) highlighting P1 (red) versus P3 (green) type motifs.

benzene−benzene contacts are also observed; however, π−π stacking plays a minor role even in the largest aggregates explored comprising 125 molecules (Figure 2b). As a consequence, the aggregates mainly reflect an agglomerate of hydrogen-bonded dimers, which tend to arrange in favor of hydrogen-bonded sheets in the course of aggregate growth. However, this ordering is observed only partially (i.e., locally for bulk-like domains), while the outer boundaries of the aggregate are characterized by benzene groups at the surface to avoid undercoordinated hydrogen bond acceptors and donors. Indeed, up to 125 molecules, the overall aggregate shape is close to spherical and thus mainly controlled by surface tension. On the other hand, benzamide aggregation from aqueous solution is found to be driven by two roughly equally strong interactions. Figure 3a illustrates the morphogenesis of a representative growth run. In particular, for the dimer, but also the binding of trimers, tetramers, etc., we clearly observe an interplay of hydrophobic segregation of the benzene moieties and hydrogen bonding between the amide groups. As a consequence, hydrogen-bonded benzamide dimer motifs are observed at much later stages of aggregate growth. The largest 2974



Article

Figure 4. (a) Motif recognition of hydrogen-bonded dimers in aggregates grown from the vapor (red) and in water (blue) as a function of aggregate size. In aqueous solution, benzamide−benzamide hydrogen bonding competes with hydrophobic segregation. Aggregate ordering in favor of dimer motifs thus occurs at later stages, once a desolvated aggregate core is developing. See also Figures 1 and 2. (b) Motif recognition of sheets comprising two parallel hydrogen-bonded dimers, shown for aggregates grown from the vapor (red) and in water (blue) as a function of aggregate size. See also Figures 1 and 2.

■

CONCLUSIONS Employing the Kawska−Zahn approach, we investigated benzamide molecular crystal nucleation as a series of molecule-by-molecule association events. Our simulation procedure starts with the very early stage of dimer association, followed by aggregation of up to 125 molecules. Within the growing nuclei, we clearly observe the self-organization in favor of either P1- or P3-type polymorphs, depending on the solvent environment. In agreement with the available experimental experience, our study thus shows the suitability of our simulation approach to rationalize the mechanisms of molecular

surface/interface effects diminishes and the quantitative P1/P3 motif ratio surely changes (to more and more pronounced favoring). Our findings should, therefore, be considered from a more qualitative point of view, i.e., the selectivity of P1- versus P3-type trimers as a consequence of different driving forces to solute association in polar and nonpolar solutions. This constitutes a bias in forming aggregates, which ordering was found to increase continuously. Even if aggregate growth at later stages should involve an abrupt disorder−order transformation, we still argue that the observed motifs reflect nucleation seeds giving rise to similar P1 versus P3 selectivity. 2975



Article

crystal aggregation and the nucleation of competing polymorphic structures. While the investigated benzamide aggregates undergo rapid relaxation, followed by vibrations and rotation of the whole aggregate, we point out that the case of prenucleation clusters is also accessible to our method. This was recently elaborated by de Yoreo and co-workers for calcium carbonate aggregation in water,21 who identified liquid−crystalline-type aggregate structures that were sampled from parallel-replica molecular dynamics simulations.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We gratefully acknowledge funding from DFG within the grant ZA-420/4. REFERENCES

(1) Ectors, P.; Zahn, D. Phys. Chem. Chem. Phys. 2013, 15, 9219− 9222. (2) Gavezzotti, A. Acc. Chem. Res. 1994, 27, 309−314. (3) (a) Schön, C.; Jansen, M. Angew. Chem. 1996, 108, 1358−1377; (b) Angew. Chem., Int. Ed. Engl. 1996, 35, 1286−1304. (4) Price, S. L. Acc. Chem. Res. 2009, 42, 117−126. (5) TenWolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. J. Chem. Phys. 1996, 104, 9932−9947. (6) Gavezzotti, A.; Filippini, G.; Kroon, J.; van Eijck, B. P.; Klewinghaus, P. Chem.Eur. J. 1997, 3, 893−899. (7) Vekilov, P. G. Cryst. Growth Des. 2004, 4, 671−685. (8) Piana, S.; Reyhani, M.; Gale, J. D. Nature 2005, 438, 70−73. (9) Lechner, W.; Dellago, C.; Bolhuis, P. G. Phys. Rev. Lett. 2011, 106, 085701. (10) Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. J. Am. Chem. Soc. 2012, 134, 17221−17223. (11) Anwar, J.; Zahn, D. Angew. Chem., Int. Ed. 2011, 50, 1996−2013. (12) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollmann, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (13) Gao, Q.; Jeffrey, G. A.; Ruble, J. R.; McMullan, R. K. Acta Crystallogr., Sect. B 1991, 47, 742−745. (14) Thun, J.; Seyfarth, L.; Senker, J.; Dinnebier, R. E.; Breu, J. Angew. Chem., Int. Ed. 2007, 46, 6729−6731. (15) David, W. I. F.; Shankland, K.; Pulham, C. R.; Blagden, N.; Davey, R. J.; Song, M. Angew. Chem., Int. Ed. 2005, 44, 7032−7035. Blagden, N.; Davey, R.; Dent, G.; Song, M.; David, W. I. F.; Pulham, C. R.; Shankland, K. Cryst. Growth Des. 2005, 5, 2218−2224. (16) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollmann, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (17) Butterhof, C.; Martin, T.; Ectors, P.; Zahn, D.; Niemietz, P.; Senker, J.; Näther, C.; Breu, J. Cryst. Growth Des. 2012, 12, 5365− 5372. (18) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (19) Kawska, A.; Brickmann, J.; Kniep, R.; Hochrein, O.; Zahn, D. J. Chem. Phys. 2006, 124, 024513. (20) Duchstein, P.; Hochrein, O.; Zahn, D. Z. Anorg. Allg. Chem. 2009, 635, 649−652. (21) Wallace, A. F.; Hedges, L. O.; Fernandez-Martinez, A.; Raiteri, P.; Gale, J. D.; Waychunas, G. A.; Whitelam, S.; Banfield, J. F.; De Yoreo, J. J. Science 2013, 341, 885−889.

2976


Nucleation Mechanisms of a Polymorphic Molecular Crystal: Solvent

Recommend Documents