Interplay of π-stacking and inter-stacking interactions in two-component crystals of neutral closed-shell aromatic compounds: periodic DFT study.

This paper bridges the gap between high-level ab initio computations of gas-phase models of 1 : 1 arene-arene complexes and calculations of the two-component (binary) organic crystals using atom-atom potentials. The studied crystals consist of electron-rich and electron-deficient compounds, which form infinite stacks (columns) of heterodimers. The sublimation enthalpy of crystals has been evaluated by DFT periodic calculations, while intermolecular interactions have been characterized by Bader analysis of the periodic electronic density. The consideration of aromatic compounds without a dipole moment makes it possible to reveal the contribution of quadrupole-quadrupole interactions to the π-stacking energy. These interactions are significant for heterodimers formed by arenes with more than 2 rings, with absolute values of the traceless quadrupole moment (Q zz) larger than 10 D Å. The further aggregation of neighboring stacks is due to the C-H⋯F interactions in arene/perfluoroarene crystals. In crystals consisting of arene and an electron-deficient compound such as pyromellitic dianhydride, aggregation occurs due to the C-H⋯O interactions. The C-H⋯F and C-H⋯O inter-stacking interactions make the main contribution to the sublimation enthalpy, which exceeds 150 kJ mol-1 for the two-component crystals formed by arenes with more than 2 rings. This journal is © The Royal Society of Chemistry.

Introduction

Next to conventional and nonconventional hydrogen bonds, the π–π interactions (π-stacking) are the second most studied non-covalent interactions in supramolecular and solid-state chemistry. The big interest in these interactions is mainly raised by their applications in molecular material design,[1] organic semiconductors,[2] biosciences,[3] polymer solar cells[4] and crystallography.[5] There is much debate about the nature of π-stacking interactions.[6] Different classes of tools were developed to translate the results of quantum computations into intuitively useful concepts, which appeal to chemists.[7] Symmetry-adapted perturbation theory (SAPT)[8] is usually used to analyze the energy component of π-stacking in non-periodic systems consisting of molecules with one aromatic ring.[9,10] This analysis dissects the attractions between two monomers into various fundamental physical components, such as electrostatics, London dispersion forces, induction/polarization forces, and exchange–repulsion terms. The geometries of π-stacking are controlled by electrostatic interactions,[11] but the major energetic contribution comes from other factors (dispersion and repulsion).[12-14] The quantum theory of atoms in molecules (Bader analysis)[15] and the non-covalent interaction index[16] visualize and identify π-stacking. These tools use non-periodic[17] or periodic[18] electron densities. Both approaches typically consider gas phase dimers[19] or a pair of neighbouring molecules extracted from the crystals.[20,21] However, a mismatch may arise between the results gained by these two methods, especially when π-stacking determines dimer stabilization.[17,22] The ongoing debate is brought about by a number of reasons. First, the use of the high-level ab initio wave function methods in combination with a large basis set is required to describe π-stacking in non-periodic systems.[23] DFT methods are used in periodic (solid-state) computations of molecular crystals.[24,25] Second, different decomposition schemes (Morokuma–Kitaura,[26] SAPT,[12] SAPT0,[27] the extended transition-state method and the natural orbitals for chemical valence theory (ETS-NOCV)[28]) are applied to analyze the energy component of π-stacking in non-periodic systems, while the PIXEL approach and the classical atom–atom Coulomb–London–Pauli (AA-CLP) model are usually used for molecular crystals.[29,30] Third, for a long time the so-called general arene–arene geometries, namely, T-shaped (edge-to-face), parallel stacked (perfectly stacked) and parallel-displaced (parallel offset or slipped-parallel) structures, were mainly considered in theoretical papers.[9,10,31-33] To the experimental end, the 1 : 1 arene–arene complex existing in the solid-state is more common than the perfectly parallel complexes.[34-36] Such orientation of neighbouring molecules is called “slipped-stacked”[34] or “tilted”.[12] Fourth, only π-stacking interactions are considered in non-periodic calculations, while inter-stacking interactions play a structure-forming role.[29] Finally, a variety of systems was considered, namely: neutral closed-shell aromatic and heteroaromatic compounds,[30] charge–transfer complexes,[6] systems consisting of a charged molecule[37] or a molecule in the electronically excited state,[38] and metal complexes with aromatic ligands.[39] Obviously, the electrostatic contribution to the π-stacking energy in such systems is caused by interactions of various types (Coulomb, dipole–dipole, etc.).[40]

