Computational modelling of solvent effects in a prolific solvatomorphic porous organic cage

A computational approach has been developed to assess the effect of solvent stabilisation on the predicted crystal structures of a porous organic cage.

Single crystal X-ray crystallography Table S1. Single crystal X-ray refinement details for 2(CC1)•7.75(CH2Cl2), 2(CC1)•11.32(CHCl3) and 2(CC1)•10(CCl4).  Refinement notes and crystal structure plots: Figure S1. Displacement ellipsoid plot from the single crystal structure, 2(CC1)•11.32(CHCl3), only ordered CHCl3 molecules, that are located in the CC1 cavity, are shown for clarity. Ellipsoids displayed at 50 % probability level. Figure S2. Displacement ellipsoid plot from the single crystal structure, 2(CC1)•10(CCl4), CC1 and CCl4 molecules shown in entirety. Ellipsoids displayed at 50 % probability level; disordered CCl4 omitted for clarity. Figure S3. Orientation of ordered CCl4 in CC1 cavity in the single crystal structure, 2(CC1)•10(CCl4). Refinement notes for, 2(CC1)•7.75(CH2Cl2): One CH2Cl2 molecules, positioned in a CC1 cavity, was severely disordered and refined with bond distance restraints (DFIX in SHELX) and constrained displacement parameters (EADP in SHELX). For a second disordered CH2Cl2 molecule the C-atom was refined with restrained displacement parameters (ISOR in SHELX). For a displacement ellipsoid plot, see Figure S5. Due to disorder, high angle X-ray diffraction was weak, despite using synchrotron radiation. A suitable resolution limit of 1.0 Å resolution limit was applied during refinement. One CHCl3, located in the CC1 cavity, was modelled over two positions, for each CHCl3 molecule, the Catoms was refined isotropically. For a displacement ellipsoid plot, see Figure S6, and for responses to the A-and B-checkCIF alerts, see below: _THETM01_ALERT_3_A Problem: The value of sine(theta_max)/wavelength is less than 0.550 Calculated sin(theta_max)/wavelength = 0.4999. Response: Due to disorder, X-ray diffraction was weak, despite using synchrotron radiation. A suitable resolution limit of 1.0 ang. was applied during refinement.

_PLAT089_ALERT_3_B
Problem: Large U3/U1 Ratio for Average U(i,j) Tensor .... 4.8 Note Response: Due to disorder, X-ray diffraction was weak, despite using synchrotron radiation. A suitable resolution limit of 1.0 ang. was applied during refinement and the number of unique reflections is lower than expected.

Refinement notes for, CC1•2(o-xylene)•CHCl3:
During refinement of the structure a 1.0 Å resolution limit was applied, and a group rigid-bond restraint was used (RIGU in SHELX). One CH2Cl2 molecule, located in the CC1 cavity, was modelled over two positions, 1,2 and 1,3 bond distance restraints were used during refinement (DFIX and DANG in SHELX). For a displacement ellipsoid plot, see Figure S7, and for responses to the A-and B-checkCIF alerts, see below: _THETM01_ALERT_3_A Problem: The value of sine(theta_max)/wavelength is less than 0.550 Calculated sin(theta_max)/wavelength = 0.4999 Response: Due to disorder, X-ray diffraction was weak, despite using synchrotron radiation. A suitable resolution limit of 1.0 ang. was applied during refinement.

_PLAT089_ALERT_3_B
Problem: Large U3/U1 Ratio for Average U(i,j) Tensor .... 4.8 Note Response: Due to disorder, X-ray diffraction was weak, despite using synchrotron radiation. A suitable resolution limit of 1.0 ang. was applied during refinement and the number of unique reflections is lower than expected.

