Anion Recognition by Pyrylium Cations and Thio-, Seleno- and Telluro- Analogues: A Combined Theoretical and Cambridge Structural Database Study

Pyrylium salts are a very important class of organic molecules containing a trivalent oxygen atom in a six-membered aromatic ring. In this manuscript, we report a theoretical study of pyrylium salts and their thio-, seleno- and telluro- analogues by means of DFT calculations. For this purpose, unsubstituted 2,4,6-trimethyl and 2,4,6-triphenyl cations and anions with different morphologies were chosen (Cl–, NO3– and BF4–). The complexes were characterized by means of natural bond orbital and “atoms-in-molecules” theories, and the physical nature of the interactions has been analyzed by means of symmetry-adapted perturbation theory calculations. Our results indicate the presence of anion-π interactions and chalcogen bonds based on both σ- and π-hole interactions and the existence of very favorable σ-complexes, especially for unsubstituted cations. The electrostatic component is dominant in the interactions, although the induction contributions are important, particularly for chloride complexes. The geometrical features of the complexes have been compared with experimental data retrieved from the Cambridge Structural Database.


Introduction
Noncovalent interactions have a constitutive role in the science of intermolecular relationships. In particular, those involving aromatic rings play a vital role in chemistry and biology [1], which OPEN ACCESS becomes prominent in drug-receptor interactions, crystal engineering and protein folding [2]. In recent years, the contact between an anion and the region above the plane of an electron-deficient aromatic ring has also been recognized as a noncovalent bonding interaction. The nature of this interaction, designated an "anion-π bond" [3], has been described by numerous theoretical studies, which demonstrate that it is energetically favorable [3][4][5][6][7][8], in addition to several experimental investigations [9][10][11][12].
While anion-π interactions continue to gain attention as their role in chemical processes is being increasingly recognized [13][14][15][16][17][18], studies of cationic oxygenated aromatic rings involved in anion-π bonding were absent in the literature until very recently [19], where the authors report an experimental and theoretical study for pyrylium-based tri-arylated series of compounds with a tetrafluoroborate anion. Pyrylium salts are a very important class of cationic organic molecules containing a trivalent oxygen atom in a six-membered aromatic ring [20,21]. They can act as intermediates for a variety of syntheses [22] and have been exploited to design sensors for anions [23], amines, amino acids and chameleon labels for quantifying proteins [24], as well as a kind of ionic liquid crystal material, which can be used in the display industry [25,26].
In this work, we present a systematic computational study on unsubstituted, trimethyl and triphenyl pyrylium cations (Scheme 1) interacting with Cl -, NO3and BF4 -. We chose a monoatomic, a trigonal planar and a tetrahedral anion to account for different shapes and to analyze the effects that this would have in the formation of the complexes. Moreover, we also included in our study the thio-, seleno-and telluro-analogues of the pyrylium cations previously mentioned to analyze the effect of the chalcogen atom in such complexes (Scheme 1). The variation of the chalcogen atom comes also from the idea of studying the existence of chalcogen bonds. These noncovalent contacts belong to the σ-hole interaction family, as tetrel, pnicogen, halogen and aerogen bonds, which has attracted considerable attention in recent years, and they are recognized by the scientific community as powerful tools in supramolecular chemistry, crystal engineering and biochemistry [27][28][29][30][31][32]. We have also carried out natural bond orbital (NBO), atoms in molecules (AIM) and symmetry-adapted perturbation theory (SAPT) calculations to characterize the complexes and to analyze the physical nature of the interactions. Scheme 1. Molecules included in our study and the numbering of C atoms.

