Host–guest binding selectivity of ethylated pillar[5]arene (EtP5A) towards octane, 1,7-octadiene, and 1,7-octadiyne: a computational investigation

Host–guest binding selectivity of the perethylated pillar[5]arene (EtP5A) macrocycles with aliphatic modified hydrocarbons, i.e., octane, 1,7-octadiene, and 1,7-octadiyne as guests, has been investigated computationally employing molecular docking simulations. Density functional theory (DFT) investigations were also performed on these host–guest complexes using the dispersion-corrected approach BLYP-D3(BJ)/TZP/COSMO calculations as implemented in the ADF program and two dispersion-corrected density functionals, ωB97XD and B97D, along with the 6-311G* basis set, coupled with the PCM solvation model as implemented in the Gaussian software. We performed analysis of the frontier molecular orbitals (FMO) and natural bond orbitals (NBO), energy decomposition analysis (EDA), and noncovalent interaction (NCI-RDG) analysis. The study sheds light on the structures and binding energetics of EtP5A with the above-mentioned guests as well as on the physicochemical nature of the noncovalent interactions involved in these host–guest inclusion complexes. Based on the docking simulations, the EtP5A host revealed slightly better binding ability in the complex with the alkyne guest than with the octane and alkene, as corroborated by the EDA analysis. The results showed that the complexation of EtP5A with the hydrocarbons is mainly governed by the interplay of electrostatic interactions and dispersive noncovalent interactions. These results agree well with NCI-RDG and NBO analysis showing that host–guest binding interactions result predominantly from electrostatic C-H···π and van der Waals interactions, the H-bonding being weak or not observed. The results obtained using different computational methods were found to be in good agreement and complementary.

In 2010, Li et al. [40] reported the change of the complexation selectivity of alkylated pillar [5]arenes (P5As) towards the neutral 1,4-bis(imidazol-1-yl)butane guest with the increase of alkyl chain length and showed this inclusion complex to exhibit potent host-guest interactions with the largest association constant, ca. 10 −4 M −1 . Furthermore, EtP5A host-guest complexes with neutral guests including n-alkanes [23][24][25][26], haloalkanes [41][42][43], or unsaturated aliphatic hydrocarbons [44] have been investigated using adsorptive separation experiments. It has been shown that selective host-guest binding is profoundly influenced by noncovalent interactions [23,45]. Interestingly, Venkataramanan et al. [26]. reported a computational study on the nature of host-guest interactions between the hexane and pillar [5]arene and its quinone modified pillararenes using the dispersion-corrected DFT and wave function methods. The authors showed that the introduction of quinone in pillararene prompted flexibility in structure and the increase of the electrophilicity and the binding ability of pillararenes towards linear guest molecules.
In particular, the selectivity of EtP5A towards neutral aliphatic and linear hydrocarbon guests such as octane, 1,7-octadiene, and 1,7-octadiyne has been investigated experimentally by Hu et al. [44]. The authors concluded that the host-guest binding strength increased in accordance with the electronegativity of the guest terminal carbon atom: alkyne > alkene > alkane. At the theoretical level, much effort has been made to achieve deeper understanding of the selective binding of hydrocarbons by the pillar [5]arene host and rationalizing the underlying interactions within its inclusion complexes [10,15,22,27,28,30,33,34].
Motivated by these studies, we report for the first time a detailed theoretical investigation of the selectivity of perethylated pillar [5]arenes towards octane, 1,7-octadiene, and 1,7-octadiyne using the molecular docking approach combined with DFT calculations. We aim to explain the host-guest binding properties and the driving forces responsible for the EtP5A selective complexation towards such aliphatic guests. The study is focused on the nature of interactions within the inclusion complexes formed by the EtP5A and the three neutral hydrocarbon molecules. To thoroughly address this question, we have performed the energy decomposition analysis (EDA) study, and molecular electrostatic potential (MEP) and atom-in-molecule (AIM) studies were carried out. The results are further ascertained with the noncovalent interaction index-reduced density gradient (NCI-RDG) study. Frontier molecular orbitals (FMOs) and natural bonding orbital (NBO) analyses were performed as well providing other relevant descriptors. The conclusions drawn from the study might shed light on improving the host-guest binding and delivery properties of pillar[n]arenes.

