Electrostatics does not dictate the slip-stacked arrangement of aromatic π–π interactions

Abstract

Benzene dimer has long been an archetype for π-stacking. According to the Hunter–Sanders model, quadrupolar electrostatics favors an edge-to-face CH⋯π geometry but competes with London dispersion that favors cofacial π-stacking, with a compromise “slip-stacked” structure emerging as the minimum-energy geometry. This model is based on classical electrostatics, however, and neglects charge penetration. A fully quantum-mechanical analysis, presented here, demonstrates that electrostatics actually exerts very little influence on the conformational landscape of (C₆H₆)₂. Electrostatics also cannot explain the slip-stacked arrangement of C₆H₆⋯C₆F₆, where the sign of the quadrupolar interaction is reversed. Instead, the slip-stacked geometry emerges in both systems due to competition between dispersion and Pauli repulsion, with electrostatics as an ambivalent spectator. This revised interpretation helps to rationalize the persistence of offset π-stacking in larger polycyclic aromatic hydrocarbons and across the highly varied electrostatic environments that characterize π–π interactions in proteins.

---

Introduction

Understanding the factors that govern geometric preferences of π–π interactions is of vital importance in crystal engineering,¹ with implications as well for understanding protein structure^2–4 and biological recognition, including drug design.⁵ Benzene dimer has long served as an archetype for understanding the geometric preferences of aromatic π–π interactions,^6–9 though its emblematic status has occasionally been questioned.¹⁰

Relevant conformations of (C₆H₆)₂ include a cofacial “sandwich” geometry (Fig. 1a), a “slip-stacked” or parallel-displaced geometry (Fig. 1b), and a T-shaped isomer characterized by a CH⋯π interaction (Fig. 1c). For (C₆H₆)₂ in the gas phase, the cofacial π-stacking arrangement is an energetic saddle point along a sliding coordinate that leads to the slip-stacked structure,⁷ which is essentially iso-energetic with the T-shaped geometry.^6–9 The latter is also a saddle point,⁸ with the slightly tilted structure depicted in Fig. 1d emerging as the global minimum.⁸ However, the potential surface is extremely flat along the pendular coordinate that connects the canonical T-shaped and the tilted isomer, which differ by only about 0.2 kcal mol⁻¹.^8,9

Fig. 1 Stationary points on the (C₆H₆)₂ potential energy surface: (a) cofacial π-stacked configuration, (b) “slip-stacked” or “offset-stacked” configuration, (c) canonical T-shaped geometry, and (d) tilted T-shaped geometry. For each structure, a schematic view of the quadrupolar charge distribution is provided.

In a landmark paper,¹¹ Hunter and Sanders developed a simple model to explain these conformational preferences, based on benzene's sizable quadrupole moment.^11–13 As illustrated in Fig. 1, the Hunter–Sanders model suggests that cofacial π-stacking requires quadrupolar electrostatic repulsion to be compensated by attractive dispersion. In contrast, the T-shaped isomer is stabilized by quadrupolar electrostatics but sacrifices some of the favorable dispersion interaction, because the π clouds of the two monomers are farther apart as compared to the cofacial geometry. The adversarial relationship between electrostatics and dispersion has since emerged as the principle paradigm for interpreting conformational preferences in aromatic systems,^13–21 and in particular is used to explain the emergence of the slip-stacked configuration in (C₆H₆)₂.¹¹

Since the seminal paper by Hunter and Sanders,¹¹ additional analysis has begun to erode the simple electrostatic picture proffered by this model. For example, a thorough survey of stationary points on the (C₆H₆)₂ potential surface reveals that electrostatic interactions are actually attractive in the sandwich geometry.⁸ The need for a quantum-mechanical treatment of electrostatics is underscored by studies of substituted benzene dimers by Sherrill and co-workers^22–26 that illustrate the importance of charge-penetration effects (missing in the quadrupolar electrostatic picture) and demonstrate that electrostatic attraction increases upon substitution, regardless of whether the substituent is electron-withdrawing or electron-donating. Work by these authors,^22,26 as well as by Wheeler and Houk,^27–32 demonstrates that the Hunter–Sanders model does not always adequately describe the effects of substitution on aromatic π-stacking energies, and draws attention to the ways in which electrostatic potential maps can sometimes be misleading.^28–32 However, none of this work directly contests the role of electrostatics in determining geometries of π-stacked complexes. Whereas Grimme notes that the Hunter–Sanders model “overemphasizes” quadrupolar electrostatic interactions,¹⁰ we will suggest that this model simply gets the molecular physics wrong, and that electrostatics in competition with dispersion is a flawed framework for understanding π-stacking interactions.

