Co-Crystal Structures of Furosemide:Urea and Carbamazepine:Indomethacin Determined from Powder X-Ray Diffraction Data

: Co-crystallization is a promising approach to improving both the solubility and the dissolution rate of active pharmaceutical ingredients. Crystal structure determination from powder di ﬀ raction data plays an important role in determining co-crystal structures, especially those generated by mechanochemical means. Here, two new structures of pharmaceutical interest are reported: a 1:1 co-crystal of furosemide with urea formed by liquid-assisted grinding and a second polymorphic form of a 1:1 co-crystal of carbamazepine with indomethacin, formed by solvent evaporation. Energy minimization using dispersion-corrected density functional theory was used in ﬁnalizing both structures. In the case of carbamazepine:indomethacin, this energy minimization step was essential in obtaining a satisfactory ﬁnal Rietveld reﬁnement.


Introduction
In pharmaceutical science, poor solubility and a low dissolution rate in an aqueous environment can limit the usefulness of an active pharmaceutical ingredient (API). Increasingly, new chemical entities have poor aqueous solubility, and this has prompted researchers to develop strategies to overcome this limitation. Although salt formation is a widely used approach to improve solubility, it is not applicable to the many APIs that are neutral molecules. Attempts to use amorphous APIs are hindered by issues of thermodynamic instability. Crystal engineering [1] offers an alternative approach, whereby one seeks to "tune" the physicochemical properties of an API by forming a multicomponent crystal of the API with a complementary molecule called a coformer. Hundreds of coformers classified as "generally regarded as safe" are available for use in pharmaceutical formulations.
Although co-crystals have been widely studied, the definition of a co-crystal varies from study to study. The United States Food & Drug Administration states that pharmaceutical co-crystals are "crystalline materials composed of two or more different molecules, typically an active pharmaceutical ingredient and a co-crystal former ('coformers'), in the same crystal lattice", and that "from a physical chemistry and regulatory perspective, co-crystals can be viewed as a special case of solvates and hydrates, wherein the second component, the coformer, is not a solvent (including water), and is typically nonvolatile" [2]. The two (or more) molecules in the co-crystal interact by nonionic, noncovalent interactions such as hydrogen bonding, π-π stacking and van der Waals interactions, differentiating them from salts that involve ionic interactions.
Although co-crystals possess many potential academic and industrial advantages, screening for co-crystals is still challenging. With the high number of coformers available, predictive models for co-crystallization, such as the knowledge-based method developed by Wood et al. [3], are highly valuable. The actual preparation of co-crystals can be carried out by a number of different methods [4]. For example, solution-based methods, including evaporation, cooling, antisolvent, and slurry crystallization, can lead to the growth of single crystals suitable for analysis by X-ray diffraction. However, controlling the factors which affect co-crystallization using these methods is critical and the failure to obtain single crystals is not uncommon. Mechanochemical methods, including neat (dry) co-grinding and liquid-assisted grinding, offer an alternative method of manufacture [5]. The addition of a small amount of a solvent in liquid-assisted grinding increases the rate of co-crystallization by a process called solution-mediated phase transformation [6]. However, mechanochemical methods inherently produce only polycrystalline products, and as such, crystal structure determination from powder diffraction data methods [7][8][9][10] has an important role to play in the structural characterization of these products.
The very nature of co-crystals (i.e., multicomponent systems with generally large numbers of atoms in the asymmetric unit) means that they present a significant challenge to powder-based approaches, and the limited information content of a powder X-ray diffraction (PXRD) pattern dictates that particular attention has to be paid to the fine detail of the refined crystal structure, e.g., the positions of those hydrogen atoms involved in hydrogen bonding within the crystal structure.
This work reports for the first time the crystal structures of two co-crystals of pharmaceutical interest, generated as part of a program of work on the application of structure determination from powder diffraction data methods to co-crystals. One involves the use of liquid-assisted grinding to form a 1:1 co-crystal (hereafter referred to as FUR:URE) of the practically water-insoluble diuretic furosemide with the widely used coformer urea. The existence of a 1:1 form has been previously reported [11], but its crystal structure has not been determined until now. The other uses solvent crystallization to obtain a second polymorphic form of a 1:1 carbamazepine:indomethacin co-crystal (hereafter referred to as CBZ:IND form II). Of particular interest, given the complexity of the systems being studied, is the use of dispersion-corrected density functional theory (DFT-D) not only as a means of crystal structure validation, but also as a valuable part of the structure refinement process.