Preliminary Studies
We chose different aromatic systems to interact with chloride, namely pyrylium, thiopyrylium and selenopyrylium (1, 2 and 3, respectively, in Scheme 1). We carried out Density Functional Theory (DFT) calculations at the BP86-D3/def2-TZVPD, B97-D3/def2-TZVPD levels and resolution of the identity (RI) second-order Møller Plesset (MP2) calculations at the RI-MP2/aug-cc-pVTZ (AVTZ) level. Geometry optimizations for all complexes. The geometries of the analyzed complexes are shown in Figure 1, and the geometric and energetic results of these calculations are gathered in Table 1. From the energetic point of view, regardless of the computational method used, all interaction energies are negative, indicating that the formation of the complexes between the aromatic systems and Clis favorable, as expected from the interaction of two ions of opposite sign. The largest binding energy corresponds to the complex of 3, which decreases for 2·Cland then again when the ring is 1.
Comparison of the DFT energy values yields larger binding energies for the BP86 than for the B97 results. The B97 interaction energies are in better agreement than the BP86 ones with the MP2 results, with energy differences as large as 2.7 kcal·mol −1 . Furthermore, the agreement is even better when B97 and CCSD(T)/CBS (CCSD(T): coupled-cluster singles doubles and non-iterative triples correction; CBS: complete basis set) interaction energies are compared.  From the geometrical point of view, a comparison has also been made among the different levels of theory, with two points worth being mentioned. On the one hand, the equilibrium distances obtained at the B97 level are generally slightly longer (≥0.07 Å) than the ones obtained at the BP86 level (see Table 1). For instance, for the 1· Clcomplex, a 0.07 Å difference in length is observed. On the other hand, the BP86 equilibrium distances are in good agreement (better than the B97 ones) with the MP2 ones, where the former are slightly longer than the latter.
From the results gathered in Table 1, it can be deduced that despite the small differences in equilibrium distances, the B97 method nicely reproduces all of the interactions. Therefore, the comparison of B97 and CCSD(T) results validates the use of the RI-B97-D3/TZVP level of theory to study the interaction of our systems with Cl -. From now on, this is the computational methodology that will be used throughout the manuscript.
We have also computed the molecular electrostatic potential (MEP) surfaces of the pyrylium cations and analogues (except from tellurium compounds), which are represented in Figure 2. From the analysis of the MEPs, several considerations can be made. The pyrylium Compound 1 has the highest electrostatic potential (ESP = 140 kcal·mol −1 ) located at the center of the molecule, although a little bit displaced towards the location of the oxygen atom. The ESPs associated with in-plane interactions are 5 kcal·mol −1 lower in energy. If we change the X atom, we observe how the ESP located above the molecular plane diminishes (to 125 and 119 kcal·mol −1 for 2 and 3, respectively). Interestingly, the highest ESP for 2 and 3 is located at both sides of the X atom (137 and 141 kcal·mol −1 , respectively). Therefore, from electrostatic considerations, the most favorable anion complex of 1 would be that where the anion is positioned above the ring, whereas the most favorable complexes of 2 and 3 would be those where the anion is located in the molecular plane and close to the X atom. However, for π interactions, polarization effects can be important, and the prediction based only on electrostatic effects may not be adequate. As inferred from Figure 2, the trimethyl Derivatives 5-7 have less positive values of ESP than their respective unsubstituted counterparts, and likewise, the triphenyl Analogues 9-11 have less positive ESP values than their respective trimethyl derivatives.

