Correcting π-delocalisation errors in conformational energies using density-corrected DFT, with application to crystal polymorphs

Abstract

We consider several molecules characterised by π-electron conjugation whose extent changes along a flexible torsional coordinate, and which represent the monomer units of polymorphic molecular crystals. Delocalisation error in density functional theory (DFT) adversely impacts conformational energies in these species, overstabilising the conformation that maximises conjugation length and leading to incorrect relative energies for the corresponding crystal polymorphs. We demonstrate that density-corrected (DC-)DFT, in which a DFT exchange-correlation functional is evaluated using a Hartree–Fock density, significantly reduces these conformational energy errors. When DC-DFT monomer energies are used as a low-cost intramolecular correction to a periodic DFT calculation of the molecular crystal, the resulting relative polymorph energies are within 1 kJ/mol of benchmark results.

GRAPHICAL ABSTRACT

Keywords:

• Crystal structure prediction
• density functional theory
• self-interaction error
• delocalisation error
• polymorphism

1. Introduction

Crystal packing plays an important but difficult-to-predict role in controlling physical properties of organic molecular crystals, including solubility, shelf stability, and bioavailability [1–3]. Polymorphism [1–4], or the potential to form multiple distinct crystal packing motifs, is critically important to the pharmaceutical industry [3,5], as highlighted by the case of the HIV drug ritonavir, where the appearance of a previously unknown and low-solubility polymorph made it impossible to produce the original polymorph and forced a drug recall [6]. Other pharmaceutical examples of such ‘disappearing polymorphs’ are known [7–9], and more recently it has been suggested that the thermodynamically stable crystal form has not been found experimentally for 15–45% of small-molecule drugs currently available on the market [10].

Over the past decade, ab initio crystal structure prediction has moved from the category of ‘impossible’, as it was famously characterised 30 years ago [11,12], to merely ‘extremely challenging’ [13–23]. Successful crystal structure prediction, including quantitative ranking of low-energy polymorphs, is most feasible for small-molecule species with limited conformational flexibility [24–26]. Quantitative prediction of thermodynamically stable polymorphs comes with stringent accuracy requirements, as an estimated 50% of polymorph pairs differ by less than 2 kJ/mol and around 95% of the conformers differ by less than 8 kJ/mol [4,27,28].

Dispersion-corrected density functional theory (DFT+D) had made enormous strides towards making crystal structures predictable from first principles [22,29–34], but this success is mitigated in certain cases by overstabilisation of delocalised π-electron systems [26,28,35,36], which is a consequence of self-interaction error (SIE) [37]. This issue is thought to impact polymorph prediction in well-known examples such as the anti-cancer drug axitinib [26], the organic semiconductor molecule rubrene [38], the well-studied ROY molecule [26,28,35,39,40], molecule X [26,41] from the third blind test of crystal structure prediction [42], as well as other systems [26,36,43]. Structures for several of these molecules are shown in Figure 1. Each of these molecules exhibits conformational polymorphism whereby changes in intramolecular conformation access different intermolecular crystal packing motifs. Notably, in each of these molecules, the extent of π-electron delocalisation changes as a function of one or more intramolecular torsional coordinates. To the extent that delocalisation error may overstabilise the more highly-conjugated conformations, we expect that DFT will produce erroneous conformational energy profiles for these flexible monomers. Those errors propagate into errors in the relative energies for the crystal polymorphs formed from these monomers [26,28,35,39]. Given the small energy differences between polymorphs, even modest conformational energy errors can have a qualitative impact on polymorph stability ordering.

Figure 1. Structures of the four molecules investigated in this work, indicating the torsional angles that modify the conjugation length. In rubrene, it is a twisting of the C${}_1$–C${}_2$–C${}_4$–C${}_3$ dihedral angle that modulates the conjugation.

In principle, these issues can be addressed using more advanced electronic structure models, but application of such methods to periodic crystals of pharmaceutical-sized molecules is generally cost-prohibitive. A more pragmatic approach combines a conventional periodic DFT calculation using a generalised gradient approximation (GGA) with an intramolecular correction based on the difference between the baseline DFT functional and a more advanced model, computed for a gas-phase molecule [26]. For example, using the B86bPBE+XDM functional [44–46], a GGA that includes a dispersion correction based on the exchange dipole moment (XDM) [46], in combination with dispersion-corrected second-order Møller–Plesset perturbation theory (MP2D) [47,48], the resulting polymorph stability rankings are consistent with experiment for crystalline ROY, acetamidobenzamide, and axitinib [26,40].