There are numerous two-component (binary) crystals composed of two organic compounds unable to form conventional hydrogen bonding.[5,6,29,30,34] The association modes mostly include stacking of flat systems, one of them usually being an aromatic hydrocarbon, e.g. arene.[30] It is an electron-rich molecule characterized by Qzz < 0, where Qzz is the z-component of the traceless quadrupole moment.[26,27] Another compound is an electron-deficient molecule with Qzz > 0, e.g. perfluoroarene. If organic crystals consist of neutral closed-shell aromatic compounds without a dipole moment, electrostatic contribution to the π-stacking energy is associated with quadrupole–quadrupole interactions. The parallel-displaced structures of the perfluoroarene–arene heterodimer are consistent with a simplified model. In this model the interaction is favoured by electrostatic attraction when the charge distributions are approximated by central quadrupoles of the opposite phase.[41] The considered crystals include hexafluorobenzene/arene and perfluoronaphthalene/arene families in which the arene varies from benzene to pyrene (Table 1 in ref. 30). This made it possible to quantify the influence of the arene size (number of rings) on the energy of the quadrupole–quadrupole interactions.

Essential features of the present theoretical study are the following. (i) Checking the level of approximation used in the periodic (solid-state) DFT calculations by comparing the calculated values of the enthalpy of sublimation with the published data. (ii) Identification and quantification of the pattern of intermolecular (non-covalent) interactions using Bader analysis of periodic electron density. (iii) Identification of the role of quadrupole–quadrupole interactions from “first principles”, that is, without using the wave functions of the crystals under consideration. As a result, the physical nature of the electrostatic interaction among flat-ring moieties of two aromatic molecules in the 1 : 1 stoichiometry is revealed and the interplay of π-stacking and inter-stacking interactions in the two-component organic crystals without conventional hydrogen bonds is clarified.

Methods

Crystals under study

In order to reveal the contribution of the quadrupole–quadrupole interactions to the π-stacking energy, a family of the two-component (binary) crystals 1 : 1 were considered. These crystals are made up of neutral closed-shell aromatic compounds without a dipole moment. One component is an electron-rich aromatic molecule (Scheme 1) which is characterized by Qzz < 0 (Table S1†). Another component is an electron-deficient aromatic (electron-poor) molecule (Scheme 2) with Qzz > 0 (Tables S1 and S2†). The stacking motif of these crystals carries over to asymmetric units made of a heterodimer.[29] Charge-transfer complexes, for example, 1,3,5-trinitrobenzene cocrystals,[6] are beyond the scope of the study.

Scheme 1

The electron-rich molecules under study. Molecule designations (the Cambridge Structural Database Refcodes) are given in parentheses.

Scheme 2

The electron-deficient molecules under study. Molecule designations (the Cambridge Structural Database Refcodes) are given in parentheses.

Periodic DFT computations

Periodic DFT computations were conducted using the CRYSTAL17 software.[42] The PBE functional[43] was employed with an all-electron Gaussian-type localized orbital basis set 6-31(F+)G**. Diffuse functions, with the exponent factor equal to 0.1076, were added to the fluorine atoms.[44] The exponent is larger than the critical value (0.06 (ref. 45)) below which problems in the SCF procedure arise. Dispersion interactions were taken into account by using the semi-empirical D3 scheme.[46] PBE-D3/6-31(F+)G** provides an adequate description of the metric and electron-density parameters of F⋯H interactions in crystals of fluorinated compounds.[47] The space groups and the unit cell parameters of the considered crystals obtained in the experimental studies were fixed for the calculations. This is a common approximation for the calculating of molecular crystals.[25,47-50] Further details of periodic DFT computations are given in ESI.†