Complexes with Cl -
First, we are going to analyze the complexes of Compounds 1-12 that are formed with chloride. After exploring the potential energy surface of the different complexes, four main structural configurations were found. Type a and b complexes are the result of a nucleophilic attack of the anion to C1 and C3 atoms, respectively, of the heteroaromatic ring, giving rise to σ-complexes. Thus, according to our nomenclature, 1a will be the σ-complex of 1 with Cl -. In Type c complexes, the anion is located above the molecular plane, either interacting with C1 and C5 or the chalcogen atom. Lastly, Type d complexes are based on hydrogen bonding interactions. Several geometrical features of the complexes with chloride are shown in Figure 3, and their corresponding binding energies are gathered in Table 2.
From the inspection of the results, several points arise. For either a and b type complexes and within every family of compounds (unsubstituted, trimethyl and triphenyl substituted), the Cl-C distance is shortened when going down the chalcogen group. This result is consistent with the electrophilic character of the heteroaromatic ring, which is increased when going from O to Te atoms. Moreover, the Cl-C distance is increasingly lengthened when trimethyl and triphenyl substitution is considered, in agreement with the electrostatic results retrieved from MEP calculations ( Figure 2). Second, a complexes have shorter Cl-C distances than b complex. Third, 12a·Cland 12b·Cladducts could not be located. All attempts to find such complexes led to 12c·Cl -. For c complexes, we observe a slightly different placement of the anion if X = O (1c) or X = S, Se, Te (2c-4c) are considered. In Complex 1c, the anion is located above the molecule and between C1 and C5, forcing the out-of-plane displacement of the oxygen atom opposite of Cl -. In fact, this complex is not a minimum on the potential energy surface, but a TS corresponding to the chloride transfer from the σ-complex on C1 (1a·Cl -) to its degenerate σcomplex on C5 [33]. However, for 2c-4c complexes, Clis parallel displaced towards the X atom, giving rise to a ClXC3 angle of ca. 107°. Moreover, for such complexes, the X atoms are pulled from the molecular plane towards the anion yielding Cl -···X distances (2.50, 2.57 and 2.63 Å for 2c, 3c and 4c, respectively) that are remarkably much shorter than the sum of the van der Waals radii (3.91, 4.01 and 4.17 Å, respectively) [34,35]. In fact, all intermolecular contacts reported throughout the manuscript are substantially shorter than the sum of the corresponding van der Waals radii. Analogously to a and b complexes, the Cl -···X distance in c complexes is increasingly lengthened when trimethyl and triphenyl substitution is considered, which could be due, on the one hand, to the decrease in electrophilic character of the heteroaromatic ring for the trimethyl and triphenyl systems ( Figure 2) and, on the other hand, to the presence of CH···Clhydrogen bonds that slightly lower the net charge on the anion. These hydrogen bonds appear to be particularly important for the triphenyl derivatives, since their respective hydrogen bond distances are short (2.65-2.87 Å).  Every attempt to obtain exclusively hydrogen bond-based complexes was unsuccessful for 1-8 compounds, and the geometry optimizations always ended up yielding one of the geometries mentioned above. However, hydrogen bond complexes (Type d) were obtained for 9-12, and their geometries are shown in Figure 3. All of these Type d complexes have the same recognition pattern, where the anion is interacting with two hydrogen atoms in ortho to two phenyl rings, with the hydrogen atom bonded to C2. The CH···Cldistance oscillates from 2.28-2.54 Å for X = O to 2.44-2.74 Å for X = Te. Table 2. Interaction energies (∆E, in kcal·mol −1 ), electron density (ρBCP), its Laplacian (∇ 2 ρBCP) and the total electron energy density (HBCP) for the bond critical points (BCPs) of Cl -···C or Cl -···X contacts of complexes of Compounds 1-12 with Cland the number of imaginary frequencies (Nimag).
The interaction energies of all complexes are large and negative, as shown in Table 2, ranging from −116.0 kcal·mol -1 for 1a·Clto −77.5 kcal·mol -1 for 10d·Cl -, values than can be rationalized from the electrostatic point of view. The largest binding energies are generally found for a and b σ-complexes: For the unsubstituted compounds, we usually obtain the largest binding energies for a complexes. Analogously, the same trend is also observed for the trisubstituted derivatives, i.e., a complexes have larger ∆E values than b complexes, with c complexes having the smallest binding energies. However, we observe that the energy gap between a and c complexes gets smaller when going from 1 (∆∆E = 4.8 kcal·mol −1 ) to 3 (∆∆E = 0.9 kcal·mol −1 ). In fact, for Compound 4, its c complex is more favorable than its a complex by 5.2 kcal·mol −1 , leading to the conclusion that the more electrophilic the X atom, the more favorable the c complex. A similar behavior is observed for complexes of 5-8, with ∆∆E values between a and c complexes of 6.2, 0.2 and −3.7 kcal·mol −1 for 1, 3 and 4, respectively. Strikingly, for Compounds 10-12, their c complexes are the most favorable (neither 9c·Clnor 12a·Clcould be located). For instance, the energy gap between a and c chloride complexes is −5.6 and −6.8 kcal·mol −1 for 10 and 11, respectively.
To investigate the origin and nature of these interactions, AIM, NBO and SAPT studies have been carried out for selected complexes. We have examined all possible intermolecular interactions between occupied (donor) Lewis-type NBOs and vacant (acceptor) non-Lewis NBOs and estimated their energetic importance by second-order perturbation theory. According to the NBO analysis, the interaction in c complexes of X = S, Se and Te is primarily based on a charge donation from the lone pairs of Clto the vacant π* orbital of the X-C1 or X-C5 bond, as derived from the calculated second-order orbital perturbation energies listed in Table 3. As a matter of fact, these charge-transfer contributions are very large, with increasing energies as we move from the lightest (44.9 kcal·mol −1 for 2c·Cl -) to the heaviest (62.3 kcal·mol −1 for 4c·Cl -) chalcogens for the unsubstituted complexes. The same trend is observed for the trimethyl derivatives, although with smaller energy contributions than those for their respective unsubstituted cases. Moreover, there is an excellent linear correlation (r 2 = 0.969) between the binding energies (∆E in Table 2) and the second-order orbital perturbation energies (E(2) in Table 3) mentioned above for 2c-4c·Cland 6c-8c·Clcomplexes, indicating that these Lp(Cl -)→π*(X-C1/C5) donor-acceptor contributions are very relevant ( Figure 4). It is worth mentioning the existence of other charge-transfer contributions of up to 3.7 kcal·mol −1 that arise from the charge donation from Cllone pairs to empty d orbitals of S, Se and Te ( Table 3). All of these results suggest the existence of a chalcogen bond through π-hole interactions. Differently from 2c-4c·Cland 6c-8c·Clcomplexes, there is no charge transfer from Clto either O-C1 or O-C5 empty π* orbitals in 1c·Cland 5c·Cl -. Instead, modest energetic contributions are obtained from Lp(Cl -)→π*(C1-C2) and Lp(Cl -)→π*(C4-C5) donor-acceptor interactions, leading to anion-π interactions. The absence of charge transfer processes in 1c·Cland 5c·Clindicates that induction forces are not important in their binding mechanism, which will be primarily based on electrostatics, as we will see later in the SAPT analysis. As already noted, the geometries of X = O complexes are different from the X = S, Se and Te ones (X is moving away with respect to the anion for the former and getting closer to the anion for the latter), giving rise to a different binding mechanism, where induction forces are supposed to play an important role in the latter compounds. Table 3. Natural bond orbital (NBO) second-order orbital perturbation energies (E(2), in kcal·mol −1 ) derived from charge donation from an occupied (donor) to an empty (acceptor) orbital for Complexes 1c-8c·Cl -.  In line with the NBO results, our AIM calculations for 2c-3c·Cl -, 6c-7c·Cland 10c-11c·Clcomplexes show the existence of one BCP connecting the X atom with Cl -( Figure 5). The values of the electron densities at these BCPs range between 0.061 au for 2c·Cland 0.041 au for 10c·Cl -( Table 2). At these BCPs, the values of the Laplacian are positive, and the total energy densities are negative (between −0.010 and −0.003 au), suggesting that these interactions have a certain covalent character [36]. In addition, the nature of the chemical bond in 1a-7a·Cland 1b-7b·Clσ-complexes has been analyzed by means of the AIM approach. From the inspection of the results (Table 2), the values of the electron density at the BCP originating from the Cl -···C1/C3 contact are large (between 0.14 and 0.06 au) with values of the Laplacian oscillating between negative (1a-3a, 7a, 1b and 3b) and positive (5a-6a, 2b and 5b-7b). However, the total energy density for all cases is negative (between −0.071 and −0.011 au; Table 2), indicating that there is a significant sharing of electrons for this contact.

