Theoretical Study on the Complexes of Benzene with Isoelectronic Nitrogen-Containing Heterocycles

The π–π interactions between benzene and the aromatic nitrogen heterocycles pyridine, pyrimidine, 1,3,5-triazine, 1,2,3-triazine, 1,2,4,5-tetrazine, and 1,2,3,4,5-pentazine are systematically investigated. The T-shaped structures of all complexes studied exhibit a contraction of the C—H bond accompanied by a rather large blue shift (40–52 cm−1) of its stretching frequency, and they are almost isoenergetic with the corresponding displaced-parallel structures at reliable levels of theory. With increasing number of nitrogen atoms in the heterocycle, the geometries, frequencies, energies, percentage of s character at C, and the electron density in the C—H σ antibonding orbital of the complexes all increase or decrease systematically. Decomposition analysis of the total binding energy showed that for all the complexes, the dispersion energy is the dominant attractive contribution, and a rather large attraction originating from electrostatic contribution is compensated by its exchange counterpart.


Introduction
Interactions between p systems have been the focus of attention for a long time because they play key roles in many fields such as supramolecular chemistry, drug design, biochemistry, crystal engineering, and many other new cross-disciplines associated with molecular science. [1][2][3][4] As a prototype of p-p interactions, the benzene dimer has been studied both theoretically and experimentally, [5][6][7][8] and these studies have greatly improved our understanding of the fundamental physics of p-p interactions. In many instances, however, not only the aromatic benzene ring but also N-containing heterocycles are involved in p-p interactions, such as p-p stacking in metal complexes with aromatic nitrogen-containing ligands [9] and nucleobase stacking in nucleic acids, [4] three nucleobases of which, namely, cytosine, thymine, and uracil, are pyrimidine derivatives. It is therefore important to systematically study p-p interactions between benzene and aromatic nitrogen heterocycles.
The very electronegative heteroatom nitrogen withdraws electron density from the aromatic heterocycle. With increasing number of nitrogen atoms in the ring, the p electron density in the ring decreases and the molecular electrostatic potential (MEP) becomes more positive. As can be seen in Figure 1, the MEP value is negative above the center of the benzene ring but positive in the aromatic nitrogen heterocycles pyridine, pyrimidine, 1,3,5-triazine, 1,2,3-triazine, 1,2,4,5-tetrazine, and 1,2,3,4,5-pentazine. The negative MEP of benzene is responsible for the existence of the T-shaped benzene dimer. The absence of the negative MEP above the center in the other compounds suggests that their dimers might not be T-shaped like the benzene dimer. Furthermore, lower p-electron density in the ring means less p-electron repulsion between p-donor and p-acceptor ring systems. Hence, the interactions between the benzene molecule and pyridine, pyrimidine, 1,3,5-triazine, 1,2,3-triazine, 1,2,4,5-tetrazine, and 1,2,3,4,5-pentazine should be stronger than that in the benzene dimer. Ugozzoli and Massera investigated the intermolecular potential of the ben-zene···1,3,5-triazine complex in detail. [10] In contrast to the results of Tsuzuki et al. for the benzene dimer, calculated at the same level of theory, [11] the interaction between benzene and 1,3,5-triazine is indeed much stronger. Note that this can be at The p-p interactions between benzene and the aromatic nitrogen heterocycles pyridine, pyrimidine, 1,3,5-triazine, 1,2,3-triazine, 1,2,4,5-tetrazine, and 1,2,3,4,5-pentazine are systematically investigated. The T-shaped structures of all complexes studied exhibit a contraction of the CÀH bond accompanied by a rather large blue shift (40-52 cm À1 ) of its stretching frequency, and they are almost isoenergetic with the corresponding displaced-parallel structures at reliable levels of theory. With increasing number of nitrogen atoms in the heterocycle, the geometries, frequencies, energies, percentage of s character at C, and the electron density in the CÀH s antibonding orbital of the complexes all increase or decrease systematically. Decomposition analysis of the total binding energy showed that for all the complexes, the dispersion energy is the dominant attractive contribution, and a rather large attraction originating from electrostatic contribution is compensated by its exchange counterpart.
least partially ascribed to more attractive dispersion energy, since the polarizability of 1,3,5-triazine is higher than that of benzene. Geerlings et al. ascertained that substituted benzenes with electron-withdrawing or electron-donating groups bind more strongly to pyridine than unsubstituted benzene, [12][13] and this is also inconsistent with the Hunter-Sanders rules, which state that the electron-donating groups will destabilize the p-p interaction. [1] The calculations mentioned mainly focus on the different structures and energies of only one complex rather than studying the interactions of benzene with different aromatic nitrogen-containing heterocycles systematically. Several interesting questions regarding the complexes formed by benzene with different aromatic nitrogen-containing heterocycles remain unanswered: 1) relative stability of T-shaped and parallel-displaced structures; 2) spectral shift of the CÀH stretching frequency on complex formation (the T-shaped benzene dimer exhibits a blue shift of the CÀH stretching frequency); [14,15] and 3) decomposition of the interaction energy. The last point is very important because it can help us further understand the nature of the p-p interaction in the complexes.
To answer the above questions, we select six complexes formed between benzene and its isoelectronic nitrogen-containing heterocycles. Both T-shaped and stacked structures are considered ( Figure 2). Note that our aim is not the evaluation of accurate absolute interaction energies but merely the relative values for different complexes. For comparison, the T-shaped and stacked structures of the benzene dimer are also included.

