Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods

Abstract

Molecular crystal structure prediction is increasingly being applied to study the solid form landscapes of larger, more flexible pharmaceutical molecules. Despite many successes in crystal structure prediction, van der Waals-inclusive density functional theory (DFT) methods exhibit serious failures predicting the polymorph stabilities for a number of systems exhibiting conformational polymorphism, where changes in intramolecular conformation lead to different intermolecular crystal packings. Here, the stabilities of the conformational polymorphs of o-acetamidobenzamide, ROY, and oxalyl dihydrazide are examined in detail. DFT functionals that have previously been very successful in crystal structure prediction perform poorly in all three systems, due primarily to the poor intramolecular conformational energies, but also due to the intermolecular description in oxalyl dihydrazide. In all three cases, a fragment-based dispersion-corrected second-order Møller–Plesset perturbation theory (MP2D) treatment of the crystals overcomes these difficulties and predicts conformational polymorph stabilities in good agreement with experiment. These results highlight the need for methods which go beyond current-generation DFT functionals to make crystal polymorph stability predictions truly reliable.

---

1 Introduction

Crystal packing influences the physical properties of organic crystals. The occurrence of multiple crystalline packing motifs, or polymorphs, of a pharmaceutical can impact its solubility, bioavailability, shelf-life/stability, and tabletting properties, for example. The importance of polymorphism to the pharmaceutical industry is highlighted by examples such as ritonavir^1,2 and rotigotine,³ where the late-stage appearance of more stable, less soluble crystal forms forced product recalls and reformulations. It was recently suggested that the thermodynamically stable crystal form has not been realized experimentally for ∼15–45% of pharmaceutical molecules,⁴ raising speculation that more such examples may occur in the future. Moreover, solid form patents play an important role in the commercial life cycle of a drug, as evidenced by the recent legal wrangling over a new polymorph of Celgene's blockbuster drug Revlimid that was discovered by generic drug manufacturer Natco.⁵

The ability to predict the molecular crystal energy landscape, which is the set of possible low-energy crystal structures for a given compound, would be a tremendous boon to the pharmaceutical industry and others. Crystal structure prediction has long been challenging⁶ due to the complexity of the search space, the small energy differences that separate polymorphs, and the complexities of crystallization kinetics. The accuracy requirements for predicting the crystal energy landscape are severe: surveys suggest that about half of all polymorph pairs are separated by less than 2 kJ mol⁻¹ in lattice energy, and around 95% are separated by less than 8 kJ mol⁻¹.^7–9

The advent of high-quality, dispersion-corrected density functional theory (DFT) models^10–13 has enabled tremendous progress in the energy ranking aspects of crystal structure prediction, as evidenced by results from the recent blind tests^14–17 and other studies.^18–30 Increasingly, DFT is being called on to explore pharmaceutical crystal energy landscapes as a complement to experimental solid form screening.^31–39 Computational prediction of a highly stable, unrealized polymorph of galunisertib played a key role in the extensive characterization of its solid form landscape, for example.³⁸

Despite many successes of DFT-driven crystal structure prediction, close inspection of the literature also finds polymorphic crystals for which widely-used DFT models fail dramatically. Many of these difficult cases involve conformational polymorphs, in which different intramolecular conformations enable different intermolecular crystal packing motifs. For example, DFT methods invert the polymorph stability ordering of α and β o-acetamidobenzamide, with errors of 5–10 kJ mol⁻¹.⁴⁰ The prolific polymorph-former 5-methyl-2-[(2-nitrophenyl)amino]-3-thiophenecarbonitrile, nicknamed “ROY” after its colorful red-orange-yellow crystals, is another example. State-of-the-art density functional models predict the Y polymorph to be one of the least stable forms,^41,42 when it is actually the most stable one. These backwards stability rankings reflect errors approaching 10 kJ mol⁻¹. In another case, crystal structure prediction failed for two of six conformationally flexible species resulting from mechanochemical aromatic disulfide metathesis reactions, with errors exceeding 6 kJ mol⁻¹ due in large part to poor intramolecular DFT conformational energies.⁴³ Erroneous intramolecular conformational energies caused a similar failure for a recent DFT study of Molecule X from an earlier blind test of crystal structure prediction.²⁴

The large errors in the relative polymorph stabilities found for many of these examples greatly exceed the few kJ mol⁻¹ errors or less one typically finds for DFT in successful crystal structure prediction cases. Furthermore, these errors are catastrophically large compared to the small energy differences that are characteristic of polymorphism. The pharmaceutical industry trend toward developing larger, more flexible drug molecules⁴⁴ makes problems with ranking conformational polymorphs particularly concerning, since it raises the possibility that such ranking problems will become more prevalent as crystal structure prediction is applied to increasingly complicated species.

The problems with popular DFT functionals are not limited to conformational polymorphism either. Delocalization error in commonly used DFT generalized gradient approximation (GGA) functionals can cause spurious salt formation in co-crystals⁴⁵ and the substantial overbinding of crystals containing halogen bonds.⁴⁶ Many functionals erroneously predict the exothermic anthracene photodimerization reaction to be strongly endothermic,^47,48 which is problematic⁴⁹ when studying a class of interesting anthracene-based photomechanical materials.^50,51

Switching to a hybrid density functional can address some of the limitations of GGA-type functionals that cause incorrect polymorph rankings and other problems,^{20,28,45,46,52,53} but it does not rectify the incorrect ROY polymorph rankings,⁴¹ for example. Approaches based on periodic second-order Møller–Plesset perturbation theory (MP2),^54–68 the random phase approximation (RPA),^68–71 and quantum Monte Carlo^72–75 are also being developed that can improve the reliability of polymorph stability rankings. Computational cost is the fundamental challenge that inhibits applying these more accurate electronic structure methods to molecular crystals. As evident from the studies cited above, higher-level electronic structure methods have generally been applied only to small-molecule crystals at present.

Fragment-based methods provide one means of lowering the computational cost of correlated electronic structure methods by decomposing the molecular crystal into monomers, dimers, and many-body contributions.^10,76–80 These methods typically compute only the key monomer and dimer contributions at the most accurate level of theory, while the many-body contributions are approximated in some fashion. Employing coupled cluster singles, doubles, and perturbative triples (CCSD(T)) in the context of a fragment method can give very reliable results, as demonstrated by quantitative prediction of the benzene lattice energy⁸¹ or the prediction of the methanol polymorph phase diagram with ∼0.5 kJ mol⁻¹ accuracy.⁸²

Unfortunately, even with fragment methods, CCSD(T) calculations are cost-prohibitive for crystals involving pharmaceutical-sized species. MP2 is more feasible computationally, but it suffers from well-known problems in the description of van der Waals interactions⁸³ that cause it to substantially overestimate the interaction energy in π-stacking complexes⁸⁴ and to over-bind the benzene crystal by ∼10–20%,¹⁰ for example. This difficulty is overcome here by using the recently developed MP2D model, which employs a Grimme D3-like⁸⁵ dispersion correction that removes the problematic dispersion treatment inherent to MP2 and replaces it with a more reliable treatment. MP2D performs well across extensive benchmark calculations of dimer interactions, molecular conformations, and reaction energies,⁴⁸ and it correctly describes the energetics of the aforementioned anthracene photodimer.⁴⁹

The present study examines three challenging cases of conformational polymorphism in detail: ortho-acetamidobenzamide,⁴⁰ ROY,^42,86 and oxalyl dihydrazide (Fig. 1).^40,87 It demonstrates how problems in the intramolecular conformational energies and, to a lesser extent, the intermolecular interactions with well-regarded dispersion-corrected DFT functionals lead to incorrect polymorph stabilities. However, modeling these systems with fragment-based correlated wavefunction methods overcomes these difficulties, restoring the crucial balance between intra- and intermolecular interactions^88–90 that is required to predict the correct stabilities in conformational polymorphs. The results here highlight how despite considerable progress with DFT, polymorph stability ranking remains challenging, and models that can achieve higher accuracy than that of commonly used DFT approximations are needed before polymorph ranking can be considered a “solved” problem.

Fig. 1 The species whose conformational polymorphs are studied here.

2 Theory and methods

When predicting relative polymorph stabilities, it is important to recognize that those stabilities depend on temperature. Sometimes the free energy variations are large enough to produce an enantiotropic relationship where the thermodynamically preferred polymorph changes depending on the temperature. However, even in monotropic cases where one polymorph is always preferred thermodynamically, the magnitude of the enthalpy and free energy differences between two polymorphs will depend on temperature. The temperature dependence of the relative stabilities arises from both phonon contributions to the vibrational partition function and from phonon-driven thermal expansion. The molar volume of typical organic crystals expands several percent upon heating from 0 K to room temperature.⁹¹ This expansion alters the lattice energy and introduces anharmonicity into the phonons. Accounting for it is important for quantitatively comparing thermochemistry,^82,92–94 mechanical properties,⁹³ and spectroscopic observables^91,95 between theory and experiment.