Despite many efforts to advance the understanding of π–π interactions beyond the Hunter–Sanders model, there has been no attempt to revisit the physical explanation for the origins of offset π-stacking. The Hunter–Sanders model suggests that π–π interactions should always favor the offset-stacked arrangement over the cofacial geometry, and this fact has been taken to suggest that no special π-stacking interaction actually exists.³³ Note, however, that the minimum-energy geometry in C₆H₆⋯C₆F₆ is also parallel-displaced,^34–36 despite the fact that the C–F bonds reverse the polarity of the charge distribution (relative to that in C₆H₆), leading to a quadrupolar electrostatic interaction that is attractive rather than repulsive in the cofacial arrangement. Aromatic heterocycles, whose multipole moments are surely quite different from those of either C₆H₆ or C₆F₆, also adopt offset-stacked geometries.^37,38 It has also been suggested that quadrupolar repulsion cannot account for the structures adopted by larger polycyclic aromatics, and that geometries of larger graphene analogues may be better understood via Pauli repulsion models.³⁹ Below, we show that even in the archetypal case of (C₆H₆)₂, the Hunter–Sanders model fails to capture the essence of π–π interactions. We propose an alternative but equally simple explanation in terms of competition between London dispersion and Pauli repulsion, i.e., a “van der Waals” (vdW) interaction model.

Methods

We investigate potential energy surfaces for both cofacial and CH⋯π arrangements of (C₆H₆)₂, along a one-dimensional center-to-center sliding coordinate, using energy decomposition analysis (EDA) based on symmetry-adapted perturbation theory (SAPT),^40–43 the preeminent ab initio theory of intermolecular interactions.⁴⁴ Specifically, we employ an “extended” SAPT formalism that includes a variational description of polarization effects,^45–47 though we also include a δE_HF correction for polarization⁴² and a many-body dispersion contribution.^48,49 This approach has been shown to provide accurate results in dispersion-dominated complexes.^49,50 That said, the same conclusions emerge from more traditional SAPT0/jun-cc-pVDZ calculations, an approach that is known to provide accurate interaction energies for (C₆H₆)₂,^51,52 as well as from EDA based on absolutely localized molecular orbitals (ALMO-EDA).^53,54 These additional calculations are presented in Section S1 of the ESI.† All calculations were performed using Q-Chem v. 5.3.⁵⁵

Results and discussion

In the analysis that follows, we will often group together the electrostatic (“elst”) and induction (“ind”, also known as polarization) components of the interaction energy. This “elst + ind” energy (E_elst + E_ind) represents the quantum-mechanical Coulomb interaction between polarized charge distributions for the two monomers, without any multipole approximation, and fully accounts for interpenetration of the monomer charge densities. We examine this and other energy components along a parallel-displacement or “sliding” coordinate, in either a parallel or a perpendicular configuration. In the former, the distance between the two molecular planes is fixed at 3.4 Å, which is characteristic of the parallel-displaced local minimum in Fig. 1b and smaller than the 3.8 Å separation that characterizes the cofacial (sandwich) saddle point in Fig. 1a. In the perpendicular orientation, the center of one monomer is displaced at a fixed distance of 5.0 Å (representing the T-shaped saddle point in Fig. 1c) from the plane containing the other monomer.

SAPT-based EDA suggests that the elst + ind interaction energy is nearly flat along the coordinate corresponding to sliding two benzene monomers between cofacial and slip-stacked geometries; see Fig. 2. In the perpendicular configuration of (C₆H₆)₂, the sum of electrostatics and induction exhibits a slight preference (by <0.5 kcal mol⁻¹) for a parallel-displaced version of the T-shaped isomer that is not actually a local minimum of the total interaction energy. For both the parallel and perpendicular configurations, these results are precisely opposite to the Hunter–Sanders prediction of a repulsive electrostatic maximum at the sandwich geometry!

Fig. 2 Polarized Coulomb interaction E_elst + E_ind along a sliding coordinate in (C₆H₆)₂. The zero-displacement position represents either the sandwich or the T-shaped structure, as shown, with center-to-center monomer distances of 3.4 Å (parallel) and 5.0 Å (perpendicular) at zero displacement.

