Intramolecular Basis Set Superposition Error Effects on the Planarity of Benzene and Other Aromatic Molecules: a Solution to the Problem Intramolecular Basis Set Superposition Error Effects on the Planarity of Benzene and Other Aromatic Molecules: a Solution to the Problem

Recently, the surprising result that ab initio calculations on benzene and other planar arenes at correlated MP2, MP3, configuration interaction with singles and doubles ͑CISD͒, and coupled cluster with singles and doubles levels of theory using standard Pople's basis sets yield nonplanar minima has been reported. The planar optimized structures turn out to be transition states presenting one or more large imaginary frequencies, whereas single-determinant-based methods lead to the expected planar minima and no imaginary frequencies. It has been suggested that such anomalous behavior can be originated by two-electron basis set incompleteness error. In this work, we show that the reported pitfalls can be interpreted in terms of intramolecular basis set superposition error ͑BSSE͒ effects, mostly between the C–H moieties constituting the arenes. We have carried out counterpoise-corrected optimizations and frequency calculations at the Hartree–Fock, B3LYP, MP2, and CISD levels of theory with several basis sets for a number of arenes. In all cases, correcting for intramolecular BSSE fixes the anomalous behavior of the correlated methods, whereas no significant differences are observed in the single-determinant case. Consequently, all systems studied are planar at all levels of theory. The effect of different intramolecular fragment definitions and the particular case of charged species, namely, cyclopentadienyl and indenyl anions, respectively, are also discussed.


INTRODUCTION
In a recent Communication, Moran et al. 1 reported a vast number of ab initio calculations on benzene and other planar arenes at different correlated levels of theory and Pople's basis sets that yielded nonplanar minima. The planar optimized structures turned out to be transition states exhibiting one ͑or more͒ large imaginary frequency. Singledeterminant-based methods such as Hartree-Fock ͑HF͒ and density functional theory ͑DFT͒ methods ͑BLYP and B3LYP͒ lead to the expected planar minima and hence no imaginary frequencies.
These results are in line with the earlier reports on failures of electronic structure methods for the correct description of low-lying out-of-plane vibrational frequencies of benzene, 2,3 several planar arenes, [4][5][6][7][8] and other nonrigid molecules. 9 Similar pitfalls have been recently found by Shabahzian 10 for the challenging ͑B 6 C͒ 2− anion, particularly at the MP2 level of theory. Rather to any deficiency of the post-HF electronic structure methods including electron correlation, the origin of the problem has been suggested to be rooted on ͑atom-centered͒ basis set deficiencies. Simandiras et al. 9 found that the use of f-type basis functions was necessary to obtain accurate bending frequencies. Martin et al. 2 dealt in detail with the benzene case and concluded that outof-plane bending modes are pathologically basis set depen-dent, putting forward also a basis set superposition origin, based on the work of Sellers and Almlöf. 11 Jensen 12 found that imaginary frequencies can also appear at the singledeterminant level for an uncontracted double-zeta quality basis set including d-type diffuse function ͑aug-pc-1͒ for a rather narrow range of diffuse function's exponent values.
Perhaps the most important error introduced by the use of truncated atom-centered basis sets is the so-called basis set superposition error ͑BSSE͒. 13 It has been recognized for years 14 that BSSE introduces some spurious extra binding in the ab initio calculations and that its correction is essential to properly describe intermolecular interactions. In the last decade, it has been also shown that BSSE does not merely affect the interaction energy but also the topology of the potential energy surface ͑PES͒ of these systems, which translates into geometrical and vibrational effects. There are a few examples in the literature where BSSE accounts for large geometric effects, namely, the hydrogen fluoride dimer 15,16 or weak C-H¯O ͑Ref. 17͒ interactions.
There are different ways to deal with BSSE. One option is to increase the basis set quality, since in the complete basis set limit BSSE should vanish. Whereas this could be accomplished at the HF or DFT levels of theory, for post-HF methods, BSSE has been shown to converge to zero very slowly with basis set improvement. 18,19 On the other hand, there are several strategies [20][21][22][23] to correct for BSSE, one is the counterpoise ͑CP͒ method 20 which is the simplest and is most widely used. Given a system formed by N fragments, the CP correction to the energy, in its simplest form, 24 where i i ͑R ជ ͒ and i full ͑R ជ ͒ represent the energy of the ith fragment of the system calculated with its own basis set and with the full basis set of the system ͑ghost-orbital calculation͒, respectively. Explicit dependence on the position of the atoms of the systems has been included to point out that the CP correction is geometry dependent. Note also that the charge, multiplicity, and electron state of the fragment's calculations must be specified, and that the value of ␦ CP ͑R ͒ will actually depend on the choice. This does not represent practical difficulties for closed-shell interacting molecules. However, the case of radical or charged molecular complexes, where the spin density or the charge is spread over the whole system, is not so trivial. We will go back to this point when dealing with the cyclopentadienyl anion case.
When the CP correction is applied as correction term to the total energy of the system, 25 it is possible to obtain BSSE-corrected energies and any molecular property that may be stated as a derivative of the energy, which includes optimized geometries or vibrational frequencies, as implemented in GAUSSIAN 03. 26 Experience gathered on the BSSE correction of intermolecular interactions over the last years shows that the results obtained with BSSE-corrected calculations using moderate basis sets are generally close to those that one can obtain with much larger basis sets with or without correction, provided that ͑a͒ the basis set is flexible enough to describe the true physical interactions in the system, and ͑b͒ BSSE effects on the geometry are also taken into account by optimization on the BSSE-corrected PES.
Much less attention has been paid to BSSE effects in single molecules, what one can refer to as intramolecular BSSE. 11,[27][28][29][30] Noteworthy exceptions are recent studies on aromatic-backbone intramolecular interactions 31-34 on peptide models, which have put forward that the accurate determination of relative stabilities of conformers is heavily affected by BSSE. In fact, the first CP calculation by Jansen and Ros 13 was applied for a single cation ͑HCO + ͒. In single molecules, there is not a priori problem with the fact that one atom or fragment can use the basis functions centered on other atoms or parts of the molecule. This helps to naturally describe polarization or charge transfer effects. The problem arises when the use of external basis functions is merely due to a lack of flexibility of the fragment's basis set; i.e., when the basis set incompleteness error ͑BSIE͒ is strongly geometry dependent, spurious geometry changes can be induced by intramolecular BSSE effects. Unbalanced descriptions of PES can also emerge. This directly translates into poor vibrational frequencies and likely changes in the topology of the stationary points, as found by Moran et al. 1 In that revealing report, the authors showed that in those cases where a nonplanar optimized structure was found for benzene, the two-electron BSIE tends to dramatically increase for geometries away from planarity. In fact, even for the levels of theory where no imaginary frequencies are observed, BSIE still shows moderate dependence on the geometry. Thus, the two-electron BSIE diagnostic seems indeed a valuable tool to detect possible spurious geometries induced by BSSE effects, but it is not definitive yet. A recent study 35 based on the use of chemical hardness profiles is also able to detect such spurious structures, but again, it does not offer a solution to the problem.
In this work, we wish to show that the reported pitfalls can be interpreted in terms of intramolecular BSSE effects and, as such, the calculations can be corrected using BSSE correction techniques typically applied to intermolecular interactions, particularly the CP method. For this purpose, we have carried out CP-corrected optimizations and subsequent CP-corrected frequency calculations for a number of the arenes ͑see Scheme 1͒ at the HF, B3LYP, MP2, and CISD levels of theory with several basis sets. We will show that, with no exception, correcting for intramolecular BSSE fixes the anomalous behavior of the correlated methods, whereas no significant differences are observed in the single-determinant case.

