Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods† †Electronic supplementary information (ESI) available: Additional computational details, tables of energies, analysis of model convergence, & optimized crystal structures. See DOI: 10.1039/c9sc05689k

Widely used crystal structure prediction models based on density functional theory can perform poorly for conformational polymorphs, but a new model corrects those polymorph stability rankings.


S1.1 DFT Monkhorst-Pack grids and quality of the optimized structures
The planewave DFT Monkhorst-Pack grids indicated in Table S1 were used for the B86bPBE-XDM geometry optimizations and single-point energies. These choices were made by systematically increasing the grid until the geometry no longer changed. Roughly speaking, 9 k-points were used for lattice constants below 4Å, 7 k-points were used for lattice constants of 4-6Å, 5 k-points for lattice constants of 6-8Å, 3 k-points for lattice constants 8-12Å, and 1 k-point for larger lattice constants. Table S1 also indicates the rmsd15 value for the 0 K (fully relaxed) and 298 K (fixed experimental lattice parameters) B86bPBE-XDM geometry optimizations. The rmsd15 metric indicates the root-mean-square deviation in the non-hydrogen atomic positions for a cluster of 15 molecules in each crystal. 1 For o-acetamidobenzamide, the room-temperature x-ray diffraction structures ACBNZA and ACBNZA01 were chosen over the more recent low-temperature neutron diffraction structures (ACBNZA02 and ACBNZA03) 2 to enable examination of the temperature-dependent behaviors. Full B86bPBE-XDM relaxation of the neutron diffraction structures resulted in structures that were virtually identical to those used here (rmsd15 differences less than 0.04Å and energy differences of 0.1 kJ/mol or less).
An earlier study 3 employing the PBE-NP dispersion-corrected DFT functional found that the ON structure relaxed considerably upon optimization, with the key intramolecular dihedral angle θ changing by ∼25 • . This suggested that DFT might be particularly problematic for the ON polymorph. However, the B86bPBE-XDM functional used here does not appear to suffer from the same problems in reproducing the experimental crystal structures. For the fixed-cell optimization, the dihedral angle θ changes only 2.6 • , which is on par with the average 3.4 ± 2.4 • change observed for the other seven polymorphs (Table S2). Relaxing the cell fully leads to a larger 9.5 • change in that dihedral angle, but this is again on par with the average change of 8.2 ± 3.6 • change seen for the other forms. Structure overlays give rmsd15 values of 0.05Å for the fixed cell and 0.22Å for the fully relaxed cell relative to experiment, which are again compatible with values seen for the other crystals (Table S1). In other words, there is no obvious structural discrepancy between the B86bPBE-XDM and experimentally reported structures that might account for the disagreement in relative stabilities.

S1.2 Basis set dependence for the MP2 models
The HMBI calculations here rely on monomer and dimer MP2 calculations. MP2-based models exhibit slower convergence with basis set than do HF or DFT models, and the use of large basis sets is important for obtaining well-converged polymorph stability predictions.
In this work, all MP2-based methods were extrapolated to the CBS limit from aug-cc-pVTZ and aug-cc-pVQZ results. Additionally, the MP2D dispersion correction was designed for use with large basis sets, and the short-range damping functions were fitted against CBS-limit results. 4 Figures S1 and S2 show the dramatic impact of basis set completeness on the predicted polymorph stabilities for acetamidobenzamide and ROY. In acetamidobenzamide, the relative stability of the two forms changes by more than 6 kJ/mol between aug-cc-pVDZ and the CBS limit. Examination of the 1-body and 2-body contributions (not explicitly shown here) indicates that the basis set dependence is approximately evenly split between the intraand intermolecular contributions. The ROY polymorph stabilities also exhibit strong basis set dependence. For polymorph ORP, the relative stability decreases 6.9 kJ/mol between aug-cc-pVDZ and the CBS limit. The average magnitude of the basis set change between aug-cc-pVDZ and the CBS limit is 3.6 kJ/mol, which again is substantial on the energy scale of polymorphism. The basis set dependence for oxalyl dihydrazide has been examined previously, 5,6 and similar results were found to those shown here for ROY and acetamidobenzamide.

Acetamidobenzamide basis set dependence
Energy Relative to α Form (kJ/mol) Figure S1: Dependence of the polymorph stabilities for o-acetamidobenzamide on the basis set used for the MP2D 1 & 2-body contributions. In all cases, periodic HF/pob-TZVP-rev2 was used for the many-body contribution.

ROY Polymorph basis set dependence
Relative Energy (kJ/mol) Figure S2: Dependence of the ROY polymorph stabilities on the basis set used for the MP2D 1 & 2-body contributions. In all cases, periodic HF/pob-TZVP-rev2 was used for the many-body contribution.