This contradiction can be understood in terms of charge penetration that is absent in a model based on classical multipoles. SAPT uses the full electron density to compute both electrostatics and induction, and this has the effect of stabilizing the cofacial structure as compared to the slip-stacked structure; see Fig. S1b.† This stabilization is largely canceled by induction, leading to the relatively flat E_elst + E_ind potential in Fig. 2. ALMO-EDA also predicts electrostatic stabilization of the sandwich geometry (Fig. S5†). At larger intermolecular separation, where the density overlap between monomers is negligible, classical electrostatics is recovered and quadrupolar repulsion does become the dominant contribution to electrostatics (see Fig. S3†), but charge penetration is significant at intermolecular separations that typify π–π interactions.^8,24

The polarized Coulomb potential (E_elst + E_ind) alone fails to afford a meaningful saddle point at the sandwich geometry (see Fig. S1b†), and the only energy component that does have a meaningful local maximum in the sandwich configuration is the exchange term (see Fig. S2†). This suggests that it is Pauli repulsion, rather than electrostatics, that provides the driving force towards the slip-stacked arrangement. In fact, when exchange repulsion is removed from the interaction energy, the sandwich geometry emerges as the most stable one, as shown in Fig. 3. With Pauli repulsion removed, the perpendicular configuration demonstrates essentially no preference between the conventional T-shaped isomer and another perpendicular configuration in which the CH⋯π interaction is moved to the edge of the aromatic ring.

Fig. 3 Benzene dimer interaction energies without exchange repulsion, E_int − E_exch = E_elst + E_ind + E_disp. The parallel-displacement or sliding coordinate is the same as in Fig. 2.

A more meaningful picture of π–π interactions can be cultivated by considering a “vdW” interaction potential (E_vdW), which we define to be the sum of Pauli repulsion and dispersion, or in other words the total interaction energy (E_int) minus the polarized Coulomb (E_elst + E_ind) contribution:

E_vdW = E_exch + E_disp = E_int − E_elst − E_ind.

In comparison to the Coulomb + dispersion potentials in Fig. 3, the vdW potential energy scans in Fig. 4 are much more suggestive of the true stationary points in (C₆H₆)₂, in both its parallel and its perpendicular orientation. The perpendicular orientation (Fig. 4a) exhibits an energy minimum at the T-shaped configuration, not because of attractive electrostatics but rather because dispersion is maximized and Pauli repulsion is minimized. The exchange interaction drives the perpendicular orientation away from “L-shaped” geometries and towards the T-shaped one, because the former exhibits a steric clash between a hydrogen atom on one monomer and the C–C bond density on the other. As noted above, the canonical T-shaped isomer (Fig. 1c) is a saddle point, and a small tilt leads to a minimum that is 0.2 kcal mol⁻¹ lower in energy (Fig. 1d).⁸ However, the potential surface along this “pendular” coordinate is quite flat, and sliding potentials for the tilted versus T-shaped configurations differ barely at all, as shown in Fig. 4a.

Fig. 4 Total interaction energies and vdW interaction potentials (E_vdW = E_exch + E_disp), along the parallel-displacement coordinate of (C₆H₆)₂ in (a) perpendicular and (b) parallel orientations. There are two sets of curves in (a), representing the canonical T-shaped isomer (Fig. 1c) and its tilted analogue (Fig. 1d).

In the cofacial arrangement, (C₆H₆)₂ adopts a slip-stacked geometry that is also a minimum on the vdW surface (Fig. 4b). Because the electron density exhibits local maxima at the nuclei, the sandwich isomer (in which the nuclei of one monomer are situated directly above those of the other monomer) maximizes Pauli repulsion, leading this configuration to be a saddle point on the total interaction potential. The parallel-displaced arrangement reduces this repulsion somewhat, without sacrificing too much of the favorable dispersion interaction. Whereas the elst + ind energy landscape is largely featureless, the vdW landscape exhibits a clear preference for the parallel-displaced and T-shaped geometries.