Computational details
The three-dimensional EtP5A host-guest complexes including the hydrocarbon molecules, octane, 1,7-octadiene, and 1,7-octadiyne were generated from the pdb format using Auto-Dock Vina software [46], which allows to predict the binding free energy of complexes. The EtP5A host was used as a rigid receptor, and AutoDock tools [47] were used to determine the possible rotations of bonds by molecular docking. During the simulation, the ligands (guest molecules) were docked to the EtP5A with box dimensions set to the grid (22 × 24 × 24) Å 3 in volume, and the Cartesian coordinate center was set to 0.017, − 0.021, and − 0.097 Å. The docking configuration with the highest free energy score (ΔG) was used as the starting geometry conformations for the DFT calculations. The structures of all species, including the host, the guests, and their three complexes, generated by AutoDock tool docking simulations, were reoptimized in both gas and implicit solvent phases using the DFT approaches as implemented in two quantum chemistry programs: Amsterdam Density Functional ADF2021.107 program release [48,49] and the Gaussian 09 program package [50].
As implemented in the ADF program, the generalized gradient approximation (GGA) BLYP exchange-correlation density functional [51,52] was used along with the Grimme's D3-BJ empirical dispersion correction [53,54]. The D3-BJ correction was used to improve the description of noncovalent interactions, i.e., hydrogen bonding and van der Waals forces. Indeed, several previous post-Hartree-Fock (CCSD(T) and MP2) studies showed that dispersion corrections allow a better assessment of energies for noncovalent interactions [41,51]. For the ADF computations, the all-electron TZP basis set composed of triple-ζ Slater-type orbitals (STO) and augmented by polarization functions was used for all atoms. The solvent (chloroform) effects were introduced via the COSMO (continuum screening model) method [55]. This solvent was chosen due to the fact that it is often used for synthesis and characterization of such EtP5A complexes and various organic compounds [14]. To ensure that the optimized structures are the real minima, vibration analysis based on the harmonic frequency approximation of the optimized structures was performed using ADF single-point calculations at the same BLYP-D3(BJ)/TZP level of theory. To elucidate the nature of the host-guest binding properties in such complexes, the analysis of the noncovalent interactions (NCI) combined with the index-reduced density gradient (RDG), commonly named the NCI-RDG technique, was carried using the Mul-tiWfn program with larger size grids [56], and noncovalent interactions were depicted with the visual molecular dynamics (VMD) software [57].
The energy decomposition analysis (EDA) method developed by Ziegler et al. [58], based on the transition-state method developed by Morokuma [59], providing insights into the balance of the different bonding electronic or electrostatic factors at work between the isolated fragments, was used to analyze in a quantitative manner the nature of host-guest noncovalent interactions. The complexation or binding energy (∆E) between the guest and the host was calculated as follows: where E complex , E Host , and E guest are energies of the optimized complex, the host (EtP5A), and the free guest, respectively.
The inclusion phenomenon is accompanied by a geometry distortion which can be described by the strain (deformation) energy ∆E strain . The latter is the energy needed to deform the EtP5A or guest molecules during the inclusion process.
The strain energy was calculated as the energy difference between the molecule (host or the guest) within the inclusion complex (E complex ) and its free optimized components (E free ) given by Eq. (2): The energy decomposition analysis (EDA) was performed considering the binding energy of all host-guest complexes as a combination of several meaningful terms: the Pauli electrostatic or repulsion energy (ΔE Pauli ), the electrostatic energy (ΔE Elstat ), the orbital interaction energy (ΔE Orbt ), and the dispersion energy (ΔE Disp ), as given by Eq. (3): The method of noncovalent interaction-reduced density gradient (NCI-RDG) analysis provides the graphical visualization of the regions where noncovalent interactions occur in real space and was demonstrated to be capable of distinguishing hydrogen bonds, van der Waals interactions, and repulsive steric interactions through simple color patches [22,60,61].
Using the Gaussian 09 program, two DFT functionals including dispersion corrections, ωB97XD [62] and B97D [63], were used in combination with the full-electron 6-311G* split-valence basis set with one set of polarization functions on heavier atoms [64,65]. These computational approaches are further referred to as ωB97XD/6-311G* and B97D/6-311G*. Geometry optimizations and frequency calculations were performed with the implicit effects of chloroform taken into account (dielectric constants ε(CHCl 3 ) = 4.7113) using the self-consistent reaction field IEF-PCM method [66], with the UFF default model used in the Gaussian 09 package, with the electrostatic scaling factor α set to 1.0.
To ensure that the basis set superposition error (BSSE) does not change significantly these energies, the binding energies including BSSE at the ωB97XD/6-311G* level were calculated for the complexes 1-3 in the gas phase (vide infra).
In order to gain further insights into the nature of bonding within the complexes and to examine the electronic features, the natural bond orbital (NBO) [67] analysis as implemented in Gaussian 09 program was carried out at the same level of theory for all optimized structures. Molecular orbitals (MOs) for the ground state structures were calculated with the implicit solvent effects included as well using the ωB97XD/6-311G* and B97D/6-311G* optimized geometries. Also, molecular electrostatic potential (MEP) was calculated with the implicit solvent effects included using the ωB97XD/6-311G* and B97D/6-311G* optimized geometries. Complex structures and MOs were visualized using OpenGL version of Molden 5.8.2 visualization software [68]. Avogadro, version 1.1.1, was used to visualize the molecular electrostatic potential (MEP) maps [69,70].