A common approach to Elatt estimation for two-component crystals is based on the total electronic energies of the crystal (Ebulk) and isolated molecules of its components (Ei) calculated per asymmetric unit:[50]where Z is a number of molecules in the unit cell, ai is a number of symmetrically equivalent molecules of each type in the unit cell. The basis set superposition error is taken into account by Boys–Bernardi counterpoise correction.[51]Elatt calculation details are given elsewhere.[52]

According to ref. 53, the theoretical value of ΔHsub at given temperature T is calculated as:

Identification and quantification of intermolecular (non-covalent) interactions

Bader analysis of periodic (crystalline) electron density[54] was used for identification and quantification of the pattern of intermolecular (non-covalent) interactions (IMIs) in the considered crystals. Within this approach, the particular IMI is associated with the existence of the bond path (i.e., the bond critical point) between the pair of atoms belonging to different molecules. The network of the bond paths yields the comprehensive bond picture; the energetic of each specific interaction is considered totally independent of others. The effects of crystal environment, long-term electrostatic, etc. are taken into account implicitly, via the periodic electronic wave function, and are coded in the bond critical point features. Only bond critical points with the values of the electron density ρb above the threshold 0.003 a.u. were considered.[55] Details of Bader analysis of periodic electron density are given elsewhere.[56]

The energy of the particular IMI, Eint, was calculated as:[57]where Gb is the positively-defined local electronic kinetic energy density at the bond critical point of the considered IMI. Indices j and i denote the atoms belonging to different molecules.

Eqn (3) with k = 1124 yields reasonable Eint values for conventional/non-conventional H-bonds and π-stacking interactions in molecular crystals.[58,59] The accuracy of the obtained values is several kJ mol−1 (Table 1 in ref. 60). To estimate the energy of F⋯F, Br⋯Br and I⋯I interactions in crystals, the coefficient value in eqn (3) must be changed.[47,48,61]k = 338 was used to compute the energy of F⋯F interactions[47] in the present study.

Eqn (3) is used for Elatt estimating in crystalline materials (see ref. 62 for examples). In this approach, Elatt is presented as a sum of energies of all IMIs within the asymmetric unit:[63]where Eint,j,i is evaluated using eqn (3). The applicability of eqn (4) for estimating the Elatt values of the considered crystals is justified by comparing them with the results obtained by eqn (1).

Estimation of quadrupole–quadrupole interaction

A general expression for the quadrupole–quadrupole interaction between two molecules is given in ref. 64 and 65, but it is too bulky. A simpler expression was obtained for axisymmetric molecules:[66]where are unit vectors as defined in the diagram (Fig. 1).

Fig. 1

The system of coordinates used for the description of the orientation of molecules forming a heterodimer. z1 and z2 denote the axis of symmetry of two molecules (D6h is the case of benzene and C6F6; D2h for the remaining molecules). R is the distance between their centers of mass.

In the case of the face-to-face configuration, i.e. the perfectly stacked arrangement of the rings: , eqn (5) derives to the following expression:where Q1 and Q2 stand for the traceless quadrupole moment Qzz of molecules forming heterodimer.

Details of the quadrupole–quadrupole interaction calculations are given in ESI.†

SAPT analysis

Decomposition of interaction energy into particular components was carried out within the framework of SAPT.[8] First, the structures of the three dimers, (C10F8/naphthalene, C10F8/anthracene) and (C10F8)2, were taken from the crystal structures and optimized at the dispersion-corrected DFT level of theory with PBE functional,[43] empirical D3 dispersion correction[46] and def2-TZVPP basis set,[67] denoted hereafter as PBE-D3/def2-TZVPP. The optimizations were repeated with the same basis set and the Møller–Plesset MP2 perturbative calculus, MP2/def2-TZVPP, with core (1s) electrons excluded from the perturbational treatment. Resolution of identity (RI) acceleration was used for DFT and MP2 as well. These calculations were carried out with the Turbomole 6.5 software.[68] The subsequent SAPT analysis using Psi4 rev. 1.3.2 program[69] employed aug-cc-pVDZ basis set,[70] following the recommendations to use the balanced correlation-consistent basis sets for interaction energy partitioning[71] and selecting the double-zeta valence variant as computationally tractable. The perturbative SAPT scheme was expanded to the highest available level, SAPT2 + 3δMP2.[72]

Results and discussions