FRAGMENT'S DEFINITION
A key point when estimating intramolecular BSSE effects is the proper definition of the fragments forming each molecule in order to apply Eq. ͑1͒. A natural and unambiguous way to proceed could be to take atomic fragments. The problem with this approach is that the number of necessary extra calculations to obtain the CP correction would be 2N, N being the number of atoms. It would be desirable to be able to define larger fragments while maintaining the molecular symmetry of the system. In this sense, it is essential to note that in all the reported cases the problems are associated with out-of-plane bending low-lying modes. This tells us that in these cases the intramolecular BSSE does not seem to appreciably affect bonds or bond angles ͑otherwise stretching and other bending modes would have been affected͒. Taking all this into account, we can infer that proper moieties that can be used as fragments are the C-H units constituting the arenes.
In order to assess the effect of different fragment definitions, we have first determined the single-point CP-corrected energy profile along the b 2g vibrational mode of benzene at the MP2/6-311G level of theory, for which an imaginary frequency of 722i cm −1 is obtained. The uncorrected and CP-corrected results using twelve atomic fragments, six C-H fragments, as well as three ͑C-H͒ 2 fragments are displayed in Fig. 1. It can be readily seen that, even though the results obviously depend on the fragments' definition, all CP-corrected energy profiles properly describe the distortion as a vibration, contrary to the uncorrected profile.

COMPUTATIONAL DETAILS
As detailed in the Supporting Information, 36 we have used both GAUSSIAN 03 and our own code that allows us to exploit the symmetry of the molecules by minimizing the number of fragment calculations in both the CP-corrected optimization and frequency calculations. This is particularly important for the case of CISD, as in any ghost-orbital calculation symmetry must be turned off ͑fragments usually do not have the same symmetry as the supermolecular system͒.
The case of cyclopentadienyl anion is a special one because, even though the 6 electrons are fully delocalized in the ring, the extra electron must be assigned, in principle, to one of the five C-H fragments, thus breaking the symmetry of the CP-corrected calculation. Since our program allows us for a flexible definition of the CP-correction term, we have carried out the following approach in order to equally share the extra electron among five identical CH fragments: We have determined the CP correction using both neutral, ␦ CH CP ͑R ជ ͒, and negatively charged, ␦ CH − CP ͑R ជ ͒ C-H moieties and combined them as thus preserving the symmetry of the calculation and also taking into account the effect of the extra electron. This approximate CP-correction term has been used for both geometry optimization and frequency calculations. The case of indenyl anion is also problematic. Note that up to five imaginary frequencies are obtained at the MP2/6-311G level of theory, the largest one exceeding 1000i cm −1 . For this system, a CP function, which keeps the C 2v symmetry of the molecule and equally shares the electron charge among the three C-H moieties of the fivemember ring, has been used, in the spirit of Eq. ͑2͒.