Molecular docking simulations
The molecular geometries of the inclusion complexes 1, 2, and 3 including the EtP5A host and the hydrocarbon ligands as guest, i.e., octane, 1,7-octadiene and 1,7-octadiyne, respectively, were generated by AutoDock tools docking simulations and are shown on As shown in Fig. 2, all ligands are completely docked in the cavity of the host molecule in the same manner, forming the inclusion complexes with 1:1 stoichiometric ratio. The docking energy results show the binding energy (ΔG) values to be all negative suggesting that the formation of these inclusion complexes is thermodynamically favored, with the ligand 3 (1,7-octadiyne) having slightly more favored binding energy (-5.1 kcal/mol) relative to the octane and 1,7-octadiene (-4.8 kcal/mol). Analysis of interactions between EtP5A and guest molecules ( Fig. 2) shows no evidence for hydrogen bond formation, and only van der Waals and C-H···π interactions should mainly contribute to the major host-guest interactions in these inclusion complexes.

Host-guest binding energies
All geometry structures of host-guest species generated by AutoDock tool docking simulations were reoptimized using the BLYP-D3(BJ)/TZP approach. The top and side views of the DFT-optimized geometries for the host molecule (EtP5A) in solvated phase (COSMO) are shown in Fig. 3.  The BLYP-D3(BJ)/TZP binding energies reported in Table 1 were computed in both gas and implicit solvent phases. The BLYP-D3(BJ)/TZP results (Table 1) are compared to the ωB97XD/6-311G* and B97D/6-311G* results.
As compared to the docking simulation results (Fig. 2), Table 1 shows all DFT binding energies to be negative, corresponding to the bound host-guest complexes, thus indicating that the inclusion phenomenon is stabilized and thermodynamically favored. The computed BLYP-D3(BJ) gas-phase binding energies are higher than those obtained in the solvent phase, by 1.61-3.45 kcal/mol. Notably, the three DFT methods show that the complex 3, associating the EtP5A and 1,7-octadiyne molecule, is found to be the most stable one in comparison to the two others.
Furthermore, our DFT computed binding energies agree well with the experimental trend observed in previous works suggesting that host-guest binding strength increased as follows: alkyne > alkene > alkane. Indeed, as reported by Hu et al. [44], the EtP5A/1,7-octadiyne complex gives the largest association constant (K a = 82 ± 5 mol −1 L), compared to the congener EtP5A/1,7-octadiene with the constant 6.8 times lower. In the cases of the unsaturated hydrocarbon (octane), the constant was not accurately measured and was found to be too low (K a < 2 mol −1 L) [44].
Indeed, the computed binding energies of the complexes using the ωB97XD/6-311G* and B97D/6-311G* approaches steadily increase from complex 1 to complex 3, the latter being the most stable, in good agreement with the BLYP-D3(BJ)/TZP results. One can notice that the ωB97XD/6-311G* binding energies are by ca. 6.25 kcal/mol larger than the B97D/6-311G* values (Table 2). Despite these differences between the three used DFT approaches, the relative stability of the host-guest complexes is predicted to behave in the same way. Interestingly, the size of guest molecules seems to play an important role in the inclusion process. Indeed, the 1,7-octadiyne has less steric effects due to the sp hybridization of its terminal carbon, which could favor the inclusion process.
The strain energy as given by Eq.
(2) (see computational details) was computed using the BLYP-D3(BJ)/TZP approach in solvated phase and reported as relative deformation energies ∆E strain in Table 2 for both the guest and host molecules, both in free and in distorted (complexes 1-3) configurations.
From the results of Table 2, it appears that the host (EtP5A), referring to complex 3 (EtP5A/1,7-octadiyne), exhibits slightly lower strain energy ∆E strain (0.61 kcal/mol)   relative to complexes 1 and 2 associating the octane and 1,7-octadiene, respectively. Previous studies on similar stable inclusion complexes showed low strain energy for both the guest and host molecule [22,42,71]. Moreover, regarding the strain energies ∆E strain of the three guest molecules, the alkyne molecule shows the lowest deformation energy (0.70 kcal/mol) relative to the alkane and alkene, thus resulting in a better guest molecule inclusion in the cavity of the host. These results are in line with the higher association constant (K a ) between the alkyne (guest) and EtP5A (host) within complex 3, which was observed to be 6.8 times larger than the association constant measured for alkene in complex 2, while for alkane-EtP5A complex 1, the association constant is much weaker [44].