PBE-D3/6-31(F+)G** method gives a reasonable description of metric parameters of the considered two-component crystals. In particular, it reproduces the distance between the centres of mass of electron-rich and electron-deficient compounds forming a heterodimer (Table 3). It also provides an adequate description of the metric, electron–density parameters and energy of intermolecular C–H⋯F–C/C–H⋯O interactions in organic crystals.[47] The ΔHsub values of the considered crystals are in reasonable agreement with the available literature data (Table 1). The calculated ΔHsub values are increasing alongside with the growing number of condensed rings, which is in agreement with the literature.[82,83] The calculated ΔHsub values of the two-component crystals formed by the molecules with an extended π-system exceed 150 kJ mol−1, which is consistent with the results of the calculations of the crystalline arene–perhaloarene adducts[49] and two-component crystals of pharmaceutical use.[84,85]

The quadrupole–quadrupole interaction energy Uavs. the sum of Coulomb (ECoul) and polarization (Epol) energy b of stacked heterodimers extracted from the experimental crystal structures of arene–perfluoroarene crystals. R denotes for the distance between the centers of mass of the molecules forming the heterodimer

Heterodimerc	U Q 1 Q 2 , kJ mol−1	E Coul + Epol, kJ mol−1	R, Åd
C6F6/benzene	−3.2	−17	3.763 (3.763)
C6F6/naphthalene	−31.1	−24	3.435 (3.430)
C6F6/anthracene	−34.5	−27	3.619 (3.6185)
C6F6/pyrene	−38.1	−26	3.473 (3.473)
C10F8/naphthalene	−14.7	−25	3.729 (3.7285)
C10F8/anthracene	−79.2	−33	3.405 (3.405)
C10F8/pyrene	−92.5	−43	3.362 (3.3625)

Evaluated using eqn (6); Q1 and Q2 were calculated at the PBE-D3/6-31(F+)** level.

Table 1 in ref. 30.

Abbreviations used to refer electron-deficient molecules are defined in Scheme 2.

R values for heterodimers extracted from the optimized structures of the considered crystals obtained from periodic DFT are given in parentheses.

The comparison of theoretical values of the sublimation's enthalpy ΔHsuba of the considered crystalsb with the available literature data. The units are kJ mol−1

Crystal (ref.)	ΔHsub
	Eqn (2), this work	Ref. 29 c
C6F6/benzene[73]	91	90
C6F6/naphthalene[74]	121	118
C6F6/anthracene[75]	138	143
C6F6/pyrene[75]	151	145
C10F8/naphthalene[76]	147	146
C10F8/anthracene[77]	165	171
C10F8/pyrene[77]	173	173
C10F8/diphenylacetylene[78]	162	165
C14F10/diphenylacetylene[79]	174	184
PYDMAN/naphthalene[80]	189	—

The Elatt values used in the ΔHsub evaluation were calculated using eqn (1).

The abbreviations used to refer electron-deficient molecules are defined in Scheme 2.

There are no experimental determinations of ΔHsub available for the considered two-component crystals.[81] Following,[29] ΔHsub is the sum of the sublimation enthalpies of the corresponding hydrocarbons.[82,83]

We can conclude that PBE-D3 with the 6-31(F+)G** basis set (where diffuse functions are added to the fluorine atoms) provides an adequate description of intermolecular interactions in the two-component aromatic crystals without conventional hydrogen bonds.

Intermolecular interactions in considered crystals

Bader analysis of the periodic electron density is used to identify and visualize the pattern of IMIs in the considered crystals. The presence of a bond path between two nuclei belonging to different molecules means that there is a IMI between them.[86-88] The pair of nuclei connected by the bond path define the particular type of the IMI denoted as H⋯C, H⋯H, etc. (Fig. 4 in ref. 86). The presence of a short interatomic distance does not imply the existence of an IMI, that is, a bond critical point.[47] For example, in crystalline naphthalene[89] the H(4) atom forms two short distances with the carbon atoms of a neighbouring molecule (Fig. S1†). However, Bader analysis of the periodic electron density shows only one H(4)···C IMI. Our results (the number of IMIs and their characteristics) are consistent with published data[86] (Table S4 and Fig. S1†).