A possible alternative solution to the π-delocalisation problem comes in the form of density-corrected (DC-)DFT [49–55]. In this approach, an approximate density functional $\tilde{E}[\rho ]$ is evaluated using a density ${\rho }_{\mathrm{HF}}$ that is obtained at the Hartree–Fock (HF) level: (1) $E_{\text{DC-DFT}}=\tilde{E}[{\rho }_{\mathrm{HF}}].$(1) This procedure, which has also been called HF-DFT [56,57], has been shown to correct some of the worst delocalisation problems associated with GGA functionals, including those with reaction barrier heights [57–59], torsional barriers [60], halogen bonding [61], potential energy surfaces for main-group radicals [49–51], and polaron defects in metal oxide materials [59]. In the present work, we investigate the ability of DC-DFT to correct conformation energies of the π-delocalised systems in Figure 1, and examine how that impacts the predicted stabilities of the corresponding polymorphs.

2. Methods

2.1. DC-DFT

In conventional Kohn–Sham DFT, the approximate energy functional $\tilde{E}[\rho ]$ is minimised with respect to ρ. The resulting energy approximates what would be obtained via optimisation of an exact functional $E[\rho ]$, if the latter were known. Burke, Sim, and coworkers [49–55] suggest that it is pedagogically useful to separate the error $\tilde{E}-E$ into two contributions: a functional-driven part (2) $\mathrm{\Delta }{E}_{\mathrm{F}}={\tilde{E}}_{\mathrm{xc}}[\rho ]-E_{\mathrm{xc}}[{\rho }_0]$(2) that originates in the approximate nature of the exchange-correlation (XC) functional, and a density-driven error (3) $\mathrm{\Delta }{E}_{\mathrm{D}}=\tilde{E}[\tilde{\rho }]-\tilde{E}[{\rho }_0].$(3) Here, $\tilde{\rho }$ indicates the density obtained by variational minimisation of $\tilde{E}$, as in an ordinary self-consistent DFT calculation, whereas ${\rho }_0$ indicates the exact ground-state density obtained upon minimisation of the exact functional $E[\rho ]$. For cases where $|\mathrm{\Delta }{E}_{\mathrm{D}}|\gg |\mathrm{\Delta }{E}_{\mathrm{F}}|$, the overall error may be reduced by evaluating the approximate functional using an SIE-free density (such as ${\rho }_{\mathrm{HF}}$) rather than a self-consistent one [49,52–55].

With this idea in mind, the DC-DFT energy functional in Equation (1(1) $E_{\text{DC-DFT}}=\tilde{E}[{\rho }_{\mathrm{HF}}].$(1) ) can be written more explicitly as (4) $E_{\text{DC-DFT}}[\rho ]=\tilde{E}[\underset{\rho }{\arg \min }(E^{\mathrm{HF}}[\rho ])].$(4) This form makes it clear that ${\rho }_{\mathrm{HF}}$ must first be obtained from a self-consistent HF calculation, starting from the input density ρ. Following this, the energy is obtained by evaluating the approximate functional $\tilde{E}[\rho ]$ using the density ${\rho }_{\mathrm{HF}}$. A more detailed expression for the same DC-DFT energy is (5) $E_{\text{DC-DFT}}[\rho ]=E^{\mathrm{HF}}[\rho ]+({\tilde{E}}_{\mathrm{xc}}[{\rho }_{\mathrm{HF}}]-E_x^{\mathrm{HF}}[\rho ]),$(5) where $E^{\mathrm{HF}}[\rho ]$ is the HF total energy functional, $E_x^{\mathrm{HF}}[\rho ]$ is its exchange component, and ${\tilde{E}}_{\mathrm{xc}}[\rho ]$ is the XC part of $\tilde{E}[\rho ]$. The expression in Equation (5(5) $E_{\text{DC-DFT}}[\rho ]=E^{\mathrm{HF}}[\rho ]+({\tilde{E}}_{\mathrm{xc}}[{\rho }_{\mathrm{HF}}]-E_x^{\mathrm{HF}}[\rho ]),$(5) ) recognises that the functionals $E^{\mathrm{HF}}[\rho ]$ and $\tilde{E}[\rho ]$ differ only in their XC parts, whereas the kinetic energy, Hartree self-repulsion, and nuclear–electron (or ‘external’) potential terms are the same.

The DC-DFT procedure sacrifices self-consistency, which implies that the functional defined in Equations (4(4) $E_{\text{DC-DFT}}[\rho ]=\tilde{E}[\underset{\rho }{\arg \min }(E^{\mathrm{HF}}[\rho ])].$(4) ) and (5(5) $E_{\text{DC-DFT}}[\rho ]=E^{\mathrm{HF}}[\rho ]+({\tilde{E}}_{\mathrm{xc}}[{\rho }_{\mathrm{HF}}]-E_x^{\mathrm{HF}}[\rho ]),$(5) ) is not variational. This necessitates modification to the analytic gradient [57,59]. Lack of self-consistency also means that DC-DFT is not beneficial in all cases [55,62], but rather only in cases that are dominated by density-driven error. Delocalisation driven by SIE tends to fall into the category of density-driven error, as evident from the fact that DC-DFT based on GGA functionals significantly improves results for radicals [49–51,59] and barrier heights [57–59]. (Transition states are usually characterised by partially-broken bonds and thus fractional charges on different chemical moieties, a classic SIE problem [63,64].)