Frontier molecular orbital analysis
It was established that for a macromolecule to act as a controlled delivery system, it should reversibly adsorb and release a ligand (e.g., a drug) via a kinetically controllable process [9,22,66]. To understand the kinetic stability of these inclusion complexes, the E gap , the energy difference between the HOMO (highest occupied molecular orbital), and LUMO (lowest unoccupied molecular orbital) was computed using the ADF/BLYP-D3BJ method. The computed HOMO and LUMO energies (eV), E gap (eV), and the global dipole moment μ (Debye) for all the species in the solvent phase are provided in Table 3.
As shown in Table 3, the HOMO-LUMO gap for the inclusion complexes is significantly smaller than for the hydrocarbons (guests), which indicates that complex formation can be kinetically controllable, in good agreement with previous studies [9,22,66]. Moreover, complexes 1, 2, and 3 show similar energy gaps, ca. 3.15 eV, close to the host EtP5A (3.16 eV). The same trend is obtained using the ωB97XD and B97D functionals (see Tables S1 and S2), showing that all complexes have the HOMO-LUMO gaps noticeably smaller than the guest hydrocarbons, again in good agreement with the BLYP/D3(BJ) results. HOMO and LUMO of the complexes are mainly those of the host molecule as it can be seen in Fig. 5.
Moreover, it appears that the computed global dipole moment (μ) of the free guest decreases with the change of the terminal carbon atom hybridization, sp > sp 2 > sp 3 . Indeed, the higher electronegativity of the sp-hybridized terminal carbons on 1,7-octadiyne was reported, resulting in stronger binding between alkyne and EtP5A compared to the other guests [44]. For the saturated hydrocarbon (octane), very weak binding is observed as can be seen from the lower value of the association constant [44]. Interestingly, the dipole moment trend is inverted for the inclusion complexes leading to the highest dipole moment (3.236 D) for the octane-containing (sp 3 ) complex 1. It can be seen from Table 4 that the dipole moment of the host in its distorted geometry (within complexes 1-3) exhibits the same trend as the dipole moment of the complexes, increasing from 0.059 to 1.637 D when passing from the free to the inclusion structure (complex 3).
This increase of dipole moment from the host and guest to the complexes is also closely related to the enhancement of induced dipole-dipole interactions (van der Waals interactions) between the guest and host molecules which are, in cooperative multiple noncovalent interactions, essential for realizing such strong complexations. Therefore, it is likely that the latter interactions play a significant role in stabilizing the inclusion process. This point will be detailed later in the text.
The FMO plots for the complexes obtained at the BLYP/ D3(BJ) level are shown in Fig. 5. For all the complexes, the HOMO density is predominantly located on the host moiety, with no contributions from the guest molecules. This indicates that no electronic charge transfer occurs within such host-guest inclusion complexes, and the encapsulated ligand does not alter the electronic structure of the host significantly. This result suggests that the encapsulation occurs by a pure physical adsorption process [10].
The comparison with the ωB97XD/6-311G* and B97D/6-311G* computed plots of the frontier MOs (cf. Figs. S3 and S4) shows good agreement with the BLYP-D3(BJ)/TZP obtained plots (Fig. 5). Even with larger isosurface value, 0.01 au, it is clearly seen that both HOMO and LUMO of the complexes are dominated by the host moiety contributions, with very small, if any, guest contributions.
The binding energy values (Tables 2 and S1), along with the HOMO and LUMO plots, indicate that the inclusion  process is dominated by rather physical adsorption and suggests that the EtP5A can be used as a reversible host for the considered hydrocarbon molecules as reported for delivery systems [72].

