Interaction energies and stabilities of heteroatom-containing aromatic compounds on graphite

Abstract

Stacked structures based on π–π stacking interactions have been widely observed in natural and artificial self-assembled systems composed of polycyclic aromatic compounds containing carbon and heteroatoms. While previous studies have quantitatively evaluated the interaction energies (E_int) and their variations with horizontal displacement (ΔE_int) for polycyclic aromatic hydrocarbons (PAHs), such intermolecular interactions for heteroatom-containing aromatic compounds (HACs) incorporating elements such as sulphur, oxygen, or nitrogen remain largely unexplored. In this study, we evaluated E_int and ΔE_int of HACs on a graphite model surface (C₉₆H₂₄) using dispersion-corrected density functional theory calculations. The E_int values of thiophene (−10.09 kcal mol⁻¹) and pyridine (−10.31 kcal mol⁻¹) were found to be comparable to those of benzene (−10.55 kcal mol⁻¹), whereas those of furan (−8.13 kcal mol⁻¹) and pyrrole (−9.18 kcal mol⁻¹) were lower (less negative). In all cases, the E_int values were smaller than those of n-alkanes containing the same total number of carbons and heteroatoms. The dispersion force plays the dominant role in the attractive interactions between the HACs and C₉₆H₂₄. It was found that E_int increases linearly with the number of heavy (non-hydrogen) atoms of HACs and benzene-fused HACs, similar to the results obtained for PAHs. The stability against horizontal displacement of thiophene, furan, pyrrole and pyridine, as inferred from the changes in ΔE_int was also similar to that of benzene, but lower than that of n-pentane. These results indicate that the presence of heteroatoms does not significantly affect the interaction energies or stacking stabilities, implying that the introduced HACs can be treated almost identically to PAHs in terms of their interaction energies and stability of stacking with graphite.

---

Introduction

Polycyclic aromatic hydrocarbons (PAHs) composed solely of fused benzene rings form π-stacked structures which play key roles in various fields, including structural biology, materials science and molecular electronics.¹ Dispersion forces are the primary source of attractive interactions of π–π stacking,² and the size of the aromatic ring affects the magnitude of the intermolecular interaction energy (E_int).³ The E_int per carbon atom has been reported to be between approximately −1.2 and −1.4 kcal mol⁻¹, with the data obtained by theoretical calculations showing good agreement with the results of temperature-programmed desorption (TPD) experiments.³ Heteroatom-containing aromatic compounds (HACs), such as those incorporating sulphur, oxygen or nitrogen atoms, also exhibit π-stacking which affects the geometries of nano-architectures and molecular recognition of receptors.⁴ While the π–π interactions of PAHs have been investigated extensively, there have been limited studies on the intermolecular interactions between PAHs and HACs.⁵

Gómez and Martínez-Magadán^5a investigated the adsorption of dibenzothiophene on (7,7) and (10,5) carbon nanotubes (CNT) using density functional theory (DFT) calculations at the BLYP/DZVP level of theory and the calculated adsorption energies suggested that physisorption of dibenzothiophene on (7,7) CNT was preferable to that on (10,5) CNT. In an experimental study, Goering and Burghaus^5b investigated the adsorption and desorption of thiophene on single-walled CNT via TPD. They reported a binding energy of 14.3 ± 0.5 kcal mol⁻¹; however, it was pointed out that this value may be influenced by the effects of impurities.^5c,d Denis and Iribarne^5d studied the adsorption of thiophene inside and outside of single-walled CNT and on graphene by DFT calculations using van der Waals (VDW-DF) and LDA functionals. On graphene, parallel adsorption was preferred, and the E_int obtained by VDW-DF calculations was −8.9 kcal mol⁻¹. Huber et al.^5e performed dispersion-corrected density functional theory (DFT-D) calculations with the ωB97XD functional for benzene complexes with 13 different five- and six-membered aromatic heterocycles containing sulphur, oxygen or nitrogen as heteroatoms. They reported that the optimal structure had a parallel-displaced geometry and that dispersion forces strongly contributed to the overall interaction energy.

