Mechanosynthesis of Higher‐Order Cocrystals: Tuning Order, Functionality and Size in Cocrystal Design

Abstract The ability to rationally design and predictably construct crystalline solids has been the hallmark of crystal engineering research. To date, numerous examples of multicomponent crystals comprising organic molecules have been reported. However, the crystal engineering of cocrystals comprising both organic and inorganic chemical units is still poorly understood and mostly unexplored. Here, we report a new diverse set of higher‐order cocrystals (HOCs) based on the structurally versatile—yet largely unexplored—phosph(V/V)azane heterosynthon building block. The novel ternary and quaternary cocrystals reported are held together by synergistic and orthogonal intermolecular interactions. Notably, the HOCs can be readily obtained either via sequential or one‐pot mechanochemical methods. Computational modelling methods reveal that the HOCs are thermodynamically driven to form and that their mechanical properties strongly depend on the composition and intermolecular forces in the crystal, offering untapped potential for optimizing material properties.

Formation of a white precipitate was observed. The reaction was then allowed to stir overnight at room temperature. Next, the reaction mixture was filtered to remove the white solids and the filtrate was transferred to a 350 mL Schlenk flask and the solvent was removed in vacuo until the mixture was about 150 mL. 11 mmol of elemental Se was added to the solution and allowed to stir overnight at room temperature. Excess unreacted Se was removed by filtration and the filtrate was collected and transferred to a 250 mL round bottom flask. The remaining solvent was removed using a rotator evaporator and the desired pure product, compound 1 was obtained as an off-white solid. 2

Single crystal X-ray diffraction studies
Diffraction-quality crystals were obtained by slow evaporation of solvent from solution in various solvents at room temperature. Data were measured at 103 K with a Bruker Kappa diffractometer equipped with a CCD detector, employing Mo K α radiation (λ = 0.71073 Å), with the SMART suite of programs. Structural solution and refinement were carried out with SHELCTL suite of programs. The structure was solved by direct methods or Patterson maps to locate the heavy atoms, followed by difference maps for the light, non-hydrogen atoms. All non-hydrogen atoms were refined with anisotropic thermal parameters. Crystal 3 and 4 crystallises with two chemically identical but crystallographically unique molecular positions, one with the presence of disordered molecules of 18crown-6 ether and 1,4-dibromotetrafluorobenzene.

Powder X-ray diffraction studies
Powder X-ray diffraction (PXRD) data were collected at 40 keV and 15 mA on a Rigaku MiniFlex 600 Benchtop diffractometer using Cu-Kα radiation (λ = 1.5418 Å) over 2 θ range of 5.0 o -50.0 o at room temperature. Baseline correction was done using Jade5.
As is common in the field, the outcome of mechanochemical experiments were determined by performing qualitative comparisons of the PXRD patterns of the ground products with those of the starting reagents or intermediate lower-order cocrystal sets. We also compared the experimental powder patterns following ball milling with the simulated PXRD patterns from the single crystal structures of the cocrystals. In all cases, we have observed that the major Bragg reflections in the ground powder samples match those predicted from the single crystal structures. However, in the absence of Rietveld refinement, we cannot be certain that minor impurities of the starting reactants or transformation products are absent.

Attenuated Total Reflectance Infrared Spectroscopy
Fourier-Transformed Infrared Radiation (FT-IR) Spectroscopic data were obtained using PerkinElmer Spectrum 100 FT-IR Spectrometer, 8 scans from 4000 cm -1 to 650 cm -1 in the transmission mode, with a resolution of 16 cm -1 .

Mechanochemical synthesis of cocrystals
All milling experiments were conducted in a 10 mL volume stainless steel milling jar that was carefully dried before use by heating at 150 o C for 2 h and cooling in a vacuum desiccator. The milling jar was then charged with a 10 mm diameter stainless steel ball bearing, Se2P2N2 (1) and the corresponding components added and the jar and sealed milled for 30 min. Loading and sealing of the milling jar was performed in a glovebox, under an inert, dry atmosphere.

Synthesis of (Se2P2N2)·(DBTFB) (B1)
The milling jar was charged with a 10 mm diameter stainless steel ball bearing, Se2P2N2 (1) and DBTFB in a 1:1 ratio added and the jar and sealed milled for 30 min. An off-white crystalline powder was obtained.

Synthesis of (Se2P2N2)2·(DMU) (B2)
The milling jar was charged with a 10 mm diameter stainless steel ball bearing, Se2P2N2 (1) and DMU in a 1:1 ratio added and the jar and sealed milled for 30 min. An off-white crystalline powder was obtained.

Synthesis of (Se2P2N2)2·(DMU)·(DBTFB) (T1)
The milling jar was charged with a 10 mm diameter stainless steel ball bearing, Se2P2N2 (1), DMU, and DBTFB in a 2:2:1 ratio and the jar and sealed milled for 30 min. An off-white crystalline powder was obtained. Diffraction quality crystals were obtained from slow evaporation of Methanol under room temperature conditions from a mixture of 2:1:1 of 1, DMU, and DBTFB. 4

1.6
Computational methods

