Experimental Electron Density and Neutron Diffraction Studies on the Polymorphs of Sulfathiazole.

High resolution X-ray diffraction data on forms I-IV of sulfathiazole and neutron diffraction data on forms II-IV have been collected at 100 K and analyzed using the Atoms in Molecules topological approach. The molecular thermal motion as judged by the anisotropic displacement parameters (adp's) is very similar in all four forms. The adp of the thiazole sulfur atom had the greatest amplitude perpendicular to the five-membered ring, and analysis of the temperature dependence of the adps indicates that this is due to genuine thermal motion rather than a concealed disorder. A minor disorder (∼1-2%) is evident for forms I and II, but a statistical analysis reveals no deleterious effect on the derived multipole populations. The topological analysis reveals an intramolecular S-O···S interaction, which is consistently present in all experimental topologies. Analysis of the gas-phase conformation of the molecule indicates two low-energy theoretical conformers, one of which possesses the same intramolecular S-O···S interaction observed in the experimental studies and the other an S-O···H-N intermolecular interaction. These two interactions appear responsible for "locking" the molecular conformation. The lattice energies of the various polymorphs computed from the experimental multipole populations are highly dependent on the exact refinement model. They are similar in magnitude to theoretically derived lattice energies, but the relatively high estimated errors mean that this method is insufficiently accurate to allow a definitive stability order for the sulfathiazole polymorphs at 0 K to be determined.

Introduction

The polymorphism[1] of the antimicrobial drug 4-amino-N-(thiazol-2-ylidene)-benzene sulfonamide 1 (trivial name sulfathiazole) has been extensively studied and is now a classic example of this phenomenon. To date, five crystalline polymorphs of unsolvated 1 are fully characterized by single-crystal X-ray diffraction[2−10] and over 100 crystalline solvates.[11] Although compound 1 can exist as imino (1a) and amido tautomers (1b), in the crystal phase it is exclusively found as the imino tautomer 1a. The literature on the polymorphism of 1 can be confusing and contradictory, especially as different authors have used different polymorph numbering schemes. Herein we adopt the CSD[12] enumeration scheme, and to avoid any possible confusion, this scheme is detailed in Table S1 (Supporting Information), together with the standard unit cell and reduced cell constants. A detailed overview of the preparation and characterization of the polymorphs of 1 has been recently been given by Nagy and co-workers,[13] in which they emphasize the difficulties of reproducibly obtaining the pure polymorphic forms by crystallization techniques. All forms crystallize in the monoclinic space group P21/c, although for convenience the alternative setting of P21/n is used in the literature for forms IV and V.[2−10] Forms II and IV crystallize with one molecule in the asymmetric unit, and the remaining forms have two molecules.

A large number of wide-ranging studies on the polymorphism of 1 have been undertaken.[14] Even answering the vexed, if important, question of the relative thermodynamic stabilities of the polymorphs of 1 has proved a challenge,[15−18] due in part to the interconversion between the phases, but also because of their closely similar energies. It is now clear that the order of thermodynamic stability is quite temperature-dependent, and the inconsistencies in the literature regarding this matter have been pointed out recently by Croker and co-workers.[15] Their own solubility and differential scanning calorimetric measurements suggest that in the temperature range 283–323 K the stability order is I < V < IV < II < III. Above 373 K, the relative order changes and form I becomes progressively more stable. Close to the melting points, the relative order is II < III < IV < V < I, that is, almost a complete reversal of the room temperature order.[15]

In this study, we report high-resolution single crystal X-ray studies on forms I–IV and single crystal neutron diffraction studies on forms II–IV at 100 K. We were unable to obtain single crystals of form V of sufficient quality to merit similar investigations on this polymorphic form. Multipole refinements and topological analysis of the resulting density models were undertaken using the experimental X-ray structure factors and theoretical static structure factors obtained from periodic density functional theory (DFT) calculations, and lattice energies were derived from the multipole populations. These results agree reasonably well with theoretical calculations on the lattice energies and intermolecular interaction energies but are insufficiently accurate to determine the relative polymorph stabilities.

Experimental Section

X-ray Diffraction

Compound 1 was obtained from commercial sources and recrystallized according to literature procedures[2−10] to afford the various polymorphs. As discussed previously,[13] these procedures do not consistently afford the pure polymorphs, and in our hands we also found this to be the case. In some cases, we observed concomitant crystallization, and crystals had to be separated by hand and their polymorphic form confirmed by X-ray diffraction. Single crystal X-ray diffraction data were collected near 100 K on a Bruker-Nonius KappaCCD (forms I, II, and IV) or Bruker AXS APEX-II (form III) diffractometer, running under Nonius Collect[19a] or APEX-2[19b] software. The oscillation axis was either the diffractometer ω- or φ-axis. After collection of the low-order reflection data, the same scan-sets were repeated with 1/10th of the exposure time. These “rapid” images were used to record the intense low-order data more accurately, including all reflections that were overexposed in the first set of images. The unit cell dimensions used for refinement purposes were determined by postrefinement of the setting angles of a significant portion of the data set, using the Scalepack[20] or SAINT software.[19b] The KappaCCD frame images were integrated using Denzo(SMN)[20] with a sufficiently large spot size (<0.85) to account for the Kα1-α2 splitting, which becomes quite significant at θ ≈ 50°. The resultant raw intensity files from Denzo(SMN) were processed using a locally modified version of DENZOX.[21] The APEX-II images were integrated and processed using the SAINT software.[19b] An absorption correction by Gaussian quadrature[22] was then applied to the reflection data (except for form III). A second semiempirical correction[23] (without a theta-dependent correction) was then applied to remove any residual absorption anisotropy due to the mounting medium and account for other errors such as machine instabilities. The data were scaled and merged using SORTAV[24] to provide a set of unique reflections without systematic absences. All data sets were essentially complete and highly redundant. A spherical atom refinement using SHELXL-2013[25] was initially undertaken, with full-matrix least-squares on F2 and using all the unique data. All non-H atoms were allowed anisotropic thermal motion. Details of these refinements are given in Table 1. Thermal ellipsoid plots were obtained using the program ORTEP-3 for Windows.[26] All calculations were carried out using the WinGX package[26] of crystallographic programs. Details of the data collection and processing procedures are given in Table 1.

Table 1

Experimental Details for the X-ray Diffraction Studiesa

compound formula	C9H9N3O2S2	C9H9N3O2S2	C9H9N3O2S2	C9H9N3O2S2
polymorphic form	I	II	III	IV
Mr	255.31	255.31	255.31	255.31
space group	P21/c	P21/c	P21/c	P21/n
crystal system	monoclinic	monoclinic	monoclinic	monoclinic
a/Å	10.5235(2)	8.1904(2)	17.4174(6)	10.7891(2)
b/Å	12.9016(2)	8.5345(2)	8.4911(3)	8.4836(1)
c/Å	17.2177(3)	15.4497(3)	15.4952(5)	11.3978(2)
β/deg	107.834(1)	94.155(1)	112.761(2)	91.643(1)
V/Å–3	2225.32(7)	1077.11(4)	2113.17(13)	1042.82(3)
Z	8	4	8	4
Dcalc/g cm–3	1.524	1.574	1.605	1.626
λ/Å	0.71073	0.71073	0.71073	0.71073
μ/mm–1	0.466	0.482	0.491	0.497
temperature/K	105(2)	100(2)	100(2)	100(2)
crystal size/mm	0.56 × 0.35 × 0.13	0.43 × 0.24 × 0.19	0.55 × 0.30 × 0.20	0.67 × 0.21 × 0.16
θ range/deg	2.009–47.962	2.493–50.072	1.27–57.88	2.563–50.588
max sin(θ)/ λ	1.045	1.079	1.191	1.09
no. of data used for merging	456552	331551	513654	309882
no. of unique data	21249	11293	29858	11135
hkl range	–21 ≤ h ≤ 21	–17 ≤ h ≤ 20	–41 ≤ h ≤ 41	–23 ≤ h ≤ 23
	–26 ≤ k ≤ 26	–14 ≤ k ≤ 25	–20 ≤ k ≤ 20	–18 ≤ k ≤ 18
	–35 ≤ l ≤ 35	–53 ≤ l ≤ 51	–36 ≤ l ≤ 36	–24 ≤ l ≤ 24
Rint	0.041	0.035	0.049	0.032
Rσ	0.034	0.026	0.031	0.026
Spherical atom refinement
no. of data in refinement	21249	11293	29858	11135
no. of refined parameters	361	181	362	181
final R [I > 2σ(I)] (all data)	0.032(0.048)	0.026 (0.029)	0.033 (0.043)	0.025(0.029)
Rw2 [I > 2σ(I)] (all data)	0.087(0.092)	0.084 (0.083)	0.099 (0.103)	0.075(0.077)
goodness of fit S	1.03	1.139	1.098	1.08
extrema in residual map/e Å–3	1.29 → −0.69	0.79 → −0.54	0.83 → −0.54	0.62→ −0.80
max shift/esd in last cycle	2.0 × 10–3	6.0 × 10–3	4.0 × 10–3	2.0 × 10–3
Multipole refinement
no. of data in refinement	17353	10566	25834	10172
no. of refined parameters	687	495	688	495
final R [I > 3σ(I)] (all data)	0.0242 (0.0448)	0.0177 (0.0214)	0.0258 (0.0414)	0.0141 (0.0191)
Rw all data	0.0217	0.0249	0.037	0.0172
goodness of fit S	1.1094	1.6495	1.984	1.1034
extrema in residual map/e Å–3
(all data)	0.59 → −0.27	0.79 → −0.32	0.33 → −0.54	0.22 → −0.20
(truncated to 0.8 e Å–1)	0.45 → −0.14	0.54 → −0.26	0.24 → −0.38	0.16 → −0.16
max shift/esd in last cycle	1.0 × 10–5	1.0 × 10–5	4.0 × 10–5	<1.0 × 10–5

R = Σ (|Fo| – |Fc|)/Σ (Fo). Rw = {Σ(w(Fo – Fc)2)/Σ(w(Fo)2)}1/2. Rw2 = {Σ(w(Fo2 – Fc2)2)/Σ(w(Fo2)2)}1/2. Rσ = Σ [σ(Fo2)]/Σ [Fo2]. Rint = Σ{n/(n – 1}1/2| Fo2 – Fo2 (mean)|/ΣFo2 (summation carried out when more than one symmetry equivalent is averaged).

Neutron Diffraction

Neutron diffraction data were collected for forms II, III, and IV at 100 K on the SXD instrument[27] at the ISIS spallation neutron source, using the time-of-flight (TOF) Laue diffraction method. Reflection intensities were reduced to structure factors (SHELX style HKLF 2) using standard SXD procedures, as implemented in the computer program SXD2001.[28] No absorption corrections were deemed necessary. The unit cells used for refinements with SHELXL-2013[25] were taken from the X-ray determinations. Anisotropic displacement parameters were used for all atoms, including the H atoms. In the final cycles, the merged and extinction-corrected reflection data obtained from a converged SHELXL refinement was used as the input reflection file. Experimental details for the three structural determinations are given in Table 2.

Table 2

Experimental Details for the TOF Single-Crystal Neutron Diffraction Studiesa

polymorphic form	II	III	IV
space group	P21/c	P21/c	P21/n
crystal system	monoclinic	monoclinic	monoclinic
a/Å	8.1904(2)	17.4174(6)	10.7891(2)
b/Å	8.5345(2)	8.4911(3)	8.4836(1)
c/Å	15.4497(3)	15.4952(5)	11.3978(2)
β/deg	94.155(1)	112.761(2)	91.643(1)
V/Å–3	1077.11(4)	2113.17(13)	1042.82(3)
λ/Å	0.42–7.64	0.42–7.64	0.42–7.64
μ/mm–1	0.00/0.00	0.00/0.00	0.00/0.00
temperature/K	100(2)	100(2)	100(2)
crystal size/mm	2 × 2 × 8	2 × 3 × 6	2 × 2 × 2
θ range/deg	8.50–84.54	8.18–85.45	8.77–83.30
max sin(θ)/λ	1.92	2.15	1.45
no. of data used for merging	7244	8375	3667
no. of unique data	3645	3848	1686
hkl range	–17 ≤ h ≤ 20	–43 ≤ h ≤ 52	–20 ≤ h ≤ 27
	–14 ≤ k ≤ 25	–12 ≤ k ≤ 20	–20 ≤ k ≤ 13
	–53 ≤ l ≤ 51	–59≤ l ≤ 34	–31 ≤ l ≤ 28
Rint	0.216	0.181	0.143
Rσ	0.155	0.129	0.124
no. of data in refinement	3645	3848	1686
no. of refined parameters	226	451	226
final R [I > 2σ(I)] (all data)	0.061 (0.094)	0.070 (0.072)	0.072(0.072)
Rw2 [I > 2σ(I)] (all data)	0.115 (0.120)	0.153 (0.154)	0.145(0.145)
goodness of fit S	1.072	1.074	1.097
extrema in residual map/fm Å–3	2.03 → −1.80	1.94 → −1.94	2.08 → −1.40
max shift/esd in last cycle	<1.0 × 10–3	<1.0 × 10–3	<1.0 × 10–3

R = Σ(|Fo| – |Fc|)/Σ(Fo). Rw = {Σ(w(Fo – Fc)2)/Σ(w(Fo)2)}1/2. Rw2 = {Σ(w(Fo2 – Fc2)2)/Σ(w(Fo2)2)}1/2. Rσ = Σ[σ(Fo2)]/Σ [Fo2]. Rint = Σ{n/(n – 1}1/2|Fo2 – Fo2(mean)|/ΣFo2 (summation is carried out only where more than one symmetry equivalent is averaged).

Multipole Refinement

The multipole formalism of Hansen and Coppens[29] as implemented in the XD-2006 program suite[30] was used.

The function minimized in the least-squares procedure was Σw(|Fo| – k|Fc|)2, with only those reflections with I > 3σ(I) included in the refinement. The multipole expansion was truncated at the hexadecapole level for the S atoms, the octupole level for the O, N, and C atoms, and the quadrupole level for the H atoms. The importance of employing anisotropic displacement parameters (adp’s) for H atoms in multipole refinements has been emphasized recently by several authors.[31] The method of Madsen[32a] (SHADE2 program[32b]) is known to provide an excellent approximation to H-atom adp’s.[33] For those forms of 1 for which neutron diffraction data were available, models using H atom adp’s from the SHADE2 calculated or the scaled[34] experimental neutron values were carefully compared, and little difference was found. In the final refinements for forms II–IV, the scaled experimental neutron adp’s were used with fixed contributions, while for form I, the SHADE determined parameters were used in the same way. Each pseudoatom was assigned a core and spherical-valence scattering factor derived from the relativistic Dirac–Fock wave functions of Su and Coppens[35] expanded in terms of the single-ζ functions of Bunge, Barrientos and Bunge.[36] The radial fit of these functions was optimized by refinement of the expansion-contraction parameter κ. The valence deformation functions for the C, O, and H atoms used a single-ζ Slater-type radial function multiplied by the density-normalized spherical harmonics. The radial fits for the chemically distinct atoms were optimized by refinement of their expansion–contraction parameters κ, κ′. The radial functions used for the S atoms were those recommended by Dominiak and Coppens.[37] For all forms, anharmonic motion for the sulfur atoms was modeled with third- and fourth-order Gram-Charlier components. Examination of the resultant pdf’s using the XDPDF module in XD2006 indicated that this modeling was physically reasonable.

The successful deconvolution of thermal motion was judged by the Hirshfeld rigid bond criterion.[38] Scatter plots of the scale factor Fobs/Fcalc against sin(θ)/λ were flat across almost the entire resolution range, while difference Fourier maps and residual density analysis[39] showed that, for forms III and IV, essentially no unmodeled features remained in the data (Figures S1–S3, Supporting Information). For forms I and II however, some positive residuals were observed, which could be attributed to minor disorder or possibly minor twinning. For form I, there is a single residual peak close to the inversion center at (0.5, 0.5, 0) with an integrated electron density of ∼0.15 e, while for form II there are two peaks that could be attributable to the S atoms from a second orientation. There are no previous reports of disorder in any polymorph of 1, but if it is present for forms I and II, then it is certainly present at no greater than the 1–2% level. There is no evidence in the neutron diffraction study for any disorder in II, and no disorder model was used in the final refinements of the X-ray data. A careful comparison of the resultant multipole models and topological properties shows no evidence of discrepancies between forms I and II and forms III and IV which could be attributed to any putative disorder (see below).

Theoretical Studies

Gas-phase DFT calculations on 1 were undertaken with the program GAUSSIAN09[40] at both the experimental and optimized geometries, using several functionals and basis sets. We quote here the results from the optimized geometry with the M02-2X functional[41a] and using the def2-TZVP[41b] basis for all atoms. Fully periodic B3LYP[42] calculations based on the experimental crystallographic parameters were also undertaken with CRYSTAL09,[43] using standard 6-31G** basis sets for all centers. Lattice energies were calculated from the multipole models using the XD2006 suite[30] and also using the Gavezzotti atom–atom Coulomb–London–Pauli (AA-CLP)[44a] and PIXEL[44b] methods. Topological analysis on the gas-phase wave function was undertaken using the AIMALL[45a] software package, while Hirshfeld surface analysis was conducted using the CrystalExplorer[45b] program.

To obtain a benchmark interaction energy, an MP2 complete basis set (CBS) limit was estimated using results from aug-cc-pVTZ[46a,46b] and aug-cc-pVQZ[46a,46b] basis sets. The HF component of the energy was extrapolated using the formula of Karton,[46c] and the correlation energies were extrapolated using the formula of Halkier.[46d] All of these calculations included the counterpoise correction.[46e] A correction for CCSD(T) correlation was also performed at the DF-LCCSD(T0)[46f]/aug-cc-pVDZ level, comparing the correlation energy with that from DF-LMP2[46g]/aug-cc-pVDZ. The difference in correlation energy was then added to the MP2/CBS total energies, and the resultant energies are referred to as CBS(T). Counterpoise corrections were not used, as the local correlation methods are BSSE free. All calculations were undertaken using the MOLPRO program.[46h]

Results and Discussion

Molecular Structure

The molecular structure, crystal packing, and H-bonding intermolecular interactions in polymorphs I–IV of 1 are very well established,[2−10] and only salient features will be discussed here. The two independent molecules for forms I and III are denoted here as Ia/b and IIIa/b respectively. Thermal ellipsoid plots for Ia and Ib are shown in Figure 1, and plots for the other forms are in Figures S3–S7 (Supporting Information). The molecular conformations in forms II–IV are essentially identical, while Ia and Ib have a different conformation, as is clearly visible from the overlay plot shown in Figure 2 and the characterizing torsion angles given in Table 3. These differences arise from differing torsions within the two rings and also a ∼180° rotation of the NH2 group. The differing molecular conformations in form I result in a distinct H-bonding arrangement and crystal packing compared with the other forms, as was previously discussed in detail by Blagden et al.[14u] and Gelbrich et al.[2] In all forms, however, the sulfur atom of the thiazole ring lies close to one of the oxygen atoms of the sulfone group, with an S12···O11 distance in the range 2.8734(4)–2.9848(5) Å, significantly shorter than the sum (3.3 Å) of their van der Waals radii. The gas-phase structure (optimized from the experimental geometry, Figure S8, Supporting Information) has the same short S···O=S interaction (2.74 Å). This intramolecular interaction is clearly detectable in the topological analysis of the electron density in all forms, and its importance in “locking” the molecular conformation in 1 is discussed in further detail below.

Figure 1

ORTEP plots of molecules Ia (top) and Ib (bottom) showing the atomic labeling scheme, with thermal ellipsoids drawn at the 50% probability level. Anisotropic thermal parameters for H atoms were calculated using the SHADE[32] procedure.

Figure 2

Molecular best-fit overlay plot (a) forms I–IV and (b) forms II–IV. Color coding is form I – mol. 1, green, mol. 2 yellow; form II, red; form III – mol. 1, purple, mol. 2 pink; form IV, blue.

Table 3

Experimental and Theoretical Torsion Angles (°) for Sulfathiazole

form	O11–S11–C111–C116	O11–S11–N10–C11	S11–N10–C11–S12
Ia	10.12(3)	–33.56(3)	–7.83(3)
Ib	15.34(3)	–39.63(3)	0.13(3)
II	–6.16(2)	–39.43(2)	17.76(2)
IIIa	–6.86(3)	–40.04(3)	19.53(3)
IIIb	–6.17(3)	–37.08(3)	15.09(3)
IV	–7.22(2)	–37.16(2)	15.05(2)
Vaa	–17.5(3)	–17.6(3)	8.4(5)
Vba	–14.1(3)	–37.1(3)	5.5(4)
theoreticalb	21.93	5.71	–0.62

Taken from ref (7).

Conformation 1–0 (see text).

Analysis of Thermal Parameters

In principle, neutron diffraction studies should provide more accurate anisotropic displacement parameters (adp’s) than conventional X-ray studies, though in practice this is not always the case, for well-known reasons.[34] On the basis of the Hirshfeld rigid bond criterion,[38] the X-ray derived adp’s (after multipole refinement) appear slightly more reliable. For form II, the average mean-square displacement amplitudes (Δ-msda’s) are 0.0019 (0.0026) Å2 for X-ray (neutron), and for form IV they are 0.0017 (0.0020) Å2 for the atoms with anisotropic thermal parameters. Several features are consistently observed for all polymorphs, in both the X-ray and neutron refinements. The adp of the thiazole sulfur atom is consistently larger and more anisotropic than that of the sulfone sulfur atom. The elongation perpendicular to the ring is clearly visible in Figure 1 and Figures S4–S7 (Supporting Information) [mean Ueq-thiazole/Ueq-sulfone = 1.73, mean λ1/λ3 for thiazole and sulfone are 3.30, 1.64 respectively]. This feature of the thiazole S atom is also observed in some other structures containing thiazole rings, for which adp’s are available (see Figure S9, Supporting Information). To ascertain whether this was due to genuine thermal motion, rather than a disguised disorder (and nonplanarity) in the thiazole ring, variable temperature X-ray diffraction data on form II in the temperature range 100–200 K were obtained. The same crystal was used for all data sets and the resolution (sin θ/λ = 0.7 Å–1) and data processing protocols were identical. Analysis of the U tensors of both sulfur atoms showed a linear dependence of components between 200 and 100 K (Figure S10, Supporting Information). This is the expected behavior for a quantum oscillator[47] and suggests that the adp for the thiazole S atom is representative of true thermal motion rather than a convolution of static or dynamic disorder. In our case, the msda’s are not zero when extrapolated to 0 K, presumably due to systematic errors in the adp’s (possibly due to some anharmonicity).[47b] While the thiazole ring in the imino tautomer 1a is not formally aromatic, the sp2 hybridization of all three C atoms imposes planarity on the ring, and this is indeed observed experimentally and in all calculations. The contributions to the atomic displacement parameters from the computed frequencies of the internal modes are shown in Figure S11 (Supporting Information). They are consistent with an extended thermal motion of the S atom perpendicular to the ring, but it should be stressed that for heavy atoms such as sulfur, the internal modes make quite minor contributions to the observed adp, which is dominated by the low frequency modes.

Another obvious feature of the adp’s for all non-H atoms in the X-ray structures is their general similarity in all the polymorphic forms, suggesting the observed adp’s reflect molecular vibrations that are hardly differentiated by the potential field due to the crystal packing. Whitten and Spackman[31a] have compared adp’s using a quantitative similarity index S12 = 100(1 – R12), where R12 measures the overlap between the probability densities functions described by two displacement tensors U1 and U2. This procedure has been coded into the SimADP routine of WinGX,[26] in which the two molecules or fragments to be compared are rotated to minimize the discrepancy in their positional coordinates, and the orthogonalized U tensors are then compared. Since S12 does not provide a direct measure of the relative orientation of the eigenvectors of the tensors (particularly if the tensors are close to isotropic), a combined figure of merit (FOM) based on the similarity index S12, Reigval (an R value based on the magnitudes of corresponding eigenvalues) and RMSeigvec (the RMS angles between corresponding eigenvectors) is also computed. A perfect fit gives a FOM of 0.0 and a value less than 0.05 indicates a very close fit. Table 4 shows values of these comparative indices for molecules II and IV. Further comparative data are given in Tables S2–S11 (Supporting Information). These results provide quantitative confirmation that corresponding atomic adp’s in forms II, III, and IV are extremely similar—the noticeably worse fits associated with comparisons involving form I may be due in part to its different conformation.

Table 4

Comparative adp Indices of Non-H Atoms for Forms II and IV

atom	RMSeigvec	Reigval	S12	FOM
S11	16.8094	10.3505	0.4568	0.0921
S12	5.8068	8.2924	0.2728	0.0479
O11	10.0567	4.0388	0.2715	0.0479
O12	20.1793	6.9539	0.6267	0.0925
N10	36.3814	2.8178	0.2866	0.1316
N11	35.6973	4.4254	0.3143	0.1348
N15	9.0122	3.2269	0.1708	0.0414
C11	18.3842	1.6060	0.2603	0.0675
C13	5.4957	7.2897	0.2339	0.0434
C14	8.3442	5.3602	0.1624	0.0462
C111	4.4141	4.7468	0.1006	0.0309
C112	6.5681	5.4185	0.1149	0.0403
C113	4.2088	6.9566	0.1099	0.0376
C114	13.8398	6.9966	0.1574	0.0700
C115	8.9280	4.0102	0.1363	0.0436
C116	22.4103	2.9482	0.1058	0.0849

Hirshfeld Surface Analysis

The Hirshfeld surface,[48] defined as the surface where wA(r) = ρpromolecule(r)/ρprocrystal(r) = 0.5, has been proposed as a useful graphical tool for examining the differing intermolecular interactions in polymorphic systems.[49] One such example of its use in characterizing the differences in polymorphs is the case of carbamazepine.[50] This surface is, of course, just one of many possible molecular surfaces, but it has the unique property of encoding the intermolecular contacts from its very definition. The Hirshfeld surfaces for all individual molecules of 1 are shown in Figure 3, where the normalized contact distance (dnorm[49]) is color-mapped onto the surface. Red points mark intermolecular atomic contacts shorter than, while blue points mark contacts longer than the van der Waals contact. It is immediately visually obvious that the sulfathiazole molecules in forms II–IV have similar intermolecular contacts, while those in form I are different from each other and also from forms II–IV. The so-called fingerprint plots,[51] where the di and de distribution of surface points is shown as a frequency-coded scatter plot, are capable of providing useful quantitative information about the relative proportions of differing types of intermolecular contact. The fingerprint plot for form IV is shown in Figure 4, specifically broken down[52] into contributions from the C···H, H···H, N···H, O···H, and S···H contacts which together make up 91.2% of the surface. The fingerprint plots for all forms are shown in Figure S12 (Supporting Information). The differing types of H-bonding intermolecular interactions are indicated as types 1–8 and average distances given in Table S12, Supporting Information. Once again, the close similarities in the intermolecular interactions in forms II–IV compared with form I are clearly obvious from this graphical representation. In particular, although the N–H···O and N–H···N hydrogen bonding networks in the various polymorphs of 1 have been discussed in detail with reference to the different packing arrangements,[2,14u] it is clear from these fingerprint plots that the C···H and H···H contacts also make up a significant portion of the Hirshfeld surface and so may make a considerable contribution to the intermolecular interaction energies and lattice energies.

Figure 3

Hirshfeld surface plots of forms I–IV mapped with dnorm.[48]

Figure 4

Fingerprint plots of form IV showing the contributions from (a) all intermolecular contacts, (b) C···H contacts (22.4%), (c) H···H contacts (26.2%), (d) N···H contacts (11.2%), (e) O···H (22.2%), and (f) S···H contacts (9.2%). de and di represent the distances from the surface to nearest external and internal atoms, respectively. The color coding gray-blue-cyan represents increasing numbers of surface contributors at individual de/di points.

Topological Analysis of Electron Density

The electron density in the six independent experimental molecules (Ia, Ib, II, IIIa, IIIb, and IV) and the theoretical optimized structures were subjected to a Quantum Theory of Atoms in Molecules (QTAIM) analysis. The topological properties at the bond critical points (bcp’s) for a few selected bonds are given in Table 5, and a full list is given in Table S13 (Supporting Information). Since in this study we have effectively six different experimental determinations of these properties, the sample mean and standard deviation were computed and are also given in the tables. The experimental molecular graphs are all homeomorphic, and a representative plot for form IV is shown in Figure 5a. Bond paths and associated bcp’s were detected for all the conventional covalent bonds. The experimental electron densities at the bcp’s for all forms agree surprisingly well, and the theoretical value is within three standard deviations of the experimental sample mean in virtually all cases. As is commonly observed, the Laplacian values, especially for the polar covalent bonds, show slightly greater discrepancies between theory and experiment.[53] This is primarily due to differences in λ3, arising from the slightly differing positions of the bcp along the bond path. In particular for the S(11)–N(10) bond, the magnitude of the Laplacian is very sensitive to the position of the bcp, and there is a linear dependence of the Laplacian on the percentage displacement of the bcp along the bond vector (see Figure S13, Supporting Information).

Table 5

Representative Distances of Bond Critical Points from Nuclei, Density, and Laplacian of Density at bcps and Eigenvalues of Hessian at bcp’sa

bond	d1b	d2b	ρ(rb)c	∇2ρ(rb)d	λ1d	λ2d	λ3d
S(11)–O(11)	0.5714	0.8736	2.08	19.34	–12.74	–12.67	44.76
	0.5712	0.8729	2.09	19.51	–12.79	–12.68	44.98
	0.5750	0.8751	1.99	19.30	–12.23	–10.31	41.83
	0.5743	0.8745	2.06	18.06	–12.83	–11.41	42.31
	0.5748	0.8757	2.06	17.71	–12.78	–11.33	41.83
	0.5736	0.8789	2.13	15.96	–12.92	–12.54	41.43
	0.5698	0.8728	2.06	24.15	−13.36	−13.07	50.57
	0.5734	0.8751	2.07	18.31	−12.72	−11.82	42.86
	0.0015	0.0019	0.04	1.25	0.22	0.88	1.48
S(11)–N(10)	0.5947	1.0141	1.63	7.82	−7.71	−7.51	23.04
	0.5971	1.0223	1.62	6.51	−7.58	−7.37	21.46
	0.6099	0.9865	1.84	−3.00	−10.68	−7.78	15.45
	0.6249	0.9725	1.86	−8.83	−10.00	−8.82	9.99
	0.6318	0.9762	1.85	−10.58	−9.88	−8.72	8.01
	0.5951	1.0132	1.72	4.86	−8.53	−8.09	21.48
	0.6416	0.9902	1.64	−7.55	−9.86	−9.17	11.48
	0.6089	0.9975	1.75	−0.54	−9.06	−8.05	16.57
	0.0148	0.0197	0.10	7.35	1.19	0.56	5.88
C(111)–S(11)	0.8119	0.9369	1.50	–10.98	–8.83	–7.65	5.50
	0.8117	0.9361	1.51	–11.02	–8.83	–7.67	5.48
	0.7887	0.9731	1.51	–9.94	–8.58	–8.10	6.74
	0.8188	0.9447	1.52	–11.15	–9.14	–7.69	5.68
	0.8193	0.9445	1.52	–11.16	–9.16	–7.68	5.68
	0.8135	0.9517	1.53	–11.24	–8.90	–8.02	5.68
	0.8124	0.9459	1.50	−12.48	−9.60	−8.88	6.00
	0.8107	0.9478	1.51	−10.91	−8.91	−7.80	5.79
	0.0103	0.0125	0.01	0.44	0.20	0.18	0.43
C(11)–S(12)	0.8587	0.8832	1.37	–6.19	–7.34	–5.89	7.04
	0.8552	0.8822	1.38	–6.35	–7.38	–5.93	6.96
	0.8462	0.9017	1.31	–6.44	–7.17	–5.26	5.99
	0.8536	0.8934	1.29	–6.06	–6.92	–5.37	6.23
	0.8510	0.8927	1.29	–6.14	–6.94	–5.39	6.19
	0.8597	0.8868	1.39	–7.01	–7.57	–6.23	6.79
	0.8764	0.8869	1.36	−9.30	−7.84	-6.26	4.80
	0.8541	0.8900	1.34	−6.36	−7.22	-5.68	6.53
	0.0046	0.0067	0.04	0.32	0.24	0.36	0.41
C(11)–N(10)	0.5914	0.7341	2.50	–23.64	–21.27	–17.78	15.40
	0.5954	0.7342	2.49	–23.00	–21.08	–17.59	15.67
	0.6121	0.7146	2.40	–19.31	–19.64	–16.60	16.94
	0.6149	0.7121	2.47	–20.96	–20.70	–17.49	17.24
	0.6153	0.7120	2.47	–20.93	–20.72	–17.46	17.25
	0.5933	0.7350	2.40	–21.08	–19.68	–16.91	15.51
	0.4796	0.8107	2.56	−30.85	−21.56	−18.52	9.22
	0.6037	0.7237	2.45	−21.49	−20.52	−17.31	16.34
	0.0105	0.0108	0.04	1.44	0.64	0.41	0.82

First six lines correspond to experimental forms Ia, Ib, II, IIIa, IIIb, IV, respectively; the next line (italic) corresponds to reference density from wave function of optimized geometry, and the final two lines (bold) are the sample mean and sample standard deviation for the experimental data.

In units of Å.

In units of e Å–3.

In units of e Å–5.

Figure 5

Molecular graphs of (a) experimental form IV, (b) theoretical conformer 1–0, (c) theoretical conformer 1–180. Atomic color coding: H pale blue, C gray, N blue, O red, S orange; bond critical points are shown as small red spheres, and ring critical points are shown as small yellow spheres.

Perhaps the most interesting feature in the molecular graphs is the bond path associated with the S12···O11 intramolecular interaction. This feature was observed in all molecules and suggests it may be responsible for the observed molecular conformation. A representative plot of the Laplacian in the plane of this interaction is shown in Figure 6, while a complete set of plots is given in the Supporting Information (Figure S14). A local charge concentration at the O atom is linked, via the bond path, to an area of relative charge depletion on the S atom, suggesting that this interaction is primarily electrostatic in nature, that is, Sδ+···Oδ−(=S). This view is confirmed by the electrostatic potential shown in Figure 7 and S15 (Supporting Information) where a relatively positive zone on the S atom (in green) is close to the negatively charged O atom (in red). The source function[54] at the S···O bcp as a reference point is also consistent with this view, as almost all atomic basins have a noticeable influence on the electron density; that is, the source is highly delocalized (see Figure S16a, Supporting Information).

Figure 6

Plot of the Laplacian function in form IIIa through the S12–N1–O11 plane, showing the bond paths and bond (red circles) and ring (yellow circles) critical points associated with the S12···O11 intramolecular interaction.

Figure 7

Electrostatic potential (e Å–1) for form IV mapped onto the 0.5 e Å–3 electron density isosurface. The potential at +1.495 e Å–1 is shown in purple and −0.044 e Å–1 in red.

In a charge density study by Guru Row and co-workers[55] on the closely related molecule 2-(4-amino-benzosulfonimido)-5-methyl-3H-1,3,4-thiadiazole (trivial name sulfamethizole) and its salts, a closely similar intramolecular Sδ+···Oδ−(=S) interaction has been observed. Recently, Jackson et al.[56] have discussed the importance of S···O (and similar) interactions in “locking” the molecular configuration in a number of thiophene based systems. They conclude that in general this S···O interaction is not the dominating factor in determining the molecular conformation, and X–H···O interactions were shown, in certain cases, to be energetically more favorable. In view of these results, we have examined the conformation of 1 in more detail, specifically the effect of driving the S11–N10–C11–S12 torsion angle φ. Two low energy conformations were observed (Figure S8, Tables S14 and S15, Supporting Information), one where φ is close to zero (as found in the experimental structures, conformation 1–0), and one where φ is ∼180°, conformation 1–180 engendering a significant N–H···O interaction. The molecular graphs for these two conformations are shown in Figure 5b,c and demonstrate the presence of bond paths for the respective weak intramolecular interactions. The source function for conformation 1–180 (Figure S16b, Supporting Information) is consistent with the interpretation of the N–H···O interaction as an intermediate to weak H-bond.[57] The evolution of the total energy and the adjacent C111–S11–N10–C11 torsion angle χ as a function of φ is shown in Figure S17 (Supporting Information) and indicates a barrier of ∼50 kJ mol–1 to interconversion between conformers. This barrier presumably arises because of partial loss of π character in the iminium C11–N10 bond (mean experimental distance 1.327(1) Å, optimized theoretical distance 1.292 Å) during the rotation. The C11–N10 internuclear distance increases by less than 0.01 Å during the rotation, but the delocalization index[58] δ(ΩC11,ΩN10) decreases from 1.325 to 1.275, consistent with a small loss of π character.

Conformer 1–180 with the N–H···O interaction is marginally more stable by ∼6 kJ mol–1 in the gas phase but has never been observed in an experimental structure. It is interesting to speculate why this is the case, but the most likely reason is that the various intermolecular interactions involving this H atom in the crystal structures [N–H···N(imine) in form I, N–H···N(NH2) in forms II – IV, and N–H···O in form V] provide a greater stabilization than can be achieved from the involvement of this group in an intramolecular interaction. This is certainly borne out by the analysis of the H-bond contributions to the lattice energies discussed below. With regard to the idea that the intramolecular S···O (and N–H···O) interactions have a “locking” effect on the molecular conformation, some clear evidence for this is provided by the evolution of the C–S–N-C torsion χ and the associated S···O and H···O distances, shown in Figure S17 (Supporting Information). This evolution is not smooth, but significant jumps in the torsion χ are observed, which can be traced to the energetic favorability of these intramolecular interactions. As soon as the driving torsion φ approaches a value whereby the S···O (or N–H···O) interactions become sterically feasible, the torsion χ immediately adjusts to accommodate the interaction. As a result, the S···O (or H···O) distances remain quite short, ∼2.7 (1.9) Å, over a significant range of φ. The calculated barrier indicates that both conformers of the imine tautomer ought to be detectable and their interchange observable by NMR spectroscopy. However, the reported NMR spectral data[59] do not indicate any such exchange.

Finally, some mention should be made of the effect of minor disorder on derived topological parameters. As intimated in the Experimental Section, there are residual density features present for forms I and II, which could be attributed to some minor disorder. The accepted wisdom is that any disorder in a data set renders that data very suspect, at least in terms of a charge density analysis.[60] In our case, the level of disorder was so low that it could not be satisfactorily modeled by a secondary site and was effectively ignored. One way of quantifying the global effect that any such disorder has on a derived parameter P is through an R(par) value:where P(model) is the refined parameter value from the experimental data and P(theor) is that obtained from theory and the summation is over all experimental measurements of the parameter. Table 6 lists the R(par) values for the derived topological parameters Prho = ρbcp and Pdelrho = ∇2(ρbcp), where the wave function density is used as the reference theoretical density. As expected, the global agreements between experiment and theoretical parameters are much better for ρbcp than for ∇2(ρbcp), but there is no obvious worsening for forms I and II (with minor disorder), compared with forms III and IV (with no detectable disorder). Similar results are obtained for other derived topological parameters. In this study at least, we find no evidence that the proposed minor disorder has any detectable deleterious effect on the derived static density. The multipole populations obtained from the least-squares process seem relatively insensitive to errors in the structure factors from the minor disorder.

Table 6

The R(parameter) Values for R(rho) and R(delrho)

form	R(rho)	R(delrho)
Ia	0.0021	0.2029
Ib	0.0031	0.2012
II	0.0015	0.2096
IIIa	0.0085	0.1837
IIIb	0.0096	0.1757
IV	0.0056	0.2316

Interaction and Lattice Energies from Electron Densities

The difficulties in predicting whether a molecule will exhibit polymorphism (and if so, which form is the thermodynamically most stable one) are well-known and have been recently summarized by Price.[61] Since the energy differences between polymorphs may only be on the order of a few kJ mol–1, any prediction of rank-order stabilities is very challenging indeed.[61] The unit cell volumes we obtain for forms I–III at 100 K agree extremely well (<0.15% discrepancy) with the previously reported[8] 100 K cell volumes. Solely on the basis of the crystal densities at 100 K (150 K for form V) the expected relative order of thermodynamic stability is V < I < II < III < IV, but the presence of strong H-bonding in 1 is very likely to change this ordering. As stated in the Introduction, the stability ranking of the sulfathiazole polymorphs is temperature-dependent.

The lattice energy calculation (LATEN option in XD2006[30]) relies on the estimation of the total intermolecular interaction energies Eint, which is decomposed into several terms:The electrostatic term is estimated using the EP/MM method[62] and is considered to be essentially accurate. The induction term is included in the experimentally determined multipole populations (since these are determined in the crystal environment), while the exchange-repulsion and dispersion terms are approximated by atom–atom potentials of Williams and Cox.[63] Initial publications by Volkov and Coppens[64] were very encouraging regarding the comparability of this method with theoretical estimates of the interaction energy, and we were interested to use this approach to determine the relative lattice energies of the polymorphs of sulfathiazole from the charge density analysis. To reassess the accuracy of interaction energies determined from experimental multipole models, we have compared the interaction energies between the two independent molecules in the asymmetric unit of form III using the INTEREN option in XD2006 and computational methods, including a complete basis set limit calculation (see Experimental Section), which may be taken as providing a bench mark value. The results are listed in Table 7 and demonstrate the importance of including dispersion effects, since the dispersion-corrected functional B97-D gives a value much closer to the bench mark than the standard B3LYP functional. The multipole-derived value is reasonably close to the bench mark value but is clearly not as good as the dispersion-corrected DFT method.

Table 7

Gas Phase Interaction Energies between the Independent Molecules in Polymorph III

method	interaction energy (kJ mol–1)
multipole model (INTEREN/XD)	–50.7
DFT (B3LYP – def2-TZVPP)	–73.3
DFT (B97-D – def2-TZVPP)	–67.6
CBS(T)	–65.9

The lattice energies computed from the experimental multipoles are highly dependent on the exact refinement model used (Table S16, Supporting Information), in particular on the treatment of the H atom thermal motion. In Table 8, values obtained from the best multipole model are compared with several theoretical estimates, including Gavezotti’s PIXEL[44b] methodology. All these methods provide estimates of the lattice energy at 0 K, that is, only the enthalpic component. It is clear that they do not even agree on the relative order of stabilities, except that most methods indicate that form I is the least stable. Although error estimates for the derived lattice energies are not available in the XD2006 program suite, a study by Destro et al.[65] suggests a typical error in the region of 10 kJ mol–1. The discrepancies between experimental lattice energies from our multipole models and several theoretical approaches are in line with this estimate, but it is equally clear that no method is sufficiently accurate to allow a definitive stability order to be determined, at least for sulfathiazole.

Table 8

Experimental and Theoretical Estimates of Lattice Energies (kJ mol–1) for Sulfathiazole

form	experimental multipole	AA-CLP[44a]	PIXEL[44b] MP2/631G**	PIXEL[44b] B3LYP/631G**	PERIODIC-DFT[43] CRYSTAL09
I	–183	–183	–213	–210	–213
II	–205	–171	–230	–226	–225
III	–230	–182	–228	–224	–226
IV	–206	–191	–227	–224	–222

Other methods of using information on the electron density topology to obtain a relative ranking of lattice energies in polymorphs are even more approximate and hence (generally) less useful. One such method involves identifying the intermolecular interactions through their QTAIM signature—the bond path. These bond paths are generally related to well established bonding situations, such as strong and weak H-bonds, and the individual H-bond bonding energies may be estimated using the approximation of Espinosa et al.[66] Recently, Abramov[67] has employed this method to estimate the stability ranking of several polymorphic compounds, including sulfathiazole, for which the relative ranking was II ∼ III > I > V.[68] Given the approximate nature of this method, it is not surprising that the stability ranking differs from other estimations.

To investigate this methodology in more detail, the intermolecular interactions that contribute to the total lattice energy calculations provided by the PIXEL method were analyzed in more detail. The results are summarized in Tables S17 and S18 (Supporting Information). Only the intermolecular interactions of the first shell, that is, the closest neighbor H-bond interactions, were taken into consideration. The total intermolecular interaction energies of the first shell are greatest in form I, followed by III, IV, and II. However, the total interaction energy within the crystal lattice shows a different ranking stability of II > III > IV > I. In the first shell, the intermolecular interactions were limited only to the H-bonds; however, there are also repulsive interactions that contribute to the total lattice energy. In the case of form I, for example, there is a considerable repulsive interaction of 15.4 kJ mol–1 between two Ia···Ia molecules; this is also observed in forms III and IV, where values of 17 kJ mol–1 and 19.0 kJ mol–1, respectively, are found. In form II, the intermolecular repulsive interaction energies are much smaller, at 3.9 kJ mol–1. These are of comparable magnitudes to the H-bond energies, and we conclude that it is not possible to estimate the relative stability of the polymorphic forms of sulfathiazole only on the basis of the H-bond interactions.

Conclusions

The topological analysis of the electron density for the polymorphs of sulfathiazole shows an essentially identical set of properties for all forms, indicating that the potential fields due to the differing crystal packings have little effect on the examined molecular properties. An intramolecular Sδ+···Oδ−(=S) interaction between the thiazole sulfur atom and an oxygen in the sulfone group is consistently observed in all polymorphs. Two low energy conformers in the gas phase have been identified, one being very similar to that found in the crystalline phase, while the other exhibiting an N–H···O interaction of the ring NH proton with a sulfone oxygen atom. The calculations indicate that these intramolecular interactions have some “locking” influence on the conformers. The lattice energies obtained from the experimental multipole populations, although in reasonable agreement with those obtained theoretically, do not have the precision and accuracy required to rank the polymorph stabilities of sulfathiazole.