A ‘density sensitivity’ metric (S) has been suggested to determine whether a given problem is dominated by $\mathrm{\Delta }{E}_{\mathrm{F}}$ or $\mathrm{\Delta }{E}_{\mathrm{D}}$ [53–55]. This metric is defined as (6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) for an approximate functional $\tilde{E}[\rho ]$, where ${\rho }_{\mathrm{HF}}$ is the self-consistent HF density and ${\rho }_{\mathrm{LDA}}$ is the self-consistent density computed using the local density approximation (LDA). Burke and co-workers [54,55] suggest S>2 kcal/mol as an indicator of cases that are likely to be dominated by density-driven error (or S>8 kJ/mol in the units that will be used herein), although that threshold was calibrated using small molecules. For larger systems with unpaired spins, one can find examples where S ranges from 150–650 kJ/mol across a set of similar chemical systems [59]. This suggests using (7) $S_{\mathrm{atom}}=S/N_{\mathrm{atoms}}$(7) as a metric that is normalised for system size [59].

2.2. Monomer-corrected DFT

In this work, DC-DFT and other models will be used to correct intramolecular conformational energies and then the monomer correction approach [26] is employed to examine the impact of the corrected conformational energies on polymorph stabilities. In this approach, one first performs a periodic DFT calculation on the full crystal, and herein that energy ($E_{\mathrm{crystal}}^{\mathrm{DFT}}$) is computed using the B86bPBE+XDM functional, as in previous work [26,40]. Then, for each of the N molecules in the unit cell, one computes a gas-phase energy difference between the low-level DFT treatment that is used for the periodic calculation ($E_{\mathrm{molec},i}^{\mathrm{DFT}}$) and a more advanced level of theory ($E_{\mathrm{molec},i}^{\mathrm{high}}$), where the ‘high-level’ method could be DC-DFT or else a correlated wave function method. The crystal energy approximated in this way is (8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) In the present work, the high-level method will be computed with an SIE-free model, whereas B86bPBE+XDM is responsible for modelling the intermolecular interactions. The summation in Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ) can usually be simplified by exploiting symmetry equivalences of molecules within the unit cell. The geometries used for the gas-phase molecular corrections are those extracted directly from the crystal structure.

2.3. Computational details

All DC-DFT calculations were performed using the Q-Chem software [65]. For consistency with previous work, DFT-optimised crystal structures for axitinib, ROY, and molecule X, along with their intramolecular conformations, were taken directly from Ref. [26]. B86bPBE+XDM crystal energies and gas-phase monomer calculations for use in Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ) were also taken from Ref. [26]. All other gas-phase DFT single-point energy calculations were performed using the aug-cc-pVQZ basis set and the SG-2 integration grid [66]. We use the D3 dispersion correction [67] with Becke–Johnson damping, DFT+D3(BJ). Density sensitivities S in Equation (6(6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) ) were evaluated using the SCAN functional [68] for $\tilde{E}[\rho ]$, with the aug-cc-pVQZ basis set.

Additional calculations using the spin-component scaled (SCS-)MP2D method [48] are provided for comparison with the DFT results. SCS-MP2D subtracts out the uncoupled HF dispersion [69] from MP2, replacing it with coupled Kohn–Sham dispersion [70], and then the spin components of the MP2 correlation energy are empirically scaled [48]. This approach is inspired by Heßelmann's corrected MP2 (MP2C) model [71], the key difference being that the SCS-MP2D correction is framed in terms of the atom–atom D3 dispersion correction [67], rather than intermolecular perturbation theory response functions. This makes SCS-MP2D applicable to both inter- and intramolecular interactions. SCS-MP2D does not suffer from SIE, and benchmark calculations demonstrate very good performance for challenging conformational energies [48]. SCS-MP2D results reported here were extrapolated to the complete basis set (CBS) limit based on calculations using aug-cc-pVTZ and aug-cc-pVQZ, following the procedure described in Ref. [26].

Relaxed conformational energy scans are examined for several of the molecules in Figure 1, using geometries taken from the literature [28,38]. For the ROY molecule, the ∠S–C–N–C torsional angle indicated in Figure 1 is used, whereas for rubrene and fluorinated derivatives thereof, the dihedral angle ∠C${}_1$–C${}_2$–C${}_3$–C${}_4$ was used. Single-point potential energy profiles here were then computed at the SCAN+D3, SCAN0+D3, and SCS-MP2D/CBS levels of theory. For ROY, the DFT calculations use the aug-cc-pVQZ basis set as indicated above, although for rubrene and its derivatives, aug-cc-pVTZ is used instead.