Materials and Methods
Furosemide (FUR), urea (URE), carbamazepine (CBZ, form III), and indomethacin (IND, γ form) were obtained from Sigma Aldrich UK. For FUR:URE, a Retsch MM400 ball mill was used to perform 60 min of liquid-assisted grinding of a 1:1 molar ratio mixture of FUR and URE in a 50 mL grinding jar, with two 5 mm stainless steel grinding balls and a few drops of acetone. The resultant material was stored for 3 weeks in a sealed glass vial at 22 • C prior to PXRD analysis. For CBZ:IND form II, a 1:1 molar ratio of CBZ form III and γ-IND was dissolved in ethyl acetate and polycrystalline material recovered via solvent evaporation.
PXRD was carried out using a Bruker D8 Advance diffractometer equipped with a LynxEye detector and monochromatic Cu Kα 1 (λ = 1.54056 Å) radiation, operating in transmission capillary mode. Powdered samples were filled in 0.7 mm borosilicate glass capillaries, and data were initially collected over an angular range of 4 to 70 • 2θ with a step size of 0.017 • . Where data for full structure determination were required, a variable-counting time scheme [12] was used. In the case of FUR:URE, the total data collection time was 18 h (see Table 1); in the case of CBZ:IND form II, the total data collection time was scaled to 42 h. Structure indexing and structure solution were performed using the DASH software package [13][14][15]. The starting molecular models used in the DASH global optimization procedure were obtained directly from the Cambridge Structural Database (CSD) [16] and are summarized in Table 2. Solved crystal structures were subjected to rigid-body Rietveld refinement [17] using TOPAS [18].  Structure indexing and structure solution were performed using the DASH software package [13][14][15]. The starting molecular models used in the DASH global optimization procedure were obtained directly from the Cambridge Structural Database (CSD) [16] and are summarized in Table  2. Solved crystal structures were subjected to rigid-body Rietveld refinement [17] using TOPAS [18]. 1 Torsional DoF = torsional degree of freedom, i.e., a torsion angle that is free to rotate, thus altering the conformation of the molecule, and whose angular value is therefore optimized during the structure determination process. Note that ring torsions and torsions that affect only the positions of hydrogen atoms are not optimized, as the former are treated as fixed and the latter do not significantly affect the calculated powder X-ray diffraction (PXRD) data.
Periodic density functional theory with van der Waals dispersion corrections (DFT-D) was used for geometry optimization of FUR:URE and CBZ:IND form II. The Perdew-Burke-Ernzerhof functional was used with projector-augmented wave pseudopotentials and the Grimme D3 correction, as implemented in the pw.x executable of the QuantumEspresso program [19,20]. The lengths of bonds involving hydrogen atoms were normalized using Mercury CSD [21], and input files for pw.x were then created from these normalized CIFs using the cif2qe script of QuantumEspresso.  Structure indexing and structure solution were performed using the DASH software package [13][14][15]. The starting molecular models used in the DASH global optimization procedure were obtained directly from the Cambridge Structural Database (CSD) [16] and are summarized in Table  2. Solved crystal structures were subjected to rigid-body Rietveld refinement [17] using TOPAS [18]. Periodic density functional theory with van der Waals dispersion corrections (DFT-D) was used for geometry optimization of FUR:URE and CBZ:IND form II. The Perdew-Burke-Ernzerhof functional was used with projector-augmented wave pseudopotentials and the Grimme D3 correction, as implemented in the pw.x executable of the QuantumEspresso program [19,20]. The lengths of bonds involving hydrogen atoms were normalized using Mercury CSD [21], and input files for pw.x were then created from these normalized CIFs using the cif2qe script of QuantumEspresso.  Structure indexing and structure solution were performed using the DASH software package [13][14][15]. The starting molecular models used in the DASH global optimization procedure were obtained directly from the Cambridge Structural Database (CSD) [16] and are summarized in Table  2. Solved crystal structures were subjected to rigid-body Rietveld refinement [17] using TOPAS [18]. Periodic density functional theory with van der Waals dispersion corrections (DFT-D) was used for geometry optimization of FUR:URE and CBZ:IND form II. The Perdew-Burke-Ernzerhof functional was used with projector-augmented wave pseudopotentials and the Grimme D3 correction, as implemented in the pw.x executable of the QuantumEspresso program [19,20]. The lengths of bonds involving hydrogen atoms were normalized using Mercury CSD [21], and input files for pw.x were then created from these normalized CIFs using the cif2qe script of QuantumEspresso. CBMZPN01 1 Indomethacin Crystals 2020, 10, x FOR PEER REVIEW 3 of 11 Structure indexing and structure solution were performed using the DASH software package [13][14][15]. The starting molecular models used in the DASH global optimization procedure were obtained directly from the Cambridge Structural Database (CSD) [16] and are summarized in Table  2. Solved crystal structures were subjected to rigid-body Rietveld refinement [17] using TOPAS [18]. 1 Torsional DoF = torsional degree of freedom, i.e., a torsion angle that is free to rotate, thus altering the conformation of the molecule, and whose angular value is therefore optimized during the structure determination process. Note that ring torsions and torsions that affect only the positions of hydrogen atoms are not optimized, as the former are treated as fixed and the latter do not significantly affect the calculated powder X-ray diffraction (PXRD) data.
Periodic density functional theory with van der Waals dispersion corrections (DFT-D) was used for geometry optimization of FUR:URE and CBZ:IND form II. The Perdew-Burke-Ernzerhof functional was used with projector-augmented wave pseudopotentials and the Grimme D3 correction, as implemented in the pw.x executable of the QuantumEspresso program [19,20]. The lengths of bonds involving hydrogen atoms were normalized using Mercury CSD [21], and input files for pw.x were then created from these normalized CIFs using the cif2qe script of QuantumEspresso. INDMET03 5 1 Torsional DoF = torsional degree of freedom, i.e., a torsion angle that is free to rotate, thus altering the conformation of the molecule, and whose angular value is therefore optimized during the structure determination process. Note that ring torsions and torsions that affect only the positions of hydrogen atoms are not optimized, as the former are treated as fixed and the latter do not significantly affect the calculated powder X-ray diffraction (PXRD) data.
Periodic density functional theory with van der Waals dispersion corrections (DFT-D) was used for geometry optimization of FUR:URE and CBZ:IND form II. The Perdew-Burke-Ernzerhof functional was used with projector-augmented wave pseudopotentials and the Grimme D3 correction, as implemented in the pw.x executable of the QuantumEspresso program [19,20]. The lengths of bonds involving hydrogen atoms were normalized using Mercury CSD [21], and input files for pw.x were then created from these normalized CIFs using the cif2qe script of QuantumEspresso. Automatic k-point sampling was used; the kinetic energy cutoffs for plane waves and charge density were 50 and 400 Ry, respectively. The convergence thresholds on total energy and forces were set to 0.0001 and 0.001 Crystals 2020, 10, 42 4 of 10 a.u., respectively. Initial geometry optimizations were carried out with lattice parameters fixed at their crystallographic values, with subsequent variable cell geometry optimizations starting from the endpoint of the fixed-cell calculations. All calculations were carried out on a Dell Precision T7810 workstation equipped with two 2.40 GHz 8-core Intel Xeon E5-2630 v3 CPUs, running the Microsoft Windows 10 operating system and using the Windows Subsystem for Linux feature to allow the Linux-compiled, MPI-enabled pw.x executable to utilize multiple cores.