It should be noted that the cofacial saddle point on the (C₆H₆)₂ potential surface is characterized by a larger intermolecular separation (3.8 Å) as compared to the parallel-displaced minimum (3.4 Å).⁸ For calculations in parallel arrangements, the one-dimensional scans in Fig. 2–4 fix the distance between the two molecular planes at 3.4 Å, which is inside of the repulsive regime for the sandwich geometry. As such, it is worth considering whether this sliding coordinate is a valid proxy for the underlying intermolecular forces. Fundamental to our revised interpretation of these forces is that the nature of the repulsive interaction responsible for offset-stacking is primarily exchange antisymmetry rather than electrostatics, so let us next consider the repulsive forces in detail.

At an intermolecular separation of 3.4 Å, consistent with the parallel-displaced minimum-energy geometry, the vdW potential for the cofacial sandwich structure is strongly repulsive (Fig. 4b) whereas the E_elst + E_ind potential is flat (Fig. 2). Pauli repulsion therefore dominates in this close-contact regime, and there are two possible avenues by which the system may ameliorate this repulsion. It can do so either by offsetting the nuclei, leading to the parallel-displaced minimum-energy structure, or else by increasing the intermolecular separation, thereby affording the cofacial saddle point at a separation of 3.8 Å rather than 3.4 Å. This increase in the intermolecular separation reduces charge penetration, making the electrostatics slightly less favorable in the cofacial arrangement as compared to the parallel-offset geometry, but we view this difference in the electrostatics as an effect (driven by Pauli repulsion) rather than a cause. The sandwich structure sacrifices some attractive components of its interaction energy by moving to larger distance in order to minimize Pauli repulsion, whereas the parallel-offset can alleviate Pauli repulsion without significantly increasing the average internuclear separation, thereby affording a lower total energy as compared to the cofacial geometry. Examination of the fixed-separation sliding coordinate is therefore more instructive than simply considering stationary points on the full potential energy surface, because consideration of the repulsive region (e.g., the cofacial arrangement at 3.4 Å) better elucidates the fundamental forces at play. Consideration of stationary points alone might lead one to misattribute the slip-stacking phenomenon to electrostatics and charge penetration.

The C₆H₆⋯C₆F₆ system offers an interesting contrast to benzene dimer because its quadrupolar electrostatic interaction is attractive in the cofacial arrangement, which can be understood by switching the signs (δ±) in the charge distribution of one monomer in Fig. 1a. Results for benzene dimer already indicate that the Hunter–Sanders model overemphasizes electrostatics, as has been suggested previously.^8,10,25 Correspondingly, this model predicts a sandwich arrangement for C₆H₆⋯C₆F₆ as shown in Fig. S6.† SAPT calculations, however, corroborate the notion of an attractive electrostatic interaction that is most significant in the cofacial geometry, but at the same time suggest a parallel-displaced minimum that is ≈1 kcal mol⁻¹ lower in total interaction energy. Noting that the total interaction energy

E_int = E_vdW + E_elst + E_ind

is the sum of vdW and elst + ind components, which are plotted separately in Fig. 5, it is clear that the salient topographical features of the potential surface are inherited from the vdW energy, not from electrostatics. The preference for the parallel-displaced geometry in C₆H₆⋯C₆F₆ arises for the same reason that it does in (C₆H₆)₂, namely, reducing the density overlap by offsetting the nuclei. Even with its more favorable electrostatics, the energy landscape of cofacial C₆H₆⋯C₆F₆ is controlled by steric effects.

Fig. 5 Energy components for a parallel arrangement of C₆H₆⋯C₆F₆ along the cofacial sliding coordinate. The total interaction energy E_int = E_vdW + E_elst + E_ind is the sum of vdW and elst + ind components.

In the spirit of preserving a simple model potential capable of qualitative predictions, and to further emphasize the utility in rethinking π–π interactions in terms of vdW forces, we have formulated a two-component model potential. It consists of a repulsive potential proportional to the overlap of atom-centered spheres and an attractive dispersion potential, for which we use the atomic-pairwise aiD3 potential.⁴⁶ (For details, consult Section S2 of the ESI†). Parallel-displacement potentials for both (C₆H₆)₂ and C₆H₆⋯C₆F₆, as predicted by this simple vdW model, are shown in Fig. 6. The model correctly predicts that repulsion is maximized in the on-top sandwich arrangement, and that a T-shaped geometry is favored in the perpendicular edge-to-face configuration. Whereas the Hunter–Sanders model predicts that the cofacial arrangement is a local minimum in C₆H₆⋯C₆F₆, our vdW model correctly predicts that the sandwich configuration is a saddle point, in agreement with ab initio calculations. The vdW model also predicts that the tilted analogue of the T-shaped isomer (Fig. 1d) is slightly lower in energy than the canonical T-shaped structure, again in agreement with ab initio results.^8,9 The model suggests that the tilt angle adopted by the edge-to-face local minimum of (C₆H₆)₂ is driven by gains in dispersion at the expense of a small increase in Pauli repulsion.