Complex Donor Acceptor E(2)
The physical nature of the interaction between the anion and Compounds 1-3, 5-7 and 9-11 has been analyzed by means of SAPT calculations in c complexes. The energy contributions obtained from the SAPT partitioning scheme are listed in Table 4. As expected, the electrostatic component, Eel, is the largest contributor to the interaction with values ranging from −103.8-−138.1 kcal·mol −1 . The Eel values vary according to the electronegativity of the X atom, i.e., the smallest and largest contributions are found for X = O and X = Se, respectively. The positive exchange part is large and partly compensates for the electrostatic component, especially for X = S, Se complexes. The dispersion contribution is more or less constant for S and Se complexes and larger than that for O complexes. The induction-polarization contribution, Eind, however, seems to play an important role in the interaction, since large variations are detected, in line with the donor-acceptor NBO results. Thus, Eind is small for the oxygen-based 1c and 5c·Clcomplexes (ca. −12 kcal·mol −1 ), large for sulfur-based 2c, 6c and 10c·Clcomplexes (from −24 to −34 kcal·mol −1 ) and the largest for selenium-based 3c, 7c, and 11c·Clcomplexes (from −34-−48 kcal·mol −1 ). The consideration of higher-order contributions to the interaction energy, δHF (Table 4), is mandatory at least for X = S compounds to get a good agreement between SAPT interaction energies, ESAPT, and the binding energies. Table 4. DF-DFT-symmetry-adapted perturbation theory (SAPT) electrostatic, exchange, induction, dispersion and Hartree-Fock higher-order energy contributions (Eel, Eexch, Eind, Edis and δHF, respectively), SAPT interaction energy (ESAPT) and binding energies (∆E) of selected complexes. Energies in kcal·mol −1 .

Complexes with NO3 -
In this second section, we are going to analyze the complexes formed by Compounds 1-12 with nitrate as the anion. The exploration of the potential energy surface leads to the formation of five main structural configurations: Type a, b and d complexes, just as defined in the previous subsection, Type c complexes, which are very similar to the ones found for chloride, but now the anion is interacting either with the chalcogen atom or a C atom and the chalcogen atom, and the new Type e complexes, where the anion is solely interacting with C atoms above the molecular heteroaromatic ring. The geometries of the complexes can be found in Figure 6, and their respective binding energies are listed in Table 5.
Complexes a and b follow the same trend observed for chloride complexes, i.e., the O2NO-C distance shortens/lengthens with the size of the chalcogen atom/with trisubstitution, and a complexes have shorter O2NO-C distances than b complexes. However, a and b complexes are only obtained for the unsubstituted Compounds 1-4. In c·NO3complexes, different orientations of nitrate are observed. For unsubstituted compounds, the nitrate is perpendicularly placed above the heteroaromatic ring, giving rise to short O2NO -···X contacts in 2c-4c·NO3 - (Figure 6), similarly to what is observed for 2c-4c·Cl -. Type c complexes for the trimethyl derivatives fall into two subcategories, depending on whether the anion is almost perpendicularly placed above, c 1 , or not, c 2 , the ring. In 5c 1 -7c 1 complexes, the anion is slightly displaced towards the C1 half of the aromatic system, yielding three short contacts, as inferred from the distances between one O atom of nitrate with both an X and an H atom of the C1 methyl group and between another O atom of NO3with C1, as shown in Figure 6. The c 1 complex of 8 could not be found. In 5c 2 -8c 2 complexes, the anion presents to different binding modes: all three O atoms of nitrate in 5c 2 are interacting with two H atoms of C1 and C5 methyl groups, keeping the N atom sufficiently close to the O atom of the ring to expect a certain interaction. Conversely, in 6c 2 -8c 2 complexes, two O atoms of the anion are equidistant to the X atom. In addition, these two O atoms are engaged in equidistant hydrogen-bonding interactions with each O atom with one H atom of the C5 methyl group. In triphenyl c complexes, two different NO3orientations are observed: those similar to the unsubstituted ones, 9c-10c, and those similar to the trimethyl c 1 complexes, 10c 1 -12c 1 . The c complexes of 11 and 12 and the c 1 complex of 9 could not be found. The main difference between c and c1 complexes resides in the contacts between an O atom of nitrate and C1 and/or C5 atoms of the heteroaromatic system ( Figure 6). Furthermore, the intermolecular O2NO -···X contacts are shorter for c1 than for c complexes. All of these c and c 1 complexes present CH···ONO2hydrogen bonds with the H atoms in ortho to the C1 and C5 phenyl moieties. Table 5. Interaction energies (∆E, in kcal·mol −1 ), electron density (ρBCP), its Laplacian (∇ 2 ρBCP) and the total electron energy density (HBCP) for the BCPs of O2NO -···C contacts of complexes of Compounds 1-12 with NO3and the number of imaginary frequencies (Nimag).  As in Clcomplexes, hydrogen bond-based d complexes could not be found for 1-8. From the geometries included in Figure 6, we observe that 9d-12d complexes have the same recognition pattern, in which the anion is interacting with two H atoms in ortho to two phenyl rings, with the H atom at C2. The CH···ONO2distance oscillates from 2.14-2.41 Å for X = O to 2.28-2.55 Å for X = Te.