These results are compared to benchmarks obtained at the level of coupled-cluster theory with single, double, and perturbative triple excitations [CCSD(T)]. The benchmark results were taken from various sources. For molecule X, ROY, and rubrene, the CCSD(T)/CBS benchmarks were taken from Refs. [26,28,38], respectively. For axitinib, CCSD(T)/CBS results were computed here at the domain-localised pair natural orbital (DLPNO-)CCSD(T${}_1$) level [72], using Orca v. 5.0 [73]. The latter benchmarks combined MP2/CBS results with a DLPNO-CCSD(T${}_1$)/aug-cc-pVDZ correction, obtained using tight PNO settings and with $\mathtt{\text{TCutMKN}}={10}^{-4}$.

3. Results and discussion

3.1. Density sensitivities

Before analysing energy landscapes for the molecules in Figure 1, we first scrutinise the degree of density-driven error in these systems, as it has been observed that DC-DFT offers a better solution than conventional (self-consistent) DFT only when the density-driven error is large [49,52–55]. (Systems with significant static correlation also do not benefit from DC-DFT [62], which is consistent with the idea that SIE may sometimes masquerade as strong correlation [74].) We use the metric S in Equation (6(6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) ), evaluated using the SCAN functional, to assess density sensitivity. Based on small-molecule results, S>8 kJ/mol has been suggested as a threshold beyond which DC-DFT is likely to improve the results [54,55].

Table 1 provides values of S for various molecules considered in this work, computed as an average over all conformers since the variations from one conformer to the next are $\lesssim 2$ kJ/mol. (Differences in S from one conformer to the next will be considered later.) We obtain values ranging from $S=124-528$ kJ/mol although our molecules are considerably larger than the ones used to set 8 kJ/mol as the threshold value. It has also been suggested that S is not a size-extensive metric until the system is large enough for delocalisation to attenuate naturally [59]. (In an infinite conjugated polymer described by a GGA functional, it is not clear that the delocalisation length must ever terminate, as has been demonstrated for excitons described by time-dependent DFT with semilocal functionals [75–77].) For the molecules considered here, the size-normalised metric $S_{\mathrm{atom}}$ [Equation (7(7) $S_{\mathrm{atom}}=S/N_{\mathrm{atoms}}$(7) )] is as large as 7.6 kJ/mol. We will see that all of these molecules are affected by π-delocalisation error, thus the 8 kJ/mol threshold (even if applied to $S_{\mathrm{atom}}$) does not appear to be a sufficient diagnostic.

Table 1. Density sensitivities [Equations (6(6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) ) and (7(7) $S_{\mathrm{atom}}=S/N_{\mathrm{atoms}}$(7) )] for the molecules considered in this work.${}^{\mathrm{a}}$

	Sensitivity (kJ/mol)
Molecule	S	$S_{\mathrm{atom}}$
ROY${}^{\mathrm{b}}$	$126.9\pm 2.2$	$4.70\pm 0.08$
Molecule X${}^{\mathrm{b}}$	$153.0\pm 1.0$	$5.88\pm 0.04$
Axitinib${}^{\mathrm{b}}$	$137.6\pm 1.4$	$2.99\pm 0.03$
Rubrene${}^{\mathrm{c}}$	$123.9\pm 0.8$	$1.77\pm 0.01$
Half-fluoro rubrene${}^{\mathrm{c}}$	$323.9\pm 1.2$	$4.63\pm 0.02$
Perfluororubrene${}^{\mathrm{c}}$	$528.2\pm 1.2$	$7.55\pm 0.02$

${}^{\mathrm{a}}$Average across all conformers; uncertainties represent one standard deviation. ${}^{\mathrm{b}}$Computed at the SCAN/aug-cc-pVQZ level. ${}^{\mathrm{c}}$Computed at the SCAN/aug-cc-pVTZ level.

3.2. Torsional energy profiles

We next consider one-dimensional energy landscapes along some of the torsion angles suggested in Figure 1. For the ROY molecule, intramolecular conformations found in its crystal polymorphs differ primarily by a single dihedral angle that is indicated in Figure 2 [35,40,78–80]. Dihedral angles between ≈50–130${}^{\circ }$ disrupt the π conjugation between the two aromatic rings, leading to a larger band gap and polymorphs that are yellow in colour (‘Y-’). More planar conformations increase the extent of π delocalisation and shift the colour towards red (polymorphs named ‘R-’) or orange (‘O-’).

Figure 2. Relaxed potential energy scans for the ROY molecule along the ∠(S–C–N–C) dihedral angle (as indicated in the inset molecular structure), computed using the SCAN+D3 and SCAN0+D3 functionals and their DC-DFT analogues. SCS-MP2D and the CCSD(T) benchmarks are also shown. Energy profiles are computed relative to the structure with a 120${}^{\circ }$ dihedral angle, which is the minimum at the CCSD(T) level.