Results
The crystallographic data for the FUR:URE and CBZ:IND form II co-crystals, as refined against PXRD data, are reported in Table 3.

Discussion
Both structures were initially solved with relative ease using the simulated annealing (SA) approach embodied in the DASH software. In the case of FUR:URE, the ratio χ 2 SA /χ 2 Pawley was around 2, a strong indication that the structure had been solved [13]. In the case of CBZ:IND form II, this ratio was significantly higher, at around 10, indicating that the best DASH solution was still some distance from the global minimum in χ 2 space. The initial visualization of the DASH solutions using Mercury showed that both putative crystals exhibited feasible molecular conformations and packing.
For FUR:URE, our standard approach of performing a rigid-body Rietveld refinement in TOPAS, starting from the best DASH solution, allowing positional, orientational, and torsional degrees of freedom to refine, proceeded straightforwardly, giving an excellent fit to the data (R wp = 3.50%). A 4 th order spherical harmonic correction for preferred orientation was included in this refinement. Optimization of this structure using DFT-D gave a structure that, after another round of rigid body Rietveld refinement, gave an even better fit to the data (R wp = 2.84%, Figure 1). A variable-cell DFT-D optimization of this final rigid-body refined structure resulted in only a small root mean square deviation from its starting position (0.055 Å; 15 molecule overlay performed using the crystal packing similarity feature of Mercury), a value well below the upper limit suggested by van de Streek and Neumann for a correctly determined structure [22]. That this crystal structure (Figure 2) is that of the form reported previously is evidenced by the visual agreement between our PXRD data and the PXRD data shown in Figure 9 of Reference [11].
Crystals 2020, 10, x FOR PEER REVIEW 5 of 11 Rietveld refinement, gave an even better fit to the data (Rwp = 2.84%, Figure 1). A variable-cell DFT-D optimization of this final rigid-body refined structure resulted in only a small root mean square deviation from its starting position (0.055 Å; 15 molecule overlay performed using the crystal packing similarity feature of Mercury), a value well below the upper limit suggested by van de Streek and Neumann for a correctly determined structure [22]. That this crystal structure (Figure 2) is that of the form reported previously is evidenced by the visual agreement between our PXRD data and the PXRD data shown in Figure 9 of Reference [11].  In the case of CBZ:IND form II, close inspection of the solved structure gave cause for concern. Aside from the fact that a rigid-body refinement in TOPAS did not give a satisfactory fit (Rwp = 5.30%), it was noted that one of the hydrogen atoms of the carbamazepine amide group was not involved in any hydrogen bonding (Figure 3a). Furthermore, it was noted that by rotating both this group and the indomethacin carboxylic acid by 180°, the core dimer motif was maintained whilst allowing both amide hydrogen atoms to now participate in hydrogen bonding (Figure 3b).
Given the limited resolution of the PXRD data set used for the real-space structure solution and the approximately equal scattering power of both configurations shown in Figure 3, it was deemed  In the case of CBZ:IND form II, close inspection of the solved structure gave cause for concern. Aside from the fact that a rigid-body refinement in TOPAS did not give a satisfactory fit (R wp = 5.30%), it was noted that one of the hydrogen atoms of the carbamazepine amide group was not involved in any hydrogen bonding (Figure 3a). Furthermore, it was noted that by rotating both this group and the indomethacin carboxylic acid by 180 • , the core dimer motif was maintained whilst allowing both amide hydrogen atoms to now participate in hydrogen bonding (Figure 3b).
Given the limited resolution of the PXRD data set used for the real-space structure solution and the approximately equal scattering power of both configurations shown in Figure 3, it was deemed quite likely that the DASH solution had settled in a local minimum and that the correct configuration was in fact that shown in Figure 3b. Rotation of the two groups within the crystal structure was achieved using the "mode fit" functionality of Olex2 [23], but this updated structure still did not give a satisfactory fit in a rigid-body Rietveld refinement. A DFT-D geometry optimization of this updated structure gave a structure which then gave an excellent fit (R wp = 2.92%, Figure 4) in a rigid-body Rietveld refinement that included a 2nd order spherical harmonics correction for preferred orientation.     A variable-cell DFT-D optimization of this final rigid-body refined structure resulted in only a small root mean square deviation (0.064 Å). The examination of the final refined structure ( Figure 5) showed some significant changes in the core geometry of the indomethacin molecule, distinct from the purely conformational changes that result from torsion angle optimization in DASH and TOPAS. Of particular note is a change in the geometry of the chlorobenzyl group, with (a) a bond angle N1-C9-C10 opening from 116.80 • (as seen in INDMET03) to 120.06 • and (b) the disposition of Cl1 and C9 relative to the mean plane of the C10-C11-C12-C13-C14-C15 ring being altered (see Table 4).
Crystals 2020, 10, x FOR PEER REVIEW 8 of 11 Taken together, these geometry changes affect the location of a significant amount of X-ray scattering power, explaining why the initial rigid-body Rietveld refinement based on the starting INDMET03 geometry was unsuccessful. , and that form II is the more thermodynamically stable form has been confirmed by competitive slurry with temperature cycling [25].  Taken together, these geometry changes affect the location of a significant amount of X-ray scattering power, explaining why the initial rigid-body Rietveld refinement based on the starting INDMET03 geometry was unsuccessful.
The core hydrogen bonding motifs are essentially the same in both the previously reported CBZ:IND form I [24] and CBZ:IND form II, taking the general form illustrated in Figure 3b, i.e., a central R 2 2 (8) synthon and an N-H···O hydrogen bond between the amide of the CBZ and the oxygen of the 4-chlorobenzoyl group of the IND. The density of form II (1.398 g cm −3 ) is higher than that of form I (1.351 g cm −3 ), which is indicative of more efficient packing in form II. The melting point of form II (151 • C) exceeds that of form I (144 • C), and that form II is the more thermodynamically stable form has been confirmed by competitive slurry with temperature cycling [25].
The observation of polymorphism in a co-crystal system is not in itself unusual; several authors have surveyed both the literature and the CSD for examples of polymorphic co-crystal systems, and a recent survey by Mnguni et al. [26] reported 125 such systems involving APIs. However, it is less common to find polymorphism arising from the use of both solid-state and liquid-state crystallizations, as is the case with CBZ:IND.

Conclusions
Structure determination from powder diffraction data is an extremely valuable tool for elucidating co-crystal structures, and dispersion-corrected density functional theory is an increasingly valuable adjunct [27][28][29], not only for crystal structure validation but also as part of crystal structure refinement [30]. Its value in mechanochemical synthesis is very well illustrated by the case of FUR:URE, where the crystal form was reported in 2012 but not elucidated, because of the absence of a single-crystal. The rigid-body approach to Rietveld refinement has much to commend it, as it preserves core geometric features of the high-quality input models that are ideally used when solving structures using global optimization-based methods. However, if (as was the case with CBZ:IND form II) that core geometry frustrates a successful Rietveld refinement, then DFT-D offers an objective way to alter that geometry in the context of the crystal structure under study, thus avoiding the need to employ the more subjective approach of restrained Rietveld refinement. DFT-D also adds considerable value to the positioning of hydrogen atoms involved in hydrogen bonds. Finally, it should now always be used as a crystal structure validation tool, according to the guidelines laid down by van de Streek and Neumann, along with other structure checking tools such as Mogul, PLATON, and checkCIF.