∆E
Finally, e complexes are characterized by showing intermolecular contacts with the C atoms of the heteroaromatic system and were obtained only for Compounds 5-9. Thus, in 5e-8e complexes, the anion is located parallel to and above the ring showing distances of 2.8-3.0 Å between the O atoms of nitrate and the ring plane. Two hydrogen bonds are also observed with the H atoms of C3 and C5 methyl groups. In Complex 9e, the nitrate is located above the heteroaromatic ring and perpendicular to it, giving rise to three intermolecular contacts between one O atom of NO3and two C atoms (C1 and C5) and between a second NO3 -O atom and C3.
The values of the interaction energies of all complexes (Table 5)  The origin and nature of these interactions have been analyzed within the frames of NBO, AIM and SAPT theories for selected complexes. According to the NBO analysis, the interaction in 4c is based on a charge donation from the lone pairs of an O atom of NO3to the vacant π* orbital of the Te-C5 bond, as previously observed for their chloride analogues, resulting in a π-hole chalcogen bond. Results for 2c and 3c complexes show some inconsistencies and therefore are not reported. In contrast, interaction in 5c 2 -8c 2 complexes come from Lp(ONO3)→σ*(X-C5) donor-acceptor contributions, except for 5c 2 , where contributions are based on charge transfer from the anion to H atoms of the methyl groups. Moreover, for this last complex, a very small contribution is found for the Lp(Opyrylium)→π*(N-O) donor-acceptor interaction. Thus, results for 6c 2 -8c 2 suggest the existence of a σ-hole chalcogen bond. Different charge-transfer contributions appear in 5c 1 -7c 1 complexes, where Lp(O)→π*(C1-C2) donor-acceptor interactions for 6c 1 and 7c 1 and donation from O lone pairs of NO3to empty p orbitals of C1 are dominant, which will be in agreement with the formation of an anion-π interaction.
Consistent with the NBO results, AIM calculations for c·NO3complexes show the existence of one BCP connecting the X atom with ONO2 - (Figure 7). The values of the electron densities at these BCPs, comprised between 0.046 au for 3c·NO3and 0.009 au for 5c 2 ·NO3 -( Table 5) and its positive values of the Laplacian are consistent with the existence of noncovalent interactions. However, the total energy densities are negative for 2c and 3c·NO3 -(−0.003 and −0.004 au, respectively), suggesting that these interactions have a small covalent character. In addition, according to the high electron density values (ρ = 0.2 au), and negative values of the Laplacian (between −0.40 and −0.24 au) and total energy densities (between −0.24 and −0.18 au) at the BCPs connecting the anion with the ring, a·NO3and b·NO3σ-complexes have a marked covalent character.
The SAPT analysis of complexes 2c-3c·NO3and 11c 2 -12c 2 ·NO3 -( Table 4) shows similar conclusions to those obtained for chloride complexes. Indeed, the electrostatic component is the largest contributor ranging from −105.4 to −80.0 kcal·mol −1 , and the smallest and largest contributions are found for X = O and X = Se, respectively. The dispersion contribution is practically the same for all complexes and the induction component varies according to the size of the chalcogen. Additionally, again, the consideration of δHF (Table 4) is very important at least for 2c-3c compounds to reach an agreement between ESAPT and ∆E values. In contrast, exchange, induction, dispersion and δHF contributions are very similar for 5e-7e·NO3complexes, showing only small differences in the electrostatic component, which in the end could be responsible for the observed small differences in the binding energies. Figure 7. Molecular graphs of 3c, 5e, 5c 1 , 5c 2 , 7c 2 , 10c 1 and 10c 2 ·NO3 -. The BCPs, ring and cage critical points are represented by green, red and blue dots, respectively. Only bond paths are depicted.