Previous work has shown that commonly-used GGA and hybrid functionals overstabilise the red and orange polymorphs relative to the yellow ones [26,28,35,39,40,79]. This can be understood as an SIE effect stemming from over-delocalisation of the π electrons in near-planar geometries, where the system is more conjugated. As shown in Figure 2, both GGAs and hybrid functionals overstabilise the more planar conformations relative to benchmark CCSD(T) calculations. Hybrid functionals reduce the SIE and improve the conformational energy profile, but errors remain significant and are large enough to erroneously alter the polymorph stability rankings [35,39]. Torsional profiles computed using (DC-)B3LYP+D3 and (DC-)PBE0+D3 are shown in Figure S1 and are similar to the (DC-)SCAN0+D3 results in Figure 2.

Figure 3 plots torsional energy profiles of rubrene and some fluorinated derivatives thereof, along the ∠C${}_1$ –C${}_2$–C${}_3$–C${}_4$ angle that is indicated in Figure 1. (The fluorinated derivatives have been suggested as a means of tuning the non-covalent interactions and therefore the crystal structure [81–83], modulating the band gap [82,83], or to enhance stability with respect to oxidation [84].) We refer to ‘planar’ versus ‘twisted’ configurations in Figure 3 as regards the core tetracene unit. As noted in previous work [38], twist angles $\gtrsim {30}^{\circ }$ can increase conjugation, bringing the π electrons of the aromatic side chains into play, and common density functionals thus overstabilise twisted orientations relative to planar ones, even while the twisted configuration is more stable according to CCSD(T) calculations [38]. In the case of the SCAN+D3 functional, DC-DFT improves the relative conformational energies of the planar ($0^{\circ }$) and fully-twisted (${60}^{\circ }$) structures by 2–4 kJ/mol, although the twisted conformation is still somewhat overstabilised even with DC-DFT. This remains the case for the B3LYP+D3 and PBE0+D3 functionals also; see Figure S2. With the SCAN0+D3 functional, DC-DFT improves the relative energies of the fully-twisted structures, and these energies match quite well with the CCSD(T) benchmark for the two out of three derivatives; see Figure 3(a,b).

Figure 3. Relaxed one-dimensional scans of (a) rubrene, (b) half-fluorinated (F${}_{14}$) rubrene, and (c) perfluororubrene, along the dihedral coordinate ∠C${}_1$–C${}_2$–C${}_3$–C${}_4$ that is suggested in Figure 1. The (DC-)DFT energies are computed using the aug-cc-pVTZ basis set. SCS-MP2D and CCSD(T) calculations are taken from Ref. [38], and all structures are taken from that work as well.

3.3. Relative conformer and polymorph energies

The remainder of this work focuses on relative energies of ROY, molecule X, and axitinib, both as isolated molecules (considering the conformations found in the relevant crystal polymorphs) and in the crystalline phase, where the lattice energies are estimated using the monomer-correction scheme in Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ).

The energy ordering of some ROY conformers is shown in Figure 4(a), as computed using various density functionals (both self-consistent and DC-DFT) as well as SCS-MP2D and CCSD(T). According to CCSD(T) calculations [26], the correct gas-phase conformational energy ordering is (9) $\begin{array}{rl}\mathrm{YT}04& <\mathrm{Y}<\mathrm{YN}<\mathrm{R}\approx \mathrm{OP}\approx \mathrm{ON}\\ & <\mathrm{PO}13<\mathrm{ORP}<\mathrm{R}05\mathrm{a}\approx \mathrm{R}05\mathrm{b}.\end{array}$(9) R05a and R05b refer to the two symmetry-unique molecular conformations found in the R05 polymorph. In the self-consistent DFT calculations, shown in the left panel of Figure 4(a), all of the dispersion-corrected GGA and meta-GGA functionals perform very poorly, overstabilising the red and orange conformers over the yellow ones by several kJ/mol. This can be understood as an SIE artifact, as discussed above. PBE+D3 shows the worst performance, completely reversing the energetic ordering relative to CCSD(T). Hybrid functionals partially mitigate the effect of SIE and behave better than the GGA functionals.  The hybrid meta-GGA functional SCAN0+D3 emerges as the best performer amongst the set of functionals tested here, although it still does not reproduce the CCSD(T) energy ordering.

Figure 4. Energetic ordering of conformers of ROY (Figure 1), computed using both self-consistent DFT and DC-DFT for the same set of functionals: (a) isolated monomer energies, and (b) crystal polymorph energies, with the latter computed using the monomer-correction scheme of Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ). Energies are plotted relative to that of the Y conformer.