It is therefore important to consider the thermodynamic conditions under which experimental measurements were made when making theoretical predictions. Ideally, one would capture temperature effects via molecular dynamics (including nuclear quantum effects, since zero-point contributions can be significant^94,96). However, molecular dynamics simulations based on high-level electronic structure methods are very computationally expensive.⁹⁶ The quasi-harmonic approximation is often successfully used to approximate the volume-dependent contributions to the phonons and lattice energies.^{13,23,28,82,91–94,97–100} Nevertheless, quasi-harmonic calculations remain considerably more expensive than purely harmonic calculations that neglect thermal expansion.

Here, a simple approximation is employed to estimate the temperature dependence of the thermochemical stabilities and facilitate comparison with experiment. Two types of crystal structure optimizations are performed. Fully relaxed crystal structures that optimize both the atomic positions and the unit cell vectors approximate the structure at 0 K (albeit without zero-point vibrational expansion⁹⁴). Room-temperature structures are mimicked via fixed-cell optimizations that relax the atomic positions subject to the constraint of the room-temperature experimental lattice parameters. Harmonic phonons are computed separately on each set of structures, thereby approximately capturing the anharmonicity that results from the change in unit cell dimensions.

Similar fixed-cell optimizations have been used by many other authors previously to examine room-temperature crystal properties, including in earlier studies on the same polymorphic systems studied here.^40,42 Constraining the lattice parameters effectively captures the thermal expansion effects and its associated phonon anharmonicity, while relaxing the atomic positions addresses any issues in the experimental molecular geometries (hydrogen atom placement, for example) and ensures the structure is at a minimum for harmonic vibrational frequency calculations. Additional support for this approach comes from the fact that nuclear magnetic resonance chemical shift predictions performed on structures relaxed with fixed lattice parameters reproduce experimental chemical shifts better than those obtained from fully relaxed structures⁹¹ or even neutron diffraction structures.¹⁰¹

Note that the approximations used here neglect thermal/large-amplitude dynamical motions that can occur in molecular crystals. Fortunately, the structures of the systems considered here do not exhibit significant disorder and are likely amenable to static modeling treatments. The differences between quasi-harmonic and molecular dynamics models are frequently (but not always) small.^102–104 In the end, combining information from the fully-relaxed 0 K structures and fixed-cell room-temperature structures provides information regarding the topology of crystal energy landscapes that facilitates comparison with experiment.

Experimental crystal structures were obtained from the Cambridge Structure Database for o-acetamidobenzamide¹⁰⁵ (reference codes ACBNZA and ACBZNA01), ROY^106–108 (QAXMEH–QAXMEH05, QAXMEH12, and QAXMEH52), and oxalyl dihydrazide⁸⁷ (VIPKIO01–VIPKIO05). Crystal structures were optimized using periodic DFT with the B86bPBE density functional^109,110 and exchange-hole dipole moment (XDM) dispersion correction.¹¹¹ This particular combination performs well in many molecular crystal applications.^22–24,111