Complexes with BF4 -
In this third section, we are going to analyze the complexes formed by Compounds 1-12 with tetrafluoroborate as the anion. The exploration of the potential energy surface leads to the formation of four main structural configurations: Type a, c, d and e complexes, analogous to those previously described. The geometries of the complexes are shown in Figure 8, and their respective binding energies are listed in Table 6. There is one Type a complex, 1a, which exhibits a quite short F···C distance (2.37 Å) coupled with two hydrogen bonds with H atoms of C1 and C2. In c·BF4complexes, different orientations of tetrafluoroborate are observed. For unsubstituted compounds, BF4is located establishing two interactions of the F···X type and one hydrogen bond with the H atom at C2 in 2c-4c·BF4 - (Figure 8). Type c complexes for the trimethyl derivatives show two different approximations of the anion. In 5c-7c, the anion is placed above the ring with three F atoms pointing at the aromatic system. In 5c, two F atoms are facing C2 and C4, and one F atom is facing the oxygen atom. The anion is slightly displaced towards the right half of the molecule in 6c-7c with two F atoms giving rise to three contacts with C1, C3 and the X atom. In triphenyl c 1 complexes, two F atoms are directed at C2, C4 and a third F atom at X, whereas in c 2 complexes, two F atoms are directed at C1 and C5 and a third F atom at X (Figure 8), with hydrogen bonds with the H atoms in ortho to the C1 and C5 phenyl moieties. In Figure 8, we observe that 9d-12d complexes have the same recognition pattern, in which the anion is interacting with two H atoms in ortho to two phenyl rings and with the H atom at C2. The CH···FBF3distance oscillates from 2.07-2.33 Å for X = O to 2.21-2.72 Å for X = Te. Lastly, e complexes have the anion located above the aromatic ring. As seen in previous BF4complexes, the most favorable disposition of the anion is that with three F atoms directed at the aromatic system. In fact, tetrafluoroborate is symmetrically placed, facing C1, C3 and C5 in 1e and 5e-8e·BF4 -, and a little bit displaced to one side of the molecule in 2e-3e·BF4 -, as depicted in Figure 8, most likely due to some sort of interaction between an F atom and X = S, Se. Table 6. Interaction energies (∆E, in kcal·mol −1 ), electron density (ρBCP), its Laplacian (∇ 2 ρBCP) and the total electron energy density (HBCP) for the BCPs of F3BF -···C contacts of complexes of Compounds 1-12 with BF4and the number of imaginary frequencies (Nimag).  The values of the interaction energies are the smallest of all of the anion series (Table 6), ranging from −90.5 kcal·mol −1 for 4c·BF4to −70.8 kcal·mol −1 for 12d·BF4 -. The largest binding energies are found for the unsubstituted complexes, and d complexes are the least favorable ones. Substitution affects the binding energies the same way as in chloride and nitrate complexes.
The NBO analysis of BF4complexes offers interesting results. In 2c-4c complexes, the interaction is mainly based on charge donation from the lone pairs of an F atom of BF4to the empty σ* orbital of the X-C1 bond (σ-hole chalcogen bond) and from Lp(F)→σ*(C1-H) to a minor extent, except for 4c, where the latter charge donation is dominant. A different picture is observed for 5c-8c complexes with several contributions arising from Lp(F)→p(C1/C3) in 5c, Lp(F)→p(C3/C4) in 6c and Lp(F)→p(C5) in 7c (in line with the formation of anion-π interactions). The Lp(F)→π*(X-C1) donor-acceptor contribution is common for these complexes, increasing its importance with the size of the chalcogen in such a way that it becomes the only contribution in 8c. Contributions for e complexes are also based on charge transfer from lone pairs of F atoms to vacant p orbitals of C3 and C5 (and C1 for 5e), consistent with the existence of an anion-π interaction. Moreover, 6e-8e complexes have charge-transfer contributions of the Lp(F)→π*(X-C1) π-hole chalcogen bond type.
AIM calculations for all complexes are in agreement with NBO results and show the existence of one BCP connecting the X, C and H atoms with FBF3 - (Figure 9). It is worth noting that the type of contacts with X (characteristic for c complexes) are also present in 2e-3e·BF4 -. The values of the electron densities (between 0.027 and 0.004 au), Laplacian (between 0.098 and 0.020 au) and total energy densities (between 0.003 and 0.001 au) at the BCPs of all complexes (Table 6) are consistent with the existence of noncovalent interactions.  graphs of 1a, 3c, 5c, 7c, 11c 2 , 1e, 3e and 5e·BF4 -. The BCPs, ring and cage critical points are represented by green, red and blue dots, respectively. Only bond paths are depicted.
The SAPT results for 1e-3e, 5e-7e, 1a, 2c-3c and 10c 2 -11c 2 ·BF4are analogous to the SAPT analysis of nitrate complexes. In general, all contributions are smaller than those for Cland NO3 -, which is reflected in the smallest binding energies of all series of anions. The largest contribution comes from electrostatics with values between −96.3 and −73.3 kcal·mol −1 . Both dispersion and exchange contributions are kept more or less constant for all complexes, and the induction component only varies in c complexes.

Cambridge Structural Database Results
The Cambridge Structural Database (CSD) [37] is a convenient and reliable storehouse for geometric information. The utility of small molecule crystallography and the CSD in analyzing geometric parameters and nonbonding interactions is well established [38]. We have explored the CSD, searching for structures containing the pyrylium moiety 1-4 and we found 103 X-ray structures, 93 of them containing the pyrylium core. For thiopyrylium, we found eight structures and only one for both seleno-and telluro-pyrylium systems. In Table 7, we have gathered the number of structures for every family of compounds along with the occurrence of the associated anion. From the inspection of the results, we observe that the most common anion is perchlorate with 35 hits, followed by tetrafluoroborate (18 hits) and triflate (11 hits). We show in Figure 10 three pyrylium selected examples (CSD Reference Codes: ASEVOK [39], GIMPOI [40] and IWATEF) [41] and in Figure 11, three selected examples for the rest of the derivatives (CSD reference codes: NUCWAK [42], AFAMEZ [43] and AFAMAV) [43] that we have retrieved from the database, where the contacts between the anion and the cation have a strong influence in determining the crystal packing. In the ASEVOK structure, a BF4is directing one F atom to the heteroaromatic ring and a second one to an aromatic H atom. However, three O atoms of ClO4and three F atoms of PF6are facing the ring in GIMPOI and IWATEF, respectively, as previously shown in the computed e·BF4complexes. A good agreement is found between our computed complexes and the rest of the structures in Figure 11, where we observe several BF4 -···X interactions for every example, as inferred from the values of the intermolecular distances. Particularly, in AFAMEZ and AFAMAZ, the anion is engaged in F3BF -···X interactions that act as a bridge joining two seleno-and telluro-pyrylium cationic molecules, giving rise to columns ( Figure 11).