In addition to the above examples, the π-stacked structure of heteroatom-containing polycyclic aromatic compounds (hetero-PACs) can also be found in the “bee structure” of asphaltenes⁶ and in self-assemblies on highly oriented pyrolytic graphite (HOPG).⁷ Asphaltene is composed of large aromatic species such as PAHs as well as hetero-PACs, and it has been reported that disruption of the π-stacked structure of asphaltene (the “bee structure”) by mixing with a saccharide-derived compound significantly alters the physical properties and surface morphology of asphalt.⁸ Self-assemblies of HACs and hetero-PACs containing alkyl chains have been investigated by scanning tunnelling microscopy (STM) at the HOPG/solvent interface. Numerous STM studies have demonstrated that sufficiently strong adsorption and positional stability of adsorbates are essential for the stable and persistent formation of physisorbed monolayers.⁹ Further, electric properties of graphene-based molecular devices is governed by several factors such as the ordering, orientation, and adsorption stability of the constituent molecules.¹⁰ Therefore, it is essential to clarify the differences in the adsorption properties of PAHs and HACs, which are widely employed as building blocks in molecular devices on graphene surfaces. Because HACs contain heteroatoms, their electrostatic interactions are expected to differ from those of PAHs, potentially leading to different adsorption behaviours. However, the actual effect of introducing heteroatoms on electrostatic interactions has not yet been fully understood. Moreover, to the best of our knowledge, neither the energy of the intermolecular interactions between HACs and large PAHs such as graphene nor the resistance to lateral movement of HACs on graphene has been systematically investigated. Thus, the effect of introducing heteroatoms on electrostatic interactions and resulting adsorption properties has not yet been elucidated.

In this study, we performed dispersion-corrected DFT calculations to evaluate the energy of the interaction (E_int) between HACs and a graphite model surface (C₉₆H₂₄), and the changes in E_int associated with the horizontal displacement of HACs (ΔE_int). The E_int values of HACs such as thiophene, furan, pyrrole, pyridine, and their fused derivatives (benzene-fused HACs), in which one or two benzene rings are fused to the heteroaromatic core, were evaluated and ΔE_int values were obtained by migrating thiophene, furan, pyrrole, and pyridine along the x- and y-directions of C₉₆H₂₄. Additionally, the E_int and ΔE_int values of the HACs were compared with those of PAHs and n-alkanes to gain insight into the energetic and structural features of the π-stacked assemblies of HACs and hetero-PACs.

Computational details

The Gaussian 16 program¹¹ (Gaussian, Wallingford, CT, USA) was used for the DFT and ab initio calculations. The geometries of coronene, C₉₆H₂₄ and their complexes were optimised at the B3LYP/6-31G* level of theory¹² with Grimme's D3 dispersion correction¹³ (B3LYP-D3/6-31G*). First, the geometries of isolated coronene and C₉₆H₂₄ were optimised. Then, these geometries were fixed in the geometry optimisations of the complexes of coronene and C₉₆H₂₄ with HACs and benzene-fused HACs. The E_int for the coronene and C₉₆H₂₄ complexes with HACs or benzene-fused HACs were calculated at the B3LYP-D3/6-311G** level (E_B3LYP-D3) using the optimised geometries, unless otherwise noted. In the geometry optimisation of complexes and calculations of their interaction energies, the basis set superposition error (BSSE)¹⁴ was corrected using the counterpoise method.¹⁵ Molecular polarisabilities were calculated at the B3LYP/aug-cc-pVDZ level of theory using B3LYP-D3/6-31G* level-optimised geometries.

E_int calculated for the coronene complexes with HACs at the B3LYP-D3 level of theory (E_B3LYP-D3) was compared with the estimated CCSD(T) level E_int¹⁶ at the basis set limit (E_{CCSD(T)(limit)}) to confirm the accuracy of the E_B3LYP-D3. Note that E_{CCSD(T)(limit)} values have been reported to be close to the experimental E_int values in the gas phase.¹⁷ The optimised geometries and E_B3LYP-D3 values calculated for the coronene complexes with HACs are shown in Fig. S1. E_{CCSD(T)(limit)} were obtained according to the previously reported method.¹⁸ Briefly, MP2-level E_int at the basis set limit (E_MP2(limit)) was estimated from the MP2-level¹⁹ E_int obtained using the aug-cc-pVDZ and aug-cc-pVTZ basis sets, according to the method of Helgaker et al.²⁰ Then, E_{CCSD(T)(limit)} was estimated from the E_MP2(limit) values and the CCSD(T) correction term (difference between CCSD(T) and MP2 E_int values), which was obtained from CCSD(T) calculations using the 6-31G* basis set. As shown in Table S1, the values of E_B3LYP-D3 and E_{CCSD(T)(limit)} were very close to each other within an error of less than 0.6 kcal mol⁻¹, indicating that B3LYP-D3/6-311G** calculations are sufficiently reliable for E_int evaluation.