Molecular electrostatic potential (MEP) analysis
Molecular electrostatic potential analysis is a useful tool for studying the charge distribution within the host cavity and orientation of noncovalent bonds in the inclusion complexes [73]. A rationalization of noncovalent intermolecular interactions based on MEP analysis has been carried out using the BLYP-D3(BJ)/COSMO results including the implicit chloroform solvation effects. Figure 6 displays molecular electron density isosurfaces (0.001 au) of the EtP5A host overlaid, e.g., with MEP of the inclusion complex 2. Similarly, the MEP analysis has been carried out using ωB97XD/6-311G* calculations (Fig. S5). Figure 6 shows that the electrostatic potential surface of the complexes is typically non-uniform. It reveals combination of electron-rich negative regions (depicted in red) corresponding mainly to the oxygen-containing groups and electron-deficient positive parts (blue) corresponding to ends of the alkyl groups. Such positive regions can be considered as electrophilic sites, and negative regions surrounding the guest molecules in attractive interactions are considered as nucleophilic sites. Indeed, as reported by previous studies [22,42,73,74], these electrostatic properties shown on the EtP5A macromolecule favor stronger attractive host − guest interactions arising from their binding during the inclusion process (Fig. S5). The polarization of the host upon complexation is significant as it can be seen comparing Fig. 6a and b.

Energy decomposition analysis (EDA)
To shed more light on the contributions of the different components of the interaction energy ΔE bind within the considered host-guest complexes, the energy decomposition analysis (EDA) was performed considering the implicit solvent phase (chloroform) using the BLYP-D3(BJ)/TZP COSMO approach. The total binding energies as well as the decomposition terms (kcal/mol) along with their contributions in percentage (%) are reported in Table 4. The different decomposition terms are given by Eq. (4) (see Computational details): In Table 4, the percentages (%) are the weights of the different stabilizing terms, electrostatic interactions (ΔE elst ), orbital interaction (ΔE orb ), and dispersion (the London or van der Waals interactions) (ΔE disp ), excluding the Pauli repulsion (ΔE Pauli ) destabilizing term from the components of the total binding energy (ΔE bind ).
As shown in Table 4, the destabilizing Pauli term (ΔE Pauli ) is computed to be significantly higher for the complex 1 (52.96 kcal/mol) than for complexes 2 and 3. This is likely due to the increased steric hindrance caused by the terminal carbon hybridization (sp 3 ) of the guest molecule (octane). The stabilizing electrostatic (ΔE elst ) term is computed to be lower for complexes 1 and 2 than for complex 3, in relation with the higher electronegativity of the terminal carbon of the alkyne (sp). Interestingly, one can notice that the dispersion term (ΔE disp ), the dominating effect in the inclusion process, is predicted to be slightly higher in complex 1 and 2 than in 3. It is noteworthy that regarding the contribution percentage, it appears that the dispersion (ΔE disp ) term is the most important factor in the stabilizing inclusion process, depending on electronegativity of the terminal atom in this case [41]. Therefore, it appears that the use of dispersion correction is essential in estimating the correct binding energy and to determine the long bonding interactions as stated in the literature [22,41,53,62].
The electrostatic contribution (ΔE elst ) for complexes 1-3 is also predicted to be significantly high, contributing to their stability. Furthermore, this electrostatic term seems to increase with the terminal carbon hybridization change (sp > sp 2 > sp 3 ), being slightly higher for complex 3 in comparison to complexes 1 and 2 (29.0 vs. 24.9 and 26.1, respectively). The orbital interaction (ΔE orb ) term has low contribution in comparison to the other terms and suggests weak charge transfer and polarization interactions within these host-guest complexes.
It appears that electrostatic and dispersion effects play a crucial role in stabilizing such inclusion complexes, corroborating well the experimental and theoretical results [22,41,44]. The computed binding energy shows the inclusion complexes to be stable, and thus, their formation can occur at room temperature. In order to check this point, basis set superposition error (BSSE) calculations have been carried out using the Gaussian software at the ωB97XD/6-311G* level in the gas phase. The computed complexation energies with or without the BSSE corrections are reported in Table 5. From this latter table, it follows that taking into account the BSSE somewhat lowers the binding energy as usual; however, the binding trend remains the same as without BSSE. Notably, the BSSE corrections do not change significantly the energies.