Computational Methods
We carried out geometry optimizations by means of DFT calculations at the B97-D3 [44,45] level of theory. This approach is based on Becke's generalized gradient approximation (GGA) functional introduced in 1997 [46] and is explicitly parameterized by including damped atom-pairwise dispersion corrections of the form C6 · R −6 . We optimized the geometry of all compounds using restricted Kohn-Sham DFT [47]. In these calculations, the Ahlrichs augmented triple-ζ basis for all atoms (def2-TZVPD basis set, abbreviated as TZVPD) [48] was used, unless stated otherwise. Our DFT calculations were carried out using the resolution of the identity (RI) B97-D3, which uses an auxiliary fitting basis set [49] to avoid treating the complete set of two-electron repulsion integrals, thus speeding up calculations by a factor of 10. Because of the time-consuming nature of the calculations, we used the parallel RI-DFT [50] methodology. Geometries of all structures were optimized with the analytical gradient method without imposing any symmetry constraints, unless stated otherwise. The vibrational frequencies were calculated at the RI-BP86-D3/def2-TZVPD level of theory. For our preliminary calculations, we carried out geometry optimizations using the MP2 method (RI-MP2) [51][52][53] with the aug-cc-pVTZ basis set and refined the energies by using the coupled-cluster method with singles and doubles excitations and triple excitations non-interatively (CCSD(T)) with a complete basis set (CBS). The CCSD(T) technique provides reliable interaction energies only if they are combined with extended AO basis sets, and the larger the basis set, the better the interaction energies that result. Owing to the rather strong dependence of the interaction energy on the AO basis set size, it is recommended that the relevant calculations be performed at the complete basis set (CBS) limit. Different extrapolation schemes have been introduced, and the scheme of Helgaker et al. [54] has become the most widely used. Here, the HF and correlation (MP2) energies are extrapolated separately as follows: where EX and ECBS are the energies for the basis set with the largest angular momentum X and for the complete basis set, respectively. The CCSD(T)/CBS level can be attained via a separate extrapolation of the MP2 and higher-order correlation energies towards the basis set limit (Equation (3)). Here, each of the components is differently sensitive to the AO basis set: the MP2 interaction energy is the more slowly converging, and the larger the basis set used in the extrapolation, the better. In our case, we have used from the aug-cc-pVTZ (AVTZ) to the aug-cc-pVQZ (AVQZ) basis sets. The second term, called the CCSD(T) correction term (∆CCSD(T)), is determined as the difference between the CCSD(T) and MP2 interaction energies and converges much faster than the previous term. The use of such a term is possible, because the MP2 and CCSD(T) interaction energies converge with basis set size in a very similar way; consequently, its difference is much less basis set dependent, and much smaller basis sets can be applied. In our case, we have used the AVDZ basis set to compute the CCSD(T) correction. All of the theoretical calculations described herein were carried out by using TURBOMOLE Version 6.5 [55,56], apart from the CCSD(T) calculations, which were computed by using the MOLPRO program [57,58].
The partitioning of the interaction energies into the individual electrostatic, induction, dispersion and exchange-repulsion components was carried out by performing density functional theory (DFT) calculations combined with the symmetry-adapted perturbation theory (DFT-SAPT) approach [59][60][61][62][63] with the MOLPRO program. The DFT-SAPT total interaction energy (Eint) is given in terms of the first-, second-and higher-order correction interaction terms that are indicated by the superscripts in Equation (4): where E el 1 and E exch 1 are the sum of the electrostatic interaction energy and the first-order exchange energy, respectively. E ind 2 , E ind-exch 2 , E disp 2 and E disp-exch 2 denote the induction (with response) energy, the second order induction-exchange (with response) energy, the dispersion energy and the exchange-dispersion contribution, respectively. The Hartree-Fock (HF) correction for third-and higher order contributions to the intermolecular interaction potential (δHF) is not included in DFT-SAPT calculations and can be estimated as the difference between the supermolecular HF energy and the sum of electrostatic, induction and their exchange counterparts obtained at the HF level. The aug-cc-pVDZ basis set was used to compute this correction apart from H, for which we used the cc-pVDZ basis set. For brevity, we will often refer to the aug-cc-pVXZ by the shorthand notation AVXZ. Physically meaningful separation of the interaction energy may be obtained by classifying the cross terms induction-exchange E ind-exch 2 and dispersion-exchange E disp-exch 2 as a part of the induction and the dispersion, respectively [64,65]. The AVTZ basis set was used for the DF-DFT-SAPT calculations. As the auxiliary fitting basis set, the JK-fitting basis of Weigend [66] was employed. The cc-pVQZ JK-fitting basis was used for all atoms. For the intermolecular correlation terms, i.e., the dispersion and exchange-dispersion terms, the related AVTZ MP2-fitting basis of Weigend, Köhn and Hättig [67] was employed. In the DFT-SAPT calculations, the BP86 functional (the B88 exchange functional [68] in combination with the P86 gradient correction) [69] was employed using the B97-D3/def2-TZVPD optimized geometries. The electronic properties of the systems have been analyzed within the atoms in molecules (AIM) theory [70,71] by using the AIMAll [72] program. The AIM theory is based on the topological analysis of the electron density, ρ(r). The presence of a path linking two nuclei in an equilibrium structure implies that the two atoms are bonded to one another, and it is characterized by the bond critical point (BCP). Topologically, a BCP corresponds to a point in the real space where the gradient of the density, ∇ ρBCP, is zero and where the curvature of ρBCP, expressed through the three eigenvalues of the diagonalized Hessian of ρBCP, is positive for an eigenvector linking two atomic centers and negative for the two others perpendicular to it. A chemical bond thus results from the competition of the parallel expansion of ρ, which separates charges in their respective atomic basins and the perpendicular contraction of ρ toward a bond path. The dominant effect is measured by the Laplacian of ρBCP, ∇ 2 ρBCP. Values of ρBCP greater than 0.2 are typical of covalent bonds, and ∇ 2 ρBCP is generally less than zero for such interactions, reflecting the concentration of electron density along the bond path linking the bonded atoms. Thus, when ∇ 2 ρBCP < 0, charge is concentrated at the critical point, while when ∇ 2 ρBCP > 0, charge is locally depleted. Additional information of the chemical bond can be obtained by applying the local expression of the virial theorem to the critical point, 1/2∇ 2 ρBCP = 2GBCP + VBCP, where GBCP and VBCP are the kinetic and potential energy densities at the BCP, respectively [73]. Cremer and Kraka demonstrated that the sign of the total energy, HBCP = GBCP + VBCP, is an index of the amount of covalency in the chemical interactions [73,74]. HBCP is negative for interactions with significant sharing of electrons. The remaining two stable critical points occur as a consequence of particular geometrical arrangements of bond paths, and they define the remaining elements of molecular structure, i.e., rings (ring critical point, RCP) and cages (cage critical point, CCP).
The bonding features of all of the studied molecules and the net atomic charges were also analyzed by means of the natural bond orbital (NBO) and the natural population analysis (NPA) [75,76]. For such purposes, the NBO 3.1 program [77] has been used. The donor-acceptor interactions are quantified by examining all possible interactions between filled (donor) Lewis-type NBOs and empty (acceptor) non-Lewis NBOs and by evaluating their energetic importance by using the second-order perturbation theory in the NBO basis [78]. When interpreting the results of such estimations, it should be noted that this approach is only performed at the SCF level of theory (i.e., the Fock or Kohn−Sham operator is analyzed in the basis of the NBO), and only bonding interactions are considered (i.e., antibonding contributions are not covered by the NBO analysis and must be calculated separately). Because these interactions lead to the loss of occupancy from localized NBOs of an idealized Lewis structure to empty non-Lewis orbitals (and thus, to departure from an idealized Lewis structure description), they are referred to as delocalization corrections to the zeroth-order natural Lewis structure.