The studied two-component crystals contain flat and stacked heterodimers that form infinite stacks (columns). IMIs between molecules in the given stack are called as π-stacking. They are associated with IMIs of the C⋯C, C⋯F types (Fig. 2 and S2–S4†). The obtained results are in general agreement with the literature data.[17,40] IMIs between infinite stacks are called as inter-stacking. In the considered crystals, an electron-deficient molecule from a given stack is surrounded by electron-rich molecules from neighboring stacks and vice versa (Fig. S5 and S6†). Thus, the inter-stacking interactions are mostly associated with IMIs of the H⋯F type in the C6F6/arene and C10F8/arene rows. Other types of IMIs (H⋯H, F⋯F etc.) are realized very rarely in the considered crystals. It should be noted that the F⋯F interactions are in a way controversial.[90,91]

Fig. 2

The pattern of IMIs in the heterodimers, extracted from the corresponding crystals. Topology of electron density [bond critical points (3, −1) are depicted as red dots, ring critical points (3, +1) as yellow, and cage critical points (3, +3) as green].

Bader analysis of the periodic electronic density followed by eqn (3) gives a possibility to reveal relative contributions of π-stacking and inter-stacking interactions to the Elatt energy.[62,84] The Elatt computation in a two-component 1 : 1 crystal is carried out in two separate steps.[92] First, the interaction energy in the heterodimer is computed. The π-stacking energy is twice as that energy. Then, the interaction energy of the heterodimer with the molecules from neighboring stacks is calculated. To do this, the number of IMIs per heterodimer should be counted and listed. In the case of arene⋯arene, C6F6⋯C6F6etc. interactions, half of the corresponding IMIs numbers should be taken into account (eqn (4) from ref. 92). The Elatt values of the considered crystals evaluated using eqn (4) are compared with the values obtained using eqn (1) in Table 2. The lattice energies calculated in those two ways are close to each other. This justifies the results on the relative contributions of π-stacking and inter-stacking interactions to the Elatt energy of the considered crystals.

The Elatt values of the two-component crystals evaluated using eqn (4) and relative contributions of the π-stacking and inter-stacking interactions to it and the types of IMIs

Crystal	E latt, kJ mol−1a	Contribution, %
		π-stackingb	Inter-stackingb
C6F6/benzene	−98 (−86)	24 (C⋯C, C⋯F)	76 (C–H⋯F, F⋯F)
C6F6/naphthalene	−110 (−116)	28 (C⋯C, C⋯F)	72 (C–H⋯F, F⋯F, H⋯H)
C6F6/anthracene	−101 (−133)c	13 (C⋯C)	87 (C–H⋯F, F⋯F, H⋯H)
C6F6/pyrene	−129 (−146)	31 (C⋯C, C⋯F)	69 (C–H⋯F, H⋯H)
C10F8/naphthalene	−142 (−142)	24 (C⋯C, C⋯F)	76 (C–H⋯F)
C10F8/anthracene	−177 (−160)	37 (C⋯C)	63 (C–H⋯F, F⋯F, H⋯H)
C10F8/pyrene	−155 (−168)	34 (C⋯C)	66 (C–H⋯F, F⋯F, H⋯H)
C10F8/diphenylacetylene	−149 (−157)	30 (C⋯C, C⋯F)	70 (C–H⋯F, F⋯F, H⋯H)
C14F10/diphenylacetylene	−169 (−174)	28 (C⋯C, C⋯F)	72 (C–H⋯F, F⋯F, H⋯H)
PYDMAN/naphthalene	−147 (−152)	26 (C⋯C)	74 (C–H⋯O, O⋯O)

E latt values evaluated using eqn (1) are given in parenthesis.

The types of IMIs are given in parenthesis.

The accuracy of the Elatt values, evaluated using eqn (1) and (4), does not exceed ±15 kJ mol−1.[92,93] In a C6F6/anthracene crystal, these values differ by ∼ 30 kJ mol−1. It is caused by a number of bond critical points with ρb value below the threshold due to longer C–H⋯F and F⋯F distances.