Fig. 6 Potential energy scans for parallel, perpendicular, and tilted-perpendicular arrangements of (C₆H₆)₂, along with the parallel configuration of C₆H₆⋯C₆F₆, obtained from an empirical vdW model that combines an overlap-based repulsive potential with a pairwise-atomic dispersion potential.

Pauli repulsion and dispersion are both size-extensive quantities but quadrupolar electrostatics is not, and the latter will saturate in larger polycyclic aromatic systems. In view of this, we examined anthracene dimer using both the Hunter–Sanders model and the new vdW model, for which potential energy surfaces are plotted in Fig. 7 and compared to ab initio geometries. The Hunter–Sanders model (Fig. 7a) affords a flat potential for T-shaped (C₁₄H₁₀)₂ whereas the vdW model (Fig. 7b) correctly predicts three stationary points corresponding to different CH⋯π motifs. By matching dispersion with commensurate Pauli repulsion, the vdW model correctly predicts that the minimum-energy geometry in the parallel arrangement is slip-stacked, driven by the reduction in Pauli repulsion that comes from offsetting the nuclei. In contrast, the Hunter–Sanders model predicts an offset that is perpendicular to the anthracene ribbon, but no offset along the direction of the ribbon itself, at odds with ab initio results. The vdW model is therefore more faithful to ab initio geometry optimizations, while retaining the simplicity of the older Hunter–Sanders model. As such, the new model would seem to be a good starting point for constructing more sophisticated empirical force fields whose underlying physics is qualitatively sound.

Fig. 7 Potential energy scans for parallel displacement of anthracene dimer, as predicted using (a) the conventional Hunter–Sanders model,¹¹versus (b) the vdW model developed in this work. In (b), the stationary-point geometries coincide with those computed at the TPSS+D3/def2-TZVPP level, which has been recommended for geometry optimizations in large π–π systems.⁵⁶

Conclusions

Models based on quadrupolar electrostatics, as a means to explain the geometric preferences of (C₆H₆)₂, are textbook material in the field of supramolecular chemistry,^15,16 despite the abject failure of such models to explain offset-stacking in C₆H₆⋯C₆F₆, where the sign of the quadrupolar interaction is reversed as compared to benzene dimer. The electrostatics-driven picture of π–π interactions, while valid at sufficiently long range, fundamentally misrepresents the nature of these interactions at typical π-stacking and CH⋯π distances, where charge penetration is significant. More robust conclusions are reached by considering steric repulsion in competition with London dispersion, with electrostatics largely sidelined. A simple vdW interaction potential, introduced here, rationalizes the ab initio results and might potentially be developed into a quantitative parameterization for use in force fields.

Armed with a better understanding of the physics that govern π–π interactions, it may be useful to revisit questions regarding the relevance of the “π-stacking” concept,³³ and the competition between London dispersion and steric repulsion more generally.⁵⁷ Here, we have sought to develop an intuitive model with broad implications for scenarios encountered in biochemistry and crystallography. Our model posits that Pauli repulsion is the dominant force that competes with dispersion in (C₆H₆)₂ and C₆H₆⋯C₆F₆. Insofar as both forces are always present at short range, this may explain why the slip-stacked geometry (also known as offset-stacking or parallel-displaced π-stacking) emerges as a recurring motif in the π–π interactions found in protein crystal structures,^2,13 across myriad electrostatic environments, and also in π-stacked complexes involving aromatic heterocycles.^23,37,38 Electrostatics likely does become relevant once substituents are incorporated (C₆H₅X),^22,24–32 introducing bond dipoles, as well as in solid-state architectures where long-range forces are more important. As an example of the latter, mixtures of C₆F₆ with benzene and substituted benzene derivatives are known to organize into columnar π-stacks,^13,14 with polymorphism that is thought to be controlled by subtle variations in bond-dipole versus quadrupolar electrostatics.^58–60 These phenomena deserve to be considered carefully in light of a new interpretation of π–π interactions.

