Crystal Structure Optimization and Gibbs Free Energy Comparison of Five Sulfathiazole Polymorphs by the Embedded Fragment QM Method at the DFT Level

: Molecular crystal plays an important role in many ﬁelds of science and technology, but it often crystallizes in di ﬀ erent polymorphs with di ﬀ erent physical properties. To guide the experimental synthesis of candidate materials, the atomic-scale model is frequently used to predict the most stable polymorph and its structural properties. Here, we show how an ab initio method can be used to achieve a rapid and accurate prediction of sulfathiazole crystal polymorphs (an antibiotic drug), based on the Gibbs free energy calculation and Raman spectra analysis. At the atmospheric pressure and the temperature of 300 K, we demonstrate that form III (FIII) is the most stable structure of sulfathiazole. The agreement between the predicted and experimental crystal structures corresponds to the order of stability for ﬁve sulfathiazole polymorphs as FI < FV < FIV < FII < FIII, which is achieved by employing the density functional theory (DFT) calculations.


Introduction
The exploration of drug polymorphism is one of the most active fields since the 20th century [1][2][3], which enables one to understand, rationalize and engineer the bulk properties of solids. Crystal structure prediction (CSP) is the calculation of the crystal structures of solids from first principles, which is based on the compositions of compounds and can determine the relative stability of crystal structures [4][5][6]. The CSP method is widely used to search the most stable crystal structures of pharmaceutical molecules [7][8][9][10]. For example, Shtukenberg et al., identified four new coumarin polymorphs based on the CSP method. They demonstrated that within the five crystal structures, coumarin was one of the few rigid molecules showing diverse polymorphisms at ambient conditions [11].
Pioneering CSP works have successfully assisted pharmaceutical industries and researchers to find thermodynamically favorable crystal packings with advanced efficacies [12,13]. However, most of them are based on calculations of the lattice energy, which usually neglects the effects of entropies and temperatures when evaluating the stabilities of crystal polymorphs. In this study, we use the Gibbs free energy instead of the lattice energy to calculate the polymorphisms of crystals. The Gibbs free energy provides complementary chemical insights at the atomistic level by calculating the relative energies of the molecular crystals and taking into account the entropy and temperature effects of the molecular crystals [14,15].
Sulfathiazole (STZ) (N -2-thiazolylsulphanilamide), an organosulfur compound discovered in 1939, is one of the most potent sulphonamides and is a typical example of this family of bacteriostatic drugs [16,17]. It occurs in more than one crystalline modification and five polymorphs have been reported (FI, FV, FIV, FII and FIII), but confusion arises over the stabilities of the polymorphs and their structures [18,19]. For example, Munroe et al. [20] investigated the relative stabilities of the five polymorphs of STZ based on the solution and solid-state methods, where FI and FV were identified as the least stable polymorphs below 50 • C using a combination of solubility measurements and isothermal suspension equilibration. Since sulfathiazole was reported in the literature as a classic example of a drug with five polymorphs and the differentiation among the forms was far from trivial, we chose it as a suitable system to perform the crystal structure prediction of pharmaceutical molecules. In this work, we proposed the first principles method to calculate the Gibbs free energy and compare the stabilities of five polymorphisms of STZ. To treat the macromolecular systems like molecular crystals of drug molecules, we adopted the embedded fragment quantum mechanical (QM) method [21][22][23] to calculate the total energy, where the total electronic energy per unit cell could be obtained by taking a proper combination of the energies of monomers and dimers that were embedded in the electrostatic field of the rest of the crystalline environment. To get stable structures at ambient pressure, we used the quasi-Newton algorithm [24] for the crystal structure optimization by minimizing the enthalpy of the unit cell. The crystal structure optimization was performed at the DFT level using ωB97XD/6-31G*, and the point charges of the embedding field were determined at the same level. The BFGS procedure [25] was used to update the Hessian matrix and 0.001 Hartree/Bohr was set as the convergence criterion for the maximum gradient. Furthermore, the Raman spectra of five polymorphs of STZ were calculated and compared with the experiment, which could be utilized as fingerprints to discriminate different polymorphs. A theoretical calculation shows that form III of STZ is the most stable phase under ambient temperature and pressure, and the stability order of five STZ polymorphs is FI < FV < FIV < FII < FIII, which matches the observed result very well. Figure 1 shows the molecular structure, with the formula of C 9 H 9 N 3 O 2 S 2 . The lattice parameters and the reference codes [26][27][28] of five polymorphs in CCDC are given in Table 1, where forms I, II and III have the same space group of P2 1 /c and forms IV and V have the same space group of P2 1 /n, respectively. polymorphs based on the CSP method. They demonstrated that within the five crystal structures, coumarin was one of the few rigid molecules showing diverse polymorphisms at ambient conditions [11]. Pioneering CSP works have successfully assisted pharmaceutical industries and researchers to find thermodynamically favorable crystal packings with advanced efficacies [12,13]. However, most of them are based on calculations of the lattice energy, which usually neglects the effects of entropies and temperatures when evaluating the stabilities of crystal polymorphs. In this study, we use the Gibbs free energy instead of the lattice energy to calculate the polymorphisms of crystals. The Gibbs free energy provides complementary chemical insights at the atomistic level by calculating the relative energies of the molecular crystals and taking into account the entropy and temperature effects of the molecular crystals [14,15]. Sulfathiazole (STZ) (N'-2-thiazolylsulphanilamide), an organosulfur compound discovered in 1939, is one of the most potent sulphonamides and is a typical example of this family of bacteriostatic drugs [16,17]. It occurs in more than one crystalline modification and five polymorphs have been reported (FI, FV, FIV, FII and FIII), but confusion arises over the stabilities of the polymorphs and their structures [18,19]. For example, Munroe et al. [20] investigated the relative stabilities of the five polymorphs of STZ based on the solution and solid-state methods, where FI and FV were identified as the least stable polymorphs below 50 °C using a combination of solubility measurements and isothermal suspension equilibration. Since sulfathiazole was reported in the literature as a classic example of a drug with five polymorphs and the differentiation among the forms was far from trivial, we chose it as a suitable system to perform the crystal structure prediction of pharmaceutical molecules. In this work, we proposed the first principles method to calculate the Gibbs free energy and compare the stabilities of five polymorphisms of STZ. To treat the macromolecular systems like molecular crystals of drug molecules, we adopted the embedded fragment quantum mechanical (QM) method [21][22][23] to calculate the total energy, where the total electronic energy per unit cell could be obtained by taking a proper combination of the energies of monomers and dimers that were embedded in the electrostatic field of the rest of the crystalline environment. To get stable structures at ambient pressure, we used the quasi-Newton algorithm [24] for the crystal structure optimization by minimizing the enthalpy of the unit cell. The crystal structure optimization was performed at the DFT level using ωB97XD/6-31G*, and the point charges of the embedding field were determined at the same level. The BFGS procedure [25] was used to update the Hessian matrix and 0.001 Hartree/Bohr was set as the convergence criterion for the maximum gradient. Furthermore, the Raman spectra of five polymorphs of STZ were calculated and compared with the experiment, which could be utilized as fingerprints to discriminate different polymorphs. A theoretical calculation shows that form III of STZ is the most stable phase under ambient temperature and pressure, and the stability order of five STZ polymorphs is FI < FV < FIV < FII < FIII, which matches the observed result very well. Figure 1 shows the molecular structure, with the formula of C9H9N3O2S2. The lattice parameters and the reference codes [26][27][28] of five polymorphs in CCDC are given in Table 1, where forms I, II and III have the same space group of P21/c and forms IV and V have the same space group of P21/n, respectively.   According to that the lower the free energy, the more stable the structure, the Gibbs free energy is calculated to estimate the stabilities of different forms of STZ. The Gibbs free energy of the unit cell G unit is calculated by

Methods
where H unit is the enthalpy, U v is the zero-point vibrational energy, T is the temperature and S is the entropy of the unit cell. When applying external pressure P, the enthalpy H unit of the unit cell is given by, where U int is the internal energy and V unit is the unit cell volume. Based on the embedded fragment QM approach [29][30][31][32], the internal energy U int of the unit cell can be obtained by, where n is the three-integer index of cells, E µ(n) is the energy of the µth molecule in the nth unit cell, and E µ(0)ν(n) is the energy of a dimer for the µth molecule in the 0th unit cell and the νth molecule in the nth unit cell. The first part of Equation (3) sums over all single molecular energies in the 0th unit cell. The second part calculates two-body QM interactions that have a shorter distance than a given cutoff distance d cut in the central unit cell. The third part calculates the QM interactions between two molecules that have a shorter distance than d cut in the 0th unit cell and the nth unit cell, respectively (d cut was set to 4 Å in this study) [33][34][35][36]. The first three terms are calculated using ωB97XD/6-31G* [37] in the electrostatic field which is represented by the electrostatic potential charges of the rest crystal fitted at the ωB97XD/6-31G* level. The interactions between two molecules of which the distance is larger than d cut are approximately treated as pairwise charge-charge Coulomb interactions. The last term E LR accounts for the long-range electrostatic interactions. Here, we choose a 3 × 3 × 3 supercell for QM calculation of the molecular crystal (namely, i = 1 in Equation (3)), an 11 × 11 × 11 supercell for the background charges and a 41 × 41 × 41 supercell for long-range electrostatic interactions, respectively.
With harmonic approximation, the zero-point vibrational energy U v and the entropy S v can be calculated by Equations (4) and (5): where k 0 is the Boltzmann constant, ω nk is the frequency in the nth phonon branch with the wave vector k, and β = 1/k 0 T. The product over k must be taken over all K that includes the evenly spaced grid points of k in the inverted lattice. In this study, the k-grid of 21 × 21 × 21 has been used (K = 9261).

Results and Discussion
Considering the large size and quantity of molecules in the supercell, we adopted a fragment-based QM method at the ωB97XD/6-31G* level for the calculation of Gibbs free energy and the vibrational spectrum. The experimental crystal data of five polymorphs of STZ were selected from the CCDC reference codes, as shown in Table 1. The five polymorphs of STZ are named as FI, FII, FIII, FIV, and FV, respectively. Figure 2 shows the crystal structures of five polymorphs of STZ (FII, FIV, FI, FIII, and FV) in a unit cell, where FII and FIV each have four molecules per unit cell and FI, FIII, FV each have eight molecules per unit cell. Figure 2 shows that different polymorphs have four molecular arrangements. The lattice constant comparison of five polymorphs of STZ is shown in Table 2, where the optimized lattice constants are in good agreement with the experimental data. The average lattice constant deviation between the calculation and experiment is about 0.1 to 0.3 Å.
where is the Boltzmann constant, is the frequency in the nth phonon branch with the wave vector k, and 1/ . The product over k must be taken over all K that includes the evenly spaced grid points of k in the inverted lattice. In this study, the k-grid of 21 × 21 × 21 has been used (K = 9261).

Results and Discussion
Considering the large size and quantity of molecules in the supercell, we adopted a fragmentbased QM method at the ωB97XD/6-31G* level for the calculation of Gibbs free energy and the vibrational spectrum. The experimental crystal data of five polymorphs of STZ were selected from the CCDC reference codes, as shown in Table 1. The five polymorphs of STZ are named as FI, FII, FIII, FIV, and FV, respectively. Figure 2 shows the crystal structures of five polymorphs of STZ (FII, FIV, FI, FIII, and FV) in a unit cell, where FII and FIV each have four molecules per unit cell and FI, FIII, FV each have eight molecules per unit cell. Figure 2 shows that different polymorphs have four molecular arrangements. The lattice constant comparison of five polymorphs of STZ is shown in Table 2, where the optimized lattice constants are in good agreement with the experimental data. The average lattice constant deviation between the calculation and experiment is about 0.1 to 0.3 Å.    [26][27][28] and calculated lattice parameters of five forms of STZ. The calculation is based on the density functional theory (DFT) using ωB97XD/6-31G*. The Exp. represents the experimental data and Cal. represents the calculated data.

Gibbs Free Energy Calculations
The Gibbs free energy differences of five polymorphs of STZ were calculated using ωB97XD/6-31G* with the embedded fragment method. Figure 3 shows the comparison of the calculated Gibbs free energy differences of five polymorphs of STZ under standard atmospheric pressure and 300 K. The values of Gibbs free energy differences per molecule of STZ are given in Table S1 of the Supplementary Materials. Figure 3 also shows the relative stability of five polymorphs of STZ measured from the experiment. At 300 K, the calculated stability order is FI < FV < FIV < FII < FIII, which matches the experimental result [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] very well. Here, the polymorph stability is correlated with free energy, with the most stable polymorph having the lowest free energy. Therefore, the most stable form of STZ is FIII at 300 K.

Gibbs Free Energy Calculations
The Gibbs free energy differences of five polymorphs of STZ were calculated using ωB97XD/6-31G* with the embedded fragment method. Figure 3 shows the comparison of the calculated Gibbs free energy differences of five polymorphs of STZ under standard atmospheric pressure and 300 K. The values of Gibbs free energy differences per molecule of STZ are given in Table S1 of the Supplementary Materials. Figure 3 also shows the relative stability of five polymorphs of STZ measured from the experiment. At 300 K, the calculated stability order is FI < FV < FIV < FII < FIII, which matches the experimental result [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] very well. Here, the polymorph stability is correlated with free energy, with the most stable polymorph having the lowest free energy. Therefore, the most stable form of STZ is FIII at 300 K.

Vibrational Spectroscopy
Different polymorphs of crystal structures have the same formula, different molecular arrangements, different lattice constants, and possess different vibrational spectra [39][40][41]. The vibrational spectrum is treated as a fingerprint to distinguish different polymorphs of molecular crystals. We calculated the full range of Raman and IR spectra of five forms of STZ, which are shown in Figures S1 and S2 of the Supplementary Materials. As shown in Figures S1 and S2, the Raman and IR spectra of the five STZ forms are mainly distributed in two ranges, from 0 to 2000 cm −1 and from Figure 3. The calculated (ellipse) Gibbs free energy differences of five polymorphs of STZ at 300 K and under standard atmospheric pressure. The calculated energy differences per molecule of STZ are for FI, FV, FIV, FII with reference to FIII, respectively. The experimental [20] relative stability of five polymorphs of STZ is shown in the rectangle.

Vibrational Spectroscopy
Different polymorphs of crystal structures have the same formula, different molecular arrangements, different lattice constants, and possess different vibrational spectra [39][40][41]. The vibrational spectrum is treated as a fingerprint to distinguish different polymorphs of molecular crystals. We calculated the full range of Raman and IR spectra of five forms of STZ, which are shown in Figures S1 and S2 of the Supplementary Materials. As shown in Figures S1 and S2, the Raman and IR spectra of the five STZ forms are mainly distributed in two ranges, from 0 to 2000 cm −1 and from 3000 cm −1 to 4000 cm −1 , respectively. However, the spectral differences for the five STZ polymorphs can be found in the low frequency and high-frequency ranges. Based on the experimental work proposed by Munroe et al. [20], where they measured the Raman spectra of five STZ polymorphs in the laboratory from 1100 to 1500 cm −1 , we picked the calculated corresponding Raman spectra and compare them with the experiment (see Figure 4). From 1100 cm −1 to 1500 cm −1 in Figure 4, FI and FII both have six distinct Raman peaks, while FIII, FIV and FV have five distinct Raman peaks. In Figure 4, we added the peak positions for the calculated Raman spectra, which display a slight difference in peak position as compared to the experiment. Such differences mainly arise from the non-harmonic effect in the proposed work, which was not taken into account in this study and makes the calculated value slightly larger than the experimental results. In general, the calculated Raman spectra are consistent with the experimental results of Munroe et al. [20], which establishes the correctness of our calculation. 3000 cm −1 to 4000 cm −1 , respectively. However, the spectral differences for the five STZ polymorphs can be found in the low frequency and high-frequency ranges. Based on the experimental work proposed by Munroe et al. [20], where they measured the Raman spectra of five STZ polymorphs in the laboratory from 1100 to 1500 cm −1 , we picked the calculated corresponding Raman spectra and compare them with the experiment (see Figure 4). From 1100 cm −1 to 1500 cm −1 in Figure 4, FI and FII both have six distinct Raman peaks, while FIII, FIV and FV have five distinct Raman peaks. In Figure  4, we added the peak positions for the calculated Raman spectra, which display a slight difference in peak position as compared to the experiment. Such differences mainly arise from the non-harmonic effect in the proposed work, which was not taken into account in this study and makes the calculated value slightly larger than the experimental results. In general, the calculated Raman spectra are consistent with the experimental results of Munroe et al. [20], which establishes the correctness of our calculation.

Conclusions
In this work, we studied the stabilities of five polymorphs of the STZ crystal using the embedded fragment QM method at the DFT level. The lattice parameters and Raman spectra of five polymorphs of STZ are reproduced accurately by calculations as compared to the experiment, along with the stability determination of FI < FV < FIV < FII < FIII. The calculated Gibbs free energy of the STZ crystal demonstrates that form FIII is the most stable structure of STZ under standard atmospheric pressure and 300 K. Furthermore, the characteristic peaks in Raman spectra offer an effective way to discriminate different forms of polymorphism. The proposed theoretical approaches for stability determination of different polymorphs of the molecular crystal as well as the accurate Gibbs free energy calculation provide an efficient platform for pharmaceutical structure design.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4352/9/5/256/s1, Figure S1: Calculated Raman spectra of five forms of STZ under standard atmospheric pressure. Figure S2: Calculated IR spectra of five forms of STZ under standard atmospheric pressure. Table S1: The Gibbs free energy differences of five polymorphs of STZ. The calculated energy differences per molecule of STZ is for FI, FV, FIV, FII, with reference to FIII, respectively.