Density Functional Theory and Density Functional Tight Binding Studies of Thiamine Hydrochloride Hydrates

Thiamine hydrochloride (THCL), also known as vitamin B1, is an active pharmaceutical ingredient (API), present on the list of essential medicines developed by the WHO, which proves its importance for public health. THCL is highly hygroscopic and can occur in the form of hydrates with varying degrees of hydration, depending on the air humidity. Although experimental characterization of the THCL hydrates has been described in the literature, the questions raised in previously published works suggest that additional research and in-depth analysis of THCL dehydration behavior are still needed. Therefore, the main aim of this study was to characterize, by means of quantum chemical calculations, the behavior of thiamine hydrates and explain the previously obtained results, including changes in the NMR spectra, at the molecular level. To achieve this goal, a series of DFT (CASTEP) and DFTB (DFTB+) calculations under periodic boundary conditions have been performed, including molecular dynamics simulations and GIPAW NMR calculations. The obtained results explain the differences in the relative stability of the studied forms and changes in the spectra observed for the samples of various degrees of hydration. This work highlights the application of periodic DFT calculations in the analysis of various solid forms of APIs.


Introduction
Most of the chemical compounds exist in a variety of crystal forms depending on the temperature, pressure, humidity, and solvents used during the crystallization process, which is the definition of the phenomenon known as polymorphism [1].Polymorphic forms of active pharmaceutical ingredients (APIs) may differ in certain crucial characteristics, including solubility in water, dissolution rate, melting point, stability, tabletability, and others, which may ultimately affect the drug's stability and bioavailability [2][3][4].
Hydrates, a subtype of solid solvates, are a distinct set of structures, although they share certain similarities with polymorphs.Once water molecules are incorporated into a compound's crystal lattice, they modify the intricate H-bonding network [5].Therefore, hydrates display a distinct structure when compared to corresponding anhydrates [6].As a result, hydrates could differ from their anhydrous counterparts in terms of their physical and chemical characteristics [7].Hydrates were once referred to as pseudo-polymorphs because of their similarities to polymorphs [8].This term, however, is no longer considered to be the proper one when referring to hydrates [9].
Among the solid API solvates, hydrates are particularly interesting for several reasons.First of all, the water molecule stands out because of its special properties, which include its tiny size and ability to form interactions as both a donor and acceptor of H-bonding, occasionally concurrently; this makes it a crucial "building material" in the field of crystal engineering [10].In addition, unlike the majority of other organic solvents, water is a safe and non-toxic substance.Finally, due to the air's moisture content, spontaneous hydration may happen at any point during the manufacture or storage of an API, resulting in the development or phase transition of hydrates [11].
Non-stoichiometric hydrates, sometimes referred to as variable hydrates, are a unique type of hydrates in which the water content is constantly changeable as a function of the atmospheric water vapor pressure [12].The water molecules in the non-stoichiometric hydrate lattice may occupy any or all of the predetermined places [13].Such hydrate lattices frequently have channels or connected routes that allow water to easily enter or exit the lattice, causing the water content to continuously fluctuate as a function of the water vapor pressure.
The stability of both the API and the drug product can be affected by changes in hydration status.A formed drug may experience a chemical breakdown of one or more formulation components as a result of released water.Three outcomes are possible when water is removed from a drug's crystal lattice: (1) the lattice structure is unaffected by the water removal, (2) the lattice collapses, or (3) the lattice packing is changed upon dehydration.Due to vacancies within the lattice, the product phase is frequently more reactive and less stable when the lattice structure is preserved following dehydration (outcome 1 above).Similar to outcome 1, the noncrystalline product phase is metastable when dehydration causes a collapsed lattice.Properties like the rate of dissolution can be impacted by changes in the physical form (outcome 3).Therefore, the process of dehydration and the subsequent release of water into the formulation might have a significant impact on overall drug stability [14,15].
As previously stated, the stability, solubility, bioavailability, and formulatability of solid-state APIs are influenced by their physical and chemical characteristics, which are the result of how molecules are arranged in the solid state.Therefore, it is intriguing from a strictly scientific perspective, as well as having significant practical significance in the pharmaceutical industry, to be able to precisely foresee and explain such qualities utilizing molecular modeling tools [16].
When it concerns solid APIs, molecular modeling techniques are typically used to forecast their physicochemical and structural characteristics, explain experimental findings, or predict the conditions necessary to produce novel forms of solid pharmaceutics in order to reduce the number of experiments or improve the experimental setup.Additionally, calculated properties like Nuclear Magnetic Resonance (NMR) shielding constants can make solid-state analysis much easier [17].
However, molecular modeling techniques that simulate a single molecule in vacuum or in solution were proven to be inadequate and erroneous as the distinctive properties of each solid form result from short-and long-distance intermolecular interactions [18].While those "single molecules" approaches are often utilized effectively in other areas of pharmaceutical sciences, such as to examine drug-biomolecule interactions or predict the formation of complexes, alternative types of calculations should be used to study solidstate pharmaceutics.Those kinds of calculations, DFT and DFTB, under periodic boundary conditions, were performed in the present study.
Thiamine (Figure 1), also known as vitamin B 1 and thiamin, is one of the most important vitamins and must be supplied with food [19].In a properly balanced diet in healthy people, the effects of its deficiency are rarely observed because of its rich content in wholegrain products and an artificial enrichment of some food products with B-group vitamins.On the other hand, vitamin B 1 is sensitive to UV radiation and high temperatures, especially in alkaline pH [20].In addition, some drugs, e.g., chemotherapeutics, may inactivate thiamine [21], while others may reduce its absorption by competing for the transporter, such as metformin, used as the first choice drug in the treatment of diabetes [22].All these factors may contribute to the occurrence of deficiencies of this vitamin, even in developed countries.Thiamine in a dose of 50 mg of thiamine hydrochloride (THCL) in tablets is on the List of Essential Medicines developed by the World Health Organization (WHO), which proves its importance for public health [23].The object of this study, THCL, is hygroscopic and can occur in the form of hydrates with varying degrees of hydration, depending on the humidity.The anhydrous hydrochloride absorbs water and transforms into a non-stoichiometric hydrate (NSH) with a heterogeneous water content in the crystal, up to a compound with an equimolar ratio.In addition, NSH at room temperature and air humidity >53% can transform into a thermodynamically more stable form, i.e., hemihydrate (HH) [24].Such transformation can start already at the stage of tablet production by wet granulation, which is used in the production of drugs from poorly tableting substances, such as NSH.The use of water in the granulation process causes the NSH to dissolve and then recrystallize into the more stable HH.As a result of an incomplete transformation, thiamine hydrochloride can be present in both forms: NSH and HH.During storage, the remaining part of NSH may be transformed into HH, while the addition of some excipients may slow down the transformation process [25].The conversion of NSH into HH changes the properties of the tablets, e.g., an increase in hardness and an increase in the disintegration time of the tablet, which results in lower bioavailability of the drug [26].
Although experimental characterization of thiamine hydrates has been described in the literature [14,24,25,[27][28][29], the questions raised in previously published works suggest that additional research and in-depth analysis of NSH and HH dehydration behavior are still needed.For example, Chakravarty et al. [28] speculated whether the respective physical stabilities of NSH and HH can be explained by the water-binding environment in the NSH and HH lattices.In their other work, they explicitly stated that "Molecular modeling studies will be necessary to further understand the observed dehydration-induced SSNMR spectral changes."[27].
Inspired by those previously raised questions, we have decided to address the topic of thiamine dehydration using a variety of molecular modeling methods.Therefore, the main aim of this study was to characterize, by means of quantum chemical calculations, the behavior of thiamine hydrates and explain the previously obtained NMR spectroscopic results at the molecular level.

Crystal Structure Analysis
As stated in the introduction, THCL exists in a variety of forms, differing in the degree of hydration.
One of them, NSH, is characterized by variable water content in the crystal lattice up to an equimolar ratio.In extreme cases, this can lead to the formation of a monohydrate.The crystal structure of NSH in the form of monohydrate was deposited in CCDC under refcodes THIAMC12 and THIAMC14.The authors of THIAMC12 also deposited the structure of NSH after complete dehydration (refcode UNEXOA), which resulted in the formation of an anhydrous version of NSH.
Another solvate form of THCL is a stoichiometric hemihydrate (HH), whose structure has been deposited under refcode WUWJAA.This hemihydrate is a thermodynamically stable form.Its forced dehydration results in the collapse of the crystal lattice and The object of this study, THCL, is hygroscopic and can occur in the form of hydrates with varying degrees of hydration, depending on the humidity.The anhydrous hydrochloride absorbs water and transforms into a non-stoichiometric hydrate (NSH) with a heterogeneous water content in the crystal, up to a compound with an equimolar ratio.In addition, NSH at room temperature and air humidity >53% can transform into a thermodynamically more stable form, i.e., hemihydrate (HH) [24].Such transformation can start already at the stage of tablet production by wet granulation, which is used in the production of drugs from poorly tableting substances, such as NSH.The use of water in the granulation process causes the NSH to dissolve and then recrystallize into the more stable HH.As a result of an incomplete transformation, thiamine hydrochloride can be present in both forms: NSH and HH.During storage, the remaining part of NSH may be transformed into HH, while the addition of some excipients may slow down the transformation process [25].The conversion of NSH into HH changes the properties of the tablets, e.g., an increase in hardness and an increase in the disintegration time of the tablet, which results in lower bioavailability of the drug [26].
Although experimental characterization of thiamine hydrates has been described in the literature [14,24,25,[27][28][29], the questions raised in previously published works suggest that additional research and in-depth analysis of NSH and HH dehydration behavior are still needed.For example, Chakravarty et al. [28] speculated whether the respective physical stabilities of NSH and HH can be explained by the water-binding environment in the NSH and HH lattices.In their other work, they explicitly stated that "Molecular modeling studies will be necessary to further understand the observed dehydration-induced SSNMR spectral changes."[27].
Inspired by those previously raised questions, we have decided to address the topic of thiamine dehydration using a variety of molecular modeling methods.Therefore, the main aim of this study was to characterize, by means of quantum chemical calculations, the behavior of thiamine hydrates and explain the previously obtained NMR spectroscopic results at the molecular level.

Crystal Structure Analysis
As stated in the introduction, THCL exists in a variety of forms, differing in the degree of hydration.
One of them, NSH, is characterized by variable water content in the crystal lattice up to an equimolar ratio.In extreme cases, this can lead to the formation of a monohydrate.The crystal structure of NSH in the form of monohydrate was deposited in CCDC under refcodes THIAMC12 and THIAMC14.The authors of THIAMC12 also deposited the structure of NSH after complete dehydration (refcode UNEXOA), which resulted in the formation of an anhydrous version of NSH.
Another solvate form of THCL is a stoichiometric hemihydrate (HH), whose structure has been deposited under refcode WUWJAA.This hemihydrate is a thermodynamically stable form.Its forced dehydration results in the collapse of the crystal lattice and sample amorphization; this is the reason why the crystal structure of the dehydrated form of HH has never been recorded before.
Chosen crystallographic information on the structures that we will be referring to in this study is presented in Table 1.
Table 1.Crystallographic information concerning the chosen structures of thiamine hydrochloride.2, there are four crystallographically equivalent water molecules in the unit cell of the monohydrate form of NSH (THIAMC12), labeled as UL (upper left), BL (bottom left), UR (upper right), and BR (bottom right).As our aim was to simulate the gradual dehydration of this form, we decided to remove the water molecules one by one.Due to the symmetry of the system, the choice of one of the four water molecules to remove had no influence on the results of the calculations.Therefore, the UL molecule was removed, resulting in the formation of the 0.75 hydrate, named MH3W, due to the three water molecules left in the unit cell.However, the number of possible different versions of unit cells increased significantly when creating the hemihydrate form of NSH.While the number of permutations describing the possible options of removing two out of four water molecules present in the unit cell equals six, among those six structures, three crystallographically equivalent pairs existed.Finally, the structure of the 0.25 hydrate was obtained by leaving one water molecule in the unit cell.Again, due to the symmetry of the system, the choice of the molecule to be left had no influence on the results of the calculations.To clarify the process of structures' generation, Table 2 was created.As presented in Figure 3, there are four crystallographically equivalent water molecules in the unit cell of hemihydrate (HH) form of THCL (WUWJAA), labeled as UL (upper left), BL (bottom left), UR (upper right), and BR (bottom right).As our aim was to simulate the gradual dehydration of this form as well, we have decided to remove the water molecules by one, similarly to NSH.Again, due to the symmetry of the system, the choice of one of the present four water molecules to remove had no influence on the results.Therefore, the UL molecule was removed, resulting in the formation of the 0.375 hydrate, named HH3W, due to the three water molecules left in the unit cell.However, the number of possible versions of unit cells increased significantly when creating the 0.25 hydrate from HH. Similarly to the NSH, the number of permutations describing the possible options of removing two out of four water molecules present in the unit cell equals six; among those six structures, three crystallographically equivalent pairs existed.Finally, the structure of the 0.125 hydrate was obtained by leaving one water molecule in the unit cell.Again, due to the symmetry of the system, the choice of the molecule to be left had no influence on the results of the calculations.To clarify the process of structures' generation, Table 3 was created.The results of the geometry optimization of structures described in Tables 2 and 3 are presented in Tables 4 and 5.For more facile comparison, chosen experimental structure

Crystal Structure Optimization
The results of the geometry optimization of structures described in Tables 2 and 3 are presented in Tables 4 and 5.For more facile comparison, chosen experimental structures (EXP) were also included in Tables 4 and 5.For better clarity, the results of optimization are also presented in Figures 4 and 5.
Table 4. Optimized unit cell dimensions of the MH structures of various hydration ratios, compared with the experimental ones.

Code Energy [kcal/mol]
Relative Energy  A good agreement has been observed between the corresponding experimental and computational results for both NSH and HH.What is interesting is that the difference between the theoretical value of the volume of MH4W and EXP I was found to be significantly lower than between EXP I and EXP II (3 Å 3 versus 33 Å 3 ); this is probably due to the temperature at which the SCXRD measurements were performed, namely 296 K and 173 K for EXP I and EXP II, respectively (Table 1).The structure recorded at a lower temperature is closer to the DFT-optimized one, as during the geometry optimization, the temperature was not included in the calculations.Also, as described in the introduction, the unit cell dimensions of NSH changed only slightly upon dehydration, with the exception of b length.
Even smaller changes in the unit cell dimensions have been observed for HH.In this case, the predictability of calculations was harder to assess due to the presence of solely one experimental structure of this form.According to the experimental results, any attempt to dehydrate the HH results in the collapse of the crystal lattice and sample amorphization.However, the results of geometry optimization from this work suggest that such structures could potentially exist; this can be explained by the fact that during geometry optimization, the influence of temperature is neglected, and the dynamic stability cannot be confirmed [30][31][32].Therefore, we have conducted the ab initio molecular dynamics simulations (aiMD), which are described in detail in Section 2.6.However, before the analysis of aiMD results, some energetic aspects should be discussed, which is carried out in the next Section 2.4.Table 5. Optimized unit cell dimensions of the HH structures of various hydration ratios, compared with the experimental ones.For MH2W, the arithmetic mean of the results (MH2W BLUR, MH2W ULBL, and MH2W ULUR) was shown.

Energetic Considerations
A comparison of the structures with the same THCL: water ratio is presented in Table 6.
The results of the calculations are in very good agreement with the experimental observations.They show that the HH4W is more stable than any of the modeled hemihydrates based on the structure of NSH, MH2W BLUR, MH2W ULBL, and MH2W ULUR by a few kcal/mol, which was reported previously; this also explains the conversion of MH2W to HH4W upon dehydration of MH4W.The positive (+1.688 kcal/mol) value of the ∆(HH-MH) for the anhydrous "0W" structures is in agreement with the experimental findings, revealing that it is possible to obtain the crystalline MH0W while the complete dehydration of HH0W results in the complete crystal structure destruction and sample amorphization.A good agreement has been observed between the corresponding experimental and computational results for both NSH and HH.What is interesting is that the difference between the theoretical value of the volume of MH4W and EXP I was found to be significantly lower than between EXP I and EXP II (3 Å 3 versus 33 Å 3 ); this is probably due to the temperature at which the SCXRD measurements were performed, namely 296 K and 173 K for EXP I and EXP II, respectively (Table 1).The structure recorded at a lower temperature is closer to the DFT-optimized one, as during the geometry optimization, the temperature was not included in the calculations.Also, as described in the introduction, the unit cell dimensions of NSH changed only slightly upon dehydration, with the exception of b length.
Even smaller changes in the unit cell dimensions have been observed for HH.In this case, the predictability of calculations was harder to assess due to the presence of solely one experimental structure of this form.According to the experimental results, any attempt to dehydrate the HH results in the collapse of the crystal lattice and sample amorphization.However, the results of geometry optimization from this work suggest that such structures could potentially exist; this can be explained by the fact that during  and c).For HH2W, the arithmetic mean of the results (HH2W BLUR, HH2W ULBL, and HH2W ULUR) was shown.Table 6.Energy values of the optimized structures.Due to the differences in Z, eight for HH, and four for MH, values obtained for HH have been divided by two to enable the comparison.Explanations of the abbreviations can be found in Tables 2 and 3. To compute the energy change upon dehydration, calculation of the energy of the single water molecule was necessary as the H 2 O is a product of such a reaction.While it is technically impossible to perform the calculations for nonperiodic systems in CASTEP, there is an alternative way to obtain such results, often called "molecule in the box" calculations.It consists of the creation of the unit cell, usually cubic, with a single molecule of interest inside it and unit cell lengths long enough to suppress any intermolecular interactions.Such a system is then optimized but with constrained unit cell dimensions.In this study, the cubic unit cell with equal lengths of 20 Å and a single water molecule inside it has been used to calculate the energy of H 2 O; this allowed us to determine the dehydration energies at each step, which are presented in Table 7. Analysis of the values presented in Table 7 shows that, as expected, each dehydration step is endothermic for both HH and MH.Also, dehydration of HH requires more energy at each stage of this process; this is in agreement with the experimental results showing that dehydration of MH is significantly less demanding than HH.The lowest values observed in the last step for both HH and MH are probably due to the creation of a higher symmetry structure upon complete dehydration and reduction of Z' to 1.

NMR Calculations
As stated in the introduction, NSH and HH have been extensively studied previously using solid-state NMR (ssNMR).However, for both those solvates the NMR experimentalists have found some observations that were not fully explained.Therefore, to understand, at the molecular level, the changes in the 13 C ssNMR spectra of NSH and HH occurring upon dehydration, the GIPAW calculations of NMR chemical shielding constants have been conducted, and the results are presented below.Due to the differences in the atom numbering between the previously published works, in this work, the atom numbering presented in Figure 6 has been used.

NMR Calculations of NSH
The results of the NMR calculations for the NSH of various degrees of hydration are presented in Table 8.As can be seen, the presence of water molecules has a major impact on the chemical shift (δ) value for most of the carbon atoms.For example, in C3, when the structure presents the monohydrate (MH4W), the values for all four atoms in the unit cell are close to 140.76 ppm (blue color), and the dehydration results in the increase of this value to 143.28 ppm (red color).Similarly, in MH3W, three "blue" (hydrated) values and one "red" (dehydrated) value were obtained, while for MH1W, one "blue" and three "red" values are present.Similar behavior has also been observed for C8, C1, C2, and C4.
Interestingly, major differences have been observed among the three structures presenting hemihydrate MH2W, which is well demonstrated by the SD values.For example, in C8, in the cases of MH2W BLUR and ULUR, the two distinct sets of values are observed, close to either 162.5 ppm (red) or 165 ppm (blue), which results in the high values of SD, over 1.2 ppm.On the contrary, in MH2W ULBL, chemical shifts for all four carbon atoms are averaged, and the SD equals only 0.27 ppm.Similar observations have been made for C9, C2 and C3.
Analysis of the changes in the experimental 13 C CP MAS NMR spectra resulting from dehydration of NSH (Figure 7) revealed that this process has the greatest impact on the chemical shift values of C3, C6, C5, and C4.
Molecules 2023, 28, x FOR PEER REVIEW 13 of 23 Analysis of the changes in the experimental 13 C CP MAS NMR spectra resulting from dehydration of NSH (Figure 7) revealed that this process has the greatest impact on the chemical shift values of C3, C6, C5, and C4.For C3, the increase of δ has been observed both experimentally (1.9 ppm) as well as in the calculation's results (2.5 ppm).Since in the spectrum of the 0.41 hydrate, the broadened peak of the averaged value of δ is observed instead of the two separate signals; this could indicate that MH2W ULBL is the dominant structure of the hemihydrate form, as in this one, the lowest value of SD for the calculated chemical shifts is observed (0.53 ppm vs. 1.96 ppm and 1.75 ppm for ULUR and BLUR, respectively).It should also be noted that among those three structures of MH2W, ULBL was the one with the lowest energy and, therefore, theoretically, the most stable one (Table 6).Similar observations can also be made for C4, C5, and C6 with even better accuracy between the experimental and theoretical increase of the δ upon dehydration.For example, for C4, the experimental one was 1.00 ppm, and the calculated one was 0.75 ppm.
For all of the carbon atoms, with the exception of C9, the sign of the calculated differences between the δ of the monohydrate and anhydrous forms (MH0W-MH4W) was in agreement with the experimental observations (EXPII MH0W-EXPII MH4W).
The NMR calculations also helped to properly assign the chemical shifts of C8 (C 2′) and C9 (C 4′).In the experimental spectrum (Figure 8), those two atoms were assigned to the broad peak at 163.20 ppm.According to the calculations results, the δ of C8 (163.95 ppm) is slightly higher than that of C9 (160.45 ppm).The experimentally observed averaging of this value, resulting in the overlapping of those two peaks, is probably a result of molecular dynamics.For C3, the increase of δ has been observed both experimentally (1.9 ppm) as well as in the calculation's results (2.5 ppm).Since in the spectrum of the 0.41 hydrate, the broadened peak of the averaged value of δ is observed instead of the two separate signals; this could indicate that MH2W ULBL is the dominant structure of the hemihydrate form, as in this one, the lowest value of SD for the calculated chemical shifts is observed (0.53 ppm vs. 1.96 ppm and 1.75 ppm for ULUR and BLUR, respectively).It should also be noted that among those three structures of MH2W, ULBL was the one with the lowest energy and, therefore, theoretically, the most stable one (Table 6).Similar observations can also be made for C4, C5, and C6 with even better accuracy between the experimental and theoretical increase of the δ upon dehydration.For example, for C4, the experimental one was 1.00 ppm, and the calculated one was 0.75 ppm.
For all of the carbon atoms, with the exception of C9, the sign of the calculated differences between the δ of the monohydrate and anhydrous forms (MH0W-MH4W) was in agreement with the experimental observations (EXPII MH0W-EXPII MH4W).
The NMR calculations also helped to properly assign the chemical shifts of C8 (C 2 ) and C9 (C 4 ).In the experimental spectrum (Figure 8), those two atoms were assigned to the broad peak at 163.20 ppm.According to the calculations results, the δ of C8 (163.95 ppm) is slightly higher than that of C9 (160.45 ppm).The experimentally observed averaging of this value, resulting in the overlapping of those two peaks, is probably a result of molecular dynamics.6. Explanations of the abbreviations (structure codes) can be found in Table 2. Due to Z = 4, for each carbon atom, four theoretical values have been obtained from calculations.To facilitate the analysis of the data in the table, the two-color scale was applied.The cell that holds the minimum is colored red, and the cell that holds the maximum value is colored blue.All other cells are colored proportionally.SD-standard deviation of the four values.
. 13 C CP MAS NMR spectra of 0.9 hydrate of NSH.Reprinted from [28] with permission from Elsevier.

NMR Calculations of HH
The results of the NMR calculations for the HH of various degrees of hydration are presented in Table 9.It should be noted that in the only work in which the 13 C ssNMR spectra of HH are being presented [27], the authors have not assigned any of the peaks to the particular carbon atom.In addition, the chemical shift values have not been provided either; therefore, the data presented in Table 9 in rows entitled "EXP HH4W" have been read by us directly from the spectra; this is why they have been rounded to either whole or half values, ±0.5 ppm.
As can be seen, similarly to the NSH, the presence of water molecules has a major impact on the chemical shift (δ) value for some of the carbon atoms.The authors of [27] have created a "partially dehydrated" HH2 form of HH, "HH was dried at 60 °C in a bench-top freeze dryer (Unitop 400L, Virtis, Gardiner, NY) under reduced pressure (20-60 mTorr) for 7 days.This product phase will be referred to as HH2."They then recorded the 13 C ssNMR spectra of both HH and HH2 and found some differences, which are presented in Figure 9.   13 C CP MAS NMR spectra of 0.9 hydrate of NSH.Reprinted from [28] with permission from Elsevier.

NMR Calculations of HH
The results of the NMR calculations for the HH of various degrees of hydration are presented in Table 9.It should be noted that in the only work in which the 13 C ssNMR spectra of HH are being presented [27], the authors have not assigned any of the peaks to the particular carbon atom.In addition, the chemical shift values have not been provided either; therefore, the data presented in Table 9 in rows entitled "EXP HH4W" have been read by us directly from the spectra; this is why they have been rounded to either whole or half values, ±0.5 ppm.
As can be seen, similarly to the NSH, the presence of water molecules has a major impact on the chemical shift (δ) value for some of the carbon atoms.The authors of [27] have created a "partially dehydrated" HH 2 form of HH, "HH was dried at 60 • C in a bench-top freeze dryer (Unitop 400L, Virtis, Gardiner, NY) under reduced pressure (20-60 mTorr) for 7 days.This product phase will be referred to as HH 2 ".They then recorded the 13 C ssNMR spectra of both HH and HH 2 and found some differences, which are presented in Figure 9.
Molecules 2023, 28, x FOR PEER REVIEW 14 of 23 Figure 8. 13 C CP MAS NMR spectra of 0.9 hydrate of NSH.Reprinted from [28] with permission from Elsevier.

NMR Calculations of HH
The results of the NMR calculations for the HH of various degrees of hydration are presented in Table 9.It should be noted that in the only work in which the 13 C ssNMR spectra of HH are being presented [27], the authors have not assigned any of the peaks to the particular carbon atom.In addition, the chemical shift values have not been provided either; therefore, the data presented in Table 9 in rows entitled "EXP HH4W" have been read by us directly from the spectra; this is why they have been rounded to either whole or half values, ±0.5 ppm.
As can be seen, similarly to the NSH, the presence of water molecules has a major impact on the chemical shift (δ) value for some of the carbon atoms.The authors of [27] have created a "partially dehydrated" HH2 form of HH, "HH was dried at 60 °C in a bench-top freeze dryer (Unitop 400L, Virtis, Gardiner, NY) under reduced pressure (20-60 mTorr) for 7 days.This product phase will be referred to as HH2."They then recorded the 13 C ssNMR spectra of both HH and HH2 and found some differences, which are presented in Figure 9.As noticed by the authors of [27], "(. ..) partial dehydration of HH under low pressure for 7 days (HH 2 ) resulted in a new peak at 157 ppm and the appearance of shoulders on the peaks at 148, 138, 65 ppm (. ..)".Those researchers were speculating on the origin of those extra peaks and concluded that "Molecular modeling studies will be necessary to further understand the observed dehydration-induced SSNMR spectral changes".Table 9. Calculated chemical shift values for carbon atoms, compared with the experimental ones.Atom numbering is presented in Figure 6.Explanations of the abbreviations (structure codes) can be found in Table 2. Due to Z = 8, for each carbon atom, eight theoretical values have been obtained from calculations.To facilitate the analysis of the data in the table, the two-color scale was applied.The cell that holds the minimum is colored red, and the cell that holds the maximum value is colored blue.All other cells are colored proportionally.SD-standard deviation of the four values.The results of the calculations, presented in Table 9, are in excellent agreement with the experimental observations.The largest difference of the calculated δ for the same atom, observed for HH4W and HH0W, was found for C1, and only for this atom did the difference in δ calculated for the "hydrated" and "dehydrated" structure exceed 1 ppm.Further, this was the only peak with direct splitting and evident separation observed in the experimental spectra (Figure 9).Also, the calculated δ was larger for C1 in the dehydrated structure (HH0W) than in the hydrated one (HH4W), which is also in agreement with the experimental observations.
To explain this observation made for the peak of C1, the authors of [27] have postulated two hypotheses.The first one stated that the loss of water causes changes in 13 C-14 N quadrupolar coupling.This hypothesis has been ruled out by the same authors based on the comparison of the spectra recorded at 300 MHz and 400 MHz spectrometers, as the differences in the width of signals were not significant.The second hypothesis was that the removal of water changes the physical location of the chloride ions within the crystal lattice.This hypothesis could not have been tested so far because of the lack of the experimental structure of the dehydrated form of HH.
It should be noted that in the structure of HH, two crystallographically inequivalent chloride ions can be found (Figure 10).The first one, Cl1, forms three intermolecular interactions with thiamine, two of them with H atoms of amine groups and one with the hydroxyl group.The second one, Cl2, forms one interaction with the H atom of water and one with the H atom of the thiazole ring of thiamine.
As noticed by the authors of [27], "(…) partial dehydration of HH under low pressure for 7 days (HH2) resulted in a new peak at 157 ppm and the appearance of shoulders on the peaks at 148, 138, 65 ppm (…)."Those researchers were speculating on the origin of those extra peaks and concluded that "Molecular modeling studies will be necessary to further understand the observed dehydration-induced SSNMR spectral changes".
The results of the calculations, presented in Table 9, are in excellent agreement with the experimental observations.The largest difference of the calculated δ for the same atom, observed for HH4W and HH0W, was found for C1, and only for this atom did the difference in δ calculated for the "hydrated" and "dehydrated" structure exceed 1 ppm.Further, this was the only peak with direct splitting and evident separation observed in the experimental spectra (Figure 9).Also, the calculated δ was larger for C1 in the dehydrated structure (HH0W) than in the hydrated one (HH4W), which is also in agreement with the experimental observations.
To explain this observation made for the peak of C1, the authors of [27] have postulated two hypotheses.The first one stated that the loss of water causes changes in 13 C-14 N quadrupolar coupling.This hypothesis has been ruled out by the same authors based on the comparison of the spectra recorded at 300 MHz and 400 MHz spectrometers, as the differences in the width of signals were not significant.The second hypothesis was that the removal of water changes the physical location of the chloride ions within the crystal lattice.This hypothesis could not have been tested so far because of the lack of the experimental structure of the dehydrated form of HH.
It should be noted that in the structure of HH, two crystallographically inequivalent chloride ions can be found (Figure 10).The first one, Cl1, forms three intermolecular interactions with thiamine, two of them with H atoms of amine groups and one with the hydroxyl group.The second one, Cl2, forms one interaction with the H atom of water and one with the H atom of the thiazole ring of thiamine.During the geometry optimization of the subsequently dehydrated structures of HH (Figure 11), we noticed that the RMSD of Cl1 was much higher than this of Cl2, which was caused by the fact that removing the water molecule that forms the intermolecular interaction with Cl2 causes this ion to move, and also results in the increase of the C1-H bond During the geometry optimization of the subsequently dehydrated structures of HH (Figure 11), we noticed that the RMSD of Cl1 was much higher than this of Cl2, which was caused by the fact that removing the water molecule that forms the intermolecular interaction with Cl2 causes this ion to move, and also results in the increase of the C1-H bond length; this consequently resulted in the increase of chemical shift of C1 and the presence of two peaks in the 13 C ssNMR spectrum of HH 2 (Figure 9).length; this consequently resulted in the increase of chemical shift of C1 and the presence of two peaks in the 13 C ssNMR spectrum of HH2 (Figure 9).In the case of the peak at 148 ppm, the authors of [27] have noticed a slight downfield shift, and this was also observed in the calculations results for C11 as the increase of the δ by 0.4 ppm.The change of this magnitude did not result in the peak separation, only with the increase of its width.
For the peak at 138 ppm, the experimentalists have observed the formation of another peak at the right shoulder of the one present in HH, resulting from partial dehydration; this was also present in the calculation results for C3-the upfield change by 0.5 ppm was observed as the result of the dehydration.
No major changes in the values of δ computed for the other atoms, occurring as a result of dehydration, have been observed, which is also in agreement with experimental results.
It should also be noted that the differences in the values of δ found between the different versions of HH2W, namely ULBL, BLUR, and ULUR, within the same carbon atom, are much smaller than those observed for the NSH.Also, the difference between the highest and lowest energy forms of HH2W, ULBL, and ULUR is much lower than in the case of MH2W (Table 6)-0.188versus 0.644 kcal/mol, respectively; this indicates that, unlike for MH2W, the various forms of HH2W are energetically similar and that the order in which the water molecules are removed from the crystal lattice is more random for HH than it was for NSH.

Molecular Dynamics Simulations
Despite the high accuracy in predicting the structure, energy, and NMR parameters, the ab initio CASTEP calculations were computationally too expensive to perform the Molecular Dynamics (MD) simulations using them.Therefore, for this purpose, we have chosen the DFTB semi-empirical method, which has been proven to provide accurate results at an affordable computational cost [33].
Before the MD simulations, all of the structures were optimized; the results of those calculations are presented in Tables 10 and 11.In the case of the peak at 148 ppm, the authors of [27] have noticed a slight downfield shift, and this was also observed in the calculations results for C11 as the increase of the δ by 0.4 ppm.The change of this magnitude did not result in the peak separation, only with the increase of its width.
For the peak at 138 ppm, the experimentalists have observed the formation of another peak at the right shoulder of the one present in HH, resulting from partial dehydration; this was also present in the calculation results for C3-the upfield change by 0.5 ppm was observed as the result of the dehydration.
No major changes in the values of δ computed for the other atoms, occurring as a result of dehydration, have been observed, which is also in agreement with experimental results.
It should also be noted that the differences in the values of δ found between the different versions of HH2W, namely ULBL, BLUR, and ULUR, within the same carbon atom, are much smaller than those observed for the NSH.Also, the difference between the highest and lowest energy forms of HH2W, ULBL, and ULUR is much lower than in the case of MH2W (Table 6)-0.188versus 0.644 kcal/mol, respectively; this indicates that, unlike for MH2W, the various forms of HH2W are energetically similar and that the order in which the water molecules are removed from the crystal lattice is more random for HH than it was for NSH.

Molecular Dynamics Simulations
Despite the high accuracy in predicting the structure, energy, and NMR parameters, the ab initio CASTEP calculations were computationally too expensive to perform the Molecular Dynamics (MD) simulations using them.Therefore, for this purpose, we have chosen the DFTB semi-empirical method, which has been proven to provide accurate results at an affordable computational cost [33].
Before the MD simulations, all of the structures were optimized; the results of those calculations are presented in Tables 10 and 11.
During the simulations, no significant conformational changes have been observed (Figures S1 and S2).Due to the experimentally determined very low stability of dehydrated HH and its almost instant amorphization, we have anticipated some changes either in the conformations of the molecules or in the unit cell dimensions.However, none of them have been observed; this could be possibly caused by a too-short simulation time, 20 ps, or an inadequate level of theory.However, during our previous studies [34,35], we have found that major changes accompanying polymorphic phase transitions occur during the first few picoseconds of simulations.However, in those mentioned cases, we have performed During geometry optimization, all atoms' positions and the cell parameters were optimized with no constraints.The convergence criteria were set at 5 × 10 −6 eV/atom for the energy, 1 × 10 −2 eV/Å for the interatomic forces, 2 × 10 −2 GPa for the stresses, and 5 × 10 −4 Å for the maximum displacement.The fixed basis set quality method for the cell optimization calculations and the 5 × 10 −7 eV/atom tolerance for SCF were used.

NMR Parameters Calculations
The computation of shielding was performed using the Gauge Including Projector Augmented Wave Density Functional Theory (GIPAW) method of Pickard et al. [43].To compare the theoretical and experimental data, the calculated chemical shielding constants (σiso) were converted to chemical shifts (δiso) using the following equation: δiso = (σGly + δGly) − σiso, where σGly and δGly stand for the shielding constant and the experimental chemical shift, respectively, of the glycine carbonyl carbon atom (176.50 ppm).

Periodic DFTB Calculations
The Density Functional Tight Binding (DFTB) calculations of geometry optimization and molecular dynamics simulations under periodic boundary conditions were carried out with the DFTB+ program [44] implemented in the Materials Studio 2020 software [37].The calculations utilized a library containing Slater-Koster atomic parameters, incorporating the UFF-based Lennard-Jones dispersion corrections and charge self-consistency.

Geometry Optimization
Geometry optimization was carried out using the "smart" algorithm for calculations, the "divide and conquer" method for diagonalizing the Hamiltonian (eigensolver), the Broyden charge mixing scheme, and the Methfessel-Paxton distribution function used for smearing.The number of Monkhorst-Pack k-points during sampling for a primitive cell Brillouin zone integration [42] were set to 4 × 1 × 2 (for thiamine monohydrate-based structures) and 1 × 4 × 1 (for thiamine hemihydrate based structures), respectively.The details on the structure preparation can be found in Section 2.1.
During geometry optimization, all atoms' positions and the cell parameters were optimized with no constraints.The convergence criteria were set at 1 × 10 −2 kcal/mol for the energy, 5 × 10 −2 kcal/mol/Å for the interatomic forces, 2 × 10 −2 GPa for the stresses, and 5 × 10 −4 Å for the maximum displacement and the 1 × 10 −8 kcal/mol tolerance for SCC were used.

Molecular Dynamics Simulations
Molecular dynamics (MD) simulations were run using an NPT ensemble maintained at a constant temperature of 293 K and pressure of 0.01 GPa, using a Nosé thermostat with 0.01 Q ratio and Berendsen barostat with 0.1 ps decay constant.The time step was set to 0.5 fs, and the total time of the simulation was set to 20 ps.All of the settings and electronic options were set at the same values as for geometry optimization (Section 3.2.1).No symmetry constraints were applied during the simulations.

Conclusions
In this work, two different forms of THCL hydrates, NSH and HH, have been studied using various molecular modeling methods.The first step of this work included the creation of the structures of hydrates with decreasing water: THCL molar ratio, starting either from monohydrate (NSH) or hemihydrate (HH) up to the fully dehydrated forms.After the optimization of all of the structures, we have found a good agreement between theoretically modeled and experimentally determined values describing the unit cell dimensions; this indicates that modeled structures obtained for the forms that have not been studied experimentally yet should also be accurate.Also, a comparison of the energy of the structures based on NSH and HH was in agreement with the experimental findings, indicating higher relative stability of the NSH.
GIPAW NMR calculations performed for the optimized structures allowed not only to assign all of the peaks in the experimental 13 C CP MAS NMR spectra to particular carbon atoms but also to explain, at the molecular level, all of the changes observed in the spectra occurring as a result of the experimental dehydration.Comparison of the NMR calculations results with the energy of the partially dehydrated structures enables us to predict the structure of the intermediate forms occurring between the fully hydrated forms of NSH and HH and their dehydrated counterparts.
During the molecular dynamics simulations performed at the semi-empirical QM level, we did not observe major changes in the studied structures.For all of the forms of NSH, those results were anticipated, as this form is a variable hydrate that can exist at various THCL: water ratios, depending on the air humidity.However, the lack of changes in the unit cell dimensions of the dehydrated structure of HH is intriguing, as the previous experimental works reported on the instability of this form.We aim to investigate this aspect deeply in the near future.
This work highlights the accuracy and versatility of the "solid-state DFT" calculations in the analysis of various forms of molecular solids, in particular solvates of active pharmaceutical ingredients of various degrees of hydration.

Figure 4 .
Figure 4. Experimental (exp I, exp II) and calculated unit cell dimensions of MH; for better comparison, the equal range of vertical axis-1.2Å-has been used for all three unit cell lengths (a, b, and c).For MH2W, the arithmetic mean of the results (MH2W BLUR, MH2W ULBL, and MH2W ULUR) was shown.

Figure 4 .
Figure 4. Experimental (exp I, exp II) and calculated unit cell dimensions of MH; for better comparison, the equal range of vertical axis-1.2Å-has been used for all three unit cell lengths (a, b, and c).For MH2W, the arithmetic mean of the results (MH2W BLUR, MH2W ULBL, and MH2W ULUR) was shown.

Figure 5 .
Figure 5. Experimental (exp I) and calculated unit cell dimensions of HH; for better comparison, the equal range of vertical axis-1.4Å-has been used for all three unit cell lengths (a, b, and c).For HH2W, the arithmetic mean of the results (HH2W BLUR, HH2W ULBL, and HH2W ULUR) was shown.

Figure 5 .
Figure 5. Experimental (exp I) and calculated unit cell dimensions of HH; for better comparison, the equal range of vertical axis-1.4Å-has been used for all three unit cell lengths (a, b, and c).For HH2W, the arithmetic mean of the results (HH2W BLUR, HH2W ULBL, and HH2W ULUR) was shown.

Figure 6 .
Figure 6.Carbon atom numbering of THCL used in the NMR analysis.

Figure 7 .
Figure 7. Chosen regions of the13 C CP MAS NMR spectra of NSH at various degrees of hydration (0.90, 0.41, and 0.04).Reprinted from[28] with permission from Elsevier.

Figure 7 .
Figure 7. Chosen regions of the13 C CP MAS NMR spectra of NSH at various degrees of hydration (0.90, 0.41, and 0.04).Reprinted from[28] with permission from Elsevier.

Figure 9 .
Figure 9.The chosen region of the 13 C CP MAS NMR spectra of HH and the partially dehydrated form of HH named HH2.Asterisk indicates spinning sideband.Reprinted from [27] with permission from Elsevier.

Figure 9 .
Figure 9.The chosen region of the 13 C CP MAS NMR spectra of HH and the partially dehydrated form of HH named HH2.Asterisk indicates spinning sideband.Reprinted from [27] with permission from Elsevier.

Figure 9 .
Figure 9.The chosen region of the 13 C CP MAS NMR spectra of HH and the partially dehydrated form of HH named HH 2 .Asterisk indicates spinning sideband.Reprinted from [27] with permission from Elsevier.

Figure 10 .
Figure 10.Optimized crystal unit cell of HH4W; pink dashed lines represent the intermolecular interactions.Atom coloring: N-blue, C grey, S-yellow, H-white, and O-red.The colors of symmetry equivalent Cl atoms are either violet (Cl1) or green (Cl2).

Figure 10 .
Figure 10.Optimized crystal unit cell of HH4W; pink dashed lines represent the intermolecular interactions.Atom coloring: N-blue, C grey, S-yellow, H-white, and O-red.The colors of symmetry equivalent Cl atoms are either violet (Cl1) or green (Cl2).

Table 2 .
Structures used for calculations prepared from the

Table 2 .
Structures used for calculations prepared from the THIAMC12

Table 3 .
Structures used for calculations prepared from the

Table 3 .
Structures used for calculations prepared from the WUWJAA-HH.

Table 5 .
Optimized unit cell dimensions of the HH structures of various hydration ratios, compared with the experimental ones.

Table 7 .
Dehydration energies are defined as the energy required to remove one water molecule for each step (1-4) of this process.HH-hemihydrate, MH-monohydrate.

Table 8 .
Calculated chemical shift values for carbon atoms, compared with the experimental ones.Atom numbering is presented in Figure