Conflicts of interest

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

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Award No. DE-SC0008550. Calculations were performed at the Ohio Supercomputer Center under project no. PAA-0003.⁶¹

Notes and references

1. R. Thakuria, N. K. Nath and B. K. Saha, The nature and applications of π-π interactions: a perspective, Cryst. Growth Des., 2019, 19, 523–528.
2. G. B. McGaughey, M. Gagné and A. K. Rappé, π-stacking interactions: alive and well in proteins, J. Biol. Chem., 1998, 273, 15458–15463.
3. P. Chakrabarti and R. Bhattacharyya, Geometry of nonbonded interactions involving planar groups in proteins, Prog. Biophys. Mol. Biol., 2007, 95, 83–137.
4. S. Marsili, R. Chelli, V. Schettino and P. Procacci, Thermodynamics of stacking interactions in proteins, Phys. Chem. Chem. Phys., 2008, 10, 2673–2685.
5. E. A. Meyer, R. K. Castellano and F. Diederich, Interactions with aromatic rings in chemical and biological recognition, Angew. Chem., Int. Ed., 2003, 42, 1210–1250.
6. M. O. Sinnokrot, E. F. Valeev and C. D. Sherrill, Estimates of the ab initio limit for π-π interactions: the benzene dimer, J. Am. Chem. Soc., 2002, 124, 10887–10893.
7. M. O. Sinnokrot and C. D. Sherrill, Highly accurate coupled cluster potential energy curves for the benzene dimer: sandwich, T-shaped, and parallel-displaced configurations, J. Phys. Chem. A, 2004, 108, 10200–10207.
8. R. Podeszwa, R. Bukowski and K. Sczalewicz, Potential energy surface for the benzene dimer and perturbational analysis of π–π interactions, J. Phys. Chem. A, 2006, 110, 10345–10354.
9. O. Bludský, M. Rubeš and P. Nachtigall, Investigation of the benzene-dimer potential energy surface: DFT/CCSD(T) correction scheme, J. Chem. Phys., 2008, 128, 114102.
10. S. Grimme, Do special noncovalent π–π stacking interactions really exist?, Angew. Chem., Int. Ed., 2008, 47, 3430–3434.
11. C. A. Hunter and J. K. M. Sanders, The nature of π–π interactions, J. Am. Chem. Soc., 1990, 112, 5525–5534.
12. M. R. Battaglia, A. D. Buckingham and J. H. Williams, The electric quadrupole moments of benzene and hexafluorobenzene, Chem. Phys. Lett., 1981, 78, 421–423.
13. C. A. Hunter, K. R. Lawson, J. Perkins and C. J. Urch, Aromatic interactions, J. Chem. Soc., Perkin Trans. 2, 2001, 651–669.
14. J. H. Williams, The molecular electric quadrupole moment and solid-state architecture, Acc. Chem. Res., 1993, 26, 593–598.
15. I. Fleming, Molecular Orbitals and Organic Chemical Reactions, John Wiley and Sons, Chichester, United Kingdom, 2010.
16. D. E. Fagnani, A. Sotuyo and R. K. Castellano, π–π Interactions, in Comprehensive Supramolecular Chemistry II, Elsevier, Oxford, 2017, vol. 1, pp. 121–148.
17. E. M. Cabaleiro-Lago and J. Rodríguez-Otero, σ-σ, and π-π stacking interactions between six-membered cyclic systems. Dispersion dominates and electrostatics commands, ChemistrySelect, 2017, 2, 5157–5166.
18. D. P. Malenov, A. J. Aladić and S. D. Zarić, Stacking interactions of borazine: important stacking at large horizontal displacements and dihydrogen bonding governed by electrostatic potentials of borazine, Phys. Chem. Chem. Phys., 2019, 21, 24554–24564.
19. A. D. Glova, S. V. Larin, V. M. Nazarychev, J. M. Kenny, A. V. Lyulin and S. V. Lyulin, Toward predictive molecular dynamics simulations of asphaltenes in toluene and heptane, ACS Omega, 2019, 4, 20005–20014.
20. G. Gryn'ova and C. Corminboeuf, Steric “attraction”: not by dispersion alone, Beilstein J. Org. Chem., 2018, 14, 1482–1490.
21. K. Patkowski, Recent developments in symmetry-adapted perturbation theory, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1452.
22. M. O. Sinnokrot and C. D. Sherrill, Substituent effects in π–π interactions: sandwich and T-shaped configurations, J. Am. Chem. Soc., 2004, 126, 7690–7697.
23. E. G. Hohenstein and D. C. Sherrill, Effects of heteroatoms on aromatic π–π interactions: benzene–pyridine and pyridine dimer, J. Phys. Chem. A, 2009, 113, 878–886.
24. E. G. Hohenstein, J. Duan and D. C. Sherrill, Origin of the surprising enhancement of electrostatic energies by electron-donating substituents in substituted sandwich benzene dimers, J. Am. Chem. Soc., 2011, 133, 13244–13247.
25. C. D. Sherrill, Energy component analysis of π interactions, Acc. Chem. Res., 2013, 46, 1020–1028.
26. R. M. Parrish and C. D. Sherrill, Quantum-mechanical evaluation of the π–π versus substituent–π interactions in π stacking: direct evidence for the Wheeler–Houk picture, J. Am. Chem. Soc., 2014, 136, 17386–17389.
27. S. E. Wheeler and K. N. Houk, Origin of substituent effects in edge-to-face aryl–aryl interactions, Mol. Phys., 2009, 107, 749–760.
28. S. E. Wheeler and K. N. Houk, Substituent effects in the benzene dimer are due to direct interactions of the substituents with the unsubstituted benzene, J. Am. Chem. Soc., 2008, 130, 10854–10855.
29. S. E. Wheeler and K. N. Houk, Through-space effects of substituents dominate molecular electrostatic potentials of substituted arenes, J. Chem. Theory Comput., 2009, 5, 2301–2312.
30. S. E. Wheeler, Local nature of substituent effects in stacking interactions, J. Am. Chem. Soc., 2011, 133, 10262–10274.
31. S. E. Wheeler, Understanding substituent effects in noncovalent interactions involving aromatic rings, Acc. Chem. Res., 2013, 46, 1029–1038.
32. S. E. Wheeler and J. W. G. Bloom, Toward a more complete understanding of noncovalent interactions involving aromatic rings, J. Phys. Chem. A, 2014, 118, 6133–6147.
33. C. R. Martinez and B. L. Iverson, Rethinking the term “pi-stacking”, Chem. Sci., 2012, 3, 2191–2201.
34. B. W. Gung and J. C. Amicangelo, Substituent effects in C₆F₆–C₆H₅X stacking interactions, J. Org. Chem., 2006, 71, 9261–9270.
35. S. Tsuzuki, T. Uchimaru and M. Mikami, Intermolecular interaction between hexafluorobenzene and benzene: ab initio calculations including CCSD(T) level electron correlation correction, J. Phys. Chem. A, 2006, 110, 2027–2033.
36. W. Wang, Y. Zhang and Y.-B. Wang, Highly accurate benchmark calculations of the interaction energies in the complexes C₆H₆⋯C₆X₆ (X = F, Cl, Br, and I), Int. J. Quantum Chem., 2017, 117, e25345.
37. I. Geronimo, E. C. Lee, N. J. Singh and K. S. Kim, How different are electron-rich and electron-deficient π interactions?, J. Chem. Theory Comput., 2010, 6, 1931–1934.
38. Y. Geng, T. Takatani, E. G. Hohenstein and C. D. Sherrill, Accurately characterizing π–π interaction energies of indole–benzene complexes, J. Phys. Chem. A, 2010, 114, 3576–3582.
39. T. Janowski and P. Pulay, A benchmark comparison of σ/σ and π/π dispersion: the dimers of naphthalene and decalin, and coronene and perhydrocoronene, J. Am. Chem. Soc., 2012, 134, 17520–17525.
40. B. Jeziorski, R. Moszynski and K. Szalewicz, Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes, Chem. Rev., 1994, 94, 1887–1930.
41. K. Szalewicz, Symmetry-adapted perturbation theory of intermolecular forces, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 254–272.
42. E. G. Hohenstein and C. D. Sherrill, Wavefunction methods for noncovalent interactions, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 304–326.
43. G. Jansen, Symmetry-adapted perturbation theory based on density functional theory for noncovalent interactions, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 127–144.
44. I. G. Kaplan, Intermolecular Interactions: Physical Picture, Computational Methods and Model Potentials, John Wiley and Sons, Chichester, United Kingdom, 2006.
45. J. M. Herbert, L. D. Jacobson, K. U. Lao and M. A. Rohrdanz, Rapid computation of intermolecular interactions in molecular and ionic clusters: self-consistent polarization plus symmetry-adapted perturbation theory, Phys. Chem. Chem. Phys., 2012, 14, 7679–7699.
46. K. U. Lao and J. M. Herbert, Accurate and efficient quantum chemistry calculations of noncovalent interactions in many-body systems: the XSAPT family of methods, J. Phys. Chem. A, 2015, 119, 235–253.
47. K. U. Lao and J. M. Herbert, Atomic orbital implementation of extended symmetry-adapted perturbation theory (XSAPT) and benchmark calculations for large supramolecular complexes, J. Chem. Theory Comput., 2018, 14, 2955–2978.
48. J. Hermann, R. A. DiStasio Jr and A. Tkatchenko, First-principles models for van der Waals interactions in molecules and materials: concepts, theory, and applications, Chem. Rev., 2017, 117, 4714–4758.
49. K. Carter-Fenk, K. U. Lao, K.-Y. Liu and J. M. Herbert, Accurate and efficient ab initio calculations for supramolecular complexes: symmetry-adapted perturbation theory with many-body dispersion, J. Phys. Chem. Lett., 2019, 10, 2706–2714.
50. K.-Y. Liu, K. Carter-Fenk and J. M. Herbert, Self-consistent charge embedding at very low cost, with application to symmetry-adapted perturbation theory, J. Chem. Phys., 2019, 151, 031102.
51. E. G. Hohenstein and C. D. Sherrill, Efficient evaluation of triple excitations in symmetry-adapted perturbation theory via second-order Møller–plesset perturbation theory natural orbitals, J. Chem. Phys., 2010, 133, 104107.
52. T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno and C. D. Sherrill, Levels of symmetry adapted perturbation theory (SAPT). I. Efficiency and performance for interaction energies, J. Chem. Phys., 2014, 140, 094106.
53. P. R. Horn, Y. Mao and M. Head-Gordon, Defining the contributions of permanent electrostatics, Pauli repulsion, and dispersion in density functional theory calculations of intermolecular interaction energies, J. Chem. Phys., 2016, 144, 114107.
54. P. R. Horn, Y. Mao and M. Head-Gordon, Probing non-covalent interactions with a second generation energy decomposition analysis using absolutely localized molecular orbitals, Phys. Chem. Chem. Phys., 2016, 18, 23067–23079.
55. Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kús, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock III, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio Jr, H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, C. M. Oana, R. Olivares-Amaya, D. P. O'Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, P. A. Pieniazek, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, N. Sergueev, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, V. Vanovschi, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard III, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer III, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill and M. Head-Gordon, Advances in molecular quantum chemistry contained in the Q-Chem 4 program package, Mol. Phys., 2015, 113, 184–215.
56. S. Ehrlich, J. Moellmann and S. Grimme, Dispersion-corrected density functional theory for aromatic interactions in complex systems, Acc. Chem. Res., 2013, 46, 916–926.
57. J. P. Wagner and P. R. Schreiner, London dispersion in molecular chemistry—Reconsidering steric effects, Angew. Chem., Int. Ed., 2015, 54, 12274–12296.
58. J. K. Cockroft, R. E. Ghosh, J. J. Shephard, A. Singh and J. H. Williams, Investigation of the phase behaviour of the 1:1 adduct of mesitylene and hexafluorobenzene, CrystEngComm, 2017, 19, 1019–1023.
59. J. K. Cockroft, A. Rosu-Finsen, A. N. Fitch and J. H. Williams, The temperature dependence of C–H⋯F–C interactions in benzene:hexafluorobenzene, CrystEngComm, 2018, 20, 6677–6682.
60. J. K. Cockroft, J. G. Y. Li and J. H. Williams, Influence of methyl-substitution on the dynamics of the C–H⋯F–C interaction in binary adducts, CrystEngComm, 2019, 21, 5578–5585.
61. Ohio Supercomputer Center, http://osc.edu/ark:/19495/f5s1ph73.

---

Footnote

† Electronic supplementary information (ESI) available. Additional calculations (including SAPT0 and ALMO-EDA), and details of the vdW model potential. See DOI: 10.1039/d0sc02667k

---