DC-DFT results for ROY are shown in the middle section of Figure 4, for the same set of functionals. The DC procedure improves upon the self-consistent DFT results in all cases, coming much closer to the correct energetic ordering except in a few quasi-degenerate cases. This validates the hypothesis that π-delocalisation error drives the erratic and mostly erroneous results obtained using conventional DFT. Lowest- and highest-energy conformers are correctly predicted by DC-DFT for all functionals tested here, and the three lowest-energy conformers are correctly predicted by all of the DC-DFT methods except for DC-PBE+D3. Failure of the latter approach may be a sign of larger functional-driven error as compared to the hybrid functionals or the SCAN meta-GGA functional. Even for PBE+D3, however, the DC procedure moves the results closer to the correct energetic ordering as compared to the corresponding self-consistent DFT calculation. The DC-SCAN0+D3 and DC-PBE0+D3 methods afford correct energetic ordering for most of the conformers, as judged by comparison to CCSD(T) benchmarks. The success of DC-SCAN and DC-SCAN0 has been noted in unrelated work [85].

Having identified SCAN+D3 and SCAN0+D3 as the best of the self-consistent DFT approaches for the ROY molecule, let us return to the one-dimensional potential energy scans in Figure 2. Both DC-SCAN+D3 and DC-SCAN0+D3 perform well in comparison to benchmark CCSD(T) calculations from Ref. [26], although neither is quite as faithful to the benchmark as is SCS-MP2D. In contrast, the self-consistent SCAN+D3 and SCAN0+D3 calculations yield erroneously large rotational barriers around 40–90${}^{\circ }$, as structures with angles of 0–30${}^{\circ }$ are overstabilised due to SIE. Torsional energy profiles computed using (DC-)B3LYP+D3 and (DC-)PBE0+D3 are shown in Figure S1 and exhibit most of the same trends that are seen for the hybrid (DC-)SCAN0+D3 level of theory in Figure 2. Thus, the good performance of DC-SCAN+D3 and DC-SCAN0+D3 is not simply a matter of choosing a good functional to correct; the correction itself improves the results substantially.

Density sensitivities for ROY conformers are plotted in Figure 5(a), computed using SCAN for the approximate functional $\tilde{E}[\rho ]$ in Equation (6(6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) ). Yellow conformers (indicated by names starting with Y) exhibit smaller values of S as compared to the other conformers, by $\approx 4$ kJ/mol, a significant amount given the rather small energy differences between the various conformers. Larger SIE in near-planar geometries associated with the orange and red conformers explains the larger values of S in those cases. Application of DC-DFT thus selectively modifies the energy landscape of the conformers with the largest density sensitivity, meaning red and orange more so than yellow in this case. The result is a potential energy surface that pushes the energetic ordering in the right direction.

Figure 5. Density sensitivities [Equation (6(6) $S=|\tilde{E}[{\rho }_{\mathrm{LDA}}]-\tilde{E}[{\rho }_{\mathrm{HF}}]|$(6) )] for different conformers of (a) ROY, (b) molecule X, and (c) axitinib, computed using SCAN/aug-cc-pVQZ for the approximate functional $\tilde{E}[\rho ]$.

Up to this point, we have discussed the behaviour of an isolated molecule of ROY. Figure 4(b) reports the relative energies of the crystal polymorphs, computed using the monomer-correction scheme in Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ). As reported previously [26], the B86bPBE+XDM functional yields extremely poor relative energies for the  form Y, predicting it to be the second-least stable polymorph when in fact it is the thermodynamically most stable form. The relative order of the other lattice energies is also poor at this level of theory, relative to experimental lattice enthalpies or CCSD(T) lattice energies. Furthermore, replacing B86bPBE+XDM with another self-consistent GGA or hybrid functional does not meaningfully improve the energy rankings. In contrast, the DC-DFT models all perform much better. The close energy spacings of the ROY polymorphs make this a challenging problem for any method, but SCS-MP2D and DC-DFT with hybrid functionals all give fairly consistent qualitative rankings, even if the quantitative lattice energies differ.

We next perform a similar analysis for molecule X. The relevant conformers involve rotation of the bonds indicated with arrows in Figure 1 and alter the extent of π conjugation with the central aromatic ring. According to CCSD(T) calculations [26], the gas-phase energetic ordering of the various conformers is (10) $\begin{array}{rl}\text{van Eijck-3}& <\text{van Eijck-2}<\text{Dzyabchenko-3}\\ & <\text{Ammon-2}<\mathrm{Expt}.<\text{Day-1}.\end{array}$(10) The nomenclature in Equation (10(10) $\begin{array}{rl}\text{van Eijck-3}& <\text{van Eijck-2}<\text{Dzyabchenko-3}\\ & <\text{Ammon-2}<\mathrm{Expt}.<\text{Day-1}.\end{array}$(10) ) reflects structures submitted by different groups in the third blind test of crystal structure prediction [42], with ‘Expt.’ indicating the conformer in the experimental crystal structure. DFT calculations were previously found to reverse the order of the experimental conformer with that of conformer Day-1, and furthermore to overstabilise the van Eijck-3 conformer by about 2 kJ/mol [26].