Computational Methods
The p-p interactions in the considered complexes are mostly governed by London dispersion forces, theoretical determination of which is quite difficult. Second-order Møller-Plesset theory (MP2) has been shown to be effective and accurate in determining the equilibrium structures and binding energies of many hydrogenbonded complexes. However, previous ab initio studies indicated that MP2 tends to overestimate binding in the case of p-p stacking interactions. [11,[16][17][18][19][20] The p-p stacking interactions are properly described only if a higher level of theory, such as the coupled-cluster method with singles, doubles, and noniterative triples [CCSD(T)] is adopted. CCSD(T) calculations with extended basis sets are, however, quite demanding. Fortunately, in their study on the relative energies of the multiply substituted benzene dimers and benzene dimer, Sinnokrot and Sherrill found that the MP2/aug-cc-pVDZ method yields relative energies very similar to the more expensive CCSD(T)/aug-cc-pVTZ method. [21] Evidently, this is due to a compensation of errors. Since the main aim of the present work is to compare the binding energies of the various dimers, the MP2/augcc-pVDZ method was used for most calculations.
The reliability of the MP2/aug-cc-pVDZ method for the present dimers was demonstrated by performing high-level calculations for the benzene···1,3,5-triazine complex. The CCSD(T) binding energy at the complete basis set (CBS) was approximated as the sum of the MP2/CBS interaction energy and the CCSD(T) correction term. The CCSD(T) correction term, DCCSD(T), is defined in the present paper as the difference between CCSD(T) and MP2 binding energies in the aug-cc-pVDZ basis set. The MP2/CBS binding energy was estimated by extrapolation of the calculated binding energies with the aug-cc-pVXZ basis sets by using a fitting of the form a + b expA C H T U N G T R E N N U N G (ÀcX), where X = 2, 3, and 4 for aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ, respectively. [22,23] All the structures of the benzene···1,3,5-triazine complex were fully optimized at the MP2/ aug-cc-pVDZ and MP2/aug-cc-pVTZ theory levels. The MP2/aug-cc-pVQZ calculations use the MP2/aug-cc-pVTZ optimized geometries. The basis set superposition error (BSSE) of the binding energy was systematically eliminated by using the standard counterpoise (CP) correction method of Boys and Bernard. [24] The frozen-core approximation was applied throughout.
The total interaction energies E int were decomposed into first-order E 1 , second-order E 2 , and higher order [dðHFÞ] terms by using the DFT-SAPT perturbation treatment. [25][26][27][28][29] E 1 is the sum of the electrostatic interaction energy E 1 pol and the first-order exchange energy E 1 exch . E 2 denotes the sum of the second-order induction energy E 2 ind and the second-order dispersion energy E disp 2 . Details of the DFT-SAPT method can be found elsewhere. [25][26][27][28][29] The DFT-SAPT calculations were performed with the MP2/aug-cc-pVDZ geometries. The PBE0AC exchange-correlation functional with density fitting and aug-cc-pVDZ basis set was used systematically. The basis set dependency of the individual interaction energy terms has been studied by Jansen et al., [28] who found that the aug-cc-pVDZ basis set is a reasonable compromise between computational cost and accuracy; this basis set used provided accurate energy terms with the exception of the dispersion energy, which was underestimated by about 10-20 %. In the PBE0AC exchange-correlation potential, the incorrect asymptotic behavior of the corresponding PBE exchange potential has been corrected by the gradient-regulated asymptotic correction approach of Grüning et al. [30]
The Gaussian03 suite of programs [31] was used for the ab initio molecular orbital calculations to evaluate the geometries, harmonic frequencies, and the total interaction energies. The natural bond orbital (NBO) analysis [32] was carried out by using the MP2-optimized structure, the MP2 electron density, and the built-in subroutine of the Gaussian03 program suite. The total interaction energy (DFT-SAPT procedure) was decomposed by using the MOLPRO 2006.1 software package. [33] 2. Results and Discussion 2.1. Binding Energies of Benzene···1,3,5-Triazine Complex Table 1 summarizes the binding energies for the T-shaped (B3T-1), sandwich (B3S), and parallel-displaced (B3P-1, B3P-2) configurations of the benzene···1,3,5-triazine complex at various levels of theory.
First, the MP2 level is discussed. Evidently, passing from augcc-pVDZ to the much larger aug-cc-pVQZ basis set leads to a significant increase in binding energy. By extrapolation, we obtain the MP2/CBS binding energies, which are even larger than the aug-cc-pVQZ values, by up to 0.08 kcal mol À1 . At the MP2/CBS level, the parallel-displaced structures are most stable, followed by sandwich and T-shaped ones. As expected, the CCSD(T) correction term is small for the T-shaped structure and considerably larger for sandwich and parallel-displaced structures. This term is systematically repulsive (by 2-5 kcal mol À1 ) for stacking structures of various DNA base pairs and amino acid pairs. [34] By adding this term to the MP2/ CBS one, we obtain the CCSD(T)/CBS binding energy (Table 1, last column). Evidently, these binding energies are smaller than the corresponding MP2 ones. Moreover, energy differences between these structures become smaller; at the MP2 level it is 1.9 kcal mol À1 , while at the CCSD(T) level it is only 0.5 kcal mol À1 . The parallel-displaced structure B3P-1 remains the global minimum, but the T-shaped structure is energetical-ly close and the difference is smaller than 0.2 kcal mol À1 . The sandwich structure B3S is now the least stable structure. Table 1 shows that the binding energy of the T-shaped structure B3T-1 determined at the MP2/aug-cc-pVDZ level of theory is À3.23 kcal mol À1 , which is close to the CCSD(T)/CBS binding energy of À3.04 kcal mol À1 . Similarly, the binding energies for the T-shaped structure of the benzene dimer are À2.86 and À2.74 kcal mol À1 at the MP2/aug-cc-pVDZ and CCSD(T)/CBS levels of theory, respectively. [35] The differences between the binding energies of the T-shaped structures of the benzene dimers and the benzene···1,3,5-triazine complex at the CCSD(T)/CBS and MP2/aug-cc-pVDZ levels are thus similar (0.3 and 0.37 kcal mol À1 ). This allows us to use the much cheaper MP2/aug-cc-pVDZ method for determining relative binding energies of the T-shaped structures of the complexes formed between benzene and its isoelectronic nitrogen-containing heterocycles. The MP2/aug-cc-pVDZ binding energies of the sandwich and parallel-displaced structures of the benzene···1,3,5-triazine complex differ from the CCSD(T)/CBS values by up to 1.6 kcal mol À1 . The differences become more pronounced when the larger aug-cc-pVTZ and aug-cc-pVQZ basis sets are used. However, the MP2/aug-cc-pVDZ binding energies of the B3P-1 and B3P-2 structures differ by less than 0.2 kcal mol À1 , which is very close to the CCSD(T)/CBS value (0.2 kcal mol À1 ).
This indicates that the relative binding energies of the particular structural type are predicted at the MP2/aug-cc-pVDZ level of theory fairly accurately. This conclusion is in agreement with the finding of Sherrill et al. for the benzene dimer. [21] Table 1 further shows that the CCSD(T)/CBS binding energies are almost identical for the T-shaped structure B3T-1 and the parallel-displaced structure B3P-2. Additionally, the binding energies of the T-shaped, sandwich, and parallel-displaced structures are systematically larger than the corresponding energies for the benzene dimer. [36] Interestingly, the binding energy of B3P-1 is larger than that of B3P-2 at all levels of theory, which contradicts our chemical intuition. The electron lone pairs of one nitrogen atom in B3P-1 are above the benzene ring, and consequently stronger electron-electron repulsion is expected. This question is resolved in Section 2.3 by a detailed energy decomposition analysis.
The structure and geometry of the complexes are discussed first. The distance between subsystems in the T-shaped, sandwich, and parallel-displaced structures decreases with increasing number of nitrogen atoms. The decrease is the smallest in the T-shaped structures (0.1096 ). A larger decrease is found for sandwich structures (0.2342 ), and the largest occurs for parallel-displaced structures (R 1 decreases by 0.3033 ). The shortening of the intermolecular distance is connected with an increase in binding energy. For the T-shaped, sandwich, and parallel-displaced structures, it is 2.02, 2.54 and 3.54 kcal mol À1 , respectively. According to the positions of the nitrogen hetero atoms, the T-shaped dimers can be categorized into two groups: a) B0T, B1T, B2T-1, B3T-2, and B3T-3, and b) B2T-2, B3T-1, B4T, and B5T. In group a, there is no nitrogen atom to either side of the CÀH bond involved in the CÀH···p interaction, while in group b, two nitrogen atoms are near the CÀH bond. Evidently, the binding energy increases with increasing number of nitrogen atoms in both groups. Note that the binding energy and geometry of B3T-2 are almost the same as those of B3T-3, which indicates free rotation about the C 6 axis of benzene. Similarly, the parallel-displaced dimers can be categorized into three groups based on the number of nitrogen atoms above the benzene ring. Group a' includes B0P, B1P, B2P, and B3P-3, in which no nitrogen atom is above the benzene ring. Group b' includes only B3P-1, in which one nitrogen atom is above the benzene ring. The remaining structures are included in group c', which has two nitrogen atoms above the benzene ring. The binding energy in all three groups increases with increasing number of nitrogen atoms.
In the case of T-shaped complexes, the binding energy of B2T-1 is larger than the binding energy of B2T-2, and the binding energy of B3T-2 is larger than that of B3T-1, whereas for parallel-displaced complexes, the binding energy of B3P-3 is larger than that of B3P-1, and that of B3P-2 is the lowest one among the three complexes. It is possible to conclude that binding energy increases with increasing distance of the nitrogen atoms from the benzene ring.
The geometry and spectral shift of the CÀH bond in the T-shaped dimers were also investigated. The C 2v T-shaped structure of the benzene dimer was the first system for which the existence of a blue-shifting hydrogen bond was predicted. [14,15] Since then various types of blue-shifting hydrogen bonds have been confirmed experimentally, not only in the gas phase, but also in liquid and solid phases. The CÀH···p blue-shifting hydrogen bond in the benzene dimer was, however, not found experimentally, and the experimental study pointed to a small red shift of the CÀH stretching frequency on benzene dimerization. [8] In a recent study, we suggested an explanation for this finding (global minimum C s T-shaped structure has a red shift, while a C 2v T-shaped structure, which is only a transition structure, has a blue shift). [37] Table 2 shows [a] D %s-char is the percentage change in s character in a carbon hybrid orbital in the CÀH bond; DED is the change in electron density in the CÀH s antibonding orbital, and CT the total charge transfer between the electron donor (benzene) and the electron acceptor (N heterocycle); dm/dr CÀH is the derivative of the permanent dipole moment of the proton donor. Cartesian co-ordinates of all the stationary point geometries are given in the Supporting Information.
that all CÀH bond lengths in the T-shaped complexes are systematically contracted on complex formation. Contraction of the CÀH bond is accompanied by an increase in CÀH stretching frequency, that is, by a blue shift. However, as pointed out by McDowell and Buckingham, [38] there is no correlation between the change in bond length and the shift in vibrational frequency. The largest CÀH bond contraction was found for the benzene dimer, and the smallest for the benzene···1,2,4,5tetrazine complex; the blue shift is, however, identical for both dimers. The largest blue shifts (52 and 51 cm À1 ) were found in benzene···pyridine and benzene···pyrimidine complexes. The explanation of blue-shifting H-bonding is, unlike the case of red-shifting H-bonding, not unambiguous and several possibilities should be investigated. The first concerns the electrostatic theory of H-bonding and is based on the fact that elongation of the XÀH bond of the proton donor is mostly connected with an increase of the proton-donor dipole moment. A larger dipole moment yields a larger dipole-dipole attraction, which is the reason for the elongation of the XÀH bond and the subsequent red shift of the XÀH stretching frequency. However, a few systems exhibit the opposite, that is, elongation of the XÀH bond is connected with a decrease in dipole moment. In other words, contraction of the XÀH bond results in an increase in the dipole moment of the proton donor. Analogously to the previous case, a larger dipole moment yields stronger dipole-dipole attraction, which is the reason for the contraction of the XÀH bond and the subsequent blue shift of the XÀH stretching frequency. The dipolemoment derivatives for the proton donors considered have mostly negative values (Table 2), consistent with bond contraction and a blue shift, with the exception of pyrimidine with the CÀH proton donor located between two nitrogen atoms in the B2T-2 structure and 1,2,3,4,5-pentazine. The large positive value found for the pyrimidine monomer indicates elongation and red shifting of the CÀH stretching vibration on dimerization. However, examining Table 2 reveals the opposite: the CÀH bond in the B2T-2 dimer is contracted and is associated with the second largest blue shift. Evidently, in most cases the CÀH spectral shift is explained by dipole-moment derivatives, but there are some exceptions. The second possibility is based on the theory developed by Alabugin et al., who proposed that the XÀH bond length in the XÀH···Y H-bond is controlled by a subtle balance of hyperconjugation and rehybridization. [39] Hyperconjugation (electron-density shift from lone pairs of the proton acceptor to the s* antibonding XÀH orbital) results in elongation of the XÀH bond and consequently a red shift of the XÀH stretching frequency. Rehybridization, on the other hand, strengthens the XÀH bond and results in a blue shift of the XÀH stretching frequency. Both effects act simultaneously, and the final effect depends on the prevalent role of either contribution. Table 2 lists the values of the total charge transfer (CT) between the proton acceptor and the proton donor as well as the increase in electron density (DED) in the s* antibonding orbital of a CÀH bond. Evidently, the total CT is much larger than DED, which indicates a very small or negligible effect of hyperconjugation on the CÀH bond length. In the case of the red-shifting H-bonding, the DED and CT values are practically identical. The changes in percentage s character (D%s-char) in the carbon hybrid orbital forming the CÀH bond of the proton donor are systematically positive, that is, the s character increases on complex formation. This increase is connected with strengthening of the CÀH bond, which is followed by an increase of the CÀH stretching frequency. We found no correlation of the D%s-char values with the contraction of the CÀH bond or the blue shift. Evidently, various factors are responsible for the magnitude of the blue shift and, among the three factors investigated, only rehybridization fully agrees with the spectral shift of benzene and nitrogen-containing heterocyclic compounds on complex formation. Table 3 shows the SAPT decomposition of the total interaction energies for all the complexes. They are ordered according to the aforementioned categories. The electrostatic term E 1 pol is systematically stabilizing for all the complexes studied. In each group of complexes, the absolute value of the electrostatic term increases with increasing number of nitrogen atoms. A larger value of this term correlates with the more positive molecular electrostatic potential of the heterocyclic compound with more nitrogen atoms. The first-order interaction energy (the sum of electrostatic and exchange-repulsion terms, E 1 in Table 3) is, however, systematically repulsive due to a rather large repulsion value of the exchange-repulsion term. This term also increases with increasing number of nitrogen atoms, with the exception of B0T in group a of the T-shaped dimers. Note that lower p-electron density in the ring does not mean smaller p-electron repulsion between the aromatic rings. It is the intermolecular distance that plays the important role for the exchange-repulsion term.