RESULTS AND DISCUSSION
The lowest vibrational frequencies for benzene along with the symmetry of the vibrational mode, both uncorrected and CP corrected, are gathered in Table I. We have used three basis sets ranging from qualitatively good behavior ͑6-31+ G * ͒ to dramatic BSSE effects ͑6-311+ + G͒ and applied the same levels of theory as in Ref. 1 for comparison.
Noticeably, the CP correction fixes the anomalous behavior of the correlated methods in all cases. Consequently, benzene is found to be planar for all levels of theory. No significant differences are observed for higher frequencies or for the already qualitatively correct single-determinant calculations such as HF or B3LYP. Also, the CP-corrected optimized frequencies do not appreciably differ from the uncorrected ones ͑tables with all CP-corrected geometries and frequencies for all levels of theory are given in the Supporting Information 36 ͒. The CP-corrected frequencies are also in better agreement with both experimental 37 and best theoretical estimates 2 than the uncorrected ones, even ignoring the lowest five vibrational modes ͑see Table II͒. Noticeably, Jensen 12 has detected similar problems even for a single-determinant case, ͑BLYP, in particular͒ in combination with a double-zeta quality uncontracted basis set  ͑aug-pc-1͒, which includes a set of diffuse d functions with exponent ␣ d = 0.10. We have noticed that this basis set suffers from strong linear dependencies at the benzene geometry. In fact, the overlap matrix is not positive definite and exhibits up to 24 eigenvalues below Gaussian's default threshold of 10 −6 ͑being numerically zero͒. We have observed that the ghost-orbital calculations of the CH moieties yield unphysical delocalizations of the charge density in the neighboring atoms, which indicates that this basis set is very unbalanced upon removal of the dependent basis sets. Using the basis set in its contracted form leads to the expected planar structures, the CP correction not changing this situation. Qualitative wrong behavior is observed also in the case of naphthalene, as seen in Table III. At the modest MP2/6-31G level of theory, the out-of-plane b 2g vibrational mode is poorly described, exhibiting an imaginary frequency of 400i cm −1 . The correction of the intramolecular BSSE leads to a value of 570 cm −1 , whereas no meaningful differences are observed for the rest of low frequencies.
Finally, the results for the charged species are gathered in Table IV. It is remarkable that the CP correction fixes even the pathological case of the indenyl anion at the MP2/6-311G level. The intramolecular BSSE effects are so large in this system that the matching between uncorrected and corrected out-of-plane normal modes is dubious. In the case of cyclopentadienyl anion, the imaginary frequency associated with a degenerate e 2 Љ normal mode also disappears upon CP correction, whereas large effects on other low-frequency normal modes is also observed. Hence, the rather involved CP func-tions used for these systems to cope with the extra electron and keep overall molecular symmetry have proved to be successful.
For testing purposes, we have also applied to these systems a simpler CP function ignoring the extra charge ͑i.e., neutral C-H and C-C moieties͒. The frequencies obtained ͑not reported͒ are again real in all cases, while the differences with the frequencies reported in Table IV do not exceed 25%. Hence, the charge ͑and electronic state͒ of the fragments does not seem to be critical to correct for BSSE for these systems. This had been observed before also in the case of charged intermolecular complexes 30 and it is at the heart of the general success of the CP method. The BSSE effects are noticeable when the BSIE is strongly geometry dependent. The latter affects the description of the occupied molecular orbitals and also the virtual ones, which is the reason why post-HF methods are more prone to suffering from strong BSSE effects. Therefore, the BSSE is not directly related to the number of electrons of the system. This also explains the fact that hardness profiles, which have been successfully applied to the present problem, 35 are much more robust to level of theory and basis set effects than energy profiles ͑the chemical hardness is computed as the difference between the vertical ionization potential and electron affin-ity͒. Indeed, invoking the CP philosophy, an energy difference-based measure will perform better than a single energy value.

CONCLUSIONS
We have provided evidences that intramolecular BSSE accounts for the existence of nonplanar optimized minima structures predicted by typical electronic structure methods at the correlated level of theory. We have shown that by taking as fragments the C-H and C-C moieties, the CP-corrected optimized structures correspond to planar minima with no imaginary frequencies. We have also addressed charged systems, wherein in order to maintain the molecular symmetry more involved CP function have been used. For this purpose, we have employed a code 38 that allows for flexible definitions of the CP function and deals with molecular symmetry.