The change in E_int associated with the horizontal displacement (ΔE_int) of aromatic molecules along x- and y-axes (Fig. 1) was calculated according to the following equation:^3d,21

ΔE_int = E_int − E_int(0),

where E_int(0) is E_int at the stable position. The optimised geometries of the isolated molecules were used for calculations. The Psi4 program²² was used for the symmetry-adapted perturbation theory calculations. Each energy term was calculated by the SAPT(DFT) procedure.^23,24

Fig. 1 Optimised geometries of C₉₆H₂₄ complexes with (a) thiophene, (b) furan, (c) pyrrole and (d) pyridine, viewed from different directions. The optimisations were performed at the B3LYP-D3/6-31G* level of theory. The carbon, hydrogen, sulphur, oxygen and nitrogen atoms of adsorbed molecules are coloured in green, pink, yellow, red and blue, respectively.

Results and discussion

E_int for C₉₆H₂₄ complexes with thiophene, furan, pyrrole and pyridine

The optimised geometries of flat-on thiophene, furan, pyrrole, and pyridine adsorbed on C₉₆H₂₄ are shown in Fig. 1(a)–(d). The centres of the five- or six-membered rings were located near the positions above the carbon atom of the underlying C₉₆H₂₄. The distances of the heteroatom from the basal plane of C₉₆H₂₄ were 3.5 Å (thiophene), 3.3 Å (furan), 3.2 Å (pyrrole), 3.4 Å (pyridine), respectively. Since the average distances of carbon atoms of the adsorbed benzene and n-pentane from the basal plane were 3.4 Å and 3.5 Å, respectively,^3d,21 the distances of HACs from the basal plane of C₉₆H₂₄ are almost the same as those of benzene and n-pentane. The changes of bond lengths, bond angles, and dihedral angles during the geometry optimisations of the C₉₆H₂₄ complexes with HACs (thiophene, furan, pyrrole, and pyridine) are small (Table S2), in contrast to the adsorption of n-perfluoroalkanes, which have conformational flexibility.^18b

The calculated E_int values for the flat-on thiophene, furan, pyrrole, and pyridine complexes are summarised in Table 1. The attractive interactions in the π-stacking of PAHs are mainly attributed to dispersion forces.² The differences between E_int (corresponding to E_B3LYP-D3) and the Hartree–Fock interaction energy (E_HF) obtained using the 6-311G** basis set (E_disp), which arise mainly from the contributions of dispersion forces, are summarised in Table 1. The contributions of the dispersion correction (E_disp-corr) are also shown in Table 1.

Table 1 Interaction energies calculated for C₉₆H₂₄ complexes with HACs, benzene and n-pentane, and contributions of dispersion forces^a

Adsorbate	Formula	Eint^b	EHF^c	Ees^d	Eind^e	Edisp^f	Edisp-corr^g
a Energy in kcal mol^−1.b Interaction energy at the B3LYP-D3/6-311G** level of theory.c Interaction energy at the HF/6-311G** level of theory.d Contribution of electrostatic interactions obtained by SAPT calculations using 6-31G* basis set.e Contribution of induction interactions obtained by SAPT calculations using 6-31G* basis set.f Dispersion energy estimated according to Edisp = Eint − EHF.g Dispersion-correction term in DFT calculations.h Not determined.
Thiophene	C4H4S	−10.09	5.39	−4.98	−0.84	−15.48	−14.28
Furan	C4H4O	−8.13	4.27	−4.00	−0.71	−12.40	−11.53
Pyrrole	C4H5N	−9.18	4.47	−4.82	−1.15	−13.65	−12.59
Pyridine	C5H5N	−10.31	5.55	−4.63	−0.92	−15.86	−14.66
Benzene^3d	C6H6	−10.55	5.77	−4.86	−0.83	−16.32	−15.11
n-Pentane^3d	C5H12	−11.72	6.80	n.d.^h	n.d.^h	−18.52	−16.90
n-Hexane^3d	C6H14	−13.81	8.05	n.d.^h	n.d^h	−21.86	−19.92

---

The SAPT energy decomposition analysis was performed for the C₉₆H₂₄ complexes with HACs and benzene using the 6-31G* basis set to evaluate the contributions of the electrostatic and induction interactions as summarized in Table S3. The electrostatic (E_es) and induction (E_ind) energies obtained by the SAPT calculations are also shown in Table 1. The contributions of the dispersion interactions were discussed based on the E_disp in Table 1, which is the difference between the interaction energy obtained by the dispersion corrected DFT calculation and that obtained by HF level calculation, since the SAPT calculations using the 6-31G* basis set underestimate the dispersion energy.²⁵