S1.3 Selection of the Hartree-Fock many-body treatment basis set
The 1-body (monomer) and 2-body (dimer) terms frequently account for 90% or more of the total lattice energy in the HMBI model, while the remaining ∼10% results from the intermolecular many-body contributions. However, in the polymorphic systems considered in this paper, the proportional impact of the many-body term on the relative polymorph energies is considerably larger. Therefore, it is important to obtain a well-converged description of the many-body term. Periodic HF is relatively inexpensive, particularly when Gaussian basis sets are employed, and can describe the frequently dominant many-body polarization effects.
Periodic HF (and DFT) calculations with large Gaussian basis sets suffer from well-known issues with linear dependencies and poor convergence of the self-consistent-field equations. Adapting a Gaussian basis sets designed for molecular systems to periodic calculations typically requires eliminating the most diffuse basis functions and subsequent re-tuning of the remaining basis functions. Several triple-ζ basis sets have been developed in recent years, including the pob-TZVP, 7 pob-TZVP-rev2, 8 and mTZVP basis set. 9 All three are modified versions of the popular def2-TZVP basis intended for periodic systems.
To assess the performance of these three basis sets, 15 five-molecule clusters were extracted from the oxalyl dihydrazide polymorphs. HF calculations were performed on the entire gas-phase cluster, and then separately subtracting out the HF monomer and dimer contributions according to a many-body expansion. Consistent with the HMBI treatment of the many-body terms, no counterpoise correction is employed in either the full cluster nor the fragment calculations. The HF/def2-QZVP results are chosen as the benchmark, and these results are compared against the traditional def2-TZVP basis set and the three aforementioned triple-ζ basis sets that are suitable for periodic calculations.
As shown in Table S3, the def2-TZVP basis set gives the root-mean-square (RMS) error of 0.31 kJ/mol per monomer versus def2-QZVP. The pob-TZVP-rev2 basis revised the pob-TZVP basis in order to reduce the occurrence of basis set superposition error (BSSE). While BSSE is more prevalent in dimer interactions, the results here show that the pob-TZVP-rev2 basis gives considerably smaller errors (0.31 kJ/mol) than pob-TZVP (0.63 kJ/mol) relative to the benchmark def2-QZVP values. Finally, the mTZVP basis, which typically exhibits even less BSSE, predicts many-body contributions with errors intermediate between the two other basis sets. In the end, the pob-TZVP-rev2 basis was adopted in this study based on these benchmarks. S2 o-Acetamidobenzamide polymorphs

S2.1 Enthalpies and Gibbs Free Energies
To augment the enthalpy and free energy curves shown in the main paper, Figure S4 plots the corresponding semi-schematic G and H curves for all four methods considered in the main paper. All three MP2-based models agree that the β form is thermodynamically preferred at high temperatures, though the specific free energies differ. MP2 predicts an enantiotropic relationship, while MP2C predicts a monotropic one. For MP2D, the relationship is also enantiotropic, though the two forms are nearly degenerate at 0 K. Consistent with experiment, all three MP2-based methods predict an exothermic ∆H α→β phase transition at elevated temperatures. In contrast, B86bPBE-XDM incorrectly predicts an endothermic phase transition and that the α form is thermodynamically more stable at all temperatures. Examination of the ∆H α→β values shown in Table 1 and Figure 3 of the main manuscript suggests that MP2C overestimates the stability of the β form due largely to the neglect of intramolecular dispersion. To understand the temperature dependence of the o-acetamidobenzamide predictions better, Table S4 decomposes the temperature dependence into lattice energy contributions (due to thermal expansion) and phonon contributions. Note that the MP2-type results use B86bPBE-XDM phonon frequencies, meaning that the vibrational contributions to the enthalpy H vib and free energy G vib are identical across all methods here.
For B86bPBE-XDM, the lattice energy changes stabilize the β form relative to α by less than 1 kJ/mol over the 0-423 K temperature range, while the vibrational enthalpy destabilizes the β form by 2.3 kJ/mol. The MP2-based models all predict more substantial lattice energy-driven stabilization of the β polymorph with increasing temperature, ranging from -7.3 kJ/mol for MP2 to -4.3-4.4 for the dispersion-corrected MP2D and MP2C models. These lattice energy changes favoring the β form at higher temperatures are partially canceled by the phonon contributions that increasingly favor the α form at higher temperatures, leading to the final predicted MP2D enthalpy that agrees very well with the experimentally measured transition enthalpy, as described in the main paper. Table S4: Lattice energy differences, vibrational enthalpy contributions, and vibrational Gibbs free energy contributions to the α-β polymorph energy differences in o-acetamidobenzamide, in kJ/mol. Positive values indicate α is more stable than β. Summing the lattice energy and vibrational enthalpy contributions gives the results plotted in Figure 3 of the main paper.

B86bPBE-XDM Lattice Energy Contributions
Phonon Contrib.  Table S5 summarizes the lattice energies computed here that were used to generate Figure 5 in the main paper.  Figure S5 compares HMBI results with MP2 and MP2D. The differences in the relative polymorph energies are less than 1 kJ/mol throughout. In all cases except ORP, the dispersion correction slightly destabilizes the polymorphs relative to Y. Note that although the impact on the relative polymorph energies is small in this instance, the dispersion correction has a big impact on the overall lattice energies. For example, MP2D reduces the strength of the 2-body interaction energies by 27-31 kJ/mol, depending on the polymorph.