Single-point refinement of the electronic energies was carried out using correlated wavefunction methods via the fragment-based hybrid many-body interaction (HMBI) model.^78,112–114 HMBI partitions the total energy of the crystal into intramolecular contributions (1-body interactions), pairwise intermolecular interactions (2-body interactions), and the remaining many-body intermolecular lattice contributions. The important 1-body and short-range (SR) 2-body terms are modeled with high-level electronic structure methods (e.g. CCSD(T) or MP2-based methods here), while the longer-range (LR) 2-body and many-body contributions are modeled with periodic Hartree–Fock (HF) theory (which makes it comparable to Stoll's method of increments⁷⁶).

U^HMBI_el = E^High_1-body + E^High_{SR 2-body} + E^HF_{LR 2-body} + E^HF_many-body (1)

As noted above, MP2 suffers from problematic description of van der Waals interactions. The related and highly successful^115,116 MP2C model addresses this problem by adding a non-empirical intermolecular dispersion correction to MP2.^117,118 However, the MP2C correction is derived from intermolecular perturbation theory and does not address problems with intramolecular dispersion. MP2D⁴⁸ expresses the dispersion correction in terms of atom-centered C₆ and C₈ dispersion coefficients¹¹⁹ computed using the scheme behind Grimme's D3 dispersion correction.⁸⁵ The five global empirical parameters in MP2D were determined previously⁴⁸ on small-molecule systems that do not include the species studied here.

In the cases of oxalyl dihydrazide and o-acetamidobenzamide, enthalpies and free energies are computed for comparison with experiment. This requires evaluating harmonic phonon contributions to the enthalpy and Helmholtz vibrational free energy F_vibvia the standard statistical mechanical expressions.⁹⁴ The phonons and their thermodynamic contributions are calculated at the B86bPBE-XDM level, using either the 0 K or room-temperature crystal structures, and these are used to augment the electronic energy computed with either DFT or HMBI. For example, the Gibbs free energy at the MP2D level is estimated as,

G(T,P) = U^HMBI_el + F^DFT_vib + PV (2)

For a crystal at ambient conditions, the PV term contributes negligibly and can be ignored. This combination of DFT geometries and phonons with higher-level single-point electronic energies has been validated previously.¹¹⁶

The DFT calculations were performed using Quantum Espresso v6.3 (ref. 120) using a 50 Ry planewave cutoff and well-converged Monkhorst–Pack k-point sampling grids (ESI Section S1.1†). Core electrons were treated according to the projector augmented wave (PAW) approach using PAW potentials for H, C, N, O, and S produced with A. Dal Corso's Atomic code v6.1.¹²¹ Gas-phase monomer and dimer DFT calculations used in the energy decompositions were performed in large unit cells with a minimum of 15 Å spacing between the central monomer/dimer atoms and all periodic image atoms. Using an even larger 18 Å spacing altered the gas-phase energies by ∼0.15 kJ mol⁻¹ or less, indicating that this vacuum spacing is appropriately large to mimic the gas phase.

Harmonic DFT phonon frequencies for oxalyl dihydrazide and o-acetamidobenzamide were computed at the Γ-point using Phonopy v1.12.6-r66 (ref. 122) with the same B86bPBE-XDM functional and basis set used for the energies and geometry optimizations. To ensure equal numbers of molecules (Z = 4) in the cell for each of the five oxalyl dihydrazide polymorphs, supercells were constructed for the α, β, δ, and ε forms by doubling the cell along the shortest crystallographic axis. Using larger supercells and/or capturing phonon dispersion away from the Γ point would certainly improve the quality of the predicted thermochemistry.¹⁰⁰ Still, earlier quasi-harmonic sublimation enthalpy calculations for several small-molecule crystals agreed with experiment to within a couple kJ mol⁻¹ despite neglecting phonon dispersion.¹¹⁶ Percentage errors in the entropic contributions were considerably larger for the same species, however. Phonons were not computed for the larger ROY system for reasons of computational expense.

For the HMBI fragment calculations employing correlated wave function methods, large basis sets must be used to ensure convergence of the polymorph energetics, as demonstrated in ESI Section S1.2† and many previous studies.^{10,81,82,93,94,116,123,124} Here, MP2 and MP2D monomer and dimer energies at the complete-basis-set (CBS) limit were obtained using a development version of PSI4.¹²⁵ The correlation energy was extrapolated¹²⁶ to the complete-basis-set (CBS) limit using data from the aug-cc-pVTZ and aug-cc-pVQZ basis sets¹²⁷ and combined with HF/aug-cc-pVQZ. The MP2C dispersion corrections were obtained with Molpro 2012.1 (ref. 128) in the aug-cc-pVTZ basis set. The MP2C dispersion correction typically converges faster with basis set than the raw correlation energy.¹¹⁸ CCSD(T) results at the CBS limit were obtained by correcting MP2/CBS energies with the difference between CCSD(T) and MP2 in the aug-cc-pVDZ (ROY, oxalyl dihydrazide) or cc-pVTZ (acetamidobenzamide intramolecular contributions) basis sets. The periodic HF many-body contributions were evaluated using Crystal 17 (ref. 129) and the pob-TZVP-rev2 basis set.¹³⁰ This basis set was chosen based on cluster benchmarks described in ESI Section S1.3.† Note that due to computational expense, large-basis set calculations were performed only on the room-temperature structures of ROY, since those are more directly comparable with experiment. Smaller-basis results on the fully relaxed structures are provided in ESI Section S3.3.† As expected, large basis sets are required to converge the relative polymorph stabilities.

As part of the analysis of the ROY system, a one-dimensional, gas-phase conformational energy scan over the key S–C–N–C dihedral angle was performed. At each of 16 fixed dihedral angle values ranging 0–150° in 10° intervals, all other degrees of freedom were fully relaxed at the B3LYP-D3(BJ)/def2-TZVP level of theory. Single-point energies were then computed on these geometries using B86bPBE-XDM, MP2, MP2D, and CCSD(T).

3 Results and discussion

3.1 o-Acetamidobenzamide

o-Acetamidobenzamide has two conformational polymorphs: the α form adopts the more stable intramolecular conformation containing an intramolecular hydrogen bond between the acetamide hydrogen and the amide oxygen (Fig. 2).¹⁰⁵ The β polymorph sacrifices the intramolecular hydrogen bond to adopt a conformation that allows better intermolecular hydrogen bonding. Experimentally, the α form converts exothermically and irreversibly to β upon heating to 150 °C (423 K), with ΔH_α→β values of −1.9 kJ mol⁻¹ (ref. 40) or −2.9 kJ mol⁻¹ (ref. 105). In other words, the β form is clearly preferred over the α one at high temperatures, though the stability ordering at lower temperatures is unclear.

Fig. 2 Local hydrogen bonding environments for crystalline o-acetamidobenzamide. (a) In the α polymorph, the molecule is nearly planar and adopts an intramolecular hydrogen bond, while (b) in the β polymorph the amide and acetamide side chains rotate out of the plane to achieve better intermolecular hydrogen bonds.

Earlier calculations using force fields and several different GGA and hybrid DFT functionals all predict the α form lattice energy to be ∼5–10 kJ mol⁻¹ more stable than β.⁴⁰ New B86bPBE-XDM DFT lattice energy calculations performed here similarly favor the α form by 5.8 kJ mol⁻¹ (Table 1). The B86bPBE-XDM harmonic zero-point vibrational energy contribution for the fully-relaxed 0 K structures stabilizes the β form by 1.2 kJ mol⁻¹ relative to α, which is reasonably similar to the 2 kJ mol⁻¹ value estimated previously.⁴⁰ In other words, the difference in zero-point energy contributions between polymorphs are much too small to alter the B86bPBE-XDM polymorph stability ordering. Moreover, the DFT enthalpic preference for the α form increases from 4.6 kJ mol⁻¹ at 0 K to 5.5 kJ mol⁻¹ at room temperature (Fig. 3). As shown in ESI Table S4,† this temperature-dependence between the 0 K and room-temperature structures arises from a 0.6 kJ mol⁻¹ relative destabilization of the α form caused by the lattice energy changes that is canceled by a larger 1.6 kJ mol⁻¹ stabilization due to the vibrational enthalpy contribution.

Table 1 Stability of the β o-acetamidobenzamide polymorph relative to the α one, in kJ mol⁻¹. Positive values indicate α is more stable than β. See ESI Table S4 for additional details of the lattice energy and phonon contributions

Method	ΔEintra(0 K)	ΔEinter(0 K)	ΔElattice(0 K)	ΔH(0 K)	ΔH(298 K)	ΔH(423 K)^a	ΔG(298 K)
a Linearly extrapolated to 423 K from ΔH(0 K) and ΔH(298 K) values. b Ref. 40. c Ref. 105.
B86bPBE-XDM	58.0	−52.2	5.8	4.6	5.5	5.9	4.3
MP2/CBS + pHF	50.1	−46.8	3.2	1.9	−1.6	−3.1	−2.8
MP2C/CBS + pHF	50.1	−52.2	−2.1	−3.4	−4.9	−5.5	−6.1
MP2D/CBS + pHF	52.6	−51.2	1.4	0.2	−1.4	−2.0	−2.6
CCSD(T)/CBS	52.3
Experiment						−1.9^b, −2.9^c

---

Fig. 3 Predicted enthalpy difference between the α and β polymorphs of o-acetamidobenzamide at 0 K, room temperature, and linearly extrapolated to 423 K. Experimental values were taken from ref. 40 and 105.

To compare more directly against experiment, the transition enthalpy at the phase transition temperature is estimated here via linear extrapolation of the 0 K and room-temperature results. The resulting endothermic ΔH_α→β(423 K) value of 5.9 kJ mol⁻¹ contradicts the exothermic phase transition observed experimentally. The B86bPBE-XDM Gibbs free energies exhibit a clear preference for the α polymorph of 4.6 kJ mol⁻¹ at 0 K that decreases only slightly to 4.3 kJ mol⁻¹ at room temperature (Table 1). Linear extrapolation of ΔG_α→β to 423 K gives 4.2 kJ mol⁻¹. This strong DFT free energy preference for the α form is inconsistent with the irreversible α → β phase transition seen experimentally. Earlier attempts to rationalize the DFT lattice energy preference for the α polymorph suggested that the two polymorphs might be enantiotropically related.⁴⁰ However, the DFT harmonic free energy calculations here predict the α polymorph to be considerably more stable throughout the temperature range (i.e. monotropically related, see Fig. 4a). Taken together, both the B86bPBE-XDM enthalpies and free energies computed here and those reported previously⁴⁰ are inconsistent with experiment.

Fig. 4 Schematic relative enthalpy H and free energy G curves for the α and β polymorphs of o-acetamidobenzamide computed using (a) B86bPBE-XDM or (b) fragment-based MP2D/CBS + pHF. The MP2D model reverses the stability ordering, giving results that are consistent with experiment.

In contrast, fragment-based MP2D predicts polymorph stabilities that agree very well with experiment (Table 1). Calculations on the fully-relaxed 0 K structures suggest that the α polymorph is still more stable in lattice energy, but by a much smaller 1.4 kJ mol⁻¹. Including the B86bPBE-XDM zero-point vibrational contribution preferentially stabilizes the β form, such that the two forms become nearly degenerate (α is 0.2 kJ mol⁻¹ more stable than β). Heating further stabilizes the β form, with ΔH_α→β = −1.4 kJ mol⁻¹ at 298 K. This temperature dependence is largely driven by a 3.1 kJ mol⁻¹ stabilization of the β form arising from the lattice energies, which is partially canceled by the DFT phonon contribution (ESI Table S4†). The exothermic phase transition at higher temperatures predicted by MP2D is consistent with experiment. Linearly extrapolating ΔH_α→β to 423 K gives −2.0 kJ mol⁻¹, in excellent agreement with both experimental values of −1.9 and −2.9 kJ mol⁻¹ (Fig. 3).

Furthermore, MP2D Gibbs free energies indicate that the β form is thermodynamically preferred at elevated temperatures (Fig. 4b), which is consistent with the irreversible α → β transition seen experimentally at 423 K. Nominally, the MP2D calculations predict an enantiotropic relationship between the two forms, though the 0.2 kJ mol⁻¹ free energy difference between the two forms at 0 K is likely smaller than the inherent uncertainties in the models. The MP2D free energies suggest that the experimentally observed α → β phase transition at 423 K corresponds to a kinetically activated transformation from the metastable α form to the stable β one, rather than a true thermodynamic phase boundary.

To understand why MP2D performs well in this system while B86bPBE-XDM does not, Table 1 decomposes the lattice energy differences between the two polymorphs into their intra- and intermolecular contributions. MP2D and B86bPBE-XDM actually predict similar intermolecular energies that differ by only 0.3 kJ mol⁻¹ for the fully relaxed 0 K structure, and by 1.3 kJ mol⁻¹ for the room-temperature structures. Rather, the erroneous DFT predictions arise almost entirely from the intramolecular conformational energies: B86bPBE-XDM over-stabilizes the intramolecular hydrogen bond conformation by ∼6 kJ mol⁻¹ (11% error) compared to gas-phase CCSD(T) benchmarks. Delocalization error is known to cause GGA and hybrid functionals to overstabilize aromatic systems.^24,131,132 The intramolecular conformation in the β polymorph disrupts not only the intramolecular hydrogen bond, but also π conjugation between the aromatic ring and the amide/acetamide side chains. In contrast to B86bPBE-XDM, MP2D reproduces the CCSD(T) conformational energy difference to within a few tenths of a kJ mol⁻¹ (<1% error).

Finally, the o-acetamidobenzamide polymorphs demonstrate the importance of correcting both the intra- and intermolecular description of dispersion in MP2, as is done in MP2D. MP2 underestimates the intermolecular preference for the β phase by 4–5 kJ mol⁻¹, and it underestimates the penalty for disrupting the intramolecular hydrogen bond found in the α form by 2 kJ mol⁻¹ (Table 1). These errors cancel somewhat, but the resulting ΔH and ΔG values appear to change too rapidly with temperature (due to how the lattice energy varies with the temperature-dependent changes in crystal structure). MP2C corrects the description of the intermolecular interactions, but it does not alter the intramolecular description. This disrupts the fortuitous error cancellation found in MP2, and MP2C overestimates the stability of the β form substantially.

3.2 ROY

ROY is among the most prolific conformational polymorph formers known. Seven polymorphs have been well-characterized for years.^86,106,107 Since 2018, the structures of two more polymorphs, R05 ⁴¹ and PO13,¹⁰⁸ have been solved, and a structure for the RPL polymorph was proposed.⁴² However, as multiple recent studies have noted, predicting the energetics of these polymorphs has proved challenging.^41,42,133 Many well-regarded van der Waals-inclusive GGA and hybrid density functionals predict highly incorrect polymorph orderings, including PBE-D3, PBE-NP, optPBE-vdW, PBE + MBD, and PBE0 + MBD (Fig. 5).^41,42 Most strikingly, the DFT calculations frequently suggest that the Y polymorph is one of the least stable forms, when it is actually the most stable polymorph experimentally. Inclusion of zero-point energies or thermal contributions does not correct the rankings, either.⁴²

Fig. 5 Comparison between predicted lattice energies and experimentally measured enthalpies⁸⁶ for 8 polymorphs of ROY, all relative to form Y. The PBE-D3, optPBE-vdW, PBE + MBD, and PBE0 + MBD results were taken from ref. 41. They omit PO13 and use fully relaxed 0 K unit cells. The other results employ fixed-cell room-temperature structures.

Experimentally, the relative free energies of the ROY polymorphs were measured in the ∼40–120 °C range by eutectic melting experiments.^{106,107,134,135} Relative enthalpies were then obtained from the slopes of ΔG/T vs. 1/T plots, where T is temperature. Because the fitted enthalpies lack explicit temperature dependence, they are most valid in the elevated temperature regime under which they were measured. To facilitate comparison against the experimental data, fixed-cell room-temperature crystal structures are used for the polymorph stability calculations here, since they will mimic the structures under the experimental conditions better than fully relaxed 0 K ones.

B86bPBE-XDM calculations on the room-temperature structures predict stabilities that are similar to earlier DFT studies (Fig. 5), with the Y form being the second least-stable polymorph in terms of lattice energy (below only the YN form). Differences in the polymorph stabilities between the fixed-cell room-temperature and fully-optimized 0 K structures are modest. Refining the lattice energies of the room-temperature structures with single-point energy calculations at the fragment-based MP2D level completely transforms the crystal energy landscape. MP2D correctly predicts the Y form to be the most stable in terms of lattice energy. Furthermore, with the exception of the ON polymorph, the MP2D stability ordering qualitatively matches the experimental enthalpy data perfectly. The predicted lattice energy differences are somewhat larger than the experimental enthalpies. That discrepancy may in part be due to the omission of phonon contributions in the predictions here. The reason for the incorrect ordering of the ON polymorph is unclear. While an earlier study using a different density functional had difficulty reproducing the experimental ON crystal structure, the structure obtained here with B86bPBE-XDM agrees well with experiment (ESI Section S1.1†).

The MP2D calculations also predict the recently discovered PO13 polymorph to be among the less stable polymorphs, below only ON and ORP. Though experimental thermochemical data is not available for the PO13 form, this prediction is consistent with experimental data that indicates PO13 is less stable than the Y polymorph and that its heat of fusion is in the mid-range compared to the other forms. Predictions for the R05 form are omitted here because, unlike the other polymorphs, its cell has a net dipole which creates difficulties for the fragment-based approach (ESI Section S3.4†).

Like for o-acetamidobenzamide, the major differences between MP2D and B86bPBE-XDM for the ROY polymorph energies stem from the problematic B86bPBE-XDM intramolecular treatment of the conformational energies. Fig. 6 plots the one-dimensional conformational energy scan along a key dihedral angle associated with the different conformations found in the ROY polymorphs. For convenience, vertical lines in the figure highlight the corresponding values of this dihedral angle found in the different polymorphs (though the other degrees of freedom along this scan may differ from those found in the actual crystals). Similar conformational energy scans can be found in earlier studies,^42,133 albeit without the CCSD(T) benchmarks provided here. Compared to CCSD(T), B86bPBE-XDM dramatically overstabilizes the conformations found in the red (R) and orange (O) polymorphs relative to those occuring in the yellow (Y) forms. This explains why the earlier DFT calculations predict the yellow forms to be so much less stable than the others. In contrast, MP2D mimics the CCSD(T) conformational energy profile much more faithfully. MP2D does underestimate the stability of the conformations with the dihedral angles adopted by the red and orange forms by up to ∼1 kJ mol⁻¹ relative to CCSD(T). Nevertheless, MP2D represents a substantial improvement over B86bPBE-XDM. The intramolecular MP2D dispersion correction improves the MP2 conformational energies modestly, by up to ∼1 kJ mol⁻¹ (Fig. 6). At the same time, the MP2D dispersion correction does not alter the qualitative MP2 polymorph stability ordering in this system (ESI Section S3.2†).

Fig. 6 Comparison of the intramolecular conformational energy scan for the key intramolecular dihedral angle in ROY at several levels of theory. Vertical dotted lines indicate the experimental dihedral angles for each polymorph. Energies are relative to the 120° conformation, which corresponds to the CCSD(T) minimum.

For further insight, Fig. 7 decomposes the relative polymorph energies into their intra- and intermolecular contributions. As expected from the one-dimensional conformational energy scan, B86bPBE-XDM overstabilizes the intramolecular conformations of the orange and red polymorphs, while MP2 and MP2D give conformational energies which are quite faithful to CCSD(T). The intermolecular trends plotted in Fig. 7b are qualitatively similar between the three models. Indeed, the intermolecular energies are roughly parallel between methods for most of the polymorphs. The most notable difference is that B86bPBE-XDM appears to stabilize the intermolecular interactions of the other forms relative to Y more so than do MP2 and MP2D, which shifts all points other than Y in the B86bPBE-XDM curve down relative to the MP2-based ones in Fig. 7b. Unfortunately, coupled cluster benchmarks that could assess the quality of the two models for the intermolecular interactions are computationally infeasible.

Fig. 7 Energy decomposition of the room-temperature structure ROY polymorph lattice energies into (a) intramolecular, (b) intermolecular, and (c) total energy contributions. The energies are plotted relative to the most stable Y polymorph in each case.

Combining the intra- and intermolecular contributions (Fig. 7c), one sees once again that the MP2 and MP2D curves are in much better agreement with experiment than the B86bPBE-XDM one is. The MP2D energy of the ON polymorph is the most notable outlier relative to experiment. The fact that MP2D predicts the intramolecular conformational energy of the ON polymorph to within 0.2 kJ mol⁻¹ of CCSD(T) suggests that any problem in the predicted energy ranking arises from the intermolecular contributions.

3.3 Oxalyl dihydrazide

The five polymorphs of oxalyl dihydrazide differ in whether the crystal packing contains purely intermolecular hydrogen bonding (α form) or exhibits a mixture of both intra- and intermolecular hydrogen bonds (β, γ, δ, and ε forms), as shown in Fig. 8.⁸⁷ Five additional high-pressure polymorphs have been reported,¹³⁶ but their structures are unknown and they are not considered here. Experimentally, the α, ε, and δ forms are the most stable polymorphs, though the ranking among those three is uncertain. The α form should arguably be the most stable polymorph based on its high density,⁸⁷ but exceptions to density-based stability arguments can occur in hydrogen bonded crystals. All three forms convert endothermically to γ near 200–210 °C. This suggests that the γ form has a lower free energy at these temperatures, but that the α, ε, and δ forms have lower enthalpies.⁸⁷ This indicates an enantiotropic relationship between γ and the other three forms according to the heat-of-transition rule.¹³⁷ Finally, the β form has proved difficult to produce and characterize, and it converts readily to the α form.^87,136 Therefore, it is assumed to be the least stable form. Overall, the inferred lattice energy ranking is (from most to least stable): α, δ, ε < γ < β.

Fig. 8 Crystalline oxalyl dihydrazide exhibits (a) purely intermolecular hydrogen bonding in the α polymorph and (b) a mixture of intra- and intermolecular hydrogen bonding in the other four polymorphs (ε form shown here).

These five polymorphs have been studied theoretically by several groups,^10,40,58,123 and the results have been summarized by a couple authors.^8,10 Initial DFT calculations involving empirical dispersion corrections predicted lattice energies consistent with the aforementioned stability ordering, but the lattice energies spanned a surprisingly large range of ∼15 kJ mol⁻¹.⁴⁰ Subsequent DFT calculations employing more modern dispersion corrections narrowed the polymorph energy range modestly to ∼10–12 kJ mol⁻¹. Prior HMBI fragment-based MP2 and MP2C calculations also achieved the same stability ordering, albeit with much smaller ∼3–4 kJ mol⁻¹ energy window.¹²³ That study found that the competition between intra- and intermolecular basis set superposition error plays a substantial role in the energetics and that large basis sets are needed when atom-centered basis functions are used. Fully periodic local MP2 calculations performed two years later found a ∼10 kJ mol⁻¹ energy range for the polymorphs,⁵⁸ which are consistent with the DFT calculations, though the double-zeta basis set is probably too small to draw firm conclusions.

Fig. 9 presents new B86bPBE-XDM DFT results and HMBI-based MP2, MP2D, and MP2C single-point energy refinements of those DFT structures. For the fully relaxed 0 K structures, all four methods exhibit energy gaps in the ∼10–12 kJ mol⁻¹ energy range, consistent with the earlier DFT and periodic local MP2 calculations. The HMBI results here should be more reliable than the previously published ones,¹²³ since the geometries were optimized with a more robust B86bPBE-XDM dispersion-corrected DFT functional and because the many-body terms are evaluated with periodic HF instead of a polarizable force field. The strong, favorable polarization that is found only in the α polymorph makes the energy gap between the α and other four forms sensitive to the many-body description. The small energy range for the polymorphs predicted in ref. 123 appears to be an erroneous artifact of the polarizable force field contributions used there.

Fig. 9 Relative stabilities of the oxalyl dihydrazide polymorphs at several levels of theory using both the fully relaxed 0 K structures and the fixed-cell room-temperature structures. For the MP2 methods, the relative stabilities of the β and γ forms differ depending on which structure optimization is used.

Like earlier calculations, the relative B86bPBE-XDM lattice energies for the fully relaxed 0 K structures are consistent with the inferred experimental lattice energy stability ordering: α < δ < ε < γ < β. Using the same structures, the MP2-based methods agree that α < δ < ε < γ, though the ε form is somewhat less stable than it is with DFT. However, MP2-based methods predict that the β form is more stable than the γ one, contrary to what has been inferred experimentally. Fragment CCSD(T) calculations similarly predict the β form to be more stable. They also predict the δ form to be marginally (0.1 kJ mol⁻¹) more stable than ε. As noted earlier, the experimental stability ranking among α, δ, and ε is unclear.

Because the limited experimental knowledge of the β form is based on its instability at ambient conditions, the calculations were repeated using the fixed-cell room-temperature crystal structures. Using these structures destabilizes the lattice energy of the β polymorph relative to the α one by 4 kJ mol⁻¹ with B86bPBE-XDM, and by 2–2.5 kJ mol⁻¹ for the MP2-based methods (Fig. 9). Moreover, the MP2 models stabilize the room-temperature γ, δ, and ε structures by 1–2 kJ mol⁻¹ relative to the α one. The end result is that all methods predict the α < δ < ε < γ < β stability ordering when room-temperature structures are used.

The predicted temperature dependence translates directly to the enthalpies and free energies shown in Fig. 10. Augmenting the electronic energies with phonon contributions stabilizes the β, γ, δ, and ε polymorphs relative to α, but it does not alter the predicted stability orderings for the DFT or MP2 models (ESI Table S7 and Fig. S7†).

Fig. 10 Relative enthalpy H and free energy G curves for the five polymorphs of oxalyl dihydrazide computed using (a) B86bPBE-XDM or (b) fragment-based MP2/CBS + pHF. The most notable difference between the two curves is that the β and γ forms are monotropically related in the DFT calculations, but enantiotropically related with MP2D (the G curves intersect near 175 K). The data points indicate computed results, while the curves connecting them are schematic.

The difference between the lattice energies and free energies is interesting. First, the free energies span a much narrower range than the lattice energies—e.g. 6 kJ mol⁻¹versus 15 kJ mol⁻¹ for MP2D at 298 K—indicating the importance of entropic contributions in this system. Second, the free energy of the γ form stabilizes more rapidly with increasing temperature than do those of the α, δ, and ε forms. That is consistent with the experimental evidence for the γ form being enantiotropically related to the other three polymorphs.

Overall, the DFT and MP2-based results here are all consistent with available experimental observations. More experimental data would be needed to discriminate between the distinct DFT and MP2-based predictions for the β and γ polymorph stabilities. However, two points can be made based on currently available information. First, the experimental crystal structure of β has greater uncertainty than the other forms.⁸⁷ That could affect the β fixed-cell room-temperature structure moreso than the fully relaxed one. Notably, the β form contracts ∼8% upon full geometry optimization, compared to only ∼4% for the other four polymorphs. This bigger structural change between 0 K and room-temperature manifests in the correspondingly large lattice energy change seen for the β form. Of course, this larger structural change in the β form might also simply reflect differences in the crystal packing that increase its thermal expansivity.

Second, energy decomposition and CCSD(T) benchmarks on the fully-relaxed structures in Fig. 11 suggest that the MP2D energetics are more accurate than those from B86bPBE-XDM. MP2D and CCSD(T) both predict the intramolecular conformations found in the β and γ polymorphs to be considerably more stable than those adopted in the δ and ε polymorphs. In contrast, B86bPBE-XDM predicts erroneously small energy differences between the conformations. For the total pairwise intermolecular interactions, B86bPBE-XDM and MP2D predict fairly similar energetics, with the DFT results actually agreeing better with CCSD(T). This is due to fortuitous error cancellation: examining all the individual dimer interactions, MP2D exhibits a root-mean-square (rms) error of 0.4 kJ mol⁻¹versus CCSD(T), compared to 0.8 kJ mol⁻¹ for B86bPBE-XDM. Similarly, looking at total two-body contributions instead of relative lattice energies, the B86bPBE-XDM 2-body contributions are systematically under-bound by rms error 7.3 kJ mol⁻¹ compared to CCSD(T). In contrast, MP2D overbinds by a much smaller rms 1.6 kJ mol⁻¹. In other words, the B86bPBE-XDM 2-body terms have much larger errors, but they exhibit better systematic error cancellation here to produce the agreement seen in Fig. 11b.

Fig. 11 Energy decomposition for the oxalyl dihydrazide polymorph stabilities into (a) intramolecular (1-body), (b) pairwise intermolecular (2-body), and (c) many-body intermolecular contributions. Note, in the HMBI fragment approach, both MP2D and CCSD(T) employ the same HF many-body treatment. The same 14 kJ mol⁻¹ energy range is used in all three figures to facilitate comparisons.

The other major difference between models occurs for the many-body contributions. The many-body contributions in all three systems explored in this paper generally amount to only ∼5% of the total lattice energy (and always less than 10%), which is typical for organic molecular crystals.¹¹⁴ However, because the many-body contributions in a system like oxalyl dihydrazide vary considerably between polymorphs, they play an out-sized role in determining the relative lattice energies. As shown in Fig. 11c, the relative many-body contributions in oxalyl dihydrazide have the same magnitude as the overall lattice energy differences. Compared to B86bPBE-XDM, periodic HF predicts a considerably stronger polarization effect for the α form that effectively shifts the relative energies of the other polymorphs up. Furthermore, HF predicts the β and γ forms to have more repulsive many-body contributions than does B86bPBE-XDM. While experimental or higher-level theoretical benchmarks are not available to determine which many-body treatment is more accurate, it is clear that B86bPBE-XDM obtains the correct experimental stability ordering only by error cancellation between the intramolecular and many-body intermolecular contributions. For example, if one corrected the erroneous intramolecular B86bPBE-XDM conformational energies with CCSD(T) ones, the result would incorrectly predict that the δ form is the least stable polymorph by 2 kJ mol⁻¹. Obtaining the proper stability ordering would also require replacing the B86bPBE-XDM many-body energies with values similar to those obtained from periodic HF. This suggests that the HF many-body contributions are likely closer to the true values.

Finally, as shown in ESI Table S8,† the MP2-based methods all predict the 1-body energies with similar accuracy. However, for the 2-body interaction energies, MP2 and MP2C systematically overbind the dimers with rms errors of 8.1 kJ mol⁻¹ and 4.2 kJ mol⁻¹, respectively, both of which are several-fold larger than the 1.6 kJ mol⁻¹ error for MP2D (also mildly over-bound). In other words, the similar relative polymorph stabilities seen for the different MP2-based methods in Fig. 9 arise from systematic cancellation of the overbinding errors that occur in MP2 and MP2C.

4 Conclusions

For many years, it has been widely recognized that balancing intra- and intermolecular conformational energies is one of the primary obstacles to crystal structure prediction in conformational polymorphs. Since the widespread adoption of DFT in crystal structure prediction, however, this issue of balancing the intra- and intermolecular interactions has been given much less attention. While periodic DFT models may have reduced the prevalence of such balance issues compared to earlier force field studies, the results here clearly demonstrate that widely used density functionals have not yet solved the problem of ranking conformational polymorphs in crystal structure prediction.

This study examined three well-known and challenging examples of conformational polymorphism: o-acetamidobenzamide, ROY, and oxalyl dihydrazide. In the first two systems, a variety of dispersion-corrected DFT models predict catastrophically wrong relative polymorph stabilities. The ∼8 kJ mol⁻¹ DFT errors in relative polymorph stabilities found in these two systems greatly exceed the few kJ mol⁻¹ error often associated with DFT polymorph rankings. More importantly, given that over 95% of polymorph pairs exhibit energy differences less than 8 kJ mol⁻¹, and the majority have energy differences less than 2 kJ mol⁻¹, such errors are unacceptable in crystal structure prediction. In both systems, the problem for DFT arises largely from a poor description of the intramolecular conformational energy.

The third system, oxalyl dihydrazide, is more nuanced, in part due to greater ambiguity in the experimental data. High-quality experimental thermochemical measurements of this (and other) polymorphic systems would be valuable for assessing the performance of different models more clearly. Based on presently available data, the overall DFT energy rankings for the oxalyl dihydrazide polymorphs do appear generally consistent with experiment. However, energy decomposition and coupled cluster theory benchmarks suggest that such nominal agreement arises from substantial cancellation of errors between the intra- and intermolecular interactions.

While the present study focused on the B86bPBE-XDM functional, it is important to recognize that these issues transcend any individual density functional. They occur for both GGA and hybrid density functionals with a variety of dispersion treatments, as evidenced by Fig. 5. Furthermore, the conformational polymorph examples here raise the question: How common are such failures of DFT for crystal structure prediction? Are these three systems outliers? Will more examples be uncovered as DFT-based crystal structure prediction techniques are increasingly applied to larger, more flexible pharmaceutical molecules?

This work also demonstrates that fragment-based MP2D calculations provide a promising path forward. For ROY and o-acetamidobenzamide, the higher-level calculations not only correct the qualitative polymorph stability ordering, but they also predict relative stabilities that are in generally good agreement with experiment. For oxalyl dihydrazide, MP2D also appears to perform better than DFT based on the energy decomposition results, though more experimental data to confirm the predictions would be helpful.

The study also examined a few variants of MP2. The balanced description of both intra- and intermolecular dispersion makes MP2D superior to MP2C, which corrects intermolecular dispersion only, and to MP2, which performs poorly for both intra- and intermolecular dispersion. The intramolecular dispersion correction contributes in all three systems, but it proves critical in predicting the stability for the o-acetamidobenzamide polymorphs. Fortunately, the MP2D dispersion correction can be computed with trivial effort once MP2 results are available.

The overall computational cost of these calculations is considerably higher than DFT, with complete-basis-set MP2D single-point energies requiring approximately 900, 7000, and 21 000 central processing unit (CPU) hours per polymorph of oxalyl dihydrazide (C₂H₆N₄O₂), acetamidobenzamide (C₉H₁₀N₂O₂), and ROY (C₁₂H₉N₃O₂S). This computational cost is dominated by the evaluation of ∼40–50 dimer interactions at the large-basis MP2D limit and the N⁵ MP2 scaling with monomer/dimer size.¹⁰ Fortunately, the number of monomer and dimer fragments scales linearly with increasing number of molecules in the asymmetric unit.^79,113 Furthermore, those fragment calculations can be run independently and in parallel, making it feasible to perform crystalline calculations in much shorter amounts of wall time in modern high performance computing environments. The present study demonstrates that such correlated wavefunction calculations are feasible in species that are comparable in size to small-molecule pharmaceuticals. Efforts are currently underway to improve the accuracy and reduce the computational costs of these correlated wavefunction models further.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Funding for this work from the National Science Foundation (CHE-1665212 to G. B. and Graduate Research Fellowship Grant 1326120 to J. M.), and supercomputer time from XSEDE (TG-CHE110064 to G. B.) are gratefully acknowledged. Additional computational resources were provided by the UC Riverside high-performance computing center, which was funded by grants from the National Science Foundation (MRI-1429826) and National Institutes of Health (1S10OD016290-01A1). We thank Dr Charlene Tsay for helpful discussions.

Notes and references

1. S. R. Chemburkar, J. Bauer, K. Deming, H. Spiwek, K. Patel, J. Morris, R. Henry, S. Spanton, W. Dziki, W. Porter, J. Quick, P. Bauer, J. Donaubauer, B. A. Narayanan, M. Soldani, D. Riley and K. Mcfarland, Org. Process Res. Dev., 2000, 4, 413–417.
2. J. Bauer, S. Spanton, R. Henry, J. Quick, W. Dziki, W. Porter and J. Morris, Pharm. Res., 2001, 18, 859–866.
3. I. B. Rietveld and R. Céolin, J. Pharm. Sci., 2015, 104, 4117–4122.
4. M. A. Neumann and J. van de Streek, Faraday Discuss., 2018, 211, 441–458.
5. D. P. Konakanchi, B. Gongalla, K. B. Sikha, C. Kandaswamy, K. Sataya, B. R. Adibhatala and V. C. Nannapaneni, US Pat., 8,877,932 B2, 2014.
6. J. Maddox, Nature, 1988, 335, 201.
7. J. Nyman and G. M. Day, CrystEngComm, 2015, 17, 5154–5165.
8. A. J. Cruz-Cabeza, S. M. Reutzel-Edens and J. Bernstein, Chem. Soc. Rev., 2015, 44, 8619–8635.
9. A. Burger and R. Ramberger, Mikrochim. Acta, 1979, 72, 273–316.
10. G. J. O. Beran, Chem. Rev., 2016, 116, 5567–5613.
11. S. Grimme, A. Hansen, J. G. Brandenburg and C. Bannwarth, Chem. Rev., 2016, 116, 5105–5154.
12. J. Hermann, R. A. DiStasio and A. Tkatchenko, Chem. Rev., 2017, 117, 4714–4758.
13. J. Hoja, A. M. Reilly and A. Tkatchenko, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2017, 7, e1294.
14. G. M. Day, T. G. Cooper, A. J. Cruz-Cabeza, K. E. Hejczyk, H. L. Ammon, S. X. M. Boerrigter, J. S. Tan, R. G. Della Valle, E. Venuti, J. Jose, S. R. Gadre, G. R. Desiraju, T. S. Thakur, B. P. van Eijck, J. C. Facelli, V. E. Bazterra, M. B. Ferraro, D. W. M. Hofmann, M. A. Neumann, F. J. J. Leusen, J. Kendrick, S. L. Price, A. J. Misquitta, P. G. Karamertzanis, G. W. A. Welch, H. A. Scheraga, Y. A. Arnautova, M. U. Schmidt, J. van de Streek, A. K. Wolf and B. Schweizer, Acta Crystallogr., Sect. B: Struct. Sci., 2009, 65, 107–125.
15. M. A. Neumann, F. J. J. Leusen and J. Kendrick, Angew. Chem., Int. Ed., 2008, 47, 2427–2430.
16. D. A. Bardwell, C. S. Adjiman, Y. A. Arnautova, E. Bartashevich, S. X. M. Boerrigter, D. E. Braun, A. J. Cruz-Cabeza, G. M. Day, R. G. Della Valle, G. R. Desiraju, B. P. van Eijck, J. C. Facelli, M. B. Ferraro, D. Grillo, M. Habgood, D. W. M. Hofmann, F. Hofmann, K. V. J. Jose, P. G. Karamertzanis, A. V. Kazantsev, J. Kendrick, L. N. Kuleshova, F. J. J. Leusen, A. V. Maleev, A. J. Misquitta, S. Mohamed, R. J. Needs, M. A. Neumann, D. Nikylov, A. M. Orendt, R. Pal, C. C. Pantelides, C. J. Pickard, L. S. Price, S. L. Price, H. A. Scheraga, J. van de Streek, T. S. Thakur, S. Tiwari, E. Venuti and I. K. Zhitkov, Acta Crystallogr., Sect. B: Struct. Sci., 2011, 67, 535–551.
17. A. M. Reilly, R. I. Cooper, C. S. Adjiman, S. Bhattacharya, A. D. Boese, J. G. Brandenburg, P. J. Bygrave, R. Bylsma, J. E. Campbell, R. Car, D. H. Case, R. Chadha, J. C. Cole, K. Cosburn, H. M. Cuppen, F. Curtis, G. M. Day, R. A. DiStasio Jr, A. Dzyabchenko, B. P. van Eijck, D. M. Elking, J. A. van den Ende, J. C. Facelli, M. B. Ferraro, L. Fusti-Molnar, C.-A. Gatsiou, T. S. Gee, R. de Gelder, L. M. Ghiringhelli, H. Goto, S. Grimme, R. Guo, D. W. M. Hofmann, J. Hoja, R. K. Hylton, L. Iuzzolino, W. Jankiewicz, D. T. de Jong, J. Kendrick, N. J. J. de Klerk, H.-Y. Ko, L. N. Kuleshova, X. Li, S. Lohani, F. J. J. Leusen, A. M. Lund, J. Lv, Y. Ma, N. Marom, A. E. Masunov, P. McCabe, D. P. McMahon, H. Meekes, M. P. Metz, A. J. Misquitta, S. Mohamed, B. Monserrat, R. J. Needs, M. A. Neumann, J. Nyman, S. Obata, H. Oberhofer, A. R. Oganov, A. M. Orendt, G. I. Pagola, C. C. Pantelides, C. J. Pickard, R. Podeszwa, L. S. Price, S. L. Price, A. Pulido, M. G. Read, K. Reuter, E. Schneider, C. Schober, G. P. Shields, P. Singh, I. J. Sugden, K. Szalewicz, C. R. Taylor, A. Tkatchenko, M. E. Tuckerman, F. Vacarro, M. Vasileiadis, A. Vazquez-Mayagoitia, L. Vogt, Y. Wang, R. E. Watson, G. A. de Wijs, J. Yang, Q. Zhu and C. R. Groom, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 439–459.
18. J. Kendrick, M. D. Gourlay, M. A. Neumann and F. J. J. Leusen, CrystEngComm, 2009, 11, 2391.
19. H. C. S. Chan, J. Kendrick, M. A. Neumann and F. J. J. Leusen, CrystEngComm, 2013, 15, 3799.
20. N. Marom, R. A. DiStasio, V. Atalla, S. Levchenko, A. M. Reilly, J. R. Chelikowsky, L. Leiserowitz and A. Tkatchenko, Angew. Chem., Int. Ed., 2013, 52, 6629–6632.
21. A. M. Reilly and A. Tkatchenko, Phys. Rev. Lett., 2014, 113, 055701.
22. A. Otero-De-La-Roza, B. H. Cao, I. K. Price, J. E. Hein and E. R. Johnson, Angew. Chem., Int. Ed., 2014, 53, 7879–7882.
23. S. R. Whittleton, A. Otero-de-la Roza and E. R. Johnson, J. Chem. Theory Comput., 2017, 13, 441–450.
24. S. R. Whittleton, A. Otero-de-la Roza and E. R. Johnson, J. Chem. Theory Comput., 2017, 13, 5332–5342.
25. A. G. Shtukenberg, Q. Zhu, D. J. Carter, L. Vogt, J. Hoja, E. Schneider, H. Song, B. Pokroy, I. Polishchuk, A. Tkatchenko, A. R. Oganov, A. L. Rohl, M. E. Tuckerman and B. Kahr, Chem. Sci., 2017, 8, 4926–4940.
26. J. Yang, C. T. Hu, X. Zhu, Q. Zhu, M. D. Ward and B. Kahr, Angew. Chem., Int. Ed., 2017, 56, 10165–10169.
27. P. Zhang, G. P. F. Wood, J. Ma, M. Yang, Y. Liu, G. Sun, Y. A. Jiang, B. C. Hancock and S. Wen, Cryst. Growth Des., 2018, 18, 6891–6900.
28. J. Hoja, H.-Y. Ko, M. A. Neumann, R. Car, R. A. DiStasio and A. Tkatchenko, Sci. Adv., 2019, 5, eaau3338.
29. E. Schur, J. Bernstein, L. S. Price, R. Guo, S. L. Price, S. H. Lapidus and P. W. Stephens, Cryst. Growth Des., 2019, 19, 4884–4893.
30. L. M. LeBlanc and E. R. Johnson, CrystEngComm, 2019, 21, 5995–6009.
31. M.-A. Perrin, M. A. Neumann, H. Elmaleh and L. Zaske, Chem. Commun., 2009, 3181–3183.
32. J. Kendrick, G. A. Stephenson, M. A. Neumann and F. J. J. Leusen, Cryst. Growth Des., 2013, 13, 581–589.
33. D. E. Braun, J. A. McMahon, L. H. Koztecki, S. L. Price and S. M. Reutzel-Edens, Cryst. Growth Des., 2014, 14, 2056–2072.
34. M. A. Neumann, J. van de Streek, F. P. A. Fabbiani, P. Hidber and O. Grassmann, Nat. Commun., 2015, 6, 7793.
35. D. E. Braun, S. R. Lingireddy, M. D. Beidelschies, R. Guo, P. Müller, S. L. Price and S. M. Reutzel-Edens, Cryst. Growth Des., 2017, 17, 5349–5365.
36. G. R. Woollam, M. A. Neumann, T. Wagner and R. J. Davey, Faraday Discuss., 2018, 211, 209–234.
37. D. E. Braun, J. A. McMahon, R. M. Bhardwaj, J. Nyman, M. A. Neumann, J. Van De Streek and S. M. Reutzel-Edens, Cryst. Growth Des., 2019, 19, 2947–2962.
38. R. M. Bhardwaj, J. A. McMahon, J. Nyman, L. S. Price, S. Konar, I. D. H. Oswald, C. R. Pulham, S. L. Price and S. M. Reutzel-Edens, J. Am. Chem. Soc., 2019, 141, 13887–13897.
39. M. Mortazavi, J. Hoja, L. Aerts, L. Quéré, J. van de Streek, M. A. Neumann and A. Tkatchenko, Commun. Chem., 2019, 2, 80.
40. P. G. Karamertzanis, G. M. Day, G. W. A. Welch, J. Kendrick, F. J. J. Leusen, M. A. Neumann and S. L. Price, J. Chem. Phys., 2008, 128, 244708.
41. M. Tan, A. G. Shtukenberg, S. Zhu, W. Xu, E. Dooryhee, S. M. Nichols, M. D. Ward, B. Kahr and Q. Zhu, Faraday Discuss., 2018, 211, 477–491.
42. J. Nyman, L. Yu and S. M. Reutzel-Edens, CrystEngComm, 2019, 21, 2080–2088.
43. P. J. Bygrave, D. H. Case and G. M. Day, Faraday Discuss., 2014, 170, 41–57.
44. M. J. Bryant, S. N. Black, H. Blade, R. Docherty, A. G. Maloney and S. C. Taylor, J. Pharm. Sci., 2019, 108, 1655–1662.
45. L. M. LeBlanc, S. G. Dale, C. R. Taylor, A. D. Becke, G. M. Day and E. R. Johnson, Angew. Chem., Int. Ed., 2018, 57, 14906–14910.
46. A. Otero-de-la Roza, L. M. LeBlanc and E. R. Johnson, J. Chem. Theory Comput., 2019, 15, 4933–4944.
47. S. Grimme, C. Diedrich and M. Korth, Angew. Chem., Int. Ed., 2006, 45, 625–629.
48. J. Řezáč, C. Greenwell and G. J. O. Beran, J. Chem. Theory Comput., 2018, 14, 4711–4721.
49. G. J. O. Beran, CrystEngComm, 2019, 21, 758–764.
50. R. O. Al-Kaysi, A. M. Müller and C. J. Bardeen, J. Am. Chem. Soc., 2006, 128, 15938–15939.
51. L. Zhu, R. O. Al-Kaysi and C. J. Bardeen, J. Am. Chem. Soc., 2011, 133, 12569–12575.
52. A. M. Reilly and A. Tkatchenko, J. Chem. Phys., 2013, 139, 024705.
53. O. A. Loboda, G. A. Dolgonos and A. D. Boese, J. Chem. Phys., 2018, 149, 124104.
54. L. Maschio, D. Usvyat, F. R. Manby, S. Casassa, C. Pisani and M. Schütz, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 075101.
55. D. Usvyat, L. Maschio, F. R. Manby, S. Casassa, C. Pisani and M. Schütz, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 075102.
56. L. Maschio, D. Usvyat and B. Civalleri, CrystEngComm, 2010, 12, 2429–2435.
57. L. Maschio, J. Chem. Theory Comput., 2011, 7, 2818–2830.
58. D. Presti, A. Pedone, M. C. Menziani, B. Civalleri and L. Maschio, CrystEngComm, 2014, 16, 102–109.
59. M. Marsman, A. Grüneis, J. Paier and G. Kresse, J. Chem. Phys., 2009, 130, 184103.
60. A. Grüneis, M. Marsman and G. Kresse, J. Chem. Phys., 2010, 133, 074107.
61. S. Hirata and T. Shimazaki, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 1–7.
62. T. Shiozaki and S. Hirata, J. Chem. Phys., 2010, 132, 151101.
63. M. Katouda and S. Nagase, J. Chem. Phys., 2010, 133, 184103.
64. Y.-Y. Ohnishi and S. Hirata, J. Chem. Phys., 2010, 133, 034106.
65. D. Usvyat, J. Chem. Phys., 2013, 139, 194101.
66. M. Del Ben, J. Hutter and J. VandeVondele, J. Chem. Theory Comput., 2012, 8, 4177–4188.
67. M. Del Ben, J. Hutter and J. VandeVondele, J. Chem. Phys., 2015, 143, 102803.
68. M. Del Ben, J. VandeVondele and B. Slater, J. Phys. Chem. Lett., 2014, 5, 4122–4128.
69. D. Lu, Y. Li, D. Rocca and G. Galli, Phys. Rev. Lett., 2009, 102, 1–4.
70. M. Macher, J. Klimeš, C. Franchini and G. Kresse, J. Chem. Phys., 2014, 140, 084502.
71. J. Klimeš, J. Chem. Phys., 2016, 145, 094506.
72. K. Hongo, M. A. Watson, R. S. Sanchez-Carrera, T. Iitaka and A. Aspuru-Guzik, J. Phys. Chem. Lett., 2010, 1, 1789–1794.
73. K. Hongo, M. a. Watson, T. Iitaka, A. Aspuru-Guzik and R. Maezono, J. Chem. Theory Comput., 2015, 11, 907–917.
74. G. H. Booth, A. Grüneis, G. Kresse and A. Alavi, Nature, 2013, 493, 365–370.
75. A. Zen, J. G. Brandenburg, J. Klimeš, A. Tkatchenko, D. Alfè and A. Michaelides, Proc. Natl. Acad. Sci. U. S. A., 2018, 201715434.
76. B. Paulus, Phys. Rep., 2006, 428, 1–52.
77. S. Wen, K. Nanda, Y. Huang and G. J. O. Beran, Phys. Chem. Chem. Phys., 2012, 14, 7578–7590.
78. G. J. O. Beran, S. Wen, K. Nanda, Y. Huang and Y. Heit, Top. Curr. Chem., 2014, 345, 59–93.
79. G. J. O. Beran, J. D. Hartman and Y. N. Heit, Acc. Chem. Res., 2016, 49, 2501–2508.
80. S. Hirata, K. Gilliard, X. He, J. Li and O. Sode, Acc. Chem. Res., 2014, 47, 2721–2730.
81. J. Yang, W. Hu, D. Usvyat, D. Matthews, M. Schutz and G. K.-L. Chan, Science, 2014, 345, 640–643.
82. C. Červinka and G. J. O. Beran, Chem. Sci., 2018, 9, 4622–4629.
83. K. E. Riley, M. Pitonak, P. Jurecka and P. Hobza, Chem. Rev., 2010, 110, 5023–5063.
84. M. O. Sinnokrot and C. D. Sherrill, J. Phys. Chem. A, 2006, 110, 10656–10668.
85. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
86. L. Yu, Acc. Chem. Res., 2010, 43, 1257–1266.
87. S. Ahn, F. Guo, B. M. Kariuki and K. D. M. Harris, J. Am. Chem. Soc., 2006, 128, 8441–8452.
88. W. T. M. Mooij, B. P. van Eijck and J. Kroon, J. Am. Chem. Soc., 2000, 122, 3500–3505.
89. B. P. van Eijck, W. T. M. Mooij and J. Kroon, J. Comput. Chem., 2001, 22, 805–815.
90. S. L. Price, Int. Rev. Phys. Chem., 2008, 27, 541–568.
91. J. L. McKinley and G. J. O. Beran, J. Chem. Theory Comput., 2019, 15, 5259–5274.
92. J. Nyman and G. M. Day, Phys. Chem. Chem. Phys., 2016, 18, 31132–31143.
93. Y. N. Heit, K. D. Nanda and G. J. O. Beran, Chem. Sci., 2016, 7, 246–255.
94. Y. N. Heit and G. J. O. Beran, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 514–529.
95. W. Sontising, Y. N. Heit, J. L. McKinley and G. J. O. Beran, Chem. Sci., 2017, 8, 7374–7382.
96. H. Y. Ko, R. A. Distasio, B. Santra and R. Car, Phys. Rev. Mater., 2018, 2, 1–7.
97. C. Červinka, M. Fulem, R. P. Stoffel and R. Dronskowski, J. Phys. Chem. A, 2016, 120, 2022–2034.
98. C. Červinka and M. Fulem, J. Chem. Theory Comput., 2017, 13, 2840–2850.
99. C. Červinka and M. Fulem, Phys. Chem. Chem. Phys., 2019, 21, 18501–18515.
100. J. Hoja and A. Tkatchenko, Faraday Discuss., 2018, 211, 253–274.
101. J. K. Harper, R. Iuliucci, M. Gruber and K. Kalakewich, CrystEngComm, 2013, 15, 8693.
102. E. C. Dybeck, N. S. Abraham, N. P. Schieber and M. R. Shirts, Cryst. Growth Des., 2017, 17, 1775–1787.
103. E. C. Dybeck, D. P. McMahon, G. M. Day and M. R. Shirts, Cryst. Growth Des., 2019, 19, 5568–5580.
104. N. S. Abraham and M. R. Shirts, Cryst. Growth Des., 2019, 19, 6911–6924.
105. L. A. Errede, M. C. Etter, R. C. Williams and S. M. Darnauer, J. Chem. Soc., Perkin Trans. 2, 1981,(2), 233–238.
106. L. Yu, G. A. Stephenson, C. A. Mitchell, C. A. Bunnell, S. V. Snorek, J. J. Bowyer, T. B. Borchardt, J. G. Stowell and S. R. Byrn, J. Am. Chem. Soc., 2000, 122, 585–591.
107. S. Chen, I. A. Guzei and L. Yu, J. Am. Chem. Soc., 2005, 127, 9881–9885.
108. K. S. Gushurst, J. Nyman and S. X. M. Boerrigter, CrystEngComm, 2019, 21, 1363–1368.
109. A. D. Becke, J. Chem. Phys., 1986, 38, 7184–7187.
110. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
111. A. Otero-de-la Roza and E. R. Johnson, J. Chem. Phys., 2012, 136, 174109.
112. G. J. O. Beran, J. Chem. Phys., 2009, 130, 164115.
113. G. J. O. Beran and K. Nanda, J. Phys. Chem. Lett., 2010, 1, 3480–3487.
114. S. Wen and G. J. O. Beran, J. Chem. Theory Comput., 2011, 7, 3733–3742.
115. L. A. Burns, M. S. Marshall and C. D. Sherrill, J. Chem. Phys., 2014, 141, 234111.
116. J. L. McKinley and G. J. O. Beran, Faraday Discuss., 2018, 211, 181–207.
117. A. Hesselmann, J. Chem. Phys., 2008, 128, 144112.
118. M. Pitonak and A. Hesselmann, J. Chem. Theory Comput., 2010, 6, 168–178.
119. A. Tkatchenko, R. A. DiStasio, M. Head-Gordon and M. Scheffler, J. Chem. Phys., 2009, 131, 094106.
120. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
121. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
122. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
123. S. Wen and G. J. O. Beran, J. Chem. Theory Comput., 2012, 8, 2698–2705.
124. C. Červinka and G. J. O. Beran, Phys. Chem. Chem. Phys., 2017, 19, 29940–29953.
125. R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, J. F. Gonthier, A. M. James, H. R. McAlexander, A. Kumar, M. Saitow, X. Wang, B. P. Pritchard, P. Verma, H. F. Schaefer, K. Patkowski, R. A. King, E. F. Valeev, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, J. Chem. Theory Comput., 2017, 13, 3185–3197.
126. T. Helgaker, W. Klopper, H. Koch and J. Noga, J. Chem. Phys., 1997, 106, 9639–9646.
127. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
128. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O'Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and M. Wang, MOLPRO, version 2012.1, a package of ab initio programs, see http://www.molpro.net.
129. R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro and B. Kirtman, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1360.
130. D. Vilela Oliveira, J. Laun, M. F. Peintinger and T. Bredow, J. Comput. Chem., 2019, 40, 2364–2376.
131. H. L. Woodcock, H. F. Schaefer and P. R. Schreiner, J. Phys. Chem. A, 2002, 106, 11923–11931.
132. T. Heaton-Burgess and W. Yang, J. Chem. Phys., 2010, 132, 234113.
133. S. P. Thomas and M. A. Spackman, Aust. J. Chem., 2018, 71, 279.
134. S. Chen, H. Xi and L. Yu, J. Am. Chem. Soc., 2005, 127, 17439–17444.
135. L. Yu, J. Huang and K. J. Jones, J. Phys. Chem. B, 2005, 109, 19915–19922.
136. X. Tan, K. Wang, T. Yan, X. Li, J. Liu, K. Yang, B. Liu, G. Zou and B. Zou, J. Phys. Chem. C, 2015, 119, 10178–10188.
137. A. Burger and R. Ramberger, Mikrochim. Acta, 1979, 72, 259–271.

---

Footnote

† Electronic supplementary information (ESI) available: Additional computational details, tables of energies, analysis of model convergence, & optimized crystal structures. See DOI: 10.1039/c9sc05689k

---