Decomposition Analysis of the Total Binding Energy
The second-order induction energy is much smaller than other energy components and increases with increasing number of nitrogen atoms. This term is least favorable for the sandwich structures and most favorable for the T-shaped structures. Following expectations, the second-order dispersion term is large and also increases with increasing number of nitrogen atoms with the exception of B0T and B1T in group a of the T-shaped dimers. The dispersion energy is largest for parallel-displaced structures followed by sandwich structures. The dispersion energy is systematically larger than the corresponding electrostatic term. The second-order interaction energies (E 2 in Table 3) are attractive for all dimers, and their major part originates in dispersion energy.
When all energy terms [including d(HF)] are put together, we obtain the total binding energy E int . Table 3 shows that the order of the SAPT binding energies in each group is the same as that of the CP-corrected MP2/aug-cc-pVDZ binding energies. For the T-shaped dimers, the difference between the two binding energies is about 1 kcal mol À1 , and the agreement will be more satisfactory if the larger aug-cc-pVTZ and aug-cc-pVQZ basis sets are used for SAPT calculations. [28] However, for the sandwich and parallel-displaced dimers, the MP2/aug-cc-pVDZ binding energies differ from the SAPT values by up to 3.0 kcal mol À1 , which is very similar to the case of the benzene dimer. [20] The SAPT binding energies are very close to the corresponding CCSD(T)/aug-cc-pVDZ values (Table 3). Clearly, the SAPT binding energy is more accurate and the difference between the two energies is clearly due to the repulsive CCSD(T) correction term, which is missing in the MP2 binding energy but is at least partly included in SAPT binding energy. The SAPT binding energies of B2T-1, B4T, and B5T are almost the same as those of B2P, B4P, and B5P. By comparing the respective energy component of B3P-1 with that of B3P-2, we ascertained that the exchange-repulsion term of B3P-1 is indeed larger than that of B3P-2, which is consistent with our chemical intuition. However, the attractive induction, dispersion, and high-order terms of B3P-1 are all larger than those of B3P-2, which makes B3P-1 bind slightly more strongly than B3P-2.