Powder X-ray diffraction
High resolution synchrotron X-ray diffraction data were collected at beamline I11 at Diamond Light Source, Didcot, Oxfordshire. Data were collected using the position sensitive Mythen-II detector on a sample of CC1 loaded with p-xylene contained in a sealed 0.5 mm diameter borosilicate glass capillary. The capillary was spun during data collection to improve powder averaging.
Indexing was performed using TOPAS-Academic, 2 which suggested a monoclinic cell, with the space group P21 tentatively assigned. A C-centred monoclinic cell was also possible, but could not be distinguished after Le Bail fitting the available data. Structure solution was attempted using simulated annealing implemented in TOPAS-Academic. 2 The model consisted of the full CC1 molecule and two p-xylene guests occupancies varied during the optimisation. Two molecules were allowed to translate and rotate freely throughout the unit cell; the y-position of the third molecule was fixed, but the position and orientation otherwise allowed to vary freely.
No torsion angles were varied during the structure solution calculation. The simulated annealing calculation was run five times and the best solution selected for Rietveld refinement.
All five solutions were similar in terms of the positions of the cage and guest molecules, with small variations in the orientation of the p-xylene guests.
Rietveld refinement was performed with geometric restraints applied to all bond lengths and angles within the cage molecule and the xylene guests specified as rigid bodies with variable positions, orientations and occupancies. Hydrogen atoms were modelled at standard geometries and refined using the riding model. Hydrogen atoms were not included for the freely rotating methyl groups on the xylene molecules. There are two guest sites in the structure: an extrinsic sites in which the p-xylene guest is positioned with its methyl substituents close to the windows of two cages, and an intrinsic site inside the cage molecule. The extrinsic xylene consistently refined to close to full occupancy, so it was fixed at unity. The xylene guest inside the cage has a much lower refined occupancy of 0.471 (2). The isotropic displacement parameters for both guest molecules are large (Biso(Carbon) = 7-11) suggesting that the xylene molecules are poorly ordered at room temperature.  Figure S9. Crystal structure of CC1•1.47(p-xylene) determined from powder X-ray diffraction. The pxylene guest located inside the cage has a refined partial occupancy of 0.471 (2). However, the large isotropic displacement parameters for both guests suggest the xylene positions are poorly ordered.

Thermogravimetric Analysis
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.
Experimental solvent stabilisation values have been determined using TGA. The solvates were made using a 'vial in vial' technique. 10 mgs of activated CC1 was placed in an open vial, this was then placed in a larger vial with the appropriate solvent in, sealed and left at room temperature for 3 days. The TGA experiments were run immediately when the material was removed from the vial to minimise solvent loss before the analysis. Tonset was calculated using the in-built feature in the software 'TA Universal Analysis'. The TGA curves and calculated Tonset values are shown in Figure S15.

Crystal Structure Prediction (CSP)
CSP calculations taken from the Z'=1 full search set of ref 3 on the Td conformer of CC1 were used to determine the likely crystal packing arrangements of CC1 in space groups P1, P21, C2, P212121, P21212, C2221, P41212, R3, P , Cc, P21/c, C2/c, Pna21, Pbcn, Pbca, and Pnma. The CSP structures were taken and the same molecular geometry as used during the Monte-Carlo calculations was overlaid onto the structures. The structures were lattice energy minimised with an isotropic atom-atom potential using DMACRYS2.1.0. Electrostatic interactions were modeled using an atomic hexadecapole description of the molecular charge distribution taken from a B3LYP/6-31G** density obtained from Gaussian09. 4 Repulsion-dispersion interactions were modeled using the W99 potential. 5 An 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. and one conformer in the asymmetric unit. For these CSP searches, we limited ourselves to only using the experimentally determined space group for each of the given systems. Further to this, a two-step minimisation procedure was used in which a computationally less expensive monopole description of the electrostatics was used, which was obtained, by using MULFIT, 6,7 through fitting atomic partial charges to the molecular electrostatic potential produced from the full multipole model (up to hexadecapoles on each atom). A Lennard-Jones form of the W99 potential (where Lennard-Jones parameters were fitted to reproduce the exp-6 atom-atom interactions) was used during the first minimisation, during which a small pressure (10 MPa) was applied to help convergence. The pressure was removed and a second minimisation using the Buckingham form of the W99 potential was performed. For the P1, Z'=2 search a total of 10000 valid lattice energy minimisations were performed. For structures higher up the lattice energy landscape (C2/c Z'=1, C2/c Z'=2 and P21/c Z=2) a more extensive search of 50000 valid minimisations was tried. The results were clustered using COMPACK to remove duplicates and a further 2 minimisation cycles were performed on the clustered set with hexadecapoles and parameters matching those used to produce the Z'=1 CSP set. The results are collected in Table S3.