DFT and DC-DFT results for molecule X in the gas phase are shown in Figure 6(a). Similar to trends observed previously [26,41], each of the GGA and meta-GGA functionals that we consider overestimates the stability of the van Eijck-3 conformer relative to the experimental one. Of the conventional (self-consistent) DFT approaches, only BLYP+D3, B3LYP+D3, and SCAN0+D3 predict qualitatively correct relative energies for the closely-spaced experimental and the Day-1 conformers, which are the least-stable conformers in isolation. On the other hand, all of the DC-DFT methods examined here place those two conformers in the correct energetic order. The DC procedure also reduces the overstabilisation of conformer van Eijck-3, relative to the self-consistent DFT results.

Figure 6. Energetic ordering of conformers of molecule X (Figure 1), computed using both self-consistent DFT and DC-DFT for the same set of functionals: (a) isolated monomer energies, and (b) crystal polymorph energies, with the latter computed using the monomer-correction scheme of Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ). Energies are plotted relative to ‘Expt.’, which is the conformer found in the experimental crystal structure.

Density sensitivities for conformers of molecule X are plotted in Figure 5(b). We mainly focus on the Day-1 conformer and the experimental one since these are the relative energies that are most susceptible to switching order as one functional is swapped for another. The value of S for the Day-1 conformer is 1.7 kJ/mol larger than the value obtained for the experimental conformer, which is perhaps significant in view of the fact that these two conformers differ only by about 1 kJ/mol, according to DLPNO-CCSD(T) calculations [26]. The large sensitivity for the Day-1 conformer is consistent with its overstabilisation (relative to the experimental conformer) by functionals such as PBE and PBE0. Overstabilisation of the van Eijck-3 conformer with conventional DFT can also be explained by the larger S associated with that conformer.

Relative polymorph energies for molecule X are given in Figure 6(b). Once again, the predicted lattice energies improve significantly when DC-DFT is used instead of self-consistent DFT. At the same time, some of the energy orderings remain sensitive to the choice of functional, even when DC-DFT is used. This largely reflects the difficulty of the problem, with relative energy changes of ∼1  kJ/mol being enough to affect the ordering. In this particular example, DC-SCAN+D3 (based on semilocal SCAN rather than the hybrid SCAN0) best reproduces the CCSD(T) benchmarks.

We next discuss results for the anti-cancer drug axitinib. This drug has five known conformational polymorphs that are denoted I, IV, VI, XXV, and XLI, and many more multi-component crystal forms [86,87]. In most polymorphs, the axitinib molecule adopts an extended conformation but in the thermodynamically most stable XLI form, however, rotation about the sulphur group leads to a more folded conformation. In the process, the π conjugation between the aromatic ring and the terminal amide moiety is disrupted; see Figure 1. This causes the XLI form to be artificially destabilised by GGA functionals, leading to erroneous polymorph energy rankings [26].

Relative energies of axitinib conformers are shown in Figure 7. New DLPNO-CCSD(T${}_1$)/CBS calculations performed here afford the following energetic ordering for the isolated monomer: (11) $\mathrm{IVb}<\mathrm{XLI}<\mathrm{VI}<\mathrm{XXV}<\mathrm{IVa}<\mathrm{I}.$(11) None of the self-consistent functionals reproduces this ordering, and many of the functionals (with notable exceptions BLYP+D3 and B3LYP+D3) significantly destabilise form XLI, which is the conformer that appears in the lowest-energy crystal polymorph. DC-DFT generally shifts the conformational energies in the correct direction by stabilising XLI relative to the other conformers. However, the extent to which this is true varies by several kJ/mol depending upon the choice of functional, and for the BLYP+D3 and B3LYP+D3 functionals that perform reasonably well in conventional self-consistent calculations, the errors increase when DC-DFT is used. Other functionals do improve when DC-DFT replaces a self-consistent calculation, and the best-performing DFT approach is DC-PBE+D3, with root-mean-square (RMS) error of 0.9 kJ/mol with respect to the CCSD(T) benchmarks. However, SCS-MP2D does slightly better, with a RMS error of 0.5 kJ/mol.

Figure 7. Energetic ordering of conformers of axitinib (Figure 1), computed using both self-consistent DFT and DC-DFT for the same set of functionals: (a) isolated monomer energies, and (b) crystal polymorph energies, with the latter computed using the monomer-correction scheme of Equation (8(8) $E_{\mathrm{crystal}}=E_{\mathrm{crystal}}^{\mathrm{DFT}}+\sum _{i=1}^N(E_{\mathrm{molec},i}^{\mathrm{high}}-E_{\mathrm{molec},i}^{\mathrm{DFT}}).$(8) ). Energies are plotted relative to that of the XLI conformer.