Conclusions
We have performed MP2 and DFT-SAPT calculations to study the p-p interactions between benzene and the aromatic nitrogen heterocycles pyridine, pyrimidine, 1,3,5-triazine, 1,2,3-triazine, 1,2,4,5-tetrazine, and 1,2,3,4,5-pentazine. By analyzing the geometries, vibrational frequencies, natural bond orbitals, and binding energies, we conclude with the following remarks: 1) The CÀH bond length in the T-shaped dimer is contracted on complex formation. This contraction is accompanied by a large blue shift of the CÀH stretching frequency (between 40 and 52 cm À1 ). The dipole-moment derivatives and change in electron density in the s* antibonding orbital of the CÀH bond mostly correlate with the blue shift, while the change of rehybridization of the carbon atom in the CÀH bond upon dimer formation always coincides with a blue shift.
2) The CCSD(T)/CBS binding energies are almost the same for T-shaped complex B3T-1 and parallel-displaced complex B3P-2. Also, the SAPT binding energies of B2T-1, B4T, and B5T are almost the same as those of B2P, B4P, and B5P. The results indicate that the T-shaped structure of each complex studied here is isoenergetic with the corresponding parallel-displaced structure, which is very similar to the case for benzene dimer. [6] 3) With increasing number of nitrogen atoms in the heterocycles, the distance between the rings decreases while the binding energy increases. For the T-shaped complexes, the percentage of s character at C and the electron density in the CÀH s* antibonding orbital and their change on complex formation increase with decreasing distance between the rings. Correspondingly, the change in CÀH bond length decreases with decreasing intermolecular distance. 4) For all complexes investigated, the dispersion energy is the dominant attractive contribution. A rather large attraction originating from electrostatic contribution is compensated by its exchange counterpart. Generally, all energy components increase with increasing number of nitrogen atoms in the heterocycle. [a] Values in parentheses are interaction energies at the CCSD(T)/aug-cc-pVDZ level of theory.