In the C10F8/arene row the π-stacking energy is much larger than that the C6F6/arene row (Table 2). The growing number of condensed rings leads to an increase in π-stacking energy. This is consistent with an increase of Qzz of molecules forming the heterodimer (Tables S1 and S2†). The contribution of the inter-stacking interactions to the lattice energy is about 70 percent (Table 2). The growing number of condensed rings leads to an increase in the total energy of inter-stacking interactions due to an increase in the number of C–H⋯F interactions between adjacent stacks, cf. Fig. S5 and S6.† As a result, the Elatt energy of the two-component organic crystals, consisting of the electron-rich and electron-deficient compounds, is much larger than that of the corresponding one-component crystals.[82,83] The reason for this difference is the presence of the inter-stacking C–H⋯F interactions in the perfluoroarene-arene crystals. The energy of these interactions is much greater than for the H⋯H and F⋯F ones.[47] The latter are realized in one-component crystals constituted by an arene or perfluororene. A similar conclusion can be made for crystals consisting of arene and an electron-deficient compound such as pyromellitic dianhydride (Tables 1 and 2). Aggregation of neighbouring stacks in these crystals is mainly due to the C–H⋯OC interactions which energy is larger than the C–H⋯F–C interactions.[47]

In order to validate the employed periodic DFT method and basis set, the optimization of the cell parameters was also performed without cell volume restrictions. The symmetry of the considered crystals was kept during these computations. Relative changes in volume are rather small (Table S5†) in accordance with published data for organic crystals.[94] Stacking and C–H⋯F distances change slightly when cell parameters are relaxed during optimization. Lattice energies calculated by eqn (4) for crystals with optimized unit cell volumes turn out to be slightly lower than the Elatt values obtained for the non relaxed cell parameters. For example, Elatt equals to −97.7 (−98.0) kJ mol−1 for C6F6/benzene and to −169.2 (−177) kJ mol−1 for C10F8/anthracene. The Elatt values evaluated for the non relaxed cell parameters are shown in parentheses.

Summing up, the three-dimensional structure of the considered two-component organic crystals is determined by the interplay of π-stacking and inter-stacking interactions. The energy of both contributions depends on the number of condensed rings of molecules forming the two-component crystals. In the case of inter-stacking interactions this is due to the size of the molecules and the number of accessible contact sites. In the case of π-stacking interactions this is due to the increase of Qzz of molecules forming the heterodimer.

The role of the quadrupole–quadrupole interaction in formation of heterodimers

To identify the contribution of the quadrupole–quadrupole interaction to the dimerization energy Edim we considered heterodimers not crystals. The analysis of this interaction makes sense in either small molecules or molecules that are very far apart. In molecules that are almost on top of each other (Fig. 2), it is known that the multipolar expansion needs not converge, which implies the use of approximate schemes, e.g.eqn (5) and (6). The results obtained are semi quantitative nature; however, they reveal the nature of electrostatic interactions in heterodimers, that is, in π-stacking.

At first step, stacked heterodimers extracted from the experimental crystal structures of the corresponding two-component crystals (Table 3). Such (hetero)dimers are often used to estimate the intermolecular interaction energy in solid state.[95-97] The quadrupole–quadrupole interaction energy of the C6F6/arene and C10F8/arene heterodimers were compared with the electrostatic interaction energy[98] in these systems evaluated using the PIXEL scheme.[30] This comparison highlights trends of the quadrupole–quadrupole interaction energy. Since the R distance varies in the considered rows, the U value changes non-monotonously. However, its absolute value increases with an increase in the number of condensed rings in the molecules forming the heterodimer, i.e. a rapid growth in the Qzz value (Table S1†).