These results indicate that the nature of the interaction between HACs and C₉₆H₂₄ is very similar to that between benzene and C₉₆H₂₄. The electrostatic energies between C₉₆H₂₄ and HACs (−4.0 to −5.0 kcal mol⁻¹) are comparable to that between C₉₆H₂₄ and benzene (−4.9 kcal mol⁻¹). The induction energy between C₉₆H₂₄ and HACs was small (−0.7 to −1.2 kcal mol⁻¹). They are nearly identical to that between C₉₆H₂₄ and benzene (−0.8 kcal mol⁻¹). The dispersion interactions between C₉₆H₂₄ and HACs (−12.4 to −15.9 kcal mol⁻¹) are significantly larger (more negative) than the electrostatic and induction interactions, indicating that the attractive interaction is dominated by dispersion, as in the case of the interaction between C₉₆H₂₄ and benzene. The negligible effects of heteroatom substitution on the electrostatic and induction interactions suggest that differences in electronegativity and polarisability of the heteroatoms, and multipole moments of molecules have only a limited influence on the interaction between HACs and C₉₆H₂₄.

The E_int values for the thiophene and pyridine complexes are nearly the same as that of the benzene complex with C₉₆H₂₄, whereas those of the furan and pyrrole complexes were approximately 1 to 2 kcal mol⁻¹ smaller (less negative) than those of the benzene, pyridine and thiophene complexes. Table 1 shows that the stronger attraction in the benzene, pyridine, and thiophene complexes compared to that in the furan and pyrrole complexes arises from stronger dispersion forces. The polarisabilities calculated for benzene, pyridine, thiophene, furan, and pyrrole were 69.8, 64.3, 64.6, 49.1 and 55.1 au, respectively. Thus, benzene, pyridine, and thiophene have larger polarisabilities and, therefore, stronger dispersion forces. The differences between the molecular polarisabilities can be understood based on the differences in their composition. The polarisability of a molecule can be approximated as the sum of the polarisabilities of its atoms. Owing to the small polarizability of hydrogen atoms, their contribution to the molecular polarisability is generally small, and molecular polarizability is controlled by the contributions from the atomic polarisabilities of heavy (non-hydrogen) atoms. Therefore, the polarisabilities of benzene and pyridine, which have six heavy atoms, are larger than those of furan and pyrrole, which have five heavy atoms, and the polarisability of thiophene is greater than that of furan because thiophene contains a sulphur atom, which is more polarisable than the oxygen atom of furan.

For a given total number of heavy atoms, the E_int values of all HAC complexes were smaller than those of n-alkanes (Table 1). This suggests that the intermolecular interactions in HACs are weaker than those in n-alkanes with the same number of heavy atoms.

Edge-on adsorption structures, in which HACs were adsorbed in different directions, were also optimised and the optimised structures and their E_int values are shown in Fig. S2 and Table S4. The E_int values calculated for edge-on HACs are smaller than those for the corresponding edge-on n-alkanes with the same number of heavy atoms. In addition, the E_int values calculated for the edge-on HACs were approximately 3 to 5 kcal mol⁻¹ smaller than those for the corresponding flat-on HACs, indicating that the flat-on orientation was more favourable than the edge-on orientation, as in the cases of benzene and n-alkanes.^3d,18b Therefore, in subsequent investigations we focused only on adsorbed molecules with flat-on orientation.

E_int of for C₉₆H₂₄ complexes with benzene-fused HACs

The geometries of the C₉₆H₂₄ complexes with benzene-fused HACs were optimised as shown in Fig. 2. The distances of the heteroatoms from the basal plane of C₉₆H₂₄ in the optimised geometries of the benzothiophene, dibenzothiophene, dibenzofuran, quinoline and acridine complexes were 3.4 Å, while in those of the benzofuran, indole and carbazole complexes were 3.3 Å, respectively. The E_int, E_disp and E_disp-corr values calculated for these benzene-fused HACs complexes are summarised in Table 2. The large E_disp and E_disp-corr values indicate that the dispersion forces are the main origin of the attractive interactions in these complexes. The E_int values for complexes of the two-ring aromatic compounds are greater than those of the single-ring aromatic compounds, and are even greater for the three-ring aromatic compounds. The E_int values for the benzothiophene and quinoline complexes are close to that for naphthalene complex, while those for benzofuran and indole complexes are 1 to 2 kcal mol⁻¹ smaller. The E_int values of the dibenzothiophene and acridine complexes were close to those of the anthracene complex, whereas those of the dibenzofuran and carbazole complexes were smaller. An increase in E_int is positively correlated with an increase in E_disp (Fig. S3), indicating that the strength of the dispersion forces determines the magnitude of E_int.