Relative energies of the axitinib crystal polymorphs are shown in Figure 7(b) and the various DC-DFT approaches perform similarly except that the energy gap between the lowest-energy polymorph (XLI in all cases) and the others varies by up to 8 kJ/mol depending on the choice of functional. The ordering of the (near-degenerate) polymorphs VI and XXV is uncertain in the DC-DFT calculations, but that is true experimentally as well [88]. Otherwise, DC-DFT consistently predicts the correct energetic ordering of the polymorphs, for all functionals tested, although the energies are spread out over too large a range in some cases, by up to 4 kJ/mol.

The relatively poor performance of DC-DFT for axitinib, as compared to ROY or molecule X, is a bit puzzling given that $S_{\mathrm{atom}}$ is smallest for axitinib. That observation, combined with the fact that self-consistent BLYP+D3 and B3LYP+D3 calculations perform well for axitinib, suggests that functional-driven error may exceed density-driven error for this molecule. Given the nature of the conformers that are considered here, the π-delocalisation issues are largely confined to one end of this molecule, and the considerably larger size of axitinib (relative to either ROY or molecule X) may wash out the advantages of DC-DFT in this case. It is possible that a fragmentation approach [89], with DC-DFT applied only to the part of the molecule that is subject to significant delocalisation problems, could improve the situation, but we have not pursued such an approach here. In any case, results for axitinib highlight how density-driven error is only one factor that affects conformational energies in large, conjugated organic molecules, and cautions against assuming that DC-DFT is a panacea for such problems.

Stepping back from the finer details of each system, Figure 8 plots errors in conformational energies for all three molecules (ROY, molecule X, and axitinib), relative to CCSD(T) benchmarks, indicating the RMS error for each DFT approach and for SCS-MP2D. (The RMS errors are also summarised in Figure S3.) Relative to self-consistent DFT, the DC-DFT approach generally reduces the width of the error distributions by up to a factor of 3–4, and DC-hybrid functionals afford the tightest error distributions. At the same time, even the best of the DC-DFT models affords errors that are about twice as large as those obtained using SCS-MP2D. Taken together, these results highlight how DC-DFT can provide a useful, pragmatic approach for improving the quality of common DFT functionals in challenging systems with density-driven delocalisation error, but also suggest that higher-level method may still be preferred for high-accuracy work.

Figure 8. Summary of conformational energy errors (in kJ/mol) for isolated molecules ROY, molecule X, and axitinib, relative to CCSD(T) benchmarks. The box-and-whisker plots summarise quartiles of the error distributions, while the points show the errors for individual conformational energies.

4. Conclusion

We have used the DC-DFT technique to study the energy landscape of several conjugated organic molecules, as well as their relative polymorph energies in the crystalline phase. These are challenging problems due to the appearance of multiple low-energy conformers (and corresponding low-energy polymorphs), and self-consistent DFT calculations with standard semilocal and hybrid functionals fail to reproduce benchmark ab initio results. We attribute this to SIE, in the particular form of π-delocalisation error as the degree of conjugation is tuned by means of various intramolecular torsion angles. Application of DC-DFT, using the same common functionals, significantly mitigates the errors in most cases and usually results in correct energetic ordering of crystal polymorphs, when computed by means of a monomer-correction scheme [26] applied to a periodic DFT calculation.

The DC-SCAN+D3 and DC-SCAN0+D3 methods work reasonably well for conformers and polymorphs of the benchmark ROY and molecule X systems, although the larger molecule axitinib proves to be more difficult. For axitinib, most (but not all) of the DC-DFT approaches improve upon their conventional self-consistent analogues, but none of the DFT methods produces fully satisfactory agreement with CCSD(T) benchmarks. We speculate that the relevant π-delocalisation error in this case may be isolated within a particular aromatic moiety, and that DC-DFT applied to the entire (much larger) molecule may upset the balance of reasons why particular functionals perform well for the parts of the molecule that do not suffer from over-delocalisation of π electrons. A fragmentation strategy might help, with DC-DFT applied only to part of the molecule, but in any case this system highlights potential difficulties in applying DC-DFT to systems where SIE problems are confined to a relatively small portion of a much larger molecule. Overall, our results underscore previous work [26,40] indicating that monomer-corrected periodic GGA calculations, with an appropriately selected level of theory for the monomer corrections, is an effective and low-cost strategy for polymorphic molecular crystals.

Acknowledgments

The authors express condolences for the family of Prof. Nick Besley, who was a valued member of the Q-Chem community.

Disclosure statement

J.M.H. serves on the board of directors of Q-Chem Inc.

Additional information

Funding

Work by B.R. and J.M.H. was supported by the National Science Foundation (NSF), [grant number CHE-1955282], with calculations performed at the Ohio Supercomputer Center [90]. G.J.O.B. acknowledges support from the NSF [grant mber CHE-1955554] and supercomputer time from XSEDE [grant number CHE110064].