David P.
McMahon
a,
Andrew
Stephenson
b,
Samantha Y.
Chong
b,
Marc A.
Little
b,
James T. A.
Jones‡
b,
Andrew I.
Cooper
*b and
Graeme M.
Day
*a
aComputational Systems Chemistry, School of Chemistry, University of Southampton, SO17 1BJ, UK. E-mail: g.m.day@soton.ac.uk
bDepartment of Chemistry and Materials Innovation Factory, University of Liverpool, Crown St., Liverpool L69 7ZD, UK. E-mail: aicooper@liverpool.ac.uk
First published on 16th April 2018
Crystal structure prediction methods can enable the in silico design of functional molecular crystals, but solvent effects can have a major influence on relative lattice energies, sometimes thwarting predictions. This is particularly true for porous solids, where solvent included in the pores can have an important energetic contribution. We present a Monte Carlo solvent insertion procedure for predicting the solvent filling of porous structures from crystal structure prediction landscapes, tested using a highly solvatomorphic porous organic cage molecule, CC1. Using this method, we can understand why the predicted global energy minimum structure for CC1 is never observed from solvent crystallisation. We also explain the formation of three different solvatomorphs of CC1 from three structurally-similar chlorinated solvents. Calculated solvent stabilisation energies are found to correlate with experimental results from thermogravimetric analysis, suggesting a future computational framework for a priori materials design that factors in solvation effects.
An important factor that is usually overlooked in CSP is the effect of solvent on the relative crystal structure stabilities and, hence, the structure that is predicted to form. Alternative crystal structures are often predicted to be close in lattice energy and the interaction of solvent molecules during crystallisation can affect the stabilisation of one potential packing arrangement with respect to another. Frequently, the influence of solvent is limited to interactions with the surface of a crystal. The total free energy of a crystal has contributions from the bulk and surface energies, such that structures with the lowest surface energies can be more stable than those with the lowest bulk energy. The balance between bulk and surface energies depends on crystallite size and shape, with surface terms increasing for small crystallite sizes.14 While the bulk energy of a crystal is unaffected by the solvent environment, surface energies can be strongly influenced. Thus, different solvents can stabilise different structures and solvent stabilisation of crystal surfaces can have an important influence at the critical nucleation stage, leading to the growth of metastable phases. These effects on polymorph stabilities have been modelled for known polymorphs of inorganic15,16 and organic17 materials using explicit15,16 or continuum17 models of solvent on the dominant crystal faces.
The influence of solvent becomes especially important for porous crystals, which tend to crystallise with solvent in the pore channels. Solvated polymorphs, or “solvatomorphs”, are of intrinsic interest because solvatomorphs can have different physical properties.18–21 For porous crystals, different porous polymorphs can be isolated by desolvating solvatomorphs with different crystal structures.10,22,23 This can be used to an advantage: for example, we showed previously that 1,4-dioxane acts as a directing solvent to induce the formation of a specific highly-porous crystal packing of organic cage molecules.23 The solvent molecules in a porous molecular crystal can also be an integral part of the crystal structure, and hence the removal of solvent from solvatomorphs with retention of the crystal packing of the host structure is not always possible.24,25
It has been shown that the inclusion frameworks of organic molecules often correspond to local minima on the lattice energy surface, even in the absence of the guest solvent molecules.26 They can therefore be located using CSP methods, usually at energies significantly above the global lattice energy minimum. These structures can be identified because they contain voids of the right size to accommodate molecular guests. It is not usually clear, a priori, which solvent will stabilise which particular crystal packing.27 As such, predicting the influence of solvent on relative crystal lattice energies is of strong value because it could allow crystallisation conditions to be explored in silico to target specific polymorphs. For example, we recently used solvent stabilisation calculations to explain why a hydrogen bonded molecule, triptycene-tris(benzimidazolone), crystallises as a series of porous polymorphs and not as a dense, non-porous phase, which is the predicted global lattice energy minimum structure.10
Here, we investigate the effect of solvent on the crystallisation of an ethylenediamine-derived [4 + 6] imine cage, CC1 (Fig. 1a).22,28CC1 is a good candidate for studying solvent effects on polymorphism: unlike many imine cages, where bulky vertices control the crystal packing, CC1 is nearly spherical and it has many structures of similar energy on its predicted crystal energy landscape.12 The CC1-α′ polymorph is formed by desolvation of an ethyl acetate (EtOAc) solvate, CC1·2.5(EtOAc), while CC1-β′ is formed by desolvation of a dichloromethane (DCM) solvate, CC1·2.5(DCM).22 The desolvated crystal structures are related directly to their respective solvatomorphs and the two different polymorphs have different properties: CC1-β′ is porous to nitrogen while CC1-α′ is not.22 The phase behaviour of CC1 is not limited to these polymorphs: for example, recrystallisation from DCM/mesitylene affords CC1·4(mesitylene) (P),29 and recrystallisation from DCM/1,4-dioxane affords 2(CC1)·7(1,4-dioxane) (P21/c).23 Hence, the solid form landscape of CC1 is particularly rich and non-intuitive from an experimental perspective.
We present here the preparation and characterisation of five new solvatomorphs of CC1 (Fig. 1), which exist as either single solvates (one cocrystallised solvent) or as bisolvates (two different cocrystallised solvents). By using a Monte Carlo solvent filling procedure in conjunction with lattice energy minimisation, we show that all the solvatomorphs are solvent stabilised compared to the hypothetical global minimum on the CSP landscape. In addition, we show for three different chlorinated solvents that the experimentally obtained crystal solvatomorph is correctly predicted by solvent stabilisation calculations; that is, the observed solvatomorph is predicted to be more stabilised by the solvent used for its crystallisation than any of the alternative CC1 structures that were tested. This is further corroborated by thermogravimetric analysis (TGA) of the experimental solvates, which shows a correlation between the experimental stability of the solvatomorph and the computed solvent stabilisation energy.
Single crystal X-ray data sets were measured on a Rigaku MicroMax-007 HF rotating anode diffractometer (Mo-Kα radiation, λ = 0.71073 Å, Kappa 4-circle goniometer, Saturn 724+ detector); at beamline I19, Diamond Light Source, UK using silicon double crystal monochromated synchrotron radiation (λ = 0.6889 Å, Saturn 724+ detector);30 or at beamline 11.3.1, Advanced Light Source, Berkeley, USA, using silicon monochromated synchrotron radiation (λ = 0.7749 Å, APEX-II detector). Solvated single crystals were immersed in a protective oil, mounted on a MiTeGen loop, and flash cooled under a dry N2 gas flow. Absorption corrections, using the multi-scan method, were performed with the program SADABS.31 Structures were solved with SHELXD,32 or by direct methods using SHELXS,33 and refined by full-matrix least squares on |F|2 by SHELXL.34 Supplementary CIFs, that include structure factors, are available free of charge from the Cambridge Crystallographic Data Centre (CCDC). 2(CC1)·7.75(DCM) #1575913, 2(CC1)·11.32(CHCl3) #1575914, 2(CC1)·10(CCl4) #1575912, CC1·2(oX)·CHCl3 #1575911, and CC1·oX·DCM #1575910. For full refinement details, see the ESI (Tables S1 and S2, Fig. S1–S7†).
Powder X-ray diffraction data was collected in transmission mode using the Mythen-II position sensitive detector and capillary spinner on the I11 beamline at the Diamond Light Source (λ = 0.827157 Å). The sample was contained in a 0.5 mm diameter borosilicate glass capillary. The structure of CC1·1.47(pX) was solved using the simulated annealing routine implemented in Topas-Academic.35 For PXRD refinement details, see the ESI (Fig. S8 and S9†).
TGA was carried out using a TA Q5000IR analyser with an automated vertical overhead thermobalance. Samples were heated at a rate of 10 °C min−1. Tonset was calculated using the in-built feature in the software ‘TA Universal Analysis’. TGA curves and calculated Tonset values are in the ESI.†
Lattice energy minimisation calculations used the DMACRYS37 software in which molecular geometries are held rigid at the gas phase B3LYP/6-31G** geometry. Intermolecular interactions are described using an empirically parameterised atom–atom repulsion–dispersion model (W99 (ref. 38)) and a distributed multipole electrostatic model, using a distributed multipole analysis39 (DMA) of the B3LYP/6-31G** electron density with multipoles up to hexadecapole on each atom. Ewald summation was used to calculate the charge–charge, charge–dipole, and dipole–dipole interactions; all other intermolecular interactions were subject to a 30 Å cutoff. Full details are provided in the ESI.†
We attempted insertion of each solvent (CHCl3, DCM, CCl4, EtOAc, 1,4-dioxane, p-xylene) and solvent mixtures (CHCl3:o-xylene and DCM:o-xylene) into the CC1 framework that is observed to form with that solvent. For mixed solvent systems (bisolvates), the experimental solvent ratios were used during the solvent insertion calculations. It is a future challenge to extend our current methods to assess mixed solvent systems without experimental information on the solvent ratio in the solvated crystal structure. We also used the Monte Carlo approach to load five of the single solvents (CHCl3, DCM, CCl4, EtOAc, 1,4-dioxane) into all of the solvent frameworks and the CSP global minimum crystal structure.
An artificially high temperature of T = 5000 K was used to ensure good sampling of solvent configurations within the crystal structures. Frames sampled at 10 cycle intervals were lattice energy minimised with the same energy model used in the CSP calculations (W99 + DMA). No space group symmetry was enforced during solvent sampling or lattice energy minimisation. Solvated structures that led to large distortions of the CC1 framework were discarded. This was tested by comparison of the CC1 arrangement within the solvated structure to the original CSP structure using the COMPACK algorithm.41
The lattice energies of the solvated crystal structures were then corrected for the cost of removing N molecules of solvent from the pure phase:
(1) |
The correction accounts for the change in the kinetic energy contribution of the solvent molecules to the internal energy. We assume solvent molecule rotational motion to be free in the bulk solvent phase, and to become librational vibrations in the solvate crystals and use equipartition estimates of the energy. Econf is an energy correction between the conformations of CC1, which is added to any structure with a higher energy conformation.
In CC1·2.5(EtOAc) (C2/c), one EtOAc molecule was located in the CC1 cavity with additional disordered EtOAc molecules being found throughout the narrow 1-D channel along the b-axis.28 In 2(CC1)·7(1,4-dioxane) (P21/c), 1,4-dioxane sits in the window of two adjacent cage molecules and directs CC1 to pack window-to-window. Extrinsic pockets also exist between CC1 molecules filled with two 1,4-dioxane molecules.23,28
By using anti-solvent vapour diffusion crystallisations with ortho-xylene (oX) and para-xylene (pX), we isolated three additional solvatomorphs with different CC1 crystal structures: CC1·DCM·oX (P212121); CC1·CHCl3·2(oX) (P21/n); and CC1·1.47(pX) (P21) (Fig. 1). CC1·CHCl3·2(oX) was obtained from a mixture of chloroform and oX, crystallising in the form of thin needles (Fig. S6†). In this structure, CHCl3 occupies the cage centres and oX occupies the interstitial spaces between the cage molecules to form extended, solvent-filled 1-D channels along the crystallographic a-axis. When DCM and oX were used, CC1·DCM·oX was obtained (Fig. S7†). As with the CHCl3 analogue, the chlorinated solvent is in the cage centre. In this case, oX is in 1-D channels along the crystallographic b-axis and the resultant packing motif of CC1 is very different. In CC1·DCM·oX, the channels are formed by arrangement of 6 CC1 molecules such that the oX sits in a relatively large hexagonal-shaped channel (Fig. 1h). By contrast, in CC1·CHCl3·oX, the oX molecules sit in smaller channels surrounded by just 4 cages (Fig. 1g). This leads to a different CC1:solvent ratio; in both cases, there is one chlorinated solvent molecule per cage, but in the CHCl3·oX solvate, there are two oX molecules per cage, while in the DCM·oX solvate there is only one. We were unable to obtain suitable single crystals for the pX solvate of CC1, CC1·1.47(pX), but its structure was solved by powder X-ray diffraction (PXRD) in the monoclinic space group, P21 (Fig. S8 and S9†). There are no solvent-filled channels between cages in this solvatomorph (Fig. 1i); instead, the pX solvent sits inside the cage cavity and also in the interstitial voids, pointing into the cage window. The intrinsic site is partially occupied, with an occupancy of less than 50% in contrast to the extrinsic window site, which shows full occupancy. The total number of pX guests per cage is 1.471(2).
The two unsolvated experimental polymorphs, CC1-α′ and CC1-β′, are both present in the predicted set of structures, 6.3 and 6.2 kJ mol−1 above the global minimum in the energy landscape, respectively, which is within the usual energetic range of polymorphism43 in organic molecular crystals. For comparison, we calculated lattice energies of the various experimental solvates (Fig. 1) after artificial removal of the lattice solvent and subsequent lattice energy minimisation. All of these ‘artificial desolvates’ are located on the CSP landscape. The resulting structures are plotted on the CSP crystal energy landscape (Fig. 3, unfilled coloured circles). These results confirm that all 7 of the experimental solvate packings for CC1 can be located using CSP methods without including solvent molecules in the calculation. In the case of the high-energy Form II structure, the CC1 framework was only found after additional sampling in the observed space group (C2/c, Z′ = 1); this is a particularly high energy structure (53.8 kJ mol−1 above the global minimum) and there are many local minima at such high energies, making the location of all local minima very challenging.
The significant challenge, then, is to identify which of the many predicted structures will be stabilised by which solvent or solvent mixture. The first step, tackled here, is to understand the relevant energetics: that is, the energy range on the CSP landscape where solvates are found and the energetic stabilisation that arises from the inclusion of solvent molecules. The calculated lattice energies for the artificially desolvated CC1 structures range from 12.6 to 56.7 kJ mol−1 above the global lattice energy minimum – that is, the lowest energy solvent-free crystal packing predicted for CC1 (Fig. 3). Some of these solvate frameworks are therefore much higher in energy than typical organic polymorphs; hence, these packing arrangements of CC1 are only observed because of the energetic influence of solvent templating.
The sampling performed in the Monte Carlo approach is important in evaluating the solvent stabilisation energies because of the strong dependence of the energy on the positions and orientations of the solvent molecules in the CC1 framework. Fig. 4 shows the range of energies for the sampled configurations of DCM in Form I at all sampled solvent loadings (all sampled configurations of each solvent system are also shown in Fig. S11†). At N = 4, where the CC1 framework is fully solvated, the different solvent configurations produced by the Monte Carlo sampling cover a range of approximately 40 kJ mol−1 and only a few of these are close in energy to the lowest energy configuration.
To evaluate the energetic stabilisation, we assume ordered solvent in the voids of the CC1 structures and, so, we take the most stable calculated solvent arrangement in each of the host CC1 frameworks Forms I–VII. We find that all the structures are stabilised significantly by solvent inclusion in their voids; this stabilisation is calculated to be 29.4–133.8 kJ mol−1 for the various solvates. Importantly, these large stabilisation energies lower the lattice energy of the known solvates below the hypothetical CSP global minimum energy structure on the crystal energy landscape (Fig. 3). This demonstrates that the solvated CC1 structures are more stable than the unsolvated CSP global minimum and rationalises why the observed solvated structures of CC1 are formed from organic solvents in preference to the lowest energy possible packing of pure CC1. The extent of the calculated solvent stabilisation with respect to the CSP global minimum ranges from 13.3 (Form VIII, pX) to 77.1 (Form III, CCl4) kJ mol−1.
The data in Fig. 4 also show how the CC1 solvate system lowers its energy quickly as the first solvent molecules are added, but then starts to level off, approaching a minimum in stabilisation energy as more solvent is added (for DCM in Form I, around N = 3 and 4). Beyond N = 4, any further addition of DCM completely disrupts the CC1 Form I framework structure. For each solvatomorph, we predict the solvent loading as the composition with the lowest energy before the CC1 arrangement is disrupted.
For the most stabilised solvent structures, the Monte Carlo predicted solvent loading agrees remarkably well with the experimental composition (Table 1). For all cage-solvent pairings, the predicted ratio of CC1 to solvent is within 1 of the experimental value, which represents the closest possible agreement, given the use of integer sampling ratios used in the computational work. This suggests that it is possible to predict, a priori, the most stable solvation level for a given solvent within a predicted structure on a CSP landscape.
Solvent | N observed | N predicted | |
---|---|---|---|
Form I | DCM | 3.88 | 4 |
Form II | CHCl3 | 5.66 | 5 |
Form III | CCl4 | 5 | 5 |
Form IV | EtOAc | 2.5 | 2 |
Form V | 1,4-Dioxane | 3.5 | 3 |
Form VI | CHCl3·oX | 1:1 | 1:1 |
Form VII | DCM·oX | 1:2 | 1:2 |
Form VIII | pX | 1.47 | 1 |
To assess the predictive potential of our methods, we considered the five solvents DCM, CHCl3, CCl4, EtOAc and 1,4-dioxane, using the Monte Carlo procedure to insert each solvent into each of the artificially desolvated CC1 Forms I–V, as well as the CSP global minimum structure (CSP min). The solvated lattice energies were then calculated for the most stable loading of each solvent into each solvate framework (Table 2).
First, these results clearly rationalise why the CSP global lattice energy minimum structure is not observed from crystallisation in the presence of any of these five solvents. The global minimum structure from the CSP landscape can accommodate each of the solvents studied here, and it is energetically stabilised by most solvents (final row, Table 2), apart from oX and pX, where the lattice energy is increased relative to the neat, solvent-free CC1 structure (Table S8†). However, for all the solvents tested, there is always at least one competing solvate, including the experimentally observed solvate in each case, that is more stable (Tables 2, S8 and S9†); as such, the global minimum predicted structure is never predicted to be favoured.
For the three chlorinated solvents, DCM, CHCl3, and CCl4, we found excellent agreement between prediction and experiment; for each solvent, the observed solvate structure is the most stable of the CC1 structures after computational solvation (Table 2, green highlighted values).
For EtOAc, we find that Form I and Form IV (the observed EtOAc solvatomorph) are effectively equi-energetic, with Form I being slightly more stable by 0.4 kJ mol−1. EtOAc is the largest and most flexible of the solvents tested in this cross comparison and it is possible that the rigid-molecule approach used in our calculations performs less well in terms of distinguishing the observed solvate structure than for essentially rigid solvent molecules (i.e., DCM, CHCl3, and CCl4). Nonetheless, from these calculations, one might predict that EtOAc would lead to either Form I or Form IV, the latter of which is observed by experiment.
In the case of the 1,4-dioxane solvate, we found that none of the four alternative solvents could solvate the CC1 framework of the observed structure (Form V); for each of DCM, CHCl3, CCl4 and EtOAc, solvent insertion leads to a rearrangement of the Form V CC1 framework. The observed packing, in which 1,4-dioxane lies flat between the cages, cannot be replicated by the other solvents. The fact that no other solvent can form this solvated crystal structure agrees well with experiment – in a previous study,23 more than 40 other solvents were trialled in the laboratory, but none was found to direct the CC1 window-to-window packing observed in the 1,4-dioxane solvate, Form V. Comparing the inclusion of 1,4-dioxane within the set of CC1 frameworks (column 5; Table 2), we find that 1,4-dioxane forms a more stable solvate in the Form IV CC1 arrangement than the observed Form V. This is due to a failure of the Monte Carlo sampling to locate the most stable configuration of 1,4-dioxane within Form V; by calculating the lattice energy using the positions of 1,4-dioxane from the experimentally determined crystal structure, we obtain an energy for the Form V 1,4-dioxane solvate of −241.9 kJ mol−1. This configuration is 28.7 kJ mol−1 more stable than the most stable configuration determined by our solvent insertion methodology and it renders Form V the lowest energy of the tested pairings between 1,4-dioxane solvent and the CC1 host frameworks. This 1,4-dioxane result shows that the solvent filling approach could be successful, but illustrates that further developments in solvent configuration sampling methods are needed to make these methods more robust in the future.
Solvent insertion into predicted crystal structures shows good results in terms of qualitatively identifying those crystal structures in a CSP landscape that can be solvent stabilised and assessing the preferential solvent stabilisation between possible host frameworks. One aspect that we are as yet unable to address with the current method is the possibility of structural transformations that can occur upon desolvation of a solvated crystal structure. For example, desolvation of CC1·2.5(EtOAc) to afford CC1-α′ results in a 10% contraction of the unit cell volume and a transformation of the space group from C2/c to P21/c.22 Artificial, in silico desolvation and lattice energy minimisation of the experimental CC1·2.5(EtOAc) solvate leads to a polymorph which is 6.3 kJ mol−1 higher on the CSP landscape than the experimentally determined structure. Likewise, the DCM solvate, 2(CC1)·7.75(DCM), transforms during desolvation which results in a change from 2 to 1 molecule in the asymmetric unit (Z′ = 2 → Z′ = 1), which is related to a conformational change in one of the CC1 molecules.22
A broad correlation is observed between (Tonset − Tbp) and the calculated solvent stabilisation energies (Fig. 5). TGA shows that 2(CC1)·10(CCl4) is by far the most thermally stable of the CC1 solvates (Tonset − Tbp = 60 °C) and this solvate also has the largest calculated stabilisation energy (−133.8 kJ mol−1). This encouraging comparison extends to all three chlorinated solvents and EtOAc, where the stability orders based on TGA and the computational results are the same. The 1,4-dioxane solvate also fits this correlation, particularly when we use the energy calculated using the experimentally-determined solvent position within the Form V CC1 framework (blue star in Fig. 5), rather than the Monte Carlo minimum energy solvent configuration.
This computational approach is also shown to have promise in predicting the preferential stabilisation that different solvents provide to different porous crystal structures: for example, the calculations explain why three structurally-similar chlorinated solvents force CC1 to adopt three different solvated crystal packing motifs. In each case, the solvated structure is most stable in the correct experimental packing motif and that packing motif is stabilised the most by the solvent from which it is crystallised.
These calculations could be used in a predictive sense by solvating the crystal structures produced by CSP with a range of solvents to identify a solvent that might be used to preferentially stabilise a targeted porous packing arrangement. Such an approach would benefit from a more efficient global optimiser to find the lowest energy solvent configuration within each framework, especially given the failure in this work to locate the most favourable 1,4-dioxane configuration within the Form V structure. The assumption of ordered solvent within porous crystal structures should also be relaxed in future developments of the method; where multiple, energetically similar solvent configurations are possible, the solvent is likely to be disordered and the related entropic stabilisation should be considered. If the target is to desolvate the resulting structure to create a porous material, the methods should also be developed to assess the possible changes in the crystal packing motif that may occur upon desolvation, as seen in the conversion of CC1-α to CC1-α′.22
We hope that this study will motivate further development in these methods, which could have an important impact on the in silico prediction of solvated structures and, hence, the a priori design of new porous molecular crystals.
Footnotes |
† Electronic supplementary information (ESI) available. CCDC 1575910–1575914. For ESI and crystallographic data in CIF or other electronic format see DOI: 10.1039/c8fd00031j |
‡ Present Address: Defence Science and Technology Laboratory, Dstl, Porton Down, Salisbury, Wiltshire UK. |
This journal is © The Royal Society of Chemistry 2018 |