ROY Polymorphs
Relative Energy (kJ/mol) Figure S5: The impact of the MP2D dispersion correction on the relative polymorph energies is fairly modest in ROY, as compared to uncorrected MP2. Figure S6 plots the relative energies with DFT and MP2D for the fully relaxed 0 K structures and the fixed-cell room-temperature structures. Note that the results here compare MP2D in the aug-cc-pVTZ basis, rather than in the CBS limit like all other results presented here. This basis set was chosen for the 0 K structures to reduce computational costs. At the DFT level, the differences in the relative lattice energies for the two sets of structures are quite small. For MP2D, they are considerably larger. The reasons behind the larger MP2D differences are not entirely clear. One potential explanation could be that the optimal MP2D geometries might differ somewhat from the DFT ones, and therefore the DFT optimal structures correspond to steeper portions of the MP2D potential energy surface (and therefore one obtains larger changes in the energies with the geometry changes).

S3.4 Omission of the R05 polymorph
The R05 polymorph is omitted from the present study because it was recently discovered that the HMBI model energies are poorly convergent in polar unit cells. The treatment of long-range electrostatics is generally problematic in periodic cells with a net dipole moment. In practice, this issue is frequently managed via application of tin-foil boundary conditions.
The working HMBI expression has the following form: where the 2-body terms are computed out to a cutoff radius of 9-10Å. A smoothing function is applied across this interval to ensure smooth potential energy surfaces as the model transitions from the high-level (MP2D) to low-level (HF). The periodic HF contribution does employ tin-foil boundary conditions. In non-polar cells, one typically observes well-behaved convergence of the lattice energy as the 2-body cutoff radius is increased, as expected from the long-range decay of the interaction energies.
In polar unit cells, however, the contributions from the 2-body terms proves highly erratic with cutoff radius. For example, in γ-glycine, which has a polar unit cell, the energy relative to the non-polar α form can vary by many kJ/mol depending on the specific cutoff. For this reason, results on the R05 polymorph are omitted from the present study. Table S6 presents the relative lattice energies for the oxalyl dihydrazide polymorphs using the 0 K and room-temperature structures-the same data used to generate Figure 9 in the main paper. For MP2 and CCSD(T), the many-body treatment is computed at the periodic HF/pob-TZVP-rev2 level of theory. Table S6: Relative lattice energies (in kJ/mol) for the polymorphs of oxalyl dihydrazide using the fully relaxed 0 K structures or fixed-cell room-temperature structures.  Table S7 shows how the lattice energies, H vib and G vib of the different polymorphs change relative to the α form between the 0 K and room temperature structures. The lattice energy contributions here amount to the differences between the 0 K and 298 K results in Table S6. At the B86bPBE-XDM, the changes in the relative lattice energies are relative modest, except for the β form, which is destabilized by 3.8 kJ/mol. For the MP2-based methods, the β form is destabilized slightly less, by ∼3 kJ/mol. On the other hand, increasing temperature stabilizes the γ, δ, and forms by 1-2 kJ/mol relative to α. The temperature-dependence of the vibrational contributions to the enthalpy are relatively modest as well, with the other four forms being destabilized by 1-2 kJ/mol relative to the α form. In the end, the lattice energies and enthalpies predict the same qualitative polymorph stability orderings, with only small changes in the quantitative values, as can be seen from Figure S7. The key temperature-dependent change in the β/γ MP2D stability ordering arises from lattice energy contributions Table S7: Change in the relative lattice energies, vibrational enthalpy contributions, and vibrational Gibbs free energy contributions of oxalyl dihydrazide polymorph stabilities between the 0 K and room-temperature structures, in kJ/mol. Values are listed relative to the α polymorph.   Figure S7: Temperature dependence of the relative lattice energies (dotted lines) and enthalpies (solid lines) for the oxalyl dihydrazide polymorphs. The qualitative stability ordering is unchanged for both lattice energies and enthalpies. Table S8 shows energy decomposition data for oxalyl dihydrazide that was used to generate Figure 11 in the main paper. MP2 and MP2D give 1-body energies in much better agreement with the CCSD(T) benchmarks than does B86bPBE-XDM. For the two-body energies, B86bPBE-XDM does extremely well in the "relative" RMS error vs. CCSD(T) (i.e. as computed from the data shown in the table, where stabilities are specified relative to the α polymorph). However, if one instead computes the "absolute" RMS error based on the total 2-body interaction energies for each form (i.e. interaction energies relative to non-interacting monomers, summed over all dimers), the DFT functional performs far worse. MP2 performs even worse than B86bPBE-XDM, while MP2C does somewhat better on the absolute RMS error. MP2D reproduces CCSD(T) extremely well, with 1.6 kJ/mol errors in both the absolute and relative RMS errors. Table S8: Decomposition of the relative polymorph electronic energies (kJ/mol) for oxalyl dihydrazide 0 K structures into intramolecular (1-body), pairwise intermolecular, and many-body intermolecular contributions. RMS errors are computed against the CCSD(T) benchmarks. All MP2 and CCSD(T) models employ the same periodic HF many-body treatment.