Fig. 2 Optimised geometries of C₉₆H₂₄ complexes with (a) benzothiophene, (b) dibenzothiophene, (c) benzofuran, (d) dibenzofuran, (e) indole, (f) carbazole, (g) quinoline and (h) acridine, viewed from different directions. The optimisations were performed at the B3LYP-D3/6-31G* level of theory. The carbon, hydrogen, sulphur, oxygen and nitrogen atoms of adsorbed molecules are coloured in green, pink, yellow, red and blue, respectively.

Table 2 Interaction energies calculated for C₉₆H₂₄ complexes with benzene-fused HACs and PAHs, and contributions of dispersion forces^a

Adsorbate	Formula	Eint^b	EHF^c	Edisp^d	Edisp-corr^e
a Energy in kcal mol^−1.b Interaction energy at the B3LYP-D3/6-311G** level of theory.c Interaction energy at the HF/6-311G** level of theory.d Dispersion energy estimated by Edisp = Eint − EHF.e Dispersion-correction term in DFT calculations.
Benzothiophene	C8H6S	–16.00	8.42	–24.42	–22.51
Benzofuran	C8H6O	–13.86	7.27	–21.13	–19.60
Indole	C8H7N	–15.14	7.47	–22.61	–20.88
Quinoline	C9H7N	–16.27	8.66	–24.93	–23.10
Dibenzothiophene	C12H8S	–21.71	11.46	–33.17	–30.59
Dibenzofuran	C12H8O	–19.43	10.35	–29.78	–27.60
Carbazole	C12H9N	–20.70	10.51	–31.21	–28.84
Acridine	C13H9N	–22.27	11.78	–34.05	–31.55
Naphthalene^3d	C10H8	–16.57	8.86	–25.43	–23.56
Anthracene^3d	C14H10	–22.62	12.11	–34.73	–32.10

---

Geometry optimisation of the C₉₆H₂₄ complexes with HACs that have a higher number of fused benzenes was performed (Fig. S4) and the E_int values were calculated (Table S5). The E_int values presented in Tables 1, 2 and Table S4 are plotted as functions of the total number of heavy atoms in Fig. 3. The data points for HACs and benzene-fused HACs fall on the same line as those of the PAHs. Based on the linear regression analysis of the data, E_int per heavy atoms was calculated as −1.43 kcal mol⁻¹, which is close to the values previously reported for PAHs (between −1.2 and −1.4 kcal mol⁻¹).³ The E_int values for HACs, benzene-fused HACs, and PAHs are smaller (less negative) than those for n-alkanes with the same total number of heavy atoms. In fact, the E_int per the heavy atom for n-alkanes was found to be −1.96 kcal mol⁻¹,^3d which is approximately 1.4 times larger than those for HACs, benzene-fused HACs and PAHs. These results further confirm that the interactions of HACs, benzene-fused HACs, and PAHs with C₉₆H₂₄ are substantially weaker than those of the corresponding n-alkanes.

Fig. 3 Plots of interaction energy (E_int): orange triangles, blue circles, red crosses, pink asterisks, and green squares indicate E_int calculated for thiophene-, furan-, pyrrole-, pyridine-based HACs and PAHs,^3d respectively. The linear regression equation showing red dotted line is E_int = −1.43N − 1.99 (R² = 0.992), whereas that of blue dotted line derived from the n-alkanes plots is E_int = −1.96N − 1.76 (R² = 0.997): kcal mol⁻¹.^3d

ΔE_int associated with horizontal displacement of thiophene, furan, pyrrole and pyridine in C₉₆H₂₄ complexes

The ΔE_int values associated with the horizontal displacements of thiophene, furan, pyrrole, and pyridine in the C₉₆H₂₄ complexes were calculated by moving the adsorbed molecules along x and y-axes parallel to the C₉₆H₂₄ plane in increments of 0.2 Å. The distances between the heteroatoms and the C₉₆H₂₄ plane were maintained at their values in the optimised structures (thiophene: 3.5 Å, furan: 3.3 Å, pyrrole: 3.2 Å, pyridine: 3.4 Å). The calculated E_int and ΔE_int values are listed in Tables S6–S11, and ΔE_int is plotted as a function of the horizontal displacement in Fig. 4(a) and (b). Moving along the x-axis, ΔE_int reaches its maximum values at 1.2 Å and then decrease, whereas moving along the y-axis, ΔE_int reaches the maximum value at 1.4 to 1.8 Å. For the geometries with maximum ΔE_int values, most of the projections of the five- or six-membered rings containing carbon atoms and heteroatoms overlapped with those of C₉₆H₂₄, as shown in Fig. S5. Note that the maximum ΔE_int values have been suggested to be correlate with the barrier heights for horizontal displacement of adsorbed molecules.^3d,21