Periodic DFT-D optimization of experimental crystal structures
The equilibrium geometry for each crystal structure (B1, B2-Form I, B2-Form II, T1, T2 and Q) was computed using a periodic dispersion-corrected density functional theory (DFT-D) model as implemented in VASP, [1][2][3] which uses the projector-augmented wave (PAW) [4] method with plane-wave basis sets and PAW pseudo-potentials. In all cases, the Perdew-Burcke-Ernzerhof (PBE) [5] generalized gradient approximation (GGA) was used for the exchange-correlational functional coupled with the D3 dispersion-correction. [6] VASP input files were generated using the cif2cell 7 program. In order to aid rapid convergence, a two-stage geometry optimization protocol was employed: initially a fixed-cell geometry optimization was performed, after which the converged structure was used as input for a full geometry optimization of both the atom positions and cell parameters. In all cases, a Γ-centered Monkhorst-Pack scheme was used to generate a tight K-point mesh with maximum K-point distance set to 2 × 0.032 Å −1 . In all calculations, the cut-off energy for the planewave basis was fixed at 500 eV.
Within each self-consistent field cycle, the convergence threshold was set at 1 × 10 −5 eV and the geometry was considered converged when all forces were below 0.03 eV Å −1 . The bulk cocrystal stabilization energy (∆ ) for each cocrystal was estimated as the difference in the DFT-D energy of the cocrystal relative to the stoichiometrically weighted sum of the energies of the reference singlecomponent crystal structures. Where reference crystal structures for the components was not available, ∆ could not be estimated. In all cases, the fully relaxed (atom positions and cell parameters) DFT-D energies of the crystal structures was used to compute ∆ .

Complexation energy calculations
The resulting optimized structure following fixed-cell geometry optimization in VASP (stage 1) was used as the starting point for estimation of the thermodynamic driving force for complex formation (∆ ). This was achieved by extracting the contents of the asymmetric unit for each VASP stage-1 optimized (see previous section for details) CIF and using GAUSSIAN09 [8] to perform a single-point DFT-D energy evaluation at the B97D/6-31G(d,p) level of theory. Complexation energies were corrected for the basis set superposition error using the counterpoise method of Boys and Bernardi. [9]

Molecular electrostatic potential (MEP) surfaces
The molecular geometries for 1 (exo-exo and exo-endo conformations), DMU and DBTFB were optimized in the gas phase at the B97D/6-31G(d,p) level of theory using GAUSSIAN09. In order to define the correct geometry for the NH donors for each conformation of 1, the dihedrals defining the relative orientation of the NH donors were constrained to their experimental values with all other degrees of freedom allowed to optimize freely. For the relatively rigid DMU and DBTFB molecules, no constraints were applied during the gas phase geometry optimization. Molecular electrostatic potentials (MEPs) were calculated at the B97D/6-31G(d,p) level of theory and visualized in GaussView 5.0. [10] Population analysis was performed using the Merz-Singh-Kollman scheme. [11][12] The atomic radius for Se was assumed to be 1.90 Å. In all cases, the equilibrium geometry for the isolated molecules were

Mechanical properties of crystals
The VASP optimized (atom positions and lattice parameters) crystal structures for B1, B2 (Forms I and II) and T1 were loaded into BIOVIA Materials Studio 8.0. [13] Crystals T2 and Q contain atom types with no available force fields in Materials Studio, thus the properties of these crystals were not modelled.
The Forcite module of BIOVIA Materials Studio 8.0 was used to perform an initial geometry optimization (atom positions and lattice parameters) using the Dreiding [14] force field supplemented with atomic charges obtained using the charge equilibration (QEq) method. [15] The convergence thresholds used were an energy of 1 × 10 −3 kcal mol -1 , maximum force of 0.5 kcal mol -1 Å -1 , maximum stress of 0.5 GPa and maximum displacement of 1.5 × 10 −2 Å. Following the geometry optimization, Forcite was used once again to compute the mechanical properties of the optimized crystal structures. This was achieved using the same energy model and same convergence thresholds specified above.

Monte Carlo binding geometries
The gas-phase optimized molecular structures of 1, DBTFB and DMU were derived at the B97D/6- configuration was evaluated using the Dreiding force field [10] complemented with the Dmol 3 derived atomic charges. The 100 most stable configurations following each Monte Carlo search were retained for visualizing the binding modes between 1 and coformer molecules.

Powder X-Ray Diffraction (PXRD) Patterns
For each route of successful formation of a binary, ternary or quaternary cocrystal, the obtained PXRD pattern is stacked against its starting materials to show an absence thereof, as well as against its simulated pattern(s) obtained from the single crystal diffraction data to show a match.
For each route of unsuccessful formation of a binary, ternary or quaternary cocrystal, the obtained PXRD pattern is stacked against its starting material to show their presence and/or against the pattern of the suspected product(s).
A simple list of all the figures is shown below:   Figure S9. Comparison of powder X-ray Diffraction (PXRD) patterns of mechanochemical combination of both obtained ternary powder of (Se2P2N2)·(DBTFB)·(DMU) with (Se2P2N2)·(18-crown-6 ether)·(KBr) and a) all obtained binary and ternary powders as labelled. As the starting matrial of this route is of a binary cocrystal, the comparison of all the other binary and ternary powders show an absence of them, hence a possible conversion to our predicted quarternary product; b) simulated pattern from SC-XRD, Roman Numerals denote the matching peaks, indicating that the bulk of the powder obtained obtained corresponds to the quaternary cocrystal, Q. The DMU which from the ternary cocrystal T1, is believed to be displaced and exists in the final powder mixture as itself; corresponding peaks have been labelled.

SUPPORTING INFORMATION
24 . Figure S18. AT-IR Spectrum of mechanochemically obtained powder of (Se2P2N2 + KBr). Figure S19. AT-IR Spectrum of mechanochemically obtained powder of ternary cocrystal of (Se2P2N2)2·(DMU)·(DBTFB) from its individual components in a 2:2:1 ratio.       For expanded view see SI  Expanded Figure 4. Fragment of the crystal structure of the ternary ICC, (Se2P2N2)·(crown)·(KBr) T2, forming one-dimensional chains linked by Se … K + and Br -… H-N interactions. The crown and K + are positionally disordered over two sites. Some C-H hydrogens are omitted for clarity. For expanded version see SI.