Available Volume Calculations
The available volume was calculated using zeo++ with a channel and probe radius of 1.2 Å and 100000 Monte-Carlo samples per unit cell. 8,9

Monte-Carlo Solvent Insertion Calculations
The experimental CC1 structures were taken, expanded to P1, lattice energy minimised using DMACRYS. 10 Monte Carlo (MC) simulations were performed to determine the role of solvation in stabilising the packing. MC simulations using Towhee-7.10 11 in the NVT ensemble at 5000 K were performed on unit cell systems containing CC1:N Solvent ratio (where N=1, 2, 3 ..) with a force field derived from a combination of ab-initio, UFF 12 and W99 5 parameters with each simulation consisting of 10000 MC cycles using the Monte-Carlo moves of configurational-bias single box molecule reinsertion move (pm1boxcbswap), configurational-bias partial molecule regrowth (pmcb), intramolecular single atom translation move (pmtraat), center-of-mass molecule translation move (pmtracm) and rotation about the center-of-mass move (pmrotate). Additionally, intrabox two molecule switch based upon the center of mass positions (pm1boxcomswitch) was used for systems involving more than one solvent type. Examples of the distribution of move types can be found in Tables S3 and S4. Frames sampled at 10 cycle intervals were used as an input for minimisation with DMACRYS2.1.0 with a 30 Å cutoff for the real-space component of the Ewald summation. Each minimisation cycle consisted of a point charge minimisation using CHELPG charges derived from a Gaussian09 B3LYP/6-31G** calculation, 4 followed by a multipole minimisation with atomcentered multipoles up to the hexadecapole level, derived from a distributed multipole analysis of the B3LYP/6-31G** electron density. A secondary minimisation was performed to ensure the structure was a stable minimum. All valid structures were compared against the desolvated experimental framework using the COMPACK algorithm 13 using molecular clusters consisting of 30 molecules and an interatomic distance tolerance of 20 % in order to eliminate those structures in which there had been a significant change in the CC1 framework (RMSD30 > 0.8 Å). The energy of the solvent (ESolvent) was calculated from configurations arising from MC equilibrated, solvent boxes containing 50 molecules fixed at the experimental density. Simulations were in the NVT ensemble (at 300 K) and consisted of 100000 MC cycles with an equal mix of pmtracm and pmrotate moves with a target acceptance rate of 50 %. Frames sampled every 50 MC cycles from the last 50000 MC cycles were subject to a constant volume DMACRYS minimisation with atom-centered multipoles up to the hexadecapole level (derived from a distributed multipole analysis of the B3LYP/6-31G** electron density). Intermolecular interactions were modeled in a similar manner to previous calculations only this time using a 15 cutoff. The results of these calculations are collected in Table S5. ESolvent, in conjunction with a factor of (3/2)RT at T=300 K [to account for the change in internal energy of the solvent in the liquid and crystalline state] was used to compute a lattice energy of CC1 (including the additional stabilisation gained by CC1 through interaction with the solvent) in the solvated system according to equation (1).

ELatt=[ELatt(CC1+Solvent)+ ECC+N×(-ESolvent+(3/2)RT)]
(1) Where ECC=0.5*17.25 kJ mol -1 , which is a correction (where applicable) due to the difference in intramolecular energy of the Td and C3 CC1 conformations.      Figure S12. Overlays of the computed (blue) and experimental (red) solvated structures with the solvent represented by a sphere at the solvent centroid position. Figure S13. Z'=1 CSP set landscape from Case et al. 3 The colour bar gives the accessible volume per CC1 (Å 3 ) (as calculated using a 1.2 Å probe). Figure S14. CSP energy landscape for CC1 (grey points) also showing the calculated latticed energies of the desolvated structures (unfilled circles), after solvent stabilisation (filled circles) and the CSP global minimum structure solvated with the corresponding solvent (stars). Additionally, the experimentally determined α′ and β′ (squares) are given. Table S10. Comparison of experimentally measured thermal stability (Tonset -Tbp) and the calculated stabilisation energy (value in brackets from the lattice energy minimised crystal structure).