NCI-RDG analysis on the inclusion complexes
The NCI-RDG analysis was carried out to visualize the nature of interactions between the EtP5A (host) and the three hydrocarbon molecules using the program MultiWfn with larger size grids [56]. Figures 7 and 8 display the NCI-RDG mapping of the EtP5A and the corresponding inclusion complexes 1-3, respectively.
Typically, the blue and green color regions are used to represent the stabilizing H-bonding and van der Waals interactions, the while red color region is used for the destabilizing steric interactions [26,56,59,60]. The interaction types are determined by analyzing the sign of λ 2 which represents the eigenvalues of the Hessian of the electron density at the BCP (bond critical points). As reported by other authors [75,76], for the bonding interactions such as H-bonding, the λ 2 value will always be less than 0, while for the van der Waals type of interactions, it is λ 2 ≤ 0. The non-bonding interactions like steric repulsions feature λ 2 > 0. As shown by the RDG isosurfaces (Figs. 7 and 8), the strength and type of noncovalent interactions can be assessed from the product sign (λ 2 )ρ (ranging from − 0.5 to + 0.5 au). For the EtP5A macromolecule (Fig. 7), the NCI results (λ 2 < 0) show that the host is stabilized by hydrogen bonds between the oxygen atoms of alkoxy groups and the methylene hydrogens; the red disk area (λ 2 > 0) is also observed indicating a steric repulsion and at the middle of the benzene rings.
In the inclusion complexes (Fig. 8a-c), one can observe green areas between the C-H of terminal carbon atoms in the hydrocarbon guests and the EtP5A host, indicating the existence of van der Waals forces and electrostatic CH⋯π interactions as pointed out in the literature [26]. The H-bonding is weak as follows from the RDG map (RDG > 0.80). Destabilizing interactions, shown by the red disk area, are also observed indicating steric repulsions and seem to appear at the middle of the benzene rings. It is noteworthy when passing from Figs. 7 and 8 i.e., after the structural reorganization upon encapsulation of guests, the steric repulsions (red zone) within the host were slightly reduced by van der Waals stabilizing interactions (green zone). Indeed, one can visualize on NCI isosurfaces the green patches to be evenly distributed inside the EtP5A host molecule in the inclusion complex systems, indicating electrostatic interactions, which are responsible for stabilization of the hydrocarbon molecules inside the EtP5A unit. This phenomenon was also shown by previous studies [21,42,56,59,60]; however, it reveals to be quite similar for the three octane (sp 3 ), 1,7-octadiene (sp 2 ), and 1,7-octadiyne (sp) inclusion complexes, corroborating the EDA results.

NBO (natural bonding orbital) analysis
With the help of NBO analysis, we can gain significant insights into different types of orbital interactions in the compounds of interest [67]. In general, the NBO analysis is performed by considering the possible interactions between the filled (donor) Lewis-type NBOs and the empty (acceptor) non-Lewis-type NBOs. Then, their energetic contributions are computed using second-order perturbation theory. These donor-acceptor interactions lead to a decrease in the localized NBOs occupancy in the idealized Lewis structure and an increase in the occupancy of the empty non-Lewis orbitals. Therefore, such interactions are referred to as "delocalization" corrections to the zero-order natural Lewis-type structures. The stronger donor-acceptor interactions have higher stabilization energies [67]. The second-order stabilization energy E (2) is calculated using Eq. (5): In this equation, ε i and ε j represent off-diagonal elements, and F i.j is the diagonal elements of the NBO Fock matrix, q i is the donor orbital possession, and E (2) is the energy of stabilization. Table S3 lists the representative values for the host-guest stabilization energies along with energy differences between donor and acceptor i and j NBO orbitals and diagonal NBO Fock matrix elements for complexes 1-3 (with the stabilization energies threshold 0.1 kcal/mol). In Tables S4-S6, the more significant stabilization energy values are listed for complexes 1-3, respectively (with the stabilization energies threshold more than 0.2 kcal/mol), all of them computed at the ωB97XD/6-311G* level in the implicit CHCl 3 solvent.
Analysis of these stabilization energy values for complexes 1-3 shows the following. (i) For all three complexes, there exist quite numerous stabilizing intramolecular interactions host→guest and guest→host, with the highest interaction energy 0.99 kcal/mol, for the host→guest LP(O)→σ*(C-H) in complexes 2 and 3. (ii) Host→guest interactions are represented by both π(C-C)→σ*(C-H)-and LP(O)→σ*(C-H)type interactions, with the former due to the relatively highly polarizable π-electrons of the benzene rings and the latter due to relatively polarizable (compared to covalent bond electron clouds) lone-pair electrons of the alkoxy groups. In complexes 1 and 2, generally the LP(O)→σ*(C-H)-type interactions have somewhat higher energies than the π(C-C)→σ*(C-H)-type interactions (cf . Tables S4 and S5), whereas in complex 3, the LP(O)→σ*(C-H)-type interactions generally have somewhat lower energies than the π(C-C)→σ*(C-H) interactions (Table S6). (iii) Guest → host interactions are represented by the σ(C-H)→π*(C-C)-and σ(C-H)→π*(C-C)-type interactions in complex 1 (see Table S4); by the π(C-C)→σ*(C-H)-, σ(C-H)→π*(C-C)-, and σ(C-H)→LP(C)-type interactions in complex 2 containing the alkene (Table S5); and by the π(C-C)→σ*(C-H)-and σ(C-H)→π*(C-C)-type interactions in the complex 3 containing the alkyne (Table S6). The total energy of the guest→host interactions in the complex 3 is higher compared to other two complexes, which can explain higher binding energy in complex 3.

