Computation predictions of organic crystal structure and thermodynamics are essential for material design, crystal engineering and drug development. However, accurate computational tools for organic crystal thermodynamics calculations are lacking, and experimental data set for validation of computational methods is limited. Most crystal structure predictions and stability calculations depend solely on potential energy, which is an insufficient representation of thermodynamics. This thesis proposes and validates both absolute and relative free energy calculation of small organic compounds, thus presenting an accurate computational tool that overcomes the shortcomings of potential-energy-based models.
The solubility of organic molecules can be computed from a thermodynamic cycle that decomposes standard state solubility into the sum of solid-vapor sublimation, i.e. thermodynamic stability of the crystal, and vapor-liquid solvation free energies ΔG°solubility=ΔG°sub+ΔG°solv. Crystal polymorphs have different ΔG°sub thus different solubility, which of critical importance to the pharmaceutical industry, however, robust computational methods to predict this quantity from first principles are lacking. Over the past few decades, alchemical simulation methods to compute solvation free energy using classical force fields have become widely used. However, analogous methods for determining the free energy of the sublimation/deposition phase transition are currently limited. This thesis describes an absolute thermodynamics approach based on growth of the asymmetric unit into a crystal via alchemy (GAUCHE). GAUCHE computes deposition free energy ΔG°dep=-ΔG°sub=ΔG°Vol+ΔG°Au+ΔG°Au→UG as the sum of an entropic term to account for compressing a 1 M vapor into the molar volume of the crystal asymmetric unit (VAU) plus two simulation steps. In the first simulation step, the deposition free energy ΔG°AU for a system composed of only NAU asymmetric unit (AU) molecule(s) is computed beginning from an arbitrary conformation in vacuum. In the second simulation step, the change in free energy ΔG°AU�UG to expand the asymmetric unit degrees of freedom into a unit cell (UC) composed of NUC independent molecules is computed. This latter step accounts for the favorable free energy of removing the constraint that every symmetry mate of the asymmetric unit has an identical conformation and intermolecular interactions. The current work is based on NVT simulations, which requires knowledge of the crystal space group and unit cell parameters from experiment, but not a priori knowledge of crystalline atomic coordinates. GAUCHE was applied to 5 organic molecules whose sublimation free energy has been measured experimentally, based on the polarizable AMOEBA force field and more than a microsecond of sampling per compound. The mean unsigned and root-mean-square errors were only 1.6 and 1.7 kcal/mol, respectively, which indicates that GAUCHE is capable of accurately predicting sublimation thermodynamics.
For polymorphic systems, we propose a relative thermodynamics approach, that is similar to the second simulation step of GAUCHE, where ΔG°P1→P2 is calculated instead of ΔG°AU→UG. A relative approach reduces statistical uncertainty upon convergence compared to an absolute calculation; thus, it is more appropriate due to the thermodynamic stability difference between polymorphs are often fairly small. For our paracetamol test system, the experimental free energy difference was only 0.93 kcal/mol. Although both quantum and AMOEBA potential calculations predict the form II of paracetamol as more stable crystal form than form I, our relative free energy calculation predict the opposite stability ranking, which agrees with the experiment. Decomposition of free energy into entropy and enthalpy indicates that the favorable entropy change contributes to the greater thermodynamic stability of form I over form II. Although the exact magnitude of entropy and enthalpy changes differs across literature data as well as our data, the favorable entropic contribution is consistent. Further calculations over the temperature range from 100 to 308 K show the temperature dependence of free energy, which follows the parabolic trend observed in experiments. Our results show that our relative polymorph stability methods can accurately capture temperature dependence of free-energy-based stability ranking and overcome the limitations of potential-energy-based ranking. Thus, the importance of crystal structure predictions based on free energy can be further emphasized.