Fig. 4 Plots of changes of interaction energy (ΔE_int) associated with horizontal displacement along (a) x-, (b) y-axes and with (c) rotation for the C₉₆H₂₄ complexes with thiophene, furan, pyrrole, pyridine, benzene and n-pentane: Orange triangles, blue circles, red crosses, pink asterisks, green squares and purple rhombi indicate ΔE_int calculated for thiophene, furan, pyrrole, pyridine, benzene and n-pentane complexes, respectively.

As shown in Fig. 4(a) and (b), the trend of the change in ΔE_int and the maximum values of ΔE_int associated with the horizontal displacement in the HAC complexes were almost the same as those in the benzene complex, regardless of the direction of the horizontal displacement. Furthermore, the maximum values of ΔE_int for the complexes of HACs were smaller than those for the complexes of n-alkanes.^3d,21 This result indicates that HACs exhibit a resistance to horizontal displacement that is comparable to that of benzene, but lower than that of n-pentane.

The ΔE_int values calculated for the horizontal displacement of the HACs along the x-axis in the C₉₆H₂₄ complexes are summarised in Tables S6–S11. The plots of ΔE_int obtained by the HF calculations showed similar trends to those obtained by the B3LYP-D3 calculations, as shown in Fig. S6. The observed similarity implies that the dispersion interactions contribute little to ΔE_int. Rather, the change in ΔE_int associated with the horizontal displacement is attributed to the interactions included in the HF calculations, particularly the exchange-repulsion interactions in the configurations where the carbon atoms of the adsorbate and those of C₉₆H₂₄ nearly overlap.

In addition to horizontal displacement, variations in ΔE_int caused by rotation were investigated. Pyridine was rotated around the axis perpendicular to the C₉₆H₂₄ plane, passing through the intersection of the lines connecting the C2 and C5 positions and the C3 and C6 positions of the pyridine ring. The other five-membered molecules were rotated around the axis perpendicular to the C₉₆H₂₄ plane, passing through the midpoint between the two carbon atoms connected to the heteroatom (C2 and C5). The magnitude of ΔE_int variations with rotation was quite small (less than 0.5 kcal mol⁻¹), as shown in Table S12 and Fig. 4(c). The maximum value of ΔE_int due to the rotational motion was much smaller than that observed for the horizontal displacement. These results indicate that exchanging a single carbon atom of PAH with a heteroatom has little effect on the adsorption energy or resistance to rotational or lateral displacement, i.e., the substitution of a single carbon of PAH with a heteroatom has minimal effect on the stability of the π-stacked structure of PAC and the interaction energy.

Conclusions

Dispersion-corrected DFT calculations were performed to evaluate the E_int of HACs (thiophene, furan, pyrrole, and pyridine) and benzene-fused HACs on a graphite model surface (C₉₆H₂₄). Variations in E_int associated with horizontal displacement (ΔE_int) were also examined for the HACs. For a given total number of heavy atoms, the E_int values for HACs and benzene-fused HACs were comparable to those of PAHs and smaller than those of n-alkanes. The dispersion forces are primarily responsible for the attractive interactions between HACs and C₉₆H₂₄ as well as benzene-fused HACs and C₉₆H₂₄. The variation in ΔE_int associated with the horizontal displacement and rotation of the HACs showed trends similar to those of benzene. However, the maximum ΔE_int values of the HACs were lower than those of n-pentane. These results indicate that the introduction of a heteroatom into PAHs has minimal effects on the adsorption energy, rotational resistance, lateral displacement resistance, and ultimately the overall stability of the π-stacked structures composed of PAHs and HACs.