Conclusions
Computational investigation of the perethylated pillar [5] arene (host) selectivity towards the octane, 1,7-octadiene, and 1,7-octadiyne (guests) using the molecular docking simulations demonstrates the stabilizing binding of their host-guest inclusion complexes, with the ΔG scores being all negative, suggesting that their formation can occur at room temperature. Corroborating the experimental outcomes, the computed 1,7-octadiyne complex has slightly more favored binding energy relative to the octane and 1,7-octadiene congeners. Moreover, DFT calculations based on the energy decomposition analysis (EDA) reveal a key role of the electrostatic and dispersion energy in stabilizing the structures of the considered complexes. The stabilizing electrostatic contribution (ΔE elst ) is also predicted to be significantly higher for complex 3 in comparison to the two other complexes 1 and 2 and increases with the ligand terminal carbon hybridization (sp > sp 2 > sp 3 ). As expected, the increase of the dipole moment of the free guests is in line with the hybridization of the terminal carbon atom in these guests, sp > sp 2 > sp 3 . The increase of the dipole moment of the host upon complexation is related to the enhancement of induced dipole-dipole interactions (van der Waals interactions) between guest and host molecules which are, in cooperative multiple noncovalent interactions, essential for realizing such strong complexations. Therefore, it is likely that the latter interactions play a significant role in stabilizing the inclusion process.
The corroborating results obtained using different DFT functionals for the binding energies, HOMO-LUMO pictures, and molecular electrostatic potential (MEP) plots indicate that the inclusion process is governed by the physical adsorption, suggesting that the EtP5A can be used as a reversible (delivery) host for considered hydrocarbon molecules. Furthermore, the HOMO-LUMO energy gaps for the inclusion complexes are calculated to be lower for the host (EtP5A) than for the guests, which indicates that the complex formation can be kinetically controllable. Notably, the NCI-RDG analysis shows that host-guest binding properties are mainly driven by weak stabilizing noncovalent interactions, i.e., van der Waals forces (dispersion forces) and electrostatic C − H···π interactions that are responsible for the EtP5A slightly selective complexation towards the alkyne. This latter is favored by the higher electronegativity of the terminal carbon atom with the sp hybridization. The occurrence of H-bonding was shown to be unlikely in this case. The NBO intramolecular interaction analysis shows that, for all three complexes, there exist numerous stabilizing hosT→guest and guesT→host intramolecular interactions. The former is represented by both π(C-C)→σ*(C-H)and LP(O)→σ*(C-H)-type interactions, whereas the later are represented by different σ(C-H)→σ*(C-H)-, π(C-H) σ(C-C)-, π(C-C)→σ*(C-H)-, and σ(C-H)→LP(C)-type interactions in complexes 1-3. These guest→host interaction energy in complex 3 is found higher compared to the other complexes 1 and 2, further explaining its slightly higher binding energy.
Funding Aleksey Kuznetsov appreciates the financial support of USM and computational facilities of the Department of Chemistry, ITA, Brazil.
Data availability All data generated or analyzed during this study are included in this published article and its supplementary information files.

Conflict of interest
The authors declare no competing interests.