At the second step, SAPT is used to analyze the energy components of π-stacking in several (hetero)dimers formed by octafluoronaphthalene (C10F8) (energy components of π-stacking in heterodimer structures formed by C6F6 have been already studied[17,27,99]). The data summarized in Table 4 are grouped according to the principal terms derived from SAPT: Eel is the electrostatic energy of interaction of the nuclear charges and monomer electron densities; Eexch is the Pauli repulsion of the electron clouds; Eind includes the effects of mutual polarization (hence the name: induction term), including corrections to the Pauli exchange due to orbital relaxation; Edisp contains dispersion energy and its effect on the Pauli repulsion. These terms sum up to the total interaction energy of the dimer Edim.[8,71] Returning to Table 4, one can appreciate that MP2 optimization resulted in narrower separation of the monomers with appropriate increase in the Pauli repulsion, and this makes the SAPT Edim slightly worse (by 3 to 5 kJ mol−1) with respect to the PBE-D3 structures. However, the changes in the particular terms of the interaction energy are far more dramatic, with emphasis on the Pauli repulsion growing ca. three times. This effect is due to a deficiency of MP2 for π-interactions resulting from overestimation of the dispersion energy at the MP2 level.[72] In view of this, it seems safer to analyse the SAPT2 + 3δMP2 results for the PBE-D3 structures. The multipole–multipole interaction is represented in the Eel term, and for the studied highly-symmetric aromatic molecules the multipole expansion is dominated by quadrupoles. The C10F8 homodimer benefits significantly less (ca. four times) from the electrostatic attraction than the heterodimers. This phenomenon results from the diversity of patterns of electronic density distribution: the perfluoronaphthalene has an outer ring of negative fluorine atoms and an inner ring of positive aromatic electrons (negative −I inductive effect). In the arene, the outer ring of hydrogen atoms is positive, while the aromatic sextets gather negative charge. On the other hand, the studied molecules are of very low polarizability, which is manifested by the Eind values of only −1.1 to −4.3 kJ mol−1 (again with the homodimer gaining the least stability). The situation is slightly different with dispersion forces. In this case, what matters is the sheer extent of the electronic cloud that can form instantaneous multipoles. Thus, the best stabilization by Edisp is recorded for the C10F8/anthracene complex; the C10F8/naphthalene and (C10F8)2 dimers have similar Edisp, but the latter complex wins due to the electron-rich fluorine lone pair regions. Finally, when one looks at the ordering of the Edim values and compares it with the electrostatic term Eel, it is possible to conclude that, from the SAPT point of view, weakening of quadrupole–quadrupole interaction in the homodimer makes overall π-stacking the least efficient in the homodimer among the three studied systems.

SAPT2 + 3δMP2a decomposition of the dimerization energy Edim (kJ mol−1) for the C6F6/benzene, C10F8/naphthalene, C10F8/anthracene and (C10F8)2 dimer

(Hetero)dimer	R b, Å	E el	E exch	E ind	E disp	E dim	U Q 1 Q 2 c
C6F6/benzened	3.600	−20.9	32.0	−3.6	−31.9	−24.3	−10.7
C10F8/naphthalene	3.682	−28.5	50.2	−2.9	−69.3	−50.5	−21.6
	3.282	−67.2	140.7	−6.1	−114.6	−47.2	−39.3
C10F8/anthracene	3.473	−35.9	64.3	−4.3	−91.8	−67.7	−71.8
	3.099	−85.8	181.8	−9.0	−149.6	−62.6	−126.7
(C10F8)2	3.740	−8.30	40.9	−1.1	−72.9	−41.4	16.2
	3.358	−40.7	128.6	−1.4	−123.7	−37.2	25.9

Dimer structures were optimized at the PBE-D3/def2-TZVPP and MP2/def2-TZVPP levels. The MP2 data are given in bold.

The distance between the centers of mass of the molecules in the dimer.

The quadrupole–quadrupole interaction energy evaluated using eqn (5); Q1 and Q2 were calculated at the PBE-D3/6-31(F+)** level.

The geometry of the PD structure and SCS-SAPT0/aVDZ decomposition of the dimerization energy was borrowed from ref. 27.

According to SAPT analysis, the nature of the interaction in C10F8/naphthalene and (C10F8)2 is different. The Edisp energy is close, while the Eel energy in (C10F8)2 is on ∼20 kJ mol−1 lower, than in C10F8/naphthalene. The quadrupole–quadrupole interaction is positive in the relaxed gas-phase (C10F8)2 dimer (Table 4). To reveal the role of the quadrupole–quadrupole interaction, configurations of these dimers in gas-phase (the relaxed PBE-D3 structures) are compared with the structures in the crystalline state (Fig. 3). The mutual orientation of the molecules in the heterodimer practically does not change upon the transition from gas to crystal. The distance between the centers of mass of the molecules R in gas is a bit longer than in the crystal, 2.68 vs. 2.43 Å, respectively. Reducing the distance leads to an increase in the absolute U value of energy by 10 kJ mol−1, cf.Tables 3 and 4. The mutual orientation of the molecules in the homodimer in the gas phase is very different from the orientation in the crystal (Fig. 3). The R value in gas is much shorter than in the crystal, 3.74 vs. 5.00 Å, cf.Tables 4 and S6†.