Conclusions
In summary, we have examined the molecular recognition in pyrylium, thio-, seleno-and telluro-pyrylium salts where anions of different morphologies were chosen, as well as different substitution of the heteroaromatic ring. The binding energy for the series of anions follows the ordering ∆ECl > ∆ENO3 > ∆EBF4, and within every family of compounds, ∆EX > ∆EX-Me3 > ∆EX-Phe3. In general, for every derivative, the most favorable a and b complexes are found for X = O. Conversely, the largest binding energies in c complexes are obtained for the most electrophilic X = Te. Complexes strictly based on hydrogen-bonding interactions (d complexes) are the least favorable. For the monoatomic chloride anion, σ-complexes are generally the most favorable structures, except for triphenyl and X = Te trimethyl complexes, where the interactions of the anion with the chalcogen atom are important. The trigonal planar nitrate σ-complexes have the largest binding energies for the unsubstituted compounds. However, for the trimethyl and triphenyl derivatives, other configurations are more important, with the O atoms of nitrate either anion-π interacting with the ring C atoms (e complexes) or establishing interactions with the chalcogen atoms (c complexes) by means of σ-and π-hole chalcogen bonds. For the tetrahedral tetrafluoroborate anion, σ-complexes could not be found, and as in nitrate, e and c complexes are predominant. These complexes present multiple intermolecular contacts (anion-π interactions and σ-and π-hole chalcogen bonds) of the F···C, F···X and F···HC type, usually resulting from the interaction of three F atoms of a single BF4molecule with the cation. The interaction in all complexes is dominated by electrostatics. However, the induction term is especially important in chloride complexes. Finally, by analyzing the CSD, we have found several X-ray structures in which both anion-π interactions and chalcogen bonds are present, giving consistency to the computational study.