Adding anisotropy to the standard quasi-harmonic approximation still fails in several ways to capture organic crystal thermodynamics
AAdding anisotropy to the standardquasi-harmonic approximation still fails in severalways to capture organic crystal thermodynamics
Nathan S. Abraham and Michael R. Shirts ∗ Department of Chemical and Biological Engineering, University of Colorado Boulder,Boulder, CO 80309, USA
E-mail: [email protected]
Abstract
We evaluate the accuracy of varying thermalexpansion models for the quasi-harmonicapproximation (QHA) relative to moleculardynamics (MD) for 10 sets of enantiotropicorganic polymorphs. Relative to experi-ment we find that MD, using an off-the-shelf point charge potential gets the signof the enthalpic contributions correct for6 of the 10 pairs of polymorphs and thesign of the entropic contributions correct forall pairs. We find that anisotropic QHAprovides little improvement to the error infree energy differences from MD relative toisotropic QHA, but does a better job cap-turing the thermal expansion of the crystals.A form of entropy–enthalpy compensationallows the free energy differences of QHA todeviate less than 0.1 kcal/mol from MD formost polymorphic pairs, despite errors up to0.4 kcal/mol in the entropy and enthalpy. Deviations in the free energy of QHA andMD do not clearly correlate with molecu-lar flexibility, clarifying a previously pub-lished finding. Much of the error previouslyfound between QHA and MD for these flex-ible molecules is reduced when QHA is runfrom a lattice minimum consistent with thesame basin as MD, rather than the energy-minimized experimental crystal structure.Specifically, performing anisotropic QHAon lattice minimum quenched from low-temperature replica exchange simulationsreduced the error previously found by 0.2kcal/mol on average. However, these con-formationally flexible molecules can havemany low-temperature conformational min-ima, and the choice of an inconsistent min-ima causes free energies estimated fromQHA to deviate from MD at temperaturesas low as 10 K. We also find finite size er-rors in the polymorph free energy differ-ences using anisotropic QHA, with free en-1 a r X i v : . [ c ond - m a t . m t r l - s c i ] O c t rgy differences as large as 0.5 kcal/mol be-tween unit and supercells loosely correlatedwith differences in anisotropic thermal ex-pansion. These larger system sizes are com-putationally more accessible because ourcheaper 1D variant of anisotropic QHA,which gives free energies within within 0.02kcal/mol of the fully anisotropic approachat all temperature studied. The errors be-tween MD and experiment are 1–2 orders ofmagnitude larger than those seen betweenQHA and MD, so the quality of the forcefield used is still of primary concern, butthis study illustrates a number of other im-portant factors that must be considered toobtain quantitative organic crystal thermo-dynamics. Introduction
Solid organics can pack stably in multi-ple forms, each of which can have sig-nificantly different chemical and materialsproperties from each other. The ability of amolecule to arrange in multiple solid forms,or polymorphs, can change crystal solu-bility, charge mobility, hardness, andreactivity. Relative stability of crystalscan change with temperature and pres-sure, so we must consider entropic ef-fects and crystal lattice changes when mod-eling the crystal thermodynamics.Methods such as molecular dynamics(MD) and the quasi-harmonic approxima-tion (QHA) are important for understand-ing the relative stability of crystals. QHAestimates the thermodynamics as a func-tion of temperature and pressure by com-puting free energy due to harmonic vibra-tions around each crystal lattice minimum,and then determining the lattice geome-try that minimizes the Gibbs free energyat each temperature as the crystal is ex-panded from the initial crystal lattice min- imum. It is common to model thermal ex-pansion isotropically, where the lattice vec-tors remain proportional and angles fixed.MD provides a more accurate descriptionof the entropy and therefore free energy bygenerating the full conformational ensembleof the crystal, rather than neglecting anhar-monic motions. Additionally, QHA assumesa single lattice is representative of the crys-tal at low temperatures whereas MD sam-ples an ensemble of lattice configurations.Our previous studies show that anisotropic QHA model can often yield freeenergy differences between polymorphs thatfall within numerical error of the polymorphfree energy differences computed by MD forsmall and rigid molecules. However, free en-ergies from QHA deviates from those gen-erated by MD by > The small-molecule crystals examinedin our previous study all showed some levelof anisotropic expansion, and the error inthe QHA free energy qualitatively increasedwith the degree of anisotropy. We concludedthat the deviations of the isotropic QHAfrom MD were either due to 1) an inaccuratethermal expansion model or 2) anharmonicmotions in the crystal lattice, or 3) somecombination of the two. Without modelingQHA anisotropically, we could not defini-tively determine the source of the deviationsof QHA from MD.We recently developed a method to ef-ficiently determine the true anisotropic freeenergy minimum within the quasi-harmonicframework. This method relies on deter-mining gradient of the lattice parameterswith respect to temperature. Our gradi-ent approach identifies high temperaturefree energy minimum that are 0.01–0.23kcal/mol lower than the minimum identi-fied with an isotropic model, altering thepolymorph free energy differences of pirac-2tam and resorcinol by 0.02–0.12 kcal/molfor unit cells of 4–8 molecules. Our fullyanisotropic approach and even a simplified1D-anisotropic approach outperformed theaccuracy of current quasi-anisotropic meth-ods that only consider the anisotropy due tothe potential energy as a function of crystallattice vectors.
In this paper we re-evaluate the 10 ofthe 12 polymorph pairs previously studiedwith MD and isotropic QHA to determinewhether an anisotropic thermal expansionmodel can improve the performance of QHArelative to MD. Specifically, we: • evaluate the use of an off-the-shelfpoint charge potential in modeling theentropy, enthalpy, stability, and ther-mal expansion of organic polymorphs; • determine how well our 1D-QHA vari-ant compares to QHA using fullanisotropic expansion over a largerdata set than our initial anisotropicmethods paper; • determine if finite size effects signif-icantly contribute to the free energyranking of polymorphs; • discuss the appropriate ways to com-pare QHA and MD; • determine if anisotropic QHA can re-duce errors with respect to MD foundpreviously with isotropic QHA; • evaluate if entropy–enthalpy compen-sation plays a significant role in differ-ences between QHA and MD; and • highlight crystal behaviors found inMD that QHA cannot account for. Methods
In this study we examine the experimen-tally known polymorphs of the 10 moleculesshown in Figure 1 using QHA and MD. All10 of these molecules were examined in ourprevious study. Most MD results in thisstudy are from our previous publication, and a methodological discussion of MD canbe found there. Modifications of the MDprocedure and where these modified simula-tions are used will be discussed in the con-text of the paper. All isotropic QHA re-sults have been updated from that previousstudy and we describe those changes madehere. In our previous study we examined12 enantiotropic pairs of polymorphs, buthave excluded two molecules due to the ex-perimentally known plastic phases with dy-namic disorder at each molecular site. Thedynamic disorder was confirmed with MD,both by others and ourselves, for the rota-tion of cyclopentane and trans/gaucheisomerization of succinonitrile in thedisordered crystals. This behavior is in-herently anharmonic, meaning there is nochance that QHA can properly model them,and hence they are left out of this study.Furthermore, both plastic phase crystalshad their lattices constrained for stabilityin MD in previous studies, so they werenot modeling fully anisotropic thermal ex-pansion. Quasi-Harmonic Approxima-tion
In the quasi-harmonic approximation(QHA) the Gibbs free energy is computedby determining the lattice geometry thatminimizes the sum of the potential and theenergy of the static, harmonic lattice vibra-3 esorcinol TolbutamidePyrazinamide Paracetamol Diiodobenzene ChlorpropamideCarbamazepine Piracetam Adenine Aripiprazole
NNH N NNH OHHONN NH O NO O NH HO HN O IIN NH O S NHO O NHO Cl S NHO O NHOCl Cl N N O HN O
Figure 1: Molecules examined in this study.tions, as shown in eq 2. G ( T ) = min y f ( y, T ) (1) f ( y, T ) = min x ( U ( y, x ))+ A v ( y, T ) + P V (2) y = V, C , λ where T is the temperature, P is the pres-sure, min x ( U ( y, x )) is the lattice energy ofthe minimum energy cell geometry y and ge-ometry optimized lattice coordinates x , and A v is the Helmholtz vibrational free energyof a harmonic oscillator. We define the cellgeometry as y , which for isotropic expan-sion is the volume V , anisotropically withthe six dimensional lattice tensor C , or us-ing our modified 1D-anisotropic expansiondefined by λ .The Helmholtz free energy ( A v ) of thelattice harmonic vibrations contribute tothe entropic portion of QHA. We use theclassical version of A v for comparison to MD, shown in eq 3. A v ( V, T ) = (cid:88) k β − ln ( β ¯ hω k ( V )) (3)where β is ( k B T ) − , ¯ h is the reduced Planckconstant, and ω k is the frequency of thephonon of the k th vibrational mode of thecrystal lattice. Gr¨uneisen Parameters
There are two common methods within theQHA framework to estimate the phonon fre-quencies of the crystal lattice at each latticegeometry of interest. The first way is tocompute and diagonalize the mass-weightedHessian of every crystal structure, which iscomputationally demanding. The secondway is to use the Gr¨uneisen parameter ap-proach assumes that the changes in the fre-quencies of a particular phonon are constantas the crystal is strained in a particular di-rection. Use of the Gr¨uneisen parameterto calculate the thermodynamics of organicand inorganic material is common. This is smallerthan the experimental error in most cases,so we will exclusively use the Gr¨uneisen pa-rameter approach for QHA calculations inthis paper when calculating polymorph freeenergy differences.The standard definition of the Gr¨uneisendue to a volume change is in eq 4, which isdirectly applicable only for constant strains,such as isotropic expansion, is: γ k = − Vω k ∂ω k ∂V (4)where γ k is the Gr¨uneisen parameter for the k th vibrational frequency. Eq 4 is solved nu-merically, which requires diagonalization ofthe mass-weighted Hessian at two volumesto produce reference frequencies. The corre-sponding Gr¨uneisen parameters can be usedto solve the frequencies at all isotropic vol-umes relative to a reference point, generallythe lattice minimum structure.The Gr¨uneisen parameter for the vol-ume can be extended to a crystal placedunder any strain, which we will usefor anisotropic expansion. Eq 5 is theanisotropic version of the isotropic equa-tion(eq 4). γ k,i = − ω k ∂ω k ∂η i (cid:12)(cid:12)(cid:12)(cid:12) η j (cid:54) = η i (5)Here, the Gr¨uneisen parameter for the k th vibrational mode due to the i th strain ( η i )applied to the crystal is γ k,i . The symmet-ric 3 × η allows us to com-pute 6 sets of Gr¨uneisen parameters nu-merically. With the Gr¨uneisen parametersand reference frequencies, we can computethe frequencies of the lattice parameters atany cell geometry of interest, by integrating eq 5. A full discussion of the use of both theisotropic and anisotropic Gr¨uneisen param-eters can be found in previous work. Thermal Expansion
We previously implemented a method to de-termine the crystals thermal expansion, which can be used to solve bi-optimizationproblems, such as QHA. Eq 6 gives thegeneral formulation for determining thecrystal thermal expansion. If we have areference structure, the 0 K minimum (i.e.classical lattice minimum), we can computethe thermal expansion and numerically in-tegrate with temperature. ∂y∂T = (cid:18) ∂ G∂y (cid:19) − ∂S∂y (6) y = V, λ, C In eq 6, T is the temperature, S is the en-tropy, and G is the Gibbs free energy. Bycomputing the numerical derivatives of theentropy and Gibbs free energy we can de-termine the gradient of y , the cell geometry,with respect to temperature. Isotropic Expansion
Isotropic expansion assumes that the latticevectors remain proportional to one anotherand the lattice angles remain fixed. Dueto these constraints the solution to QHA ineq 2 requires minimization of a single vari-able, the volume ( y = V ). We computethe rate of isotropic thermal expansion us-ing eq 6 and replacing y with the crystalvolume ( V ). Anisotropic Expansion
Anisotropic expansion allows the crystal lat-tice to relax to the harmonic free energyminimum structure by removing all con-straints in isotropic expansion. In eq 2 weminimize the free energy as a function ofall six crystal lattice parameters ( y = C ),which makes the problem of minimizing the5ibbs free energy for QHA more complex.We compute the thermal expansion for allsix parameters by using eq 6 and replacing y with a array of the six lattice parameters( C ). Computing eq 6 for anisotropic expansion( C ) requires 73 lattice optimizations to de-termine a single six-dimensional gradient,which becomes increasingly expensive tocompute for the entire temperature range ofinterest. In our previous work, we foundif the ratio of anisotropic expansion waskept constant at all temperatures we couldachieve the same free energy differenceswithin 0.005 kcal/mol of full anisotropic ex-pansion. In that paper, we presented eq 7, C i ( λ ) = C i ( λ = 0) + κ i λ ( T ) (7)where the lattice parameter C i is a functionof the variable λ and κ i . In this case, λ isa single parameter describing the expansionin eq 2, which is zero at the 0 K lattice min-imum value of C i . κ i is the gradient com-puted in eq 6 for y = C at 0 K. Molecular Dynamics
A more complete approach to compute thethermodynamics of the crystals, includingthe free energy difference between poly-morphs, is molecular dynamics (MD), whichgenerates (in theory) the full configurationalensemble at a given temperature of inter-est. For our method, there are two nec-essary steps for computing the free energydifferences of polymorphs using MD: 1) de-termine ∆ G between polymorphs at a refer-ence temperature with simulations of a se-ries of non-physical intermediate states and2) determine ∆ G as a function of temper-ature for each polymorph using simulationsat range of temperatures. For the first step,we can drive each crystal along a reversible thermodynamic path to an ideal gas stateby turning off intermolecular energies todetermine the reference free energy differ-ences. Then, with simulations at intermedi-ate temperatures we can relate each temper-ature point back to the reference energy dif-ference between polymorphs using the Mul-tistate Bennett Acceptance Ratio (MBAR).A full discussion of this approach can befound in our previous work. Computational Details
We have provided as supporting infor-mation a .zip file including all computa-tional details to re-run simulations for bothMD and QHA. For specifications on howwe chose certain parameters we refer thereader back to our previous publicationsthat present the methods for MD andQHA. Molecular Dynamics
For a number of crystals, we found thatthe lattice parameters from the previousstudy did not continuously change withtemperature at low T for MD, so we re-ran the crystals with temperature replicaexchange (REMD) to better escape frommetastable states and thus allow for bet-ter convergence to the 0 K lattice energyminimum. Replica exchange parameterswere chosen to give an average exchangeprobability of ∼ . Tem-perature replica exchange was required tocompute the low-temperature ensembles offor tolbutamide, chlorpropamide, and arip-iprazole. We provide input files for bothMD and QHA all within the input files.zip,6hich contains all parameters for individualcrystals and simulation settings for GRO-MACS 2018.
Quasi-Harmonic Approximation
All QHA calculations were performed us-ing our Python based lattice dynam-ics code available on GitHub at http://github.com/shirtsgroup/Lattice_dynamics . The code currently wrapsaround a number of molecular modelingpackages. In this paper, all lattice vibra-tions, energy evaluations, and optimizationswere performed using Tinker 8.1 molecularmodeling distributed TINKER code to cor-rect for errors in the Hessian calculation.Those modifications are discussed in themain paper and supporting information ofour previous work and have since been up-dated in Tinker 8.7. Lattice structures were retrieved fromthe Cambridge Crystallographic DataCenter ( ) and a unit and supercell werelattice optimized. Supercells were createdassuring that the three lattice vectors weretwice the van der Waals cutoff, (8 ˚A). Spec-ifications on cell size can be found in theTable S3. Each molecule was parameter-ized using the OPLS-AA classical fixedpoint charge potential: U total = U bond + U angle + U dihedral + U vdw + U coulombic (8)where the exact forms of the each of the en-ergy terms are defined in the references.The crystal structures were geometryand lattice optimized to the lattice mini-mum structure. Optimization of the crys-tal structure was performed with Tin-ker’s xtalmin executable to an RMS gra-dient/atom of 10 − . The molecules centersof mass, keeping the coordinates relative tothe center of mass constant constant. Thecrystal was then geometry optimized using Tinker’s minimize executable to an RMSgradient/atom of 10 − . These values werechosen to maximize convergence and numer-ical stability for the gradient method, andthis choice is discussed fully in the support-ing information of our previous work. To determine y ( T ) across the entire tem-perature range we use a 4 th order Runge-Kutta integrator using the thermal gradi-ent approach to satisfy eq 1. For most crys-tals, we found that 3 steps of 100 K eachup to 300 K ensures that the crystal was ata locally metastable free energy minimumat all temperatures. If the structure is notat a free energy minimum with respect tobox geometry variables at 300 K, the struc-ture was re-run with 6 Runge-Kutta steps.For the unit cells only, if 6 steps could notmaintain the crystal at a free energy min-imum then 20 steps were run. All resultsshown are for the maximum temperaturestep that could be achieved while assuringthat it remained at a free energy minimumwith respect to box geometry under expan-sion (temperatures given in Table S4). Re-sults for the graphed intermediate pointsbetween Runge-Kutta steps were calculatedwith a third order spline that was fit to thelattice parameters with temperature. Details of experimental dataused for comparison
We have exhaustively collected experimen-tal results reported in literature to com-pare thermal expansion and thermodynamicstability to for MD. In previous workwe compared entropic and enthalpic con-tributions to polymorph relative stabili-ties, but have found some inconsistenciesin those reported results and reevalu-ate these calculations here. Provided withour supporting information are two files,experimental expansion.csv and experimen-7al stability.csv, containing the literatureresults we used for comparison. In theexpansion results file we provide the re-ported experimental lattice parameters, ci-tation references and links, reported tem-peratures, and reference codes if the struc-ture was submitted to the CCDC. For thestability results we report the values seen inFigure 2 along with the reported tempera-tures and literature references. Results and Discussion
The approaches presented above allow us toevaluate the effectiveness of QHA methodswith different treatments of thermal expan-sion relative to MD for polymorph thermo-dynamics, including free energy differencesbetween polymorphs, and lattice expansion.The differences between approaches for es-timating the temperature dependence aregreatest at high temperatures, so all com-parisons are shown at the maximum tem-perature at which QHA is stable.For QHA, we found that discontinuousreadjustment of the molecules within thelattice occurs frequently during expansion,causing the method to expand to a struc-ture that no longer corresponds with freeenergy minimum of the chosen expansionvariable. As the crystal expands, the freevolume around the molecules increases, al-lowing molecules to readjust into alternateminima more favorable to the lattice energythan the minima found by continuous defor-mation of lower temperature minima. Thesestructural disruptions are problematic forthe numerical stability of the gradient ap-proach and generally cause the crystal to ex-pand to a structure that is not a free energyminimum. For example, upon expansion oftolbutamide, at a certain point the alkyl tailcan move into a newly created free volume,causing the quasi-harmonic free energy to become discontinuous with temperature.We determine if the crystal is at a freeenergy minimum with respect to the geom-etry parameter y by checking if: (cid:18) ∂G∂y (cid:19) Backward < < (cid:18) ∂G∂y (cid:19) F orward (9)If our crystal is at or near a free energy min-imum, then the forward and backward nu-merical solution to eq 9 must be positiveand negative respectively. If both are thesame sign, then we must no longer be at aminimum.All results using QHA are shown atthe maximum temperature ( T max ), whichis the highest temperature to satisfy eq 9.The values of T max are reported in Ta-ble S4 and Table S5. Plots of the poly-morph free energy differences and latticeexpansions versus temperature for the 10molecules studied are provided (Figures S20– S56). Further quantification of the numer-ical and structural instability is discussed inthe work where we initially presented thegradient method. Comparisons to ExperimentalResults
We found that the OPLS-AA point chargepotential is a poor predictor of the sign ofenthalpic differences for polymorphs rela-tive to experiment, as is expected for thesimplicity of the energy function, but doesaccurately estimate the sign of the entropicdifferences. In Figure 2a and 2b the en-thalpic and entropic contributions to thepolymorph free energy differences at 300 Kfor the supercells using MD are shown rela-tive to their experimental values, which arereported at varying temperatures that gen-erally corresponding with the polymorphtransition temperature.
Thesign of the enthalpic contributions is only8 .0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 H Exp . [kcal/mol] H M D [ k c a l / m o l ] a) RMSD = 4.0 kcal/mol
T S
Exp . [kcal/mol] T S M D [ k c a l / m o l ] b) RMSD = 2.0 kcal/mol
Figure 2: 300 K polymorph differences from MD for a) enthalpy and b) entropy versus exper-imental values taken from literature.
The sign of the enthalpy is only correctfor MD for 6 of the 10 molecules, with no clear distinguishing molecular or crystallographicfeatures separating the crystals with the incorrect sign. This contrasts with the entropy,which is correct in sign for all 10 molecules. If the polymorph pairs fall in the green quad-rants than MD has the same sign as experiment, and the wrong sign in the red quadrants.9orrect for 6 of the 10 polymorph pairs.The 4 molecules where the simulation hasthe wrong enthalpic ranking are pirac-etam, pyrazinamide, resorcinol, and parac-etamol, and none of the enthalpic errors inthese molecules have obvious correlations tochemical groups present, flexibility, crystalZ and Z’ values, or relative packing. Thiscontrasts with the polymorph entropy dif-ferences, where the sign is correct for MD forall 10 molecules. The RMSD between the-ory and experiment are 4.0 and 2.0 kcal/molfor ∆ H and T ∆ S , respectively. The mostlikely large sources of deviation from exper-iment could be due to 1) the experimentalvalues coming from the polymorph transi-tion temperature, not 300 K and 2) the useof a point charge potential to model thecrystal energetics.The lattice geometries for MD differmoderately from experiment, though thevolumes match experiment more closely.We can quantitatively compare to experi-ment using exhaustive experimental resultsfrom Brandenburg and Wilson for carba-mazepine form III and paracetamol form I,respectively. The expansion of the lat-tice vectors, angles, and volume are shownfor carbamazepine form III (Figure 3) andparacetamol form I (Figure 4). Predictionaccuracy of lattice vectors varies both withcrystal and lattice vector. We find that ifthe lattice vectors are orthogonal for exper-iment, MD always gets the correspondingangle correct or within error as seen with α and β angles for carbamazepine and parac-etamol. Despite errors in the individual lat-tice parameters, the crystal volume is sim-ilar between MD and experimental values.In the case of paracetamol (Figure 4) theexperimental volume is within error of MD.The change in the lattice vectors and or-thogonal angles in MD generally agrees withexperiment despite errors in the actual ge-ometry. In Figures 3 and 4 we provide the thermal expansion of each parameter by fit-ting a linear fit to MD and experimentaldata at all temperatures. In carbamazepinethe change in the slope of experimental lat-tice parameters at low temperatures couldbe indicative of quantum behavior from thezero point energy, so for this crystal we per-form the linear fit on values above 100 K.For all of the lattice vectors we see thatthe thermal expansion between MD and ex-periment are the same order of magnitudeand sign. For the orthogonal angles there isno expansion. However, the thermal expan-sion of the β angle in both polymorphs hasthe wrong sign and in the case of paraceta-mol form I is an order of magnitude differ-ent. For carbamazepine, we see larger diver-gence between experiment and MD at lowtemperatures ( T <
100 K), which would beindicative of quantum zero-point energy ef-fects. For both polymorphs there is alsogood agreement with experiment for thevolumetric thermal expansion. Accuratelymodeled thermal expansion helps to explainwhy we also get the correct sign of the en-tropic differences for all 10 polymorph pairsin Figure 2b.
Testing the Validity of 1D-Anisotropic QHA
Anisotropic expansion is relatively ex-pensive compared to isotropic expansion,though still 2 orders of magnitude cheaperthan free energy estimation with MD. Ourprevious work shows that a 1D-variant ofanisotropic expansion can be a sufficientsubstitute to speed up the method withlittle effect on the accuracy of computedfree energy differences. Fully anisotropicexpansion requires 73 structure optimiza-tions every time the thermal gradient iscomputed. This contrasts with the 1D-approach, which requires the initial 73 op-10
100 200 3007.507.55 L a tt i c e V e c t o r [ Å ] dadT = 2.3×10 [Ang./K] dadT = 3.2×10 [Ang./K]a-Vector MDBrandenburg 2017 Hassan 2013Horstman 2015 Lisgarten 1989 Nievergelt 2017 Sovago 2016 dbdT = 6.2×10 [Ang./K] dbdT = 5.9×10 [Ang./K]b-Vector dcdT = 8.3×10 [Ang./K] dcdT = 6.6×10 [Ang./K]c-Vector A n g l e [ D e g . ] ddT = 0.0 [Deg./K] ddT = 0.0 [Deg./K]-Angle ddT = -2.9×10 [Deg./K] ddT = 0.8×10 [Deg./K]-Angle ddT = 0.0 [Deg./K] ddT = 0.0 [Deg./K]-Angle Temperature [K] V o l u m e [ Å ] dVdT = 0.17 [Ang. /K] dVdT = 0.16 [Ang. /K] Temperature [K]
Figure 3: Thermal expansion of carbamazepine form III comparing MD to experimentalresults provided in literature.
While the lattice parameters differ between theory andexperiment, the thermal expansion is similar between MD and experiment for all lattice pa-rameters and volume. At low temperatures we find the largest deviation between experimentand MD. The shaded green area is the standard deviation in the lattice parameter duringthe MD simulation. 11
100 200 30012.512.612.712.812.9 L a tt i c e V e c t o r [ Å ] dadT = 7.4×10 [Ang./K] dadT = 11.0×10 [Ang./K]a-Vector MD Haisa 1976 Shankland 1997 Wilson 1997 Wilson 2009 dbdT = 6.9×10 [Ang./K] dbdT = 4.0×10 [Ang./K]b-Vector dcdT = 0.6×10 [Ang./K] dcdT = 0.4×10 [Ang./K]c-Vector A n g l e [ D e g . ] ddT = 0.0 [Deg./K] ddT = 0.0 [Deg./K]-Angle ddT = 3.4×10 [Deg./K] ddT = -54.2×10 [Deg./K]-Angle ddT = 0.0 [Deg./K] ddT = 0.0 [Deg./K]-Angle Temperature [K] V o l u m e [ Å ] dVdT = 0.14 [Ang. /K] dVdT = 0.1 [Ang. /K] Temperature [K]
Figure 4: Thermal expansion of paracetamol form I comparing MD to experimental resultsprovided in literature.
While the lattice parameters differ between theory and experi-ment, the thermal expansion is similar between MD and experiment for all lattice parametersand volume, except for the β angle. The shaded green area is the standard deviation in thelattice parameter during the MD simulation.timizations at 0 K followed by 3 optimiza-tions at all subsequent points for numericalintegration. When we first presented thegradient approach we found that the 1D-anisotropic QHA approach computed poly-morph free energy differences within 0.01kcal/mol of full anisotropic expansion forthe two sets of polymorphs tested. For thisstudy, we ran both 1D-QHA and anisotropicQHA on unit cells of the quenched experi-mental structures (independent of MD) toevaluate the effectiveness of the 1D con-straint.The constraints applied to 1D-anisotropic approach provide numerical sta-bility, allowing eq 9 to be satisfied at 300K more frequently. At each integrationstep for expansion our program checks ifeq 9 is satisfied, verifying if the crystal isat a free energy minimum. We will onlyreport result up to the temperature ( T max ) where eq 9 is satisfied. When using fullyanisotropic expansion, only 11 of the 20crystals were able to satisfy eq 9 at all tem-peratures ( T max = 300 K). This contrastswith the 1D-approach, where 19 of the 20crystals remained at a free energy minimumup to 300 K. The only molecules where thefully anisotropic approach remained at afree energy minimum for both polymorphswere the relatively rigid resorcinol, adenine,and carbamazepine, suggesting that rigidmolecules are less likely to experience nu-merical disruptions in QHA.For the 10 molecules in this study, theRMSD in polymorph free energy differencebetween 1D-anisotropic and anisotropicQHA is 0.0066 kcal/mol, which is well belowphysically meaningful sensitivity thresh-olds. Figure 5a shows the polymorph freeenergy differences computed with 1D-QHAversus those computed with anisotropic12 G ( T max ) [kcal/mol]Aniso. QHA G ( T m a x ) [ k c a l / m o l ] D - A n i s o . Q H A RMSD = 0.0066[kcal/mol] a) Percent ExpansionAniso. QHA P e r c e n t E x p a n s i o n D - A n i s o . Q H A RMSD = 0.29% b) Figure 5: Comparison of the results of 1D- and fully-anisotropic QHA applied to the polymor-phic unit cells. In all plots, the filled-in markers are the molecules with multiple molecularconformers (tolbutamide, chlorpropamide, and aripiprazole). a) Polymorph free energy dif-ferences (∆ G ( T max )) computed with 1D and full anisotropic approaches are nearly equal.The black line is plotted as y = x and the data has a RMSD = 0.0066 kcal/mol relative tothat line. b) The percent expansion from 0 K to T max for the three lattice vectors and anglesare plotted. The RMSD in percent expansion between 1D and anisotropic QHA is 0.29%.13HA at T max . The free energies of mostpolymorphic pairs vary by less than 0.010kcal/mol between 1D-QHA and anisotropicQHA. The only exception is aripiprazole,where the two QHA methods for determin-ing ∆ G between polymorphs differ by 0.016kcal/mol. The error in aripiprazole comesprimarily from form X. The error in theindividual free energies of polymorphs be-tween 1D and fully anisotropic QHA is < .
006 kcal/mol (excluding aripiprazole formX), whereas the error in isotropic expansionis between 0.01–0.21 kcal/mol (Figure S11).Similarly, the high temperature lattice ge-ometries of the two approaches have anRMSD of 0.29% for the expansion relativeto lattice minimum structure.For the molecules studied, 1D-anisotropic QHA is thus a sufficiently ac-curate approach to model anisotropic ex-pansion for most purposes, and is muchmore efficient. The 1D approach producespolymorph free energy differences within0.02 kcal/mol of the fully anisotropic ap-proach at about 10% of the computationalcost. More importantly, the 1D approachhas greater numerical stability than fullanisotropic expansion, allowing us obtainresults at higher temperatures. The closenumerical agreement between the 1D andfully anisotropic approach up to T max showsthat the 1D approach is a realistic con-straint for these molecules. We will there-fore use the 1D approach for all furtherevaluations for anisotropic QHA. The Cell Size Has a Non-negligible Effect on PolymorphStability
We observe differences in the quasi-harmonic free energy differences betweenunit cells and supercells ranging from 0.04–0.5 kcal/mol and are loosely correlated to differences in the high temperature cell ge-ometry. In our previous study we assumedthat an energy contribution less than 0.12kcal/mol would give us 90% confidence thata re-ranking in crystal stability would notoccur, a number based on a previousmuch larger study of lattice energies versusquasi-harmonic free energies. Six of theten polymorphic pairs in Figure 6 exceedthat cutoff, showing that finite size errorsare sufficiently large to affect the stabilityranking of polymorphs. All unit cells werequenched from the experimental crystalsprior to running QHA and the supercellswere constructed from the quenched unitcells. These results used the 1D-anisotropicapproach for QHA.14 .00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25
RMSD (% h ) | G un i t G s u p e r | ( T m a x ) [ k c a l / m o l ] R = 0.79 Figure 6: Differences between 1D-QHA inunit and supercells are shown for polymorphfree energy differences against the averageRMSD in percent expansion at T max . Thereis a moderate ( R = 0 .
79) correlation be-tween the difference in unit and supercellpercent expansion and the finite size errorin free energy. The grey line is the 0.12kcal/mol threshold, where values below theline have a 90% estimated confidence thatre-ranking in crystal stability would not oc-cur. In all plots, the filled-in markers arethe molecules with multiple molecular con-formers (tolbutamide, chlorpropamide, andaripiprazole).The error in the free energy differencesbetween unit and supercells is correlated tothe difference in expansion of the lattice pa-rameters between unit and supercells. InFigure 6 we plot the deviations in the poly-morph free energy differences versus the av-erage RMSD of the percent expansions inboth crystals. We found that there wasa moderate linear correlation ( R = 0 . R = 0 .
109 in Fig-ure S13). The cell dimensions of supercellsrelative to unit cells can be found in TableS3.
Peforming QHA with a LatticeMinima Consistent with MDEnsembles
The lattice minimum structure quenchedfrom experiment does not always coincidewith the low temperature structures foundfrom MD simulation of the same model,raising the question of what it means tocompare QHA to MD. For some of thepolymorphs studied here, we have foundthat the crystal supercells clearly reorga-nized into different configurational ensem-bles when heated up to 300 K using MD andcan be quenched to a number of lattice min-imum.
This suggests than QHA can beperformed either independent or dependentof MD equilibration. QHA independent ofMD equilibration, using a minimized exper-imental structure, has the potential risk ofoccupying a different basin than the MDsimulation does, leading to potentially largedifferences in both lattice energy and config-urational entropy between the approaches.In contrast, if QHA is run from a thermallyannealed structure of the MD simulations,we have a better chance of matching thehigh-temperature MD ensemble with QHA,15s the lattice energy and vibrational modeswill be more nearly the same, giving a purertest of whether the high temperature behav-ior is indeed quasi-harmonic.To maximize the chance that we canidentify the free energy minimum at lowtemperature, we performed temperaturereplica exchange molecular dynamics forthe polymorphs of tolbutamide, chlor-propamide, and aripiprazole. For thesesix polymorphs, five frames from the equi-librated 10 K REMD simulation werequenched in order to find a lattice mini-mum for QHA consistent with the MD en-semble. All frames from the MD trajec-tory quenched to minima that had theirpotential energy within 0.02 kcal/mol, lat-tice vectors within 0.2 ˚A, and angles within1.0 ° from one another. Differences in min-ima arose from conformations of the alkyltails in both polymorphs of tolbutamide andchlorpropamide form V. In Table 1 we re-port the potential energy of the polymorphswhen quenched from experiment, the lat-tice minimum used in our 2017 paper, andthe lowest energy structure quenched fromREMD, relative to minimum found usingREMD. Quenching from REMD producesthe lowest lattice energy structure for allpolymorphs. Overlaid structures of restruc-tured crystals are provided in Figures S14–S19.When the lowest lattice energy mini-mum is chosen, the error between isotropicQHA and MD is reduced by 0.083 kcal/molfrom the error previously reported. Theerror previously reported for these three sys-tems was 0.083 – 0.397 kcal/mol, but withthe new minima is 0.014 – 0.247 kcal/mol.In Table 2 the error of QHA for polymorphfree energy differences relative to MD us-ing replica exchange is shown at T max for the minimum found from our previous studyand the minimum quenched from the 10 KREMD simulation. We exclude the latticeminimum quenched from the experimentalstructure since the lattice energy differencesare so large. For isotropic QHA, the min-imum quenched from REMD reduces theerror relative to MD by 0.15 kcal/mol fortolbutamide and chlorpropamide. We seean increase in the error for aripiprazole,which is minimal compared to the reduc-tion in error for the other two sets of poly-morphs.When using the lower lattice energyminimum, anisotropic QHA determines thestability of chlorpropamide and aripipra-zole within 0.05 kcal/mol relative to MD.The error in QHA for tolbutamide, 0.206kcal/mol, is still large, but is half the mag-nitude of what was previously reported. Despite the new minimum of aripiprazoleproducing higher error in isotropic QHA,using anisotropic expansion reduced the er-ror to 0.045 kcal/mol, implying that anisotropic constraint is inappropriate. Thesmall increase in error when switchingto anisotropic treatments seen in chlor-propamide for anisotropic QHA is mostlikely due to cancellation of error betweenpolymorphs. For these polymorphs, all fur-ther analysis of QHA will be shown for therestructured crystals found from quenchinglow-temperatures replicas of REMD simula-tions.
Comparison of thermal expan-sion and thermodynamics usingQHA and MD
Anisotropic QHA allows the crystal lat-tice to relax into geometries that are moresimilar to MD, but there are still persis-tent errors in the high temperature geome-16
Percent Expansion ofMD Vectors P e r c e n t E x p a n s i o n o f Q H A V e c t o r s a) RMSD=1.6%RMSD=1.3%
Iso. QHAAniso. QHA
Percent Expansion ofMD Angles P e r c e n t E x p a n s i o n o f Q H A A n g l e s b) RMSD=1.3%RMSD=1.3%
Iso. QHAAniso. QHA
Figure 7: Scatter plots of the percent expansion of QHA versus MD for the lattice a) vectorsand b) angles for isotropic and 1D-anisotropic QHA lattice parameters. The RMSD to MDis marginally improved for the lattice vectors and not improved for the angles when usinganisotropic expansion over isotropic expansion. A least square fit between MD anisotropicQHA shows that anisotropic expansion provides a better y = x fit (black line) to MD thanisotropic expansion. The slopes for anisotropic and isotropic expansion are a) 0.89 ± ± ± U Exp . a U Previous b U REMD a Potential energy differences of the lattice minimum quenched from the experimentalstructure and quenched from a 10 K REMD simulation; b Potential energy difference of thelattice minimum quenched from the restructured crystal found in our previous study andquenched from a 10 K REMD simulation.Table 2: Error of QHA Relative to MD at T max Previously Reported and Using a LatticeMinimum Quenched from REMDTolbutamide Chlorpropamide AripiprazoleQHA Iso. Aniso. Iso. Aniso. Iso. Aniso. δ (∆ G Prev . ) c δ (∆ G REMD ) d c Error between QHA and MD using the lattice minimum in our previous study; d Errorbetween QHA and MD using the lattice minimum quenched from the 10 K replica ofREMD.tries computed with QHA. In Figure 7aand Figure 7b we show the percent ex-pansion of QHA lattice vectors and anglesagainst MD, respectively. The black linerepresents a perfect agreement between theQHA method and MD, while the dashedlines are least squared fits between the twoQHA methods and MD. For the lattice vec-tors, there is marginal improvement in theRMSD when using anisotropic expansion(1.3%) over isotropic QHA (1.6%). Thereis no improvement to RMSD for the latticeangles (anisotropic 1.3%, isotropic 1.3%).Despite there being little-to-no reductionin the RMSD’s, the least square fit be-tween anisotropic QHA and MD is similarto the y = x line for both lattice vectorsand angles. This contrasts with isotropicQHA, which is essentially independent ofMD box vectors, which is expected for thefixed isotropic angles. Quantitatively, the slopes for anisotropic and isotropic expan-sion are 0.89 ± ± a , b , and c ) and 0.60 ± α , β and γ ), indicating that while anisotropic expan-sion may over or undershoot the changes,the directions of change are generally beingproperly predicted.Anisotropic expansion provides minimalimprovement in the error of polymorph freeenergy differences between QHA and MD.In Figure 8a we plot the Gibbs free energydifferences of both QHA methods againstthe results for MD at T max . The RMSD of∆ G for isotropic and anisotropic QHA rel-ative to MD are 0.113 and 0.079 kcal/mol,respectively. There is only a 0.034 kcal/molimprovement to the RMSD when usingan anisotropic thermal expansion model,which is comparable to the bootstrappederror in the RMSD. Anisotropic QHA in18 G ( T max ) MD [kcal/mol] G ( T m a x ) Q H A [ k c a l / m o l ] a) RMSD iso . = 0.113±0.026 RMSD aniso . = 0.079±0.025Iso. QHAAniso. QHA ( H QHA ) [kcal/mol] ( T S Q H A ) [ k c a l / m o l ] b) Iso. QHAAniso. QHA
Figure 8: a) Gibbs free energy differences between polymorphs for QHA relative to MD andb) entropy–enthalpy compensation in errors of QHA relative of MD. a) Overall, there is aslight improvement of RMSD in the Gibbs free energy when using anisotropic expansionover isotropic expansion. This improvement is only slightly larger than the bootstrappederror in the RMSD. b) Grey shaded area indicates free energy differences < .
12 kcal/molbetween QHA and MD, which for most of the graph indicates significant entropy–enthalpycompensation. The only molecule to fall outside of these bounds for anisotropic QHA istolbutamide, which has the largest error in QHA relative to MD. In all plots, the filled-inmarkers are the molecules with multiple molecular conformers (tolbutamide, chlorpropamide,and aripiprazole). 19act has a larger deviation from MD thanisotropic QHA for 3 of the 10 molecules,though the deviation from MD is small.Those molecules and their deviation ofanisotropic QHA from MD ( δ (∆ G D ) MD )are pyrazinamide (0.037 kcal/mol), resorci-nol (0.062 kcal/mol), and chlorpropamide(0.033 kcal/mol). These results demon-strate that the deviations in the free energydifferences between MD and QHA are dueto something other than an insufficiently ac-curate thermal expansion model. Summaryresults for each molecule is provided in theFigure S57.The difference between QHA and MDis minimized when a lattice minimum cor-responding the sampled MD lattice pa-rameters and anisotropic QHA is used.The largest errors for anisotropic QHA isseen for tolbutamide ( δ (∆ G D ) MD = 0.206kcal/mol), paracetamol (0.088 kcal/mol),and resorcinol (0.062 kcal/mol) with allother molecules < However, molec-ular size and flexibility are correlated witha greater number of alternative lattice min-ima. We also see no correlation betweenthe error in QHA and the number of hydro-gen bonds in the crystal ( R < . (cid:48) of 4 and all other polymorphs are ≤
2, while this is not statistically significantit could explain why the error in QHA is anoutlier. We also compared the gas and crys-tal phase torsions, but found no differencesthat distinguish tolbutamide from the restof the molecules (Torsions shown in figuresS60–S66).We see entropy–enthalpy compensationin the differences between anisotropic QHAand MD of up to 0.12 kcal/mol for allmolecules except for tolbutamide. Despite the entropy and enthalpy for QHA havingerrors up to 0.5 kcal/mol from MD, the can-cellation of those errors allows QHA to givefree energy differences within 0.09 kcal/molof the MD-calculated free energy differencesfor all systems except tolbutamide. In Fig-ure 8b the error in the entropic differences ofpolymorphs is plotted versus the error in en-thalpic differences of QHA relative to MD.We define a set of polymorphs demonstrat-ing entropy–enthalpy compensation (grayshaded area) if the differences between er-rors in entropy and enthalpy (i.e. the er-ror in the free energy difference) is < . Anharmonic mo-tions may be important for computing theenthalpy and entropy differences of poly-morphs, but for most of the systems studiedhere the error is minimal for computing thefree energy. The only molecule that fallsoutside of this region for anisotropic QHAis tolbutamide. Out of the 10 molecules, 8have larger entropic error in QHA than en-thalpic error, which is most likely due toanharmonic motions. The only two setsof polymorphs that have larger enthalpicerror in QHA are tolbutamide and chlor-propamide, which are further explained inthe next section.
Molecular Conformational En-sembles Can Still Exist LowTemperatures
Even at low temperatures, some crystalshave substantial conformational degenera-cies allowing quenched structures to settleinto lattice minima that break the symme-try of the supercell. In figure 9b we showframes of the trajectory of chlorpropamide,with each frame taken 0.1 ns apart for thelast 5 ns of the 10 K ensemble produced20
50 100 150 200 250 300
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD U x Lattice Minimum10K Simulation U x a) b) c)d) e) f) Figure 9: Snapshots of the 10 K ensemble for chlorpropamide form a) I and b) V taken fromthe last 5 ns of the replica exchange simulation, with 0.1 ns between each frame. Even at lowtemperatures, this crystal structure samples a surprisingly diverse conformational ensemble,with the greatest variability in the fluctuations around the alkyl tail. The RMSD of theentire 10 K simulation ensemble relative to the lattice minimum for form I is 0.165 ± ± U ) against the configurationalspace ( x ). f) Gibbs free energy difference of chlorpropamide form V relative to form I.21uring replica exchange. From this lowtemperature ensemble, molecules are ableto sample a number of configurations dueto the degrees of freedom accessible to thealkyl tail, and the low energy differencesbetween different alkyl rotamer configura-tions. When quenched, the tails fall intomostly independent conformations, break-ing the crystal symmetry.The thermodynamic accessibility ofmany configurations at low temperatures isa result of the existence of many configu-rationally diverse lattice minima with sim-ilar low energies. In Figure 9b the RMSFof form V is 0.69 ± ± G ( T ) plot from MD, between0 and 10, the polymorph free energy differ-ence for MD quickly changes, causing QHAto diverge by 0.06 kcal/mol at 10 K. WhileQHA can accurately model form I at lowtemperatures, it cannot capture the behav-ior of the alkyl tails in form V, which iswhy we see the immediate divergence of thetwo methods as T increases from zero. Attemperatures greater than 10 K the devia-tion of the free energy from MD does notgrow because of the entropy–enthalpy com-pensation, but we cannot necessarily con- clude that this would be the case for othercrystals that could have alternate configura-tional minimum. Both polymorphs of tolbu-tamide also sample a large number of con-formations at low temperatures, causing asimilar divergent behavior in ∆ G for QHAand MD at low temperatures (Figure S48).The concavity in the free energy of chlor-propamide at low temperatures is due todominating enthalpic changes close to 0 K,which quickly diminish relative to the en-tropic contributions as the temperature in-creases. Surprisingly, although form V hasa larger low T configurational ensemble, itinitially becomes less stable as a functionof T relative to the more constrained formI. Up to 10 K, form V expands faster thanform I (Figure 10a) causing an unfavorableenthalpic contribution to form V relative toform I (Figure 10b). Since T ∆ S is smallat low temperatures, the changes in the en-thalpic difference dominates the free energydifferences. At temperatures greater than10 K the difference in volume and enthalpybetween polymorphs is constant. The initialexpansion of form V allows the alkyl tails tohave more conformational space than formI and by 20 K the entropic contributions inform V are large enough to switch the sta-bility of the polymorphs. Conclusions
Although we have studied a relatively smallnumber of enantiotropic polymorphs in thisstudy, we find there here are persistentshortcomings of QHA that highlight theimportance of anharmonic configurationalsampling in determining accurate free en-ergy differences of crystal polymorphs. Theshortcomings identified in this study are:1) that including anisotropic expansion pro-vides minimal improvement in the error ofpolymorph free energy differences between22
20 40 60 80 100
T [K] V [ A n g . / m o l ] a) T [K] C o n t r i b u t i o n t o G [ k c a l / m o l ] b) T SU ( T ) - U ( T = 0 K ) Figure 10: (a) Volume differences of chlorpropamide form V relative to form I using MD.Initially, form V expands much faster than form I. Above 10 K, the difference in volumebetween the forms is almost constant. (b) The entropic and enthalpic contributions to theGibbs free energy difference of chlorpropamide form V relative to I. Below 10 K, the potentialenergy is the dominating free energy contribution making form V less favorable, explainedby the sudden expansion of form V seen at low temperatures.QHA and MD, 2) the number of accessiblelow-energy minima causes a rapid deviationof QHA relative to MD for temperatures aslow as 10 K, 3) some crystals reorganize intoclearly different structures when heated upto 300 K using MD, and 4) polymorph re-structuring prevents QHA from being runindependently of MD because the minimumenergy configuration is likely to be inacces-sible from an arbitrary low-energy startingpoint.We found that the OPLS-AA pointcharge potential incorrectly ranked 4 of the10 enthalpic differences of polymorphs rel-ative to experiment, but correctly rankedthe entropic differences of all 10 molecules.Despite the lattice geometries for MD dif-fering moderately from experiment, the ex-pansion of vectors and orthogonal angles forMD generally agree with experiment. TheRMSD of both enthalpic and entropic dif-ferences between MD and experiment were4.0 and 2.0 kcal/mol, respectively, which isroughly 2 orders of magnitude larger than errors seen of QHA relative to MD if prop-erly performed.Our one-dimensional approach to modelanisotropic thermal expansion predicts es-sentially identical high temperature latticegeometries and polymorph free energy dif-ferences as the full anisotropic QHA for the10 polymorph pairs studied. Specifically,the computed polymorph free energy differ-ences are within 0.02 kcal/mol at T max forour 1D- and fully-anisotropic QHA meth-ods for the 10 molecules studied. The lat-tice geometries of the two approaches re-mained essentially the same, with an RMSDof 0.29% for the expansion relative to lat-tice minimum structure. Additionally, 1D-anisotropic QHA is able to find a free energyminimum up to 300 K for all but one crystalunit cell, whereas the fully anisotropic ap-proach fails to reach those desired temper-atures for 50% of the crystals, due to thelarger likelihood of a rearrangement into analternate minimum when multiple deriva-tives are used, resulting in numerical insta-23ility of the approach.Finite size effects do not always can-cel out for polymorph free energy differ-ences; instead, differences between unit cellsand supercells can be up to 0.5 kcal/moland therefore must be accounted for whenattempting to accurately rank crystal sta-bility. We found that this error is corre-lated with the error in the lattice geometries( R = 0.79), meaning they are reduced butnot eliminated if the high temperature lat-tice geometries are the same between unitand supercells.Anisotropic QHA only marginally im-proves the error in free energy differencesfrom MD for the polymorph pairs studied.Using anisotropic QHA reduces the RMSDof ∆ G relative to MD by 0.034 kcal/mol,which is only slightly larger than the stan-dard error in the RMSD. The errors inQHA with respect to MD are due primarilyto anharmonic motions, because using ananisotropic thermal expansion model doesnot improve the deviations between QHAand MD. Differences between anisotropicQHA and MD do not appear strongly cor-related with molecular flexibility. Tolbu-tamide with both a large deviation andhigh flexibility is an exception, and thisdifference could also be due to the largeZ (cid:48) value of 4 for form II. The free energydifference in tolbutamide had an error of0.21 kcal/mol between anisotropic QHA andMD, whereas all other molecules deviated < .
09 kcal/mol.The degree of flexibility of the moleculescan result in a more significant chanceof crystal restructuring upon annealingwith replica exchange. Once a lower en-ergy lattice minimum was found, the er-ror in isotropic QHA was on average 0.08kcal/mol lower than previously reported. Anisotropic QHA from the new minimumwas able to reduce the free energy differ-ences between chlorpropamide and aripipra- zole to < .
05 kcal/mol. We were consis-tently able to find lower lattice minimumthan those quenched directly from exper-iment for more flexible molecules by run-ning replica exchange molecular dynamicsand quenching frames from the low temper-ature replica.The polymorphs of conformationallyflexible molecules exhibit differences in ∆ G of QHA relative to MD on the magnitude of0.06 kcal/mol even as low as 10 K, despitethe 0 K lattice energies being within 0.01kcal/mol at 0 K, showing that agreementbetween QHA and MD cannot be assumedeven at very low temperatures. These dif-ferences are due to either one or both of thepolymorphs expanding quickly at low tem-peratures, which the provides enough freevolume to access a diverse number of confor-mations around the alkyl tail. We note thatthe temperature at which a crystal can ac-cess these conformations is structure depen-dent, as seen with the different polymorphsof chlorpropamide. Acknowledgments
The authors thank Eric Dybeck for pro-ducing majority of the MD results usedfrom our previous paper. Final resultsfor QHA were performed on the ExtremeScience and Engineering Discovery Envi-ronment (XSEDE), which is supported byNational Science Foundation grant num-ber ACI-1548562. Specifically, it used theBridges system, which is supported by NSFaward number ACI-1445606, at the Pitts-burgh Supercomputing Center (PSC). Thiswork was also supported financially by NSFthrough the grant CBET-1351635.24 upporting Information for “Adding anisotropy to the standardquasi-harmonic approximation still fails in several ways to captureorganic crystal thermodynamics”
Systems Studied
In Table 3 we provide the systems stud-ied, their polymorphs names/numberings,and corresponding Cambridge StructureDatabase reference code (CSD Refcode),and the dimensions of the unit and super- cells relative to the symmetric cell pulledfrom the database. The supercell sizes arethe same as our previous study and the unitcells are the smaller cell size that assuresthat both crystals have the same number ofmolecules.
Piracetam
Form I Form IIICSD Refcode BISMEV03 BISMEV02Space Group P2 /n P2 /nZ / Z (cid:48) × × × × × × × × Carbamazepine
Form I Form IIICSD Refcode CBMZPN11 CBMZPN02Space Group P¯1 P2 /nZ / Z (cid:48) × × × × × × × × Adenine
Form I Form IICSD Refcode KOBFUD KOBFUD01Space Group P2 /c Fdd2Z / Z (cid:48) × × × × × × × × Pyrazinamide
Form γ Form δ CSD Refcode PYRZIN20 PYRZIN16Space Group Pc P¯1Z / Z (cid:48) × × × × × × × × Resorcinol
Form α Form β CSD Refcode RESORA03 RESORA08Space Group Pna2 Pna2 Z / Z (cid:48) × × × × × × × × Form α Form β CSD Refcode ZZZPRO03 ZZZPRO04Space Group Pbca PccnZ / Z (cid:48) × × × × × × × × Paracetamol
Form I Form IICSD Refcode HXACAN01 HXACANSpace Group P2 /a PcabZ / Z (cid:48) × × × × × × × × Aripiprazole
Form I Form XCSD Refcode MELFIT01 MELFIT05Space Group P2 P2 Z / Z (cid:48) × × × × × × × × Tolbutamide
Form I Form IICSD Refcode ZZZPUS04 ZZZPUS05Space Group Pna2 PcZ / Z (cid:48) × × × × × × × × Chlorpropamide
Form I Form VCSD Refcode BEDMIG BEDMIG04Space Group P2 Pna2 Z / Z (cid:48) × × × × × × × ×
1D vs. Full AnisotropicQHA
In Table 4 we report the free energy differ-ences of the unit cell of each system using1D and fully anisotropic QHA at the max-imum temperature that all methods could achieve T max . We also report the deviationbetween the two methods ( δ ( | ∆ G ( T max ) | ))in the far right column, which shows thatboth methods are within 0.02 kcal/mol ofeach other for all systems except for chlor-propamide.Differences in ∆ G of the anisotropicmethods are correlated with the errors in26able 4: Free energy differences for 1D and fully anisotropic QHA at T max for the unit cells. T max is the maximum temperature that both crystals and both methods can reach while stillbeing at a free energy minimum. For all systems there is < G with aripiprazole having the largest difference of 0.016 kcal/mol.∆ G ( T max ) [kcal/mol] δ ( | ∆ G ( T max ) | )System T max [K] 1D-QHA Aniso.-QHAPiracetam 255 -1.374 -1.384 0.010Carbamazepine 300 -0.160 -0.157 0.003Adenine 300 1.330 1.328 0.002Pyrazinamide 180 -0.552 -0.551 0.001Resorcinol 300 -0.578 -0.578 0.0001,4-Diiodobenzene 210 0.431 0.431 0.000Paracetamol 150 -0.648 -0.652 0.004Aripiprazole 180 -0.696 -0.712 0.016Tolbutamide 100 -3.592 -3.593 0.001Chlorpropamide 150 4.830 4.824 0.006the predicted high temperature lattice ge-ometries. Figure 11 has the absolute freeenergy deviations between 1D- and fully-anisotropic QHA versus the root-mean-square deviation (RMSD) between T max lat-tice geometries for individual polymorph.The RMSD is computed as: RM SD = (cid:118)(cid:117)(cid:117)(cid:116) (cid:88) i =1 ( C Di ( T max ) − C F ulli ( T max ))(10)where the RMSD is taken between the sixlattice parameters of the upper triangularcrystal tensor form. From Figure 5c, wecan see that as the RMSD between the1D- and fully-anisotropic T max structure be-comes larger, so does the error in the com-puted free energy of the polymorph. Theerror in thermal expansion and free energyfor isotropic QHA is reduced by using the1D-QHA variant. The differences in the freeenergy are a direct result of the crystal ex-panding to a slightly different geometry, butwe can see that this has little effect on thefree energy differences.The 1D-variant reduces the error of the isotropic QHA up to 90% for the sys-tems studied relative to the fully anisotropicapproach. The deviations between 1D-QHA and fully ansiotropic QHA are < Unit vs. Supercell 1D-QHA
In Table 5 we report the free energy differ-ences of the unit and supercell of each sys-tem using 1D-QHA at the maximum tem-perature all cell sizes could achieve T max .We also report the deviation between thetwo cell sizes ( δ ( | ∆ G ( T max ) | )) in the farright column. We see that there are sig-nificant finite size errors for all crystals thatrange between 0.03 – 0.50 kcal/mol.27 .0 0.1 0.2 0.3 0.4 0.5 RMSD ( C ( T max )) [ Å ] | G D / I s o . Q H A G A n i s o . Q H A | ( T m a x ) [ k c a l / m o l ] Figure 11: Deviations in the free energy of the constrained QHA (isotropic QHA and 1D-QHA) from anisotropic QHA versus the box vector RMSD from the unconstrained QHAstructure. All values are reported at T max . For individual free energies the 1D-variant devi-ates < T max using 1D-QHA are reported. T max is the maximum temperature that both crystals and cell sizes can reach while still beingat a free energy minimum. We find that size errors on the magnitude of 0.03 – 0.5 kcal/molfor polymorph free energy differences. ∆ G ( T max ) [kcal/mol] δ ( | ∆ G ( T max ) | )System T max [K] Unit SuperPiracetam 300 -1.364 -1.457 0.093Carbamazepine 300 -0.160 0.044 0.204Adenine 300 1.330 1.603 0.273Pyrazinamide 300 -0.756 -0.613 0.143Resorcinol 300 -0.578 -0.439 0.1391,4-Diiodobenzene 300 0.427 0.389 0.038Paracetamol 300 -1.241 -0.767 0.474Aripiprazole 300 -0.764 -0.950 0.186Tolbutamide 180 -3.602 -3.489 0.113Chlorpropamide 50 5.374 5.287 0.08728or individual polymorphs, there is notrend between the differences in free energyand the RMSD in the percent expansion upto T max for unit and supercells; the trendsonly appear for polymorph differences. Fig-ure 12 compares the absolute deviation infree energy with the RMSD in the percentexpansion for unit and supercells. The er-ror in free energy between unit and supercells is between 0 – 3 kcal/mol and does notcorrelate with the cells percent expansion at T max . RMSD (% h ) | G un i t G s u p e r | ( T m a x ) [ k c a l / m o l ] Figure 12: The deviation in the free energybetween unit and super cell is shown relativeto the RMSD in the percent expansion. Forindividual polymorphs there is no correla-tion between error in the free energy and theerror in the high temperature lattice geome-tries. The filled in markers are the systemswith multiple molecular conformers (tolbu-tamide, chlorpropamide, and aripiprazole)The finite size errors do not correlatewith the ratio of number of molecules inthe unit and supercells. We have supercells that have 6 – 24 times the numberof molecules in the unit cells. We wouldexpect the systems with similar number of molecules in the unit cell and supercell tohave smaller deviations in the free energy,but that is not the case. In Figure 13 thefree energy deviation between the unit andsupercell at 300 K is plotted against the ra-tio of molecules in the super cell relative tothe unit cell. The greatest error is seen forparacetamol, which is one of the systemswith the most similar unit and super cells.There is no correlation between the finitesize error and ratio of number of molecules( R = 0 . N supercell N unitcell | G un i t G s u p e r | [ k c a l / m o l ] R = 0.109 Figure 13: Plotting the finite size errors at T max versus the ratio of number of moleculesbetween the super- and unit cell.29 e-structuring Figure 14: Lattice minimum of chlor-propamide form I from quenching the ex-perimental structure (green), the struc-ture used in our previous study (blue), quenched from a 10 K replica exchange sim-ulation (red). The RMSD between all struc-tures is < quenched from a 10 K replica exchangesimulation (red). The RMSD between thestructure quenched from experiment andfrom our previous paper relative to the min-ima found from REMD are 1.26 and 1.07 ˚A,respectively.30igure 16: Lattice minimum of aripipra-zole form I from quenching the experimen-tal structure (green), the structure used inour previous study (blue), quenched froma 10 K replica exchange simulation (red).The RMSD between the structure from ourprevious paper relative to the minima foundfrom REMD is 0.36 ˚A.Figure 17: Lattice minimum of aripipra-zole form X from quenching the experimen-tal structure (green), the structure used inour previous study (blue), quenched froma 10 K replica exchange simulation (red).The RMSD between the structure from ourprevious paper relative to the minima foundfrom REMD is 0.51 ˚A. Figure 18: Lattice minimum of tolbutamideform I from quenching the experimentalstructure (green), the structure used in ourprevious study (blue), quenched from a10 K replica exchange simulation (red).The RMSD between the structure quenchedfrom experiment and from our previouspaper relative to the minima found fromREMD are 1.48 and 1.09 ˚A, respectively.Figure 19: Lattice minimum of tolbutamideform I from quenching the experimentalstructure (green), the structure used in ourprevious study (blue), quenched from a10 K replica exchange simulation (red).The RMSD between the structure quenchedfrom experiment and from our previouspaper relative to the minima found fromREMD are 1.35 and 1.35 ˚A, respectively.31 ree Energy Differencesand Expansion with Tem-perature The original MD data for several systemswere run with constraints that preventedthe angles from changing. While there fluc-tuations in the angles should have been con-sidered, the angles remain 90 ° across the en- tire temperature range. For the followingsystems the MD angles were fixed at 90 ° :adenine from II, resorcinol from α and β ,diiodobenzene form α and β -restructured,and paracetamol form II. For all crystalswhere the angles are fixed at 90 ° , we ran1 ns simulations at 50, 100, 200, and 300K with the angles unfixed and verified thataverage value of the angle remains at 90 ° . Piracetam
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 20: Polymorph free energy differences of piracetam form I relative to form III32
100 200 3006.7256.7506.7756.8006.8256.8506.8756.900 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit) Fabbiani 2005 Louer 1995 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 21: Piracetam form I percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Admiraal 1982 Chambrier 2011Galdecki 1983 Shipilova 2018Tilborg 2011 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 22: Piracetam form III percent expansion from 0 with experimental results. arbamazepine T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 23: Polymorph free energy differences of carbamazepine form I relative to form III34
100 200 3005.005.055.105.15 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD 1D-QHA (unit) Aniso. QHA (unit) Fernandes 2007 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 24: Carbamazepine form I percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHAMD 1D-QHA (unit)Aniso. QHA (unit) Brandenburg 2017Hassan 2013 Horstman 2015Lisgarten 1989 Nievergelt 2017Sovago 2016 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 25: Carbamazepine form III percent expansion from 0 K with experimental re-sults. denine T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 26: Polymorph free energy differences of adenine form II relative to form I36
100 200 3007.857.907.958.008.05 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD 1D-QHA (unit) Aniso. QHA (unit) Mahapatra 2008 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 27: Adenine from I percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD 1D-QHA (unit) Aniso. QHA (unit) Stolar 2016 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 28: Adenine form II percent expansion from 0 K with experimental results. yrazinamide T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 29: Polymorph free energy differences of pyrazinamide form γ relative to form δ
100 200 3007.157.207.257.30 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Castro 2010 Cherukuvada 2010Nakata 1987 Tamura 1961 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 30: Pyrazinamide form γ percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit) Nangia 2005 Ro 1972 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 31: Pyrazinamide form δ percent expansion from 0 K with experimental results. esorcinol T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 32: Polymorph free energy differences of carbamazepine form β relative to form α
100 200 30010.3010.3510.4010.4510.5010.55 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Bacon 1979 Bacon 1980Dru z bicki 2015 Fronczek 2001Robertson 1934 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 33: Resorcinol form α percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Bacon 1980 Dru z bicki 2015 Robertson 1938 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 34: Resorcinol form β percent expansion from 0 K with experimental results. ,4-Diiodobenzene For MD, form β was found to restructureat high temperature to form β -restructured.We show the results for form β usingQHA only, where the lattice minimum wasfound by directly quenching the experimen- tal crystal structure. The results for form β were used for the comparison of the 1D vari-ant and fully anisotropic QHA, as well asthe error with cell size. Form β -restructuredwas used for the comparison to MD. T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD
Figure 35: Polymorph free energy differences of diiodobenzene form β -restructured relativeto form α
100 200 30016.216.416.616.817.0 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Alcobe 1994 Hendricks 1933 Hinchliffe 1985 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 36: Diiodobenzene form α percent expansion from 0 K with experimental re-sults. L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 37: Diiodobenzene form β percent expansion from 0 K with experimental results.MD is not shown for the original structure of β because the system restructured at hightemperature. 43
100 200 30016.116.216.316.4 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 38: Diiodobenzene form β -restructured percent expansion from 0 K with experimentalresults. 44 aracetamol T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD1D-QHA (unit)Aniso. QHA (unit)
Figure 39: Polymorph free energy differences of paracetamol form I relative to form II45
100 200 30012.512.612.712.812.913.0 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Haisa 1976 Shankland 1997Wilson 1997 Wilson 2009 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 40: Paracetamol form I percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MD1D-QHA (unit) Aniso. QHA (unit)Drebushchak 2009 Haisa 1974 Nichols 1998 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 41: Paracetamol form II percent expansion from 0 K with experimental results. ripiprazole For MD, both crystal structures werefound to restructure at high temperature toform. We show the results for for I and Xusing QHA only, where the lattice minimumwas found by directly quenching the experi-mental crystal structure. These results were used for the comparison of the 1D variantand fully anisotropic QHA, as well as theerror with cell size. Form I-restructured andform X-restructured were used for the com-parison to MD. L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 42: Aripiprazole form I percent expansion from 0 K. MD is not shown for the originalstructure of I because the system restructured at high temperature.47
100 200 3008.748.768.788.808.828.848.86 L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 43: Aripiprazole form X percent expansion from 0 K. MD is not shown for the originalstructure of X because the system restructured at high temperature.
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD
Figure 44: Polymorph free energy differences of aripiprazole form X-restructured relative toform I-restructured 48
100 200 3007.67.88.08.28.48.68.8 L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD Braun 2009 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 45: Aripiprazole form X-restructured percent expansion from 0 K with experimentalresults. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD Braun 2009 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 46: Aripiprazole form I-restructured percent expansion from 0 K with experimentalresults. olbutamide For MD, form I was found to restructureat high temperature to form I-restructured.We show the results for form I usingQHA only, where the lattice minimum wasfound by directly quenching the experimen- tal crystal structure. The results for form Iwere used for the comparison of the 1D vari-ant and fully anisotropic QHA, as well asthe error with cell size. Form I-restructuredwas used for the comparison to MD. L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 47: Tolbutamide form I percent expansion from 0 K.MD is not shown for the originalstructure of I because the system restructured at high temperature.50
50 100 150 200 250 300
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD
Figure 48: Polymorph free energy differences of tolbutamide form II relative to form I-restructured. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA1D-QHA MDDonaldson 1981 Drebushchak 2016 Schreyer 2014 Thirunahari 2010 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 49: Tolbutamide form I-restructured percent expansion from 0 K with experimentalresults.
100 2009.059.109.159.20 L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 50: Tolbutamide form II percent expansion from 0 K with experimental results. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD Drebushchak 2016 Thirunahari 2010 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 51: Tolbutamide form II-restructured percent expansion from 0 K with experimentalresults. hlorpropamide For MD, form V was found to re-structure at high temperature to form V-restructured. We show the results for formV using QHA only, where the lattice min-imum was found by directly quenching theexperimental crystal structure. The results for form V were used for the comparison ofthe 1D variant and fully anisotropic QHA,as well as the error with cell size. Form V-restructured was used for the comparison toMD.
T [K] G [ k c a l / m o l ] Iso. QHA1D-QHAMD
Figure 52: Polymorph free energy differences of chlorpropamide form V-restructured relativeto form I-Restructured 53
100 200 30026.626.827.027.227.4 L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 53: Chlorpropamide form I percent expansion from 0 K with experimental results.
MD is not shown here because this is QHA using the minima directly quenched from theexperimental structure, which does not correspond with the 0 K minima found with MD. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD Drebushchak 2009 Koo 1980 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 54: Chlorpropamide form I-Restructured percent expansion from 0 K with experi-mental results.
100 200 30019.219.419.619.820.0 L a tt i c e V e c t o r [ Å ] a-Vector b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 55: Chlorpropamide form V percent expansion from 0 K. MD is not shown herebecause this is QHA using the minima directly quenched from the experimental structure,which does not correspond with the 0 K minima found with MD. L a tt i c e V e c t o r [ Å ] a-Vector Iso. QHA 1D-QHA MD Drebushchak 2008 Drebushchak 2009 b-Vector c-Vector A n g l e [ D e g . ] -Angle -Angle -Angle Temperature [K] V o l u m e [ Å ] Temperature [K]
Figure 56: Chlorpropamide form V-restructured percent expansion from 0 K with experi-mental results. ummary of Results P i r a c e t a m C a r b a m a z e p i n e A d e n i n e P y r a z i n a m i d e R e s o r c i n o l D ii o d o b e n z e n e P a r a c e t a m o l A r i p i p r a z o l e T o l b u t a m i d e C h l o r p r o p a m i d e ( G ) M D [ k c a l / m o l ] Iso. QHA Aniso. QHA
Figure 57: Deviations in the free energy differences from MD at 300 K (except for tolbutamideat 200 K, and chlorpropamide at 250 K). We report results for isotropic QHA and 1D-anisotropic QHA relative to MD. Overall, there is no evidence that a more complex thermalexpansion model reduces the error between QHA and MD. P i r a c e t a m C a r b a m a z e p i n e A d e n i n e P y r a z i n a m i d e R e s o r c i n o l D ii o d o b e n z e n e P a r a c e t a m o l A r i p i p r a z o l e T o l b u t a m i d e C h l o r p r o p a m i d e ( T S ) M D [ k c a l / m o l ] Iso. QHA Aniso. QHA
Figure 58: Deviations in the entropic contribution to the Gibbs free energy differences fromMD at 300 K (except for tolbutamide at 200 K, and chlorpropamide at 250 K). We reportresults for isotropic QHA and 1D-anisotropic QHA relative to MD. The entropy differencesbetween methods are compensated by a larger enthalpy difference, resulting in a smaller freeenergy differences. 56 i r a c e t a m C a r b a m a z e p i n e A d e n i n e P y r a z i n a m i d e R e s o r c i n o l D ii o d o b e n z e n e P a r a c e t a m o l A r i p i p r a z o l e T o l b u t a m i d e C h l o r p r o p a m i d e ( H ) M D [ k c a l / m o l ] Iso. QHA Aniso. QHA
Figure 59: Deviations in the enthalpic contribution to the Gibbs free energy differences fromMD at 300 K (except for tolbutamide at 200 K, and chlorpropamide at 250 K). We reportresults for isotropic QHA and 1D-anisotropic QHA relative to MD. The enthalpy differencesbetween methods are compensated by a larger enthalpy difference, resulting in a smaller freeenergy differences. 57 eplica Exchange
System Simulation Time [ns] Temperatures Sampled [K]Aripiprazole 20 10.0, 10.2, 10.4, 10.6, 10.8, 11.0, 11.2, 11.5,11.8, 12.1, 12.4, 12.7, 13.0, 13.3, 13.6, 13.9,14.2, 14.6, 15.0, 15.4, 15.8, 16.2, 16.6, 17.0,17.4, 17.9, 18.4, 18.9, 19.4, 19.9, 20.4, 21.0,21.6, 22.2, 22.8, 23.4, 24.1, 24.8, 25.5, 26.2,27.0, 27.8, 28.6, 29.4, 30.3, 31.2, 32.1, 33.1,34.1, 35.1, 36.2, 37.3, 38.4, 39.6, 40.8, 42.0,43.3, 44.6, 46.0, 47.4, 48.9, 50.4, 52.0, 53.6,55.3, 57.0, 58.8, 60.7, 62.6, 64.6, 66.6, 68.7,70.9, 73.2, 75.5, 77.9, 80.4, 81.1, 83.0, 85.7,88.5, 91.4, 94.4, 97.5, 100.7, 104.0, 107.4,110.9, 114.5, 118.2, 122.1, 126.1, 130.2, 134.5,138.9, 143.5, 148.2, 153.1, 158.1, 163.3, 168.7,174.3, 179.1, 180.0, 185.9, 188.1, 192.0, 200.0,204.9, 211.7, 218.7, 225.9, 233.4, 241.1, 249.1,257.4, 265.9, 274.7, 283.8, 293.2, 302.9, 313.0,323.4, 334.2, 345.3, 350.0Tolbutamide 20 10.0, 10.30, 10.6, 10.9, 11.2, 11.5, 11.8, 12.1,12.5, 12.9, 13.3, 13.7, 14.1, 14.5, 15.0, 15.5,16.0, 16.5, 17.0, 17.6, 18.2, 18.8, 19.4, 20.1,20.8, 21.5, 22.3, 23.1, 23.9, 24.8, 25.7, 26.6,27.6, 28.6, 29.7, 30.8, 31.9, 33.1, 34.3, 35.6,36.9, 38.3, 39.8, 41.3, 42.9, 44.6, 46.3, 48.1,50.0, 52.0, 54.0, 56.1, 58.3, 60.6, 63.0, 65.5,68.1, 70.8, 73.6, 76.6, 79.7, 82.9, 86.2, 89.7,93.3, 97.1, 101.0, 105.1, 109.4, 113.9, 118.6,123.4, 128.5, 133.8, 139.3, 145.0, 151.0, 157.2,163.7, 170.5, 177.5, 184.8, 192.4, 200.0, 208.7,217.3, 226.3, 235.7, 245.5, 255.7, 266.3, 277.4,288.9, 300.9, 313.4, 326.5, 340.1, 350.058hlorpropamide 30 3.00, 3.12, 3.25, 3.38, 3.52, 3.67, 3.83, 4.00,4.18, 4.37, 4.57, 4.78, 5.00, 5.24, 5.49, 5.75,6.03, 6.32, 6.63, 6.96, 7.31, 7.68, 8.07, 8.48,8.92, 9.38, 9.87, 10.39, 10.94, 11.52, 12.13,12.78, 13.46, 14.10, 14.18, 14.94, 15.75, 16.20,16.60, 17.50, 18.46, 19.47, 20.54, 21.67, 22.87,24.14, 25.48, 26.9, 28.4, 29.98, 31.65, 33.42,35.29, 37.27, 39.36, 41.58, 43.92, 45.80, 46.40,49.02, 51.79, 54.72, 57.82, 61.10, 64.57, 67.30,68.24, 72.12, 74.80, 75.10, 76.22, 80.56, 85.15,90.00, 95.13, 100.56, 106.30, 112.37, 118.79,125.58, 132.76, 140.36, 148.39, 156.89,160.90, 165.88, 175.38, 183.60, 185.43,187.90, 194.40, 196.06, 200.00, 205.30,207.30, 219.19, 231.76, 245.06, 253.20,259.13, 274.01, 289.74, 305.70, 306.38,323.98, 332.20, 333.50, 342.59, 362.27,382.20, 383.09, 390.70, 400.00Table 6: Systems that replica exchange we run on, the amount of time run for, and thetemperatures used for exchange. Chlorpropamide was run to lower temperatures to validateour findings at low temperatures and therefore required 10 ns longer of simulation time forthese low temperatures to converge conformationally.
Comparison of Gas &Crystal Torsions
For all molecules that contain 1 or more freerotating torsions, we compared the torsional distributions in the gas state with the crys-tals at 300 K. The gas states were run withREMD from 300 to 1000 K with about 0.2probability exchange between replica.59
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm IForm III
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 2
GasForm IForm III
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 3
GasForm IForm III
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 4
GasForm IForm III
Figure 60: Torsional distributions of piracetam at 300 K for the molecule in gas and bothpolymorphs.
NN NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm IForm III
Figure 61: Torsional distributions of carbamazepine at 300 K for the molecule in gas andboth polymorphs. 60
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm IForm II
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 2
GasForm IForm II
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 3
GasForm IForm II
Figure 62: Torsional distributions of paracetamol at 300 K for the molecule in gas and bothpolymorphs. 61
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm IForm V
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 2
GasForm IForm V
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 3
GasForm IForm V
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 4
GasForm IForm V
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 5
GasForm IForm V
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 6
GasForm IForm V
Figure 63: Torsional distributions of chlorpropamide at 300 K for the molecule in gas andboth polymorphs. 62
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 2
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 3
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 4
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 5
GasForm XForm I
Figure 64: Torsional distributions of aripiprazole at 300 K for the molecule in gas and bothpolymorphs for torsions 1 – 5. 63
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 6
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 7
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 8
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 9
GasForm XForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 10
GasForm XForm I
Figure 65: Torsional distributions of aripiprazole at 300 K for the molecule in gas and bothpolymorphs for torsions 6 – 10. 64
N NH O NO O NH HO HN ON NH O S NHO O NHOCl S NHO O NHOCl Cl N N O HN O
T1 T1T2T3 T1T2 T3T4T1 T1 T2 T3 T4 T5 T6 T7T1 T2 T3 T4 T5 T6T1 T2T3 T4 T5 T7T6 T8 T9 T10
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 1
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 2
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 3
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 4
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 5
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 6
GasForm IIForm I
180 120 60 0 60 120 180
Torsional Angle [Deg.]
Torsion 7
GasForm IIForm I
Figure 66: Torsional distributions of tolbutamide at 300 K for the molecule in gas and bothpolymorphs.
Error in QHA due to Hy-drogen Bonds
We found that there was little to no cor-relation between the hydrogen bonds and the error in anisotropic QHA relative toMD. We determine the number of hydro-gen bonds on a per molecule basis for thelattice minimum (Figure 67) and the per-65ent of hydrogen bonds remaining at 300 K(Figure 68) for MD for all systems to seeif there was any correlation between errorin 1D-QHA from MD ( δ (∆ G QHA )) and thenumber of hydrogen bonds. Despite hav-ing no hydrogen bonds present, diiodoben-zne has the lowest error in the QHA free en-ergy difference and we therefore left it outof both Figure 67 and 68. In both cases, alinear fit shows the R value is less than 0.1showing that there is little to no correlation. H-Bonds at 0 K [bonds/mol] ( G Q H A ) [ k c a l / m o l ] R = 0.012 Figure 67: Error in the polymorph free en-ergy differences for 1D-QHA from MD ver-sus the number of hydrogen bonds in thelattice minimum structure. There is no cor-relation between QHA error and 0 K hydro-gen bonds. % H-Bonds from 0 K (T = 300 K) ( G Q H A ) [ k c a l / m o l ] R = 0.071 Figure 68: Error in the polymorph free en-ergy differences for 1D-QHA from MD ver-sus the percent of hydrogen bonds remain-ing at 300 K relative to the lattice minimum.There is no correlation between QHA errorand 300 K hydrogen bonds.66 eferences (1) Cheuk, D.; Khamar, D.; McArdle, P.;Rasmuson, . C. Solid Forms, CrystalHabits, and Solubility of Danthron.
J. Chem. Eng. Data , , 2110–2118.(2) Maher, A.; Rasmuson, . C.; Cro-ker, D. M.; Hodnett, B. K. Solubilityof the Metastable Polymorph of Pirac-etam (Form II) in a Range of Solvents. J. Chem. Eng. Data , , 3525–3531.(3) Valle, R. G. D.; Brillante, A.; Fa-rina, L.; Venuti, E.; Masino, M.;Girlando, A. Organic Semiconduc-tors: Polymorphism, Phonon Dynam-ics and Carrier-Phonon Coupling inPentacene. Mol. Cryst. Liq. Cryst. Sci. , , 145–154.(4) Stevens, L. A.; Goetz, K. P.;Fonari, A.; Shu, Y.;Williamson, R. M.; Brdas, J.-L.;Coropceanu, V.; Jurchescu, O. D.;Collis, G. E. Temperature-MediatedPolymorphism in Molecular Crystals:The Impact on Crystal Packing andCharge Transport. Chem. Mater. , , 112–118.(5) Singhal, D.; Curatolo, W. Drug poly-morphism and dosage form design: apractical perspective. Adv. Drug Deliv.Rev , , 335–347.(6) Ghosh, M.; Banerjee, S.; Khan, S.;Sikder, N.; Sikder, A. UnderstandingMetastable Phase Transformation dur-ing Crystallization of RDX, HMX andCL-20: Experimental and DFT Stud-ies. Phys. Chem. Chem. Phys. ,23554–23571. (7) Dreger, Z. A.; Tao, Y.; Gupta, Y. M.Phase Diagram and Decompositionof 1,1-Diamino-2,2-dinitroethene Sin-gle Crystals at High Pressures andTemperatures.
J. Phys. Chem. C , 11092–11098.(8) van der Heijden, A. E. D. M.;Bouma, R. H. B. Crystallization andCharacterization of RDX, HMX, andCL-20.
Cryst. Growth Des. , ,999–1007.(9) Grzesiak, A. L.; Lang, M.; Kim, K.;Matzger, A. J. Comparison of theFour Anhydrous Polymorphs of Car-bamazepine and the Crystal Structureof Form I. J. Pharm. Sci. , ,2260–2271.(10) Yoshino, M.; Takahashi, K.;Okuda, Y.; Yoshizawa, T.;Fukushima, N.; Naoki, M. Con-tribution of Hydrogen Bonds toEquilibrium Transition of Resorci-nol. J. Phys. Chem. A , ,2775–2783.(11) Vemavarapu, C.; Mollan, M. J.; Need-ham, T. E. Crystal doping aided byrapid expansion of supercritical solu-tions. AAPS PharmSciTech , ,17.(12) Badea, E.; Blanco, I.; Della Gatta, G.Fusion and solid-to-solid transitionsof a homologous series of alkane-,-dinitriles. J. Chem. Thermodyn. , , 1392–1398.(13) Yu, L.; Huang, J.; Jones, K. J. Mea-suring Free-Energy Difference betweenCrystal Polymorphs through EutecticMelting. J. Phys. Chem. B , ,19915–19922.6714) Cherukuvada, S.; Thakuria, R.; Nan-gia, A. Pyrazinamide Polymorphs:Relative Stability and VibrationalSpectroscopy. Cryst. Growth Des. , , 3931–3941.(15) Torrisi, A.; Leech, C. K.; Shank-land, K.; David, W. I. F.; Ib-berson, R. M.; Benet-Buchholz, J.;Boese, R.; Leslie, M.; Catlow, C.R. A.; Price, S. L. Solid Phases ofCyclopentane: Combined Experimen-tal and Simulation Study. J. Phys.Chem. B , , 3746–3758.(16) Stolar, T.; Lukin, S.; Poar, J.;Rubi, M.; Day, G. M.; Biljan, I.;Jung, D. .; Horvat, G.; Uarevi, K.;Metrovi, E.; Halasz, I. Solid-StateChemistry and Polymorphism of theNucleobase Adenine. Cryst. GrowthDes. , , 3262–3270.(17) Alcob, X.; Estop, E.; Aliev, A. E.;Harris, K. D. M.; Rodrguez-Carvajal, J.; Rius, J. Temperature-Dependent Structural Properties ofp-Diiodobenzene: Neutron Diffractionand High-Resolution Solid State 13CNMR Investigations. J. Solid StateChem. , , 20–27.(18) Sacchetti, M. Thermodynamic Analy-sis of DSC Data for AcetaminophenPolymorphs. J. Therm. Anal. Calorim , , 345–350.(19) Cesaro, A.; Starec, G. Thermody-namic properties of caffeine crystalforms. J. Phys. Chem. , , 1345–1346.(20) Seryotkin, Y. V.; Drebushchak, T. N.;Boldyreva, E. V. A high-pressure poly-morph of chlorpropamide formed onhydrostatic compression of the -form in saturated ethanol solution. ActaCryst B , , 77–85.(21) Boldyreva, E. V.; Shakhtshnei-der, T. P.; Ahsbahs, H.; Sowa, H.;Uchtmann, H. Effect of High Pressureon the Polymorphs of Paracetamol. J. Therm. Anal. Calorim , ,437–452.(22) Cansell, F.; Fabre, D.; Petitet, J.Phase transitions and chemical trans-formations of benzene up to 550C and30 GPa. J. Chem. Phys. , ,7300–7304.(23) Gajda, R.; Katrusiak, A. Pressure-Promoted CHO Hydrogen Bonds inFormamide Aggregates. Cryst. GrowthDes. , , 4768–4774.(24) Dybeck, E. C.; Abraham, N. S.;Schieber, N. P.; Shirts, M. R. Cap-turing Entropic Contributions toTemperature-Mediated PolymorphicTransformations Through MolecularModeling. Cryst. Growth Des. , , 1775–1787.(25) Abraham, N. S.; Shirts, M. R. Ther-mal Gradient Approach for the Quasi-harmonic Approximation and Its Ap-plication to Improved Treatment ofAnisotropic Expansion. J. Chem. The-ory Comput. , , 5904–5919.(26) Erba, A.; Shahrokhi, M.; Mora-dian, R.; Dovesi, R. On how differ-ently the quasi-harmonic approxima-tion works for two isostructural crys-tals: Thermal properties of periclaseand lime. J. Chem. Phys. , ,044114.(27) Erba, A.; Maul, J.; Demichelis, R.;Dovesi, R. Assessing thermochemical68roperties of materials through ab ini-tio quantum-mechanical methods: thecase of -Al2O3. Phys. Chem. Chem.Phys. , , 11670–11677.(28) Erba, A.; Maul, J.; De La Pierre, M.;Dovesi, R. Structural and elasticanisotropy of crystals at high pres-sures and temperatures from quan-tum mechanical methods: The caseof Mg2SiO4 forsterite. J. Chem. Phys. , , 204502.(29) Erba, A. On combining temperatureand pressure effects on structuralproperties of crystals with standardab initio techniques. J. Chem. Phys. , , 124115.(30) Maul, J.; Santos, I. M. G.; Sam-brano, J. R.; Erba, A. Thermal proper-ties of the orthorhombic CaSnO3 per-ovskite under pressure from ab initioquasi-harmonic calculations. Theor.Chem. Acc , , 36.(31) Dybeck, E. C.; Schieber, N. P.;Shirts, M. R. Effects of a More Accu-rate Polarizable Hamiltonian on Poly-morph Free Energies Computed Effi-ciently by Reweighting Point-ChargePotentials. J. Chem. Theory Comput. , , 3491–3505.(32) Hore, S.; Dinnebier, R.; Wen, W.;Hanson, J.; Maier, J. Structure ofplastic crystalline succinonitrile: High-resolution in situ powder diffraction. Z. Anorg. Allg. Chem. , , 88–93.(33) Heit, Y. N.; Nanda, K. D.; Beran, G.J. O. Predicting finite-temperatureproperties of crystalline carbon diox-ide from first principles with quanti-tative accuracy. Chem. Sci. , ,246–255. (34) Ramrez, R.; Neuerburg, N.; Fernndez-Serra, M.-V.; Herrero, C. P. Quasi-harmonic approximation of thermody-namic properties of ice Ih, II, and III. J. Chem. Phys. , , 044502.(35) Nath, P.; Plata, J. J.; Usanmaz, D.;Al Rahal Al Orabi, R.; Fornari, M.;Nardelli, M. B.; Toher, C.; Cur-tarolo, S. High-throughput predictionof finite-temperature properties us-ing the quasi-harmonic approximation. Comput. Mater. Sci , , 82–91.(36) Huang, L.-F.; Lu, X.-Z.; Ten-nessen, E.; Rondinelli, J. M. Anefficient ab-initio quasiharmonic ap-proach for the thermodynamics ofsolids. Comput. Mater. Sci , ,84–93.(37) Choy, C. L.; Wong, S. P.; Young, K.Thermal expansion and Gr¨uneisen pa-rameters for anisotropic solids. Phys.Rev. B , , 1741–1747.(38) Gould, S.; Fernando, B.; Cherian, A.;Anderson, P.; Cruz, R. S.; Guo, E. OnDifferentiating Parameterized Argminand Argmax Problems with Ap-plication to Bi-level Optimization. arXiv:1607.05447 [cs, math] ,arXiv: 1607.05447.(39) Zhang, W.; Wu, C.; Duan, Y. Con-vergence of replica exchange moleculardynamics. J. Chem. Phys. , ,154105.(40) Robertson, M. J.; Tirado-Rives, J.;Jorgensen, W. L. Improved Peptideand Protein Torsional Energetics withthe OPLS-AA Force Field. J. Chem.Theory Comput. , , 3499–3509.(41) Jorgensen, W. L.; Maxwell, D. S.;Tirado-Rives, J. Development andTesting of the OPLS All-Atom Force69ield on Conformational Energeticsand Properties of Organic Liquids. J.Am. Chem. Soc. , , 11225–11236.(42) Hairer, E.; Nrsett, S. P.; Wanner, G. Solving Ordinary Differential Equa-tions I: Nonstiff Problems ; SpringerScience & Business Media, 2008.(43) Kimura, K.; Hirayama, F.; Uekama, K.Characterization of tolbutamide poly-morphs (burger’s forms II and IV)and polymorphic transition behavior.
J. Pharm. Sci. , , 385–391.(44) Braun, D. E.; Gelbrich, T.; Kahlen-berg, V.; Tessadri, R.; Wieser, J.;Griesser, U. J. Conformational poly-morphism in aripiprazole: Prepara-tion, stability and structure of fivemodifications. J. Pharm. Sci. , , 2010–2026.(45) van Miltenburg, J. C.; Oonk, H.A. J.; van den Berg, G. J. K.Low-Temperature Heat Capacities andDerived Thermodynamic Functionsof Para-Substituted Halogen Ben-zenes. 2. p-Bromoiodobenzene and p-Diiodobenzene. J. Chem. Eng. Data , , 84–89.(46) Drebushchak, V. A.;Drebushchak, T. N.; Chukanov, N. V.;Boldyreva, E. V. Transitions amongfive polymorphs of chlorpropamidenear the melting point. J. Therm.Anal. Calorim. , , 343–351.(47) Brandenburg, J. G.; Potticary, J.;Sparkes, H. A.; Price, S. L.;Hall, S. R. Thermal Expansion ofCarbamazepine: Systematic Crystal-lographic Measurements ChallengeQuantum Chemical Calculations. J.Phys. Chem. Lett. , , 4319–4324. (48) Wilson, C. Variable temperature studyof the crystal structure of paracetamol(p-hydroxyacetanilide), by single crys-tal neutron diffraction. Z. Kristallogr.Cryst. Mater. , , 693–701.(49) Horstman, E. M.; Goyal, S.;Pawate, A.; Lee, G.; Zhang, G. G. Z.;Gong, Y.; Kenis, P. J. A. Crystalliza-tion Optimization of PharmaceuticalSolid Forms with X-ray CompatibleMicrofluidic Platforms. Cryst. GrowthDes. , , 1201–1209.(50) Lisgarten, J. N.; Palmer, R. A.;Saldanha, J. W. Crystal and molec-ular structure of 5-carbamyl-5H-dibenzo[b,f] azepine. J. Crystallogr.Spectrosc. Res. , , 641–649.(51) Nievergelt, P. P.; Spingler, B. Growingsingle crystals of small molecules bythermal recrystallization, a viable op-tion even for minute amounts of mate-rial? CrystEngComm , , 142–147.(52) Sovago, I.; Gutmann, M. J.;Senn, H. M.; Thomas, L. H.; Wil-son, C. C.; Farrugia, L. J. Electrondensity, disorder and polymorphism:high-resolution diffraction studies ofthe highly polymorphic neuralgic drugcarbamazepine. Acta Cryst B , , 39–50.(53) El Hassan, N.; Ikni, A.; Gillet, J.-M.;Spasojevic-de Bir, A.; Ghermani, N. E.Electron Properties of CarbamazepineDrug in Form III. Cryst. Growth Des. , , 2887–2896.(54) Haisa, M.; Kashino, S.; Kawai, R.;Maeda, H. The Monoclinic Form ofp-Hydroxyacetanilide. Acta Cryst B , , 1283–1285.7055) Wilson, C. C.; Shankland, N.; Flo-rence, A. J.; Frampton, C. S. Single-crystal neutron diffraction of bio-active small molecules. Physica B , , 84–86.(56) Wilson, C. C. Neutron diffraction of p-hydroxyacetanilide (Paracetamol): li-bration or disorder of the methyl groupat 100 K. J. Mol. Struct , ,207–217.(57) Nyman, J.; Day, G. M. Static andlattice vibrational energy differencesbetween polymorphs. CrystEngComm , , 5154–5165.(58) Dybeck, E. C.; McMahon, D. P.;Day, G. M.; Shirts, M. R. Exploringthe Multi-minima Behavior of SmallMolecule Crystal Polymorphs at FiniteTemperature. Crystal Growth & De-sign , , 5568–5580.(59) Fabbiani, F. P. A.; Allan, D. R.; Par-sons, S.; Pulham, C. R. An explo-ration of the polymorphism of pirac-etam using high pressure. CrystEng-Comm , , 179–186.(60) Lou¨er, D.; Lou¨er, M.;Dzyabchenko, V. A.; Agafonov, V.;Ceolin, R. Structure of a metastablephase of piracetami from X-ray pow-der diffraction using the atom–atompotential method. Acta Cryst B , , 182–187.(61) Admiraal, G.; Eikelenboom, J. C.;Vos, A. Structures of the triclinic andmonoclinic modifications of (2-oxo-1-pyrrolidinyl)acetamide. Acta Cryst B , , 2600–2605.(62) Chambrier, M.-H.; Bouhmaida, N.;Bonhomme, F.; Leb`egue, S.; Gillet, J.-M.; Jelsch, C.; Ghermani, N. E. Elec-tron and Electrostatic Properties of Three Crystal Forms of Piracetam. Cryst. Growth Des. , , 2528–2539.(63) Galdecki, Z.; Glowka, M. L. Pol. J.Chem. , , 1307.(64) Shipilova, A. S.; Knyazev, A. V.;Gusarova, E. V.; Amosov, A. A.;Knyazeva, S. S. X-ray studies of pirac-etam in a wide range of temperatures. J. Solid State Chem. , , 51–56.(65) Tilborg, A.; Jacquemin, D.; Nor-berg, B.; Perp`ete, E.; Michaux, C.;Wouters, J. Structural study of pirac-etam polymorphs and cocrystals: crys-tallography redetermination and quan-tum mechanics calculations. ActaCryst B , , 499–507.(66) Fernandes, P.; Shankland, K.; Flo-rence, A. J.; Shankland, N.; John-ston, A. Solving Molecular CrystalStructures from X-ray Powder Diffrac-tion Data: The Challenges Posed by γ -Carbamazepine and ChlorothiazideN,N,-Dimethylformamide (1/2) Sol-vate. J. Pharm. Sci , , 1192–1202.(67) Mahapatra, S.; Nayak, S. K.;Prathapa, S. J.; Guru Row, T. N.Anhydrous Adenine: Crystallization,Structure, and Correlation with OtherNucleobases. Cryst. Growth Des. , , 1223–1225.(68) Nangia, A.; Srinivasulu, A. CSDComm. ,(69) Rø, G.; Sørum, H. The crystaland molecular structure of δ -pyrazinecarboxamide. Acta CrystB , , 1677–1684.7170) Castro, R. A. E.; Maria, T. M. R.;´Evora, A. O. L.; Feiteira, J. C.;Silva, M. R.; Beja, A. M.; Can-otilho, J.; Eus´ebio, M. E. S. A New In-sight into Pyrazinamide PolymorphicForms and their Thermodynamic Re-lationships. Cryst. Growth Des. , , 274–282.(71) Nakata, K.; Takaki, Y. Memoirs of Os-aka Kyoiku University , , 93.(72) Tamura, C.; Kuwano, H. Crystal-lographic data of carboxylic acidsand carboxyamides of picoline andpyrazine derivatives. Acta Cryst , , 693.(73) Bacon, G. E.; Lisher, E. J.; Paw-ley, G. S. A neutron powder diffrac-tion study of deuterated `a-resorcinol:a test of profile refinement using TLSconstraints. Acta Cryst B , ,1400–1403.(74) Bacon, G. E.; Lisher, E. J. A neutronpowder diffraction study of deuter-ated α - and β -resorcinol. Acta CrystB , , 1908–1916.(75) Dru˙zbicki, K.; Mikuli, E.;Pa(cid:32)lka, N.; Zalewski, S.; Ossowska-Chru´sciel, M. D. Polymorphism ofResorcinol Explored by Comple-mentary Vibrational Spectroscopy(FT-RS, THz-TDS, INS) and First-Principles Solid-State Computations(Plane-Wave DFT). J. Phys. Chem.B , , 1681–1695.(76) Fronczek, F. R. CSD Comm. ,(77) Robertson J. Monteath,; UbbelohdeA. R., A New Form of Resorcinol. I.Structure Determination by X-Rays.
Proc. Royal Soc. Lond. , ,122–135. (78) Hendricks, S. B.; Maxwell, L. R.;Mosley, V. L.; Jefferson, M. E. X-Rayand Electron Diffraction of Iodine andthe Diiodobenzenes. J. Chem. Phys. , , 549–565.(79) Hinchliffe, A.; Munn, R. W.;Pritchard, R. G.; Spicer, C. J.The structure of a solution-growncrystal of 1,4-diiodobenzene. J. Mol.Struct , , 93–96.(80) Haisa, M.; Kashino, S.; Maeda, H.The orthorhombic form of p-hydroxyacetanilide. Acta Cryst B , , 2510–2512.(81) Drebushchak, T. N.; Boldyreva, E. V.Variable temperature (100–360 K)single-crystal X-ray diffraction studyof the orthorhombic polymorph ofparacetamol (p-hydroxyacetanilide). Z. Kristallogr. Cryst. Mater. , , 506–512.(82) Nichols, G.; Frampton, C. S. Physico-chemical Characterization of the Or-thorhombic Polymorph of Paraceta-mol Crystallized from Solution. J.Pharm. Sci. , , 684–693.(83) Donaldson, J. D.; Leary, J. R.;Ross, S. D.; Thomas, M. J. K.;Smith, C. H. The structure of the or-thorhombic form of tolbutamide (1-n-butyl-3-p-toluenesulphonylurea). ActaCryst B , , 2245–2248.(84) Drebushchak, T. N.;Drebushchak, V. A.;Pankrushina, N. A.; Boldyreva, E. V.Single-crystal to single-crystal confor-mational polymorphic transformationin tolbutamide at 313 K. Relation toother polymorphic transformationsin tolbutamide and chlorpropamide. CrystEngComm , , 5736–5743.7285) Schreyer, M.; Guo, L.; Thirunahari, S.;Gao, F.; Garland, M. Simultaneous de-termination of several crystal struc-tures from powder mixtures: the com-bination of powder X-ray diffraction,band-target entropy minimization andRietveld methods. J Appl Cryst , , 659–667.(86) Thirunahari, S.; Aitipamula, S.;Chow, P. S.; Tan, R. B. H. Conforma-tional Polymorphism of Tolbutamide:A Structural, Spectroscopic, andThermodynamic Characterization ofBurger’s Forms I–IV. J. Pharm. Sci. , , 2975–2990. (87) Drebushchak, T. N.; Chesalov, Y. A.;Boldyreva, E. V. A conformationalpolymorphic transition in the high-temperature (cid:15) -form of chlorpropamideon cooling: a new (cid:15) ’-form. Acta CrystB , , 770–781.(88) Koo, C. H.; Cho, S. I.; Yeon, Y. H.The crystal and molecular structureof chlorpropamide. Arch. Pharm. Res. , , 37–49.(89) Drebushchak, T. N.; Chukanov, N. V.;Boldyreva, E. V. Two polymorphs ofchlorpropamide: the δ -form and thehigh-temperature (cid:15) -form. Acta CrystC ,64