Fig. 3

The mutual orientation of the molecules in the C10F8/naphthalene and (C10F8)2 dimers extracted from the crystal (element color) and after structure relaxation in gas phase (turquoise).

Next, we studied the applicability of eqn (5) in the description of the intermolecular interaction energy in the slipped-stacked configuration of homodimers existing in the single-component crystals of C10F8,[100]C14F10 (ref. 101) and octafluorobiphenylene C12F8 (ref. 102) (Fig. S7†).

The homodimers are characterized by small negative values of U (Table S6†). The negative values of U are due to a geometric effect arising from the relative shift of the molecules in the homodimer. This effect can be examined using (C10F8)2. The octafluoronaphthalene molecules are parallel to each other, but are strongly shifted, so that the cosines of the angles between the normals to the planes and R are less than 1: −0.654 and 0.654. Substitution of these values in eqn (5) gives: [1 − 2.14 − 2.14 + 2 − 8.55 + 6.39] = [−3.44]. As a result, this equation gives a small negative value for the quadrupole–quadrupole interaction in the (C10F8)2 dimer, which structure was extracted from the corresponding single-component crystal.[100]Eqn (5) leads to the positive U value for the gas-phase structure of the (C10F8)2 dimer (Table 4). Thus, the slipped-stacked arrangement of molecules in the single-component crystals is stabilized, in particular, due to the weak quadrupole–quadrupole interaction. This effect is considered in more detail in ref. 11 and 103.

Conclusions

PBE-D3 with the 6-31(F+)G** basis set (where diffuse functions are added to the fluorine atoms) provides an adequate description of metric parameters of the two-component crystals formed by electron-rich and electron-deficient aromatics. It also gives reasonable values of the sublimation enthalpy. We conclude that the PBE-D3/6-31(F+)G** approximation provides a grounded trade-off between the accuracy and the rate of calculations of intermolecular interactions in the two-component aromatic crystals without conventional hydrogen bonds.

The theoretical sublimation enthalpy ΔHsub of the considered two-component crystals is much larger than that of the corresponding single-component crystals. The reason for this is the presence of the intermolecular C–H⋯F interactions between neighbouring stacks (columns) in the two-component crystals. The ΔHsub values of the crystals formed by molecules with an extended π-system exceed 150 kJ mol−1.

The increase in the number of condensed rings in the molecules forming the heterodimer causes a rapid growth in the energy of the quadrupole–quadrupole interactions. The contribution of these interactions to the π-stacking energy dominates in arenes with more than 2 rings. The absolute value of the traceless quadrupole moment Qzz in these compounds exceeds 10 D Å. The electron-deficient heteroaromatic compounds (pyromellitic dianhydride, tetracyanoquinodimethane etc.) are characterized by large positive Qzz values (>10 D Å). All these facts point to a greater influence of the quadrupole–quadrupole interaction in the heterodimers formed by these electron-rich and electron-deficient compounds.

The results of this work allow us to propose the following scenario for the formation of two-component crystals formed by electron-rich and electron-deficient aromatics. First, heterodimers are formed. The quadrupole–quadrupole interactions determine the orientation of molecules in the heterodimer and make a significant contribution to the energy of π-stacking interactions. Then stacks (columns) are formed. Further aggregation of neighbouring stacks (columns) is due to the C–H⋯F interactions in the case of arene/perfluoroarene and the C–H⋯O interactions in the case of arene/electron-deficient molecule such as pyromellitic dianhydride. The inter-stacking C–H⋯F and C–H⋯O interactions make the main contribution to the sublimation enthalpy of the two-component aromatic crystals without conventional hydrogen bonds.

Conflicts of interest

There are no conflicts to declare.