SAPT energy decomposition analysis was performed to reveal the contributions of each energy term. These calculations indicate that the nature of the interaction between HACs and C₉₆H₂₄ are very similar to that of the interaction between benzene and C₉₆H₂₄. The magnitude of electrostatic and induction interactions between C₉₆H₂₄ and HACs is nearly identical to that between C₉₆H₂₄ and benzene, suggesting that our results can also be applied to adsorption of PAHs and HACs on other polycyclic aromatic hydrocarbons such as carbon nanotubes, fullerene and larger PAHs including graphene. The interactions between C₉₆H₂₄ and HACs are stronger than the heteroatom-mediated interactions (hydrogen bonds and halogen bonds),²⁶ suggesting that the influence of these heteroatom-mediated interactions are limited when they compete with the interactions of HACs with large PAHs. The interactions between C₉₆H₂₄ and benzene-fused HACs are much stronger. Therefore, heteroatom-mediated interactions will be even more limited in the interactions of polycyclic HACs with graphene. Our results are expected to provide guidance for the fabrication, decomposition and structural modification of functional materials based on the π–π interactions of PAHs and hetero-PACs.

Author contributions

Yoshihiro Kikkawa: conceptualisation, methodology, data curation, formal analysis, investigation, visualisation, writing – original draft, writing – review, and editing. Seiji Tsuzuki: conceptualisation, methodology, data curation, formal analysis, investigation, visualisation, writing the original draft, writing – review, and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: estimated CCSD(T) level intermolecular interaction energies at the basis set limit; small molecular geometry changes upon adsorption onto C₉₆H₂₄; intermolecular interaction energies of C₉₆H₂₄-HACs and C₉₆H₂₄-benzene complexes based on SAPT calculations; optimised geometries and intermolecular interaction energies of edge-on oriented thiophene, furan, pyrrole, and pyridine on C₉₆H₂₄; optimised geometries and intermolecular interaction energies of C₉₆H₂₄ complexes with HACs with a higher number of fused benzenes; displaced geometries of C₉₆H₂₄ complexes with thiophene, furan, pyrrole and pyridine and change in interaction energies associated with horizontal displacement. See DOI: https://doi.org/10.1039/d5cp04566e.

Acknowledgements

This work was partially supported by a JSPS KAKENHI grant (JP 23K26395) to Y. K. and a JST CREST grant (JPMJCR18J2) to S. T. Part of the computation was performed using the Research Center for Computational Science, Okazaki, Japan (Project: 25-IMS-C094).

References

1. For examples: (a) R. Thakuria, N. K. Nath and B. K. Saha, Cryst. Growth Des., 2019, 19, 523–528; (b) S. Grimme, Angew. Chem., Int. Ed., 2008, 47, 3430; (c) T. Chen, M. Li and J. Liu, Cryst. Growth Des., 2018, 18, 2765.
2. (a) S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami and K. Tanabe, J. Am. Chem. Soc., 2002, 124, 104–112; (b) S. Tsuzuki, K. Honda, T. Uchimaru and M. Mikami, J. Chem. Phys., 2004, 120, 647–659; (c) C. D. Sherrill, Acc. Chem. Res., 2013, 46, 1020–1028.
3. (a) R. Zacharia, H. Ulbricht and T. Hertel, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 155406; (b) J. Björk, F. Hanke, C.-A. Palma, P. Samori, M. Cecchini and M. Persson, J. Phys. Chem. Lett., 2010, 1, 3407–3412; (c) J. Weippert, J. Hauns, J. Bachmann, A. Böttcher, X. Yao, B. Yang, A. Narita, K. Müllen and M. M. Kappes, J. Chem. Phys., 2018, 149, 194701; (d) Y. Kikkawa and S. Tsuzuki, Phys. Chem. Chem. Phys., 2025, 27, 7421–7428.
4. E. A. Meyer, R. K. Castellano and F. Diederich, Angew. Chem., Int. Ed., 2003, 42, 1210–1250.
5. (a) B. Gómez and J. M. Martínez-Magadán, J. Phys. Chem. B, 2005, 109, 14868–14875; (b) J. Goering and U. Burghaus, Chem. Phys. Lett., 2007, 447, 121–126; (c) M. Komarneni, A. Sand, J. Goering and U. Burghaus, Chem. Phys. Lett., 2009, 473, 131–134; (d) P. A. Denis and F. Iribarne, THEOCHEM, 2010, 957, 114–119; (e) R. G. Huber, M. A. Margreiter, J. E. Fuchs, S. von Grafenstein, C. S. Tautermann, K. R. Liedl and T. Fox, J. Chem. Inf. Model., 2014, 54, 1371–1379; (f) E. N. Voloshina, D. Mollenhauer, L. Chiappisi and B. Paulus, Chem. Phys. Lett., 2011, 510, 220–223.
6. (a) S. Dos Santos, M. N. Partl and L. D. Poulikakos, Constr. Build. Mater., 2014, 71, 618–627; (b) X. Ji, Y. Hou, H. Zou, B. Chen and Y. Jiang, Constr. Build. Mater., 2020, 242, 118025; (c) Y. Gong, J. Xu, R. Chang and E. Yan, Constr. Build. Mater., 2021, 273, 121758; (d) B. Schuler, G. Meyer, D. Peña, O. C. Mullins and L. Gross, J. Am. Chem. Soc., 2015, 137, 9870–9876; (e) S. Ok, J. Samuel, D. Bahzad, M. A. Safa, M.-A. Hejazi and L. Trabzon, Energy Fuels, 2024, 38, 10421–10444.
7. (a) E. Mena-Osteritz, Adv. Mater., 2002, 14, 609; (b) M. Linares, L. Scifo, R. Demadrille, P. Brocorens, D. Beljonne, R. Lazzaroni and B. Grevin, J. Phys. Chem. C, 2008, 112, 6850–6859; (c) Y. Wu, J. Li, B. Zha, X. Miao, L. Ying and W. Deng, Surf. Interface Anal., 2017, 49, 735–739; (d) C. Chen, S. Zhang, B. Tu, T. Meng, J. Li, Y. Qian, P. Li, B. Liu, W. Duan, H. Xu, F. Zhao, Y. Peng, J. Li and Q. Zeng, Langmuir, 2020, 36, 3879–3886.
8. (a) R. Iwaura and S. Komba, ACS Sustainable Chem. Eng., 2022, 10, 7447–7452; (b) R. Iwaura, Y. Kikkawa, Y. Kawashima, S. Komba, M. I. Kumano-Kuramochi, M. Ohnuma and I. Sasaki, Sustainable Mater. Technol., 2024, 41, e01025.
9. (a) J. A. A. W. Elemans, Adv. Funct. Mater., 2016, 26, 8932–8951; (b) J. Teyssandier, K. S. Mali and S. De Feyter, ChemistryOpen, 2020, 9, 225–241; (c) L. Verstraete and S. De Feyter, Chem. Soc. Rev., 2021, 50, 5884–5897; (d) S. Liu, Y. Norikane and Y. Kikkawa, Beilstein J. Nanotechnol., 2023, 14, 872–892.
10. (a) L. Kong, A. Enders, T. S. Rahman and P. A. Dowben, J. Phys.: Condens. Matter, 2014, 26, 2443001; (b) V. Skrypnychuk, N. Boulanger, V. Yu, M. Hilke, S. C. B. Mannsfeld, M. F. Toney and D. R. Barbero, Adv. Funct. Mater., 2015, 25, 664–670; (c) J. Cervenka, A. Budi, N. Dontschuk, A. Stacey, A. Tadich, K. J. Rietwyk, A. Schenk, M. T. Edmonds, Y. Yin, N. Medhekar, M. Kalba and C. I. Pakes, Nanoscale, 2015, 7, 1471–1478.
11. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian16, Revision C.01, Gaussian, Inc., Wallingford, CT, 2016.
12. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
13. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
14. B. J. Ransil, J. Chem. Phys., 1961, 34, 2109–2118.
15. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
16. J. A. Pople, M. Head-Gordon and K. Raghavachari, J. Chem. Phys., 1987, 87, 5968–5975.
17. T. H. Dunning Jr, J. Phys. Chem. A, 2000, 104, 9062–9080.
18. (a) S. Tsuzuki and T. Uchimaru, Phys. Chem. Chem. Phys., 2020, 22, 22508–22519; (b) Y. Kikkawa and S. Tsuzuki, Phys. Chem. Chem. Phys., 2023, 25, 11331–11337.
19. (a) C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622; (b) M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
20. T. Helgaker, W. Klopper, H. Koch and J. Noga, J. Chem. Phys., 1997, 106, 9639–9646.
21. Y. Kikkawa and S. Tsuzuki, Phys. Chem. Chem. Phys., 2024, 26, 24314–24321.
22. J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill and T. D. Crawford, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 556–565.
23. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887–1930.
24. K. Szalewicz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 254–272.
25. S. Tsuzuki, R. Ono, S. Inoue, S. Matsuoka and T. Hasegawa, Commun. Chem., 2024, 7, 253.
26. S. Tsuzuki, A. Wakisaka, T. Ono and T. Sonoda, Chem. – Eur. J., 2012, 18, 951–960.

---