Potential Energy Surface for the Benzene Dimer and Perturbational

It is a model system in studies of π···π interactions1 that play key roles in molecular .... It may be noted that choosing a different set of wei...
31 downloads 0 Views 235KB Size
Download PDF
J. Phys. Chem. A 2006, 110, 10345-10354

10345

Potential Energy Surface for the Benzene Dimer and Perturbational Analysis of π-π Interactions Rafał Podeszwa,* Robert Bukowski, and Krzysztof Szalewicz Department of Physics and Astronomy, UniVersity of Delaware, Newark, Delaware 19716 ReceiVed: June 29, 2006; In Final Form: July 10, 2006

We present a complete 6-dimensional potential energy surface for the benzene dimer obtained using symmetryadapted perturbation theory (SAPT) of intermolecular interactions based on Kohn-Sham’s description of monomers. Ab initio calculations were performed for 491 dimer geometries in a triple-ζ-quality basis set supplemented by bond functions. An accurate analytic fit to the ab initio results has been developed and low-energy stationary points on the potential energy surface have been found. We have determined that there are three minima on the surface. Two of them, the tilted T-shape and the parallel-displaced, are nearly isoenergetic with interaction energies of -2.77 and -2.74 kcal/mol, respectively. The third minimum, a twisted edge-to-edge conformation, is significantly less attractive, with the interaction energy equal to -1.82 kcal/mol. Both the T-shape and sandwich geometries, sometimes assumed to be minima, are shown to be only saddle points. The potential energy surface is extremely flat between the two lowest minima, the barrier being only 0.10 kcal/mol above the global minimum. The second-virial coefficient obtained with the new potential agrees well with experimental results over a wide range of temperatures. The SAPT approach rigorously decomposes the interaction energy into physical components. The relative importance of these components has been analyzed.

I. Introduction The benzene dimer is a system of great interest. It is a model system in studies of π‚‚‚π interactions1 that play key roles in molecular recognition,2 protein folding,3 stacking of DNA bases,4 and intercalation of drugs into DNA.5 Although the system has been studied extensively both experimentally6-11 and theoretically,12-29 the features of the benzene dimer potential energy surface (PES) and the nature of the π‚‚‚π interactions are still far from being completely understood. The difficulties arise from the relative weakness of the interaction and the shallowness of PES. Because benzene lacks a permanent dipole moment, the major part of the binding force comes from the dispersion interactions, which require expensive calculations at a correlated level of theory. Furthermore, it has been shown that for the benzene dimer one has to use a high level of theory, such as the coupled cluster method with singles, doubles, and noniterative triples [CCSD(T)], because a low-level approach, such as the second-order perturbation theory based on the Møller-Plesset partitioning of the Hamiltonian (MP2), significantly overestimates the binding energy.16,19,21,23,24 Also, the dispersion energy converges slowly with basis set and therefore large basis sets have to be used. Due to the high costs, the results of previous ab initio calculations have been limited to a few points on the PES. Only one complete modern ab initio potential exists for the benzene dimer.25 However, it was computed using the supermolecular MP2 and a small, double-ζ size basis sets both not adequate for this system. Density functional theory (DFT) methods are much less timeconsuming than the CCSD(T) method; however, the existing versions of DFT, when applied within the supermolecular approach, fail to reproduce the dispersion interaction, an important part of the van der Waals force.30 This problem results from the local nature of the current density functionals, which

cannot describe long-range dispersion interaction. Recently, a few DFT-based methods were developed to overcome this deficiency and applied to the benzene dimer20,27,28,31 but although the results show large improvements over standard DFT functionals, the accuracy is not yet comparable with that of CCSD(T) results. Another approach to calculations of interaction energies was developed by Misquitta et al.32,33 and independently by Hesselmann and Jansen,34 following ideas of Williams and Chabalowski.35 This approach is based on symmetry-adapted perturbation theory (SAPT)36 but utilizes the description of the interacting monomers in terms of Kohn-Sham (KS) orbitals, orbital energies, and frequency-dependent density susceptibility (FDDS) functions. The DFT-based SAPT will be called SAPT(DFT). This method can be shown to be potentially exact for all major components of the interaction energy (asymptotically for exchange interactions) in the sense that these components would be exact if the DFT description of the monomers were exact.32,33,37 Applications to a number of small dimers have shown that SAPT(DFT) provides surprisingly accurate results, sometimes more accurate than the standard SAPT at the currently programmed level.37,38 The method has been recently applied to the benzene dimer for selected configurations39,40 and the results were very close to the much more expensive CCSD(T) results. If density-fitting (DF) (also known as the resolution of identity) technique41 is employed, the DF-SAPT(DFT) method33,39,42,43 becomes relatively inexpensive. In fact, most of the time in a SAPT(DFT) calculation is spent in the DFT calculations for the monomers.39,43 In the present work, because of this high efficiency of the method, we were able perform calculations for 491 dimer geometries with a relatively large triple-ζ-quality basis set. Such calculations would not be possible with the CCSD(T) approach.

10.1021/jp064095o CCC: $33.50 © 2006 American Chemical Society Published on Web 08/08/2006

10346 J. Phys. Chem. A, Vol. 110, No. 34, 2006

Podeszwa et al.

II. Method In the SAPT(DFT) method, the interaction energy (up to the second order) is expressed as38 (1) (2) ESAPT(DFT) ) E(1) int elst(KS) + Eexch(KS) + Eind(CKS) + (2) (2) (CKS) + E(2) ˜ exch-disp (CKS) (1) E ˜ exch-ind disp(CKS) + E

For definitions of the terms, see ref 38. All the components in eq 1 have a clear physical interpretation and correspond to the electrostatic (Coulomb interactions of charge densities of the unperturbed monomers), exchange (effect of Pauli repulsion or, equivalently, of antisymmetrization of the unperturbed wave functions of the monomers), induction (interactions of induced multipole moments with permanent moments of the partner), exchange-induction (effect of antisymmetrization of induction wave functions), dispersion (interaction of instantaneous multipole moments), and exchange-dispersion (effect of antisymmetrization of dispersion wave functions) interactions, respectively. Note that, as in ref 40, we have not used the correction for the terms of order higher than second, which can be extracted from supermolecular Hartree-Fock calculations but is not recommended for nearly nonpolar monomers. For SAPT(DFT) calculations, we used our density-fitting implementation described in detail in ref 43. With density-fitting, the calculations for benzene dimer in a triple-ζ-quality basis set are sped up by more than an order of magnitude with insignificant losses of accuracy compared to the SAPT(DFT) approach without density fitting.39,42,43 The DFT calculations for the monomers were performed using the DALTON44 program. For all calculations, we used the PBE0 DFT functional45,46 with the Fermi-Amaldi asymptotic correction and with the Tozer-Handy splicing scheme.47 With this functional, very accurate results were obtained for small systems.32,33,38,40 To compute the asymptotic correction, we used the experimental ionization potential of 0.3397 a.u.48 Vibrationally averaged monomer geometry of RCC ) 1.3965 Å and RCH ) 1.085 Å has been taken from ref 49. We used the aug-cc-pVDZ and aug-cc-pVTZ basis sets of Dunning et al.50 In all cases, the bases were extended with a set of 36 midbond functions, consisting of three s and three p shells with exponents (0.9, 0.3, 0.1), and of two d and two f shells with exponents (0.6, 0.2). The midbond functions were placed at a position rmb calculated as a weighted average of atom-atom midpoints:51

rmb )

∑∑ a∈A b∈B

wab(ra + rb)

(2)

where a and b index the atoms of monomers A and B, respectively, rx is the position of atom x, and all vectors are relative to the origin of the coordinate system. The weights wab ) rab-6/∑abrab-6, where rab is the distance between the atoms a and b, are motivated by the 1/R6 decay of the dispersion energy, the component requiring the use of bond functions.52 It may be noted that choosing a different set of weights, wab ) mamb/ MAMB, with mx being the mass of atom x and MX being the total mass of monomer X, would correspond to placing the midbond functions midway between the centers of mass of the monomers. Although the latter choice can be successfully employed in calculations for smaller monomers, it is inappropriate when the spatial extent of a monomer is comparable to the distance between the centers of mass for typical dimer configurations. The majority of our calculations were performed with the monomer-centered “plus” basis set (MC+BS) scheme,52 where

orbitals of a given monomer were expanded in terms of this monomer’s own basis, the midbond functions, and the isotropic part of the basis of the other monomer (i.e., the basis without the polarization functions). The density-fitting approximation was accomplished using auxiliary basis sets of ref 53 (corresponding to the main basis set used, i.e., aug-cc-pVTZ) for all SAPT terms except electrostatics, where the cc-pVTZ auxiliary basis set of ref 54 was used and certain calculations required for the electrostatic component, described in ref 43, were performed in quadruple precision. As was the case with the main bases, the auxiliary ones were extended with a set of midbond functions, containing five each of the uncontracted spd shells with exponents (1.8, 1.2, 0.6, 0.4, 0.2), four f shells with exponents (1.5, 0.9, 0.5, 0.3), and three g shells with exponents (1.5, 0.9, 0.3), chosen to approximately reproduce the products of midbond functions. Only the dimer-centered “plus” basis set (DC+BS) types52 were used for auxiliary bases (even if the main basis set was of MC+BS type). III. Fit of Potential Energy Surface The set of 491 dimer geometries was based on a regular grid constructed in the space of 6 dimer coordinates: two Euler angles of monomer A, three Euler angles of monomer B, and the distance R between the centers of mass of the monomers. The grid was based on 128 symmetry-unique angular configurations with several values of R for each configuration, chosen to cover the regions of the potential well, the repulsive wall, and the tail. This regular grid was then extended to 491 geometries by appending energies computed around various characteristic points on a potential surface obtained from preliminary fits. The coordinates of all the points (and their precise definition) and the interaction energies are given in the Supporting Information.55 A Fortran program for converting the internal coordinates into Cartesian can be found in ref 56. The potential surface was fitted to a site-site formula

V)

∑ ∑uab(rab)

(3)

a∈A b∈B

where the summation runs over sites a of monomer A (atomic and off-atomic; see below), sites b of monomer B, and rab denotes the distance between two such sites. The function uab given by 2

uab ) (1 +

m R -β r aab ∑ m rab)e m)1 ab

ab ab

+

qaqb

f1(δab 1 ,rab)

+

rab

∑

n)6,8,10

Cab n

fn(δab n ,rab)

rnab

(4)

may be considered a generalization of the popular Buckinghamtype potential, with an exponential component, the electrostatic component involving charges qa, qb, and the part involving the Cab n coefficients, responsible for modeling of the long-range dispersion and induction interactions, which decay as inverse powers of the distance. To alleviate the divergent character of the latter two terms at short intermolecular distances, these terms are multiplied by the Tang-Toennies damping functions57

fn(δ,r) ) 1 - e-δr

n

(δr)m

m)0

m!

∑

(5)

which are close to 1 for large r but continuously go to zero when r decreases.

Potential Energy Surface for the Benzene Dimer The sites involved in the summation in eq 3 are all the C and H atoms, as well as 13 additional off-atomic sites on each monomer. One set of six off-atomic sites are located on each of the C-H bonds, 0.752 214 Å away from the C atom. Another set of six sites were placed on each bisector between the C atoms, 1.451 29 Å from the geometric center of the molecule. The last off-atomic site is in the geometric center of the molecule. Thus, there were 5 symmetry-distinct sites per monomer. The positions of the off-atomic sites (except for the central one) were roughly optimized while the site charges were fit to molecular multipole moments (vide infra). It should be mentioned that not all types of sites carry all types of interactions implied by formula (4). Specifically, the charge of the central site is set equal to zero, which implies that the damping parameters δab 1 with either a or b corresponding to the central site need not be considered. Likewise, only the C and H atomic sites participate in modeling of the long-range dispersion and induction interactions, so that the only nonvanishing Cab n parameters correspond to ab ) CC, HH, or CH, and similarly for δab n , n ) 6, 8, 10. Overall, the fit is determined by 92 parameters (4 charges qa, 15 Rab parameters and 15 βab parameters, 9 asymptotic coefficients Cab n , 19 damping paramab eters δab n , and 30 polynomial coefficients am ). The process of adjusting the parameters in eq 4 consists of three stages. First, the charges qx are obtained by fitting to the set of multipole moments of the benzene molecule, calculated at the level consistent with the interaction energy calculation, i.e., using the PBE0 Kohn-Sham density generated in the monomer-centered part of the basis set (i.e., the aug-cc-pVTZ basis). All multipoles through l ) 6 were used in the fit, with extra weight put on the total charge (i.e., the multipole l ) 0) to keep the molecule electrically neutral. The positions of the noncentral off-atomic sites were also optimized at this stage. In the second stage, the asymptotic parameters Cab n were obtained from a linear least-squares fit to the sum of the asymptotic dispersion and induction energies, computed ab initio on a set of dimer geometries obtained from the original set of 491 points by adding 3 Å to the R coordinate of each geometry. To compute the ab initio dispersion and induction energies in the asymptotic region, we used the van der Waals asymptotic constants through R-12 obtained from the Kohn-Sham multipole moments and the coupled Kohn-Sham (CKS) static and dynamic polarizabilities, both generated in the aug-cc-pVTZ basis set using the PBE0 functional. The calculation of the asymptotic constants was performed with the help of programs DISPER and INDUCT from the POLCOR suite by Wormer and Hettema.58 The asymptotic parameters are given in the Supporting Information.55 The computed C6 coefficient equal to -1726 a.u. is in an excellent agreement with the experimental value of -1723 a.u.59 In the third and final stage of the fitting process, all the remaining parameters in eq 4 were adjusted by least-squares fitting to the SAPT(DFT) results at 491 geometries. During this nonlinear fit, the site charges and the Cab n coefficients, obtained independently in the previous stages, were kept constant. The space of nonlinear parameters Rab, βab, and δab n , n ) 1, 6, 8, 10, was explored using a variant of the Powell minimization routine. For each set of nonlinear parameter values visited by this routine, the coefficients aab m , m ) 1, 2 were obtained by means of a call to a linear least-squares fit. The weight assigned to configuration i was dependent on the interaction energy Ei at this configuration and equal to exp((E0 - Ei)/(kcal/mol)) if Ei < E0 and (E0/Ei)2 otherwise, with the parameter E0 chosen as 3 kcal/mol. Thus, the low-energy regions of the potential surface were weighted

J. Phys. Chem. A, Vol. 110, No. 34, 2006 10347

Figure 1. Minimum structures of the benzene dimer.

Figure 2. Saddle point structures S1-S3 of the benzene dimer.

much stronger than the repulsive regions. As a result of such a scheme, the unweighted RMSE of the fit calculated for important points with Ei < 0 was 0.02 kcal/mol, much smaller than the overall RMSE of 0.15 kcal/mol. The parameters of the fit are given in the Supporting Information.55 A Fortran program implementing the fit can be found in ref 56. Our choice of parameters and the fitting procedure ensures that the fit behaves correctly at large intermolecular distances. The second term on the rhs of eq 4 [or, more strictly, the sum of such terms, as implied by eq 3] represents then the electrostatic component of the interaction energy, and the third term, the sum of induction and dispersion components. For smaller R, starting from about twice the radial minimum distance, there is no simple correspondence between the individual terms of the fitting function and the interaction energy components. For such distances, only the sum of all the terms is meaningful. IV. Characteristic Points on the Potential Energy Surface The fitted potential energy surface of the dimer was explored with a simple implementation of the eigenvector-following local optimization method60 using randomly selected configurations as starting points. This procedure resulted in three minima, six saddle points of index 1, and two structures of index 2 and 3, shown in Figures 1-4. The list of higher-index stationary points given here is far from completesit only includes characteristic “sandwich” structures often studied in the literature. Geometrical parameters of all stationary points and the corresponding interaction energies obtained directly from the fit and from ab initio SAPT(DFT) calculations at the optimized geometries have been collected in Table 1 and compared to high-level ab initio data from literature or computed in this work. The results obtained from our fit at literature geometries are also given.

10348 J. Phys. Chem. A, Vol. 110, No. 34, 2006

Podeszwa et al. TABLE 1: Stationary Points on the Fitted SAPT(DFT) Potential Energy Surfacea structure symmetry

Figure 3. Saddle point structures S4-S6 of the benzene dimer.

Figure 4. Saddle point structures S7 and S8 of the benzene dimer.

Furthermore, we have calculated MP2 and CCSD(T) energies using our optimized geometries. In Table 1, we present CCSD(T) results obtained from a sum of the MP2 aug-cc-pVTZ interaction energy and ∆CCSD(T) ) ECCSD(T) - EMP2 int int obtained using the aug-cc-pVDZ basis set. Both basis sets where supplemented with midbond functions described in section II. The MP2 and CCSD(T) results were obtained in the frozen-core approximation with the standard Boys-Bernardi counterpoise correction using the MOLPRO suite of programs.61 It should be noted that due to the size of the benzene dimer, exhaustive explorations of the potential surface with ab initio methods are not feasible and only limited-dimensionality optimizations have been performed in the literature using supermolecular DFT28 and MP215,23 approaches, as well as methods that include higher-order electron correlation effects, up to the CCSD(T) level.18,21,23,24 In view of the fact that electron correlation beyond MP2 is crucial for a proper description of energetics of the benzene dimer, it is somewhat surprising that optimal structures derived from low-level methods are generally in quite good agreement with higher-level calculations. However, apart from 1- or 2-dimensional cross-sections for a handful of characteristic dimer configurations, the shape of the potential surface has not been known. The present work provides the first high-level ab initio potential that allows thorough exploration of the interaction landscape and full characterization of stationary points. The global minimum on the fitted surface is the paralleldisplaced configuration M1, nearly isoenergetic with the tilted T-shape minimum M2. In fact, the energetic sequence of these minima is reversed if the actual ab initio energies are considered instead of the fit results. The M2 structure has not been reported in the literature; however, a more tilted structure (slipped edge to plane) has been found as the global minimum on the empirical (fitted to the properties of bulk benzene) dimer potential surface of Jorgensen and Severance.62 The latter structure has been found to be almost isoenergetic with the T-shape structure in

M1 M1c M1d M1e M2 M3 S1 S2 S3 S3c S3d S3e S4 S4f S5 S6 S7 S8 S8c S8d S8e

C2h C2h C2h C2h Cs C2V Cs Cs C2V C2V C2V C2V Cs Cs C2V Cs C6V D6h D6h D6h D6h

R

θ1

θ2

Efit int

3.962 3.9 3.9 3.8 4.954 6.104 4.960 3.973 4.982 5.0 5.0 4.894 4.453 4.49 5.026 5.908 3.793 3.807 3.8 3.9 3.7

60.96 63 66 65 99.63

60.96 63 66 65 11.71

-2.764 -2.688 -2.625 -2.628 -2.763 -1.815 -2.745 -2.734 -2.689 -2.688 -2.688 -2.689 -2.642 -2.572 -2.427 -1.763 -1.868 -1.835 -1.835 -1.802 -1.777

98.59 10.90 63.16 58.88

66.89 35.97 69 39 29.83 19.09

Ecalc int

b ECCSD(T) int

-2.742 -2.699 -2.48 -2.63 -2.78 -2.771 -2.795 -1.816 -1.805 -2.754 -2.782 -2.725 -2.688 -2.662 -2.683 -2.46 -2.61 -2.74 -2.675 -2.611 -2.39 -2.422 -2.376 -1.760 -1.771 -1.857 -1.679 -1.850 -1.682 -1.48 -1.70 -1.81

a R, θ , θ are the geometrical parameters (in Å and deg) shown in 1 2 Figures 1-4. Energies are in kcal/mol. Efit int is the energy obtained from our fit at a given geometry. Ecalc int is the calculated SAPT(DFT) energy. ECCSD(T) is the CCSD(T) energy from our calculations or a int literature CCSD(T) value. Structures M1-M3 are minima, S1-S6 saddle points of index 1, and S7 and S8 saddle points of index 2 and 3, respectively b This work, except where indicated. Computed as MP2/ aug-cc-pVTZ interaction energy plus ∆CCSD(T) ) ECCSD(T) - EMP2 int int obtained using aug-cc-pVDZ. Both basis sets were supplemented with midbond functions described in text. c Tsuzuki et al.21 estimated CCSD(T) results calculated using models AIMI III (energies) and AIMI II (1- or 2-dimensional geometry optimizations). d Sinnokrot and Sherrill24 CCSD(T) results computed as MP2/aug-cc-pVQZ results plus ∆CCSD(T)/aug-cc-pVDZ*, where the asterisk denotes that some functions were removed from the original basis sets. e Sinnokrot et al.23 estimated complete basis set (CBS) limit CCSD(T) energies computed at geometries optimized at the MP2/aug-cc-pVTZ level. f Tsuzuki et al.22 geometry and energy at the CCSD(T) level from the AIMI model II. The geometry of ref 22 is slightly different from our geometry: the bottom monomer of S4 in Figure 3 should be rotated by 30° around its 6-fold axis.

early studies of Hobza et al.13 The geometric parameters of M1 agree quite well with the results of the crude optimizations at the CCSD(T)21,24 and MP223 levels, although the intermonomer distance from the SAPT(DFT) potential is slightly larger and the monomers are tilted somewhat less with respect to the intermolecular axis. The present (fit) result at M1 is 0.15 kcal/ mol above the estimated CBS-limit value of CCSD(T) interaction energy of ref 23, which contains the MP2 component calculated using the explicitly correlated MP2-R12 approach and ∆CCSD(T) computed in aug-cc-pVDZ basis set. Note that this comparison is made at a dimer geometry slightly different from our M1 geometry (and at slightly different monomer geometries), obtained from an MP2/aug-cc-pVTZ optimization in ref 23. The 0.15 kcal/mol difference is consistent with our estimates of the accuracy of SAPT(DFT) results reported below. The SAPT(DFT) interaction energies for the structures near M1 (and for most other structures) are very close to the CCSD(T) results of ref 24. The Tsuzuki et al. CCSD(T) calculations from refs 21 and 22 give significantly smaller binding energies. The structures M1 and M2 are separated by the saddle point S4 with the energy of only 0.10 kcal/mol above M2. A structure similar to S4 (differing from S4 by the 30° rotation of the bottom molecule around its 6-fold axis) was studied by Tsuzuki et al. in ref 22 at the CCSD(T) level, and a similar conclusion was

Potential Energy Surface for the Benzene Dimer reached about the height of the interconversion barrier between the M1 and T-shaped configurations (the latter approximately representing the M2 minimum). Such a low barrier suggests that the rovibrational ground-state wave function has a large amplitude in both minima. This effect may be responsible for ambiguities in the experimental assessment of the dimer geometry.6-10 Although the microwave spectra of the benzene dimer were interpreted by Arunan and Gutowsky8 as indicating the T-shape minimum structure, the possibility that other structures, including the parallel-displaced structure, are energetically close was not ruled out (see also a recent discussion in ref 63). In view of our results, the experimental findings should perhaps be reinterpreted to see if they can be reconciled with the existence of two nearly isoenergetic minima structures. Let us note that in the benzene crystal structure64 molecular pairs similar to both M1 and M2 structures are present,65 despite the fact that the cooperative two-body and pairwise nonadditive forces may in general lead to crystal structures very different from those of the dimer. Of course, the high symmetry of the benzene molecule implies the existence of multiple equivalent minima of the form M1 and M2. Tunneling between the equivalent M1 minima occurs via the saddle point S2, which is only 0.05 kcal/mol above that for M1. Therefore, this tunneling may be thought of as an essentially free relative in-plane rotation of the two molecules in the M1 configuration. Saddle points S1 and the T-shaped S3, with ab initio calculated barrier heights 0.02 and 0.11 kcal/ mol above the M2 energy, respectively, separate different equivalent M2 minima. An additional M2-M2 transition mechanism, corresponding to the S5 saddle point, is also possible, with a somewhat higher barrier of 0.35 kcal/mol. The geometries and barrier heights of the saddle points S1, S3, and S5 suggest very high mobility of the monomer, constituting the “stem”′ of the “T” in the near-T-shaped M2 configuration. In addition to the almost free rotation around the intermolecular axis, this monomer performs large-amplitude librations around its own 6-fold axis. Besides the M1 and M2 minima, we have also found another minimum structure M3, which was suggested in early MP2 studies of Hobza et al.13 The M3 minimum is about 1 kcal/mol less attractive than the other two minima, mostly due to large interatomic distances, diminishing the dispersion interactions. The configuration M3 can be easily converted into M2 via the saddle point S6, only 0.06 kcal/mol above M3. Other characteristic dimer configurations are the “sandwich” configurations S7 and S8, which turn out to be stationary points of indexes 2 and 3, respectively. Although the interatomic distances and the geometry of both these structures are somewhat reminiscent of M1, energetically S7 and S8 are about 1 kcal/mol above M1 and very close to the M3 minimum. The main reason for this difference is the electrostatic interaction, repulsive in the S7 and S8 structures (except for small R), but attractive in the M1 configuration. To determine which structures should be observed, the effects of vibrational zero-point energy (ZPE) should be taken into account. Using our PES, we have found the harmonic ZPEs equal to 0.319, 0.315, and 0.313 kcal/mol for M1, M2, and M3 structures, respectively. The harmonic frequencies of the structures are given in the Supporting Information.55 The resulting dissociation energies equal to 2.42 and 2.46 kcal/mol for M1 and M2 structures, respectively, are in excellent agreement with the experimental energy of 2.4 ( 0.4 kcal/mol of Grover et al.,7 but significantly above the 1.8 ( 0.2 kcal/ mol result of Krause et al.11

J. Phys. Chem. A, Vol. 110, No. 34, 2006 10349 TABLE 2: SAPT(DFT) Results Compared to the CCSD(T) and MP2 Results for the M1, M2 and S8 Structures (Dimer Geometries Optimized in the Present Work)a SCF/aVDZ+b SCF/aVTZ+b SCF/aVQZ+b MP2/aVDZ+b MP2/aVTZ+b MP2/aVQZ+b CCSD(T)/aVDZ+b CCSD(T)/aVTZ+b ∆CCSD(T)/aVDZ+b ∆CCSD(T)/aVTZ+b CCSD(T)b CCSD(T) best estimate (see text) SAPT(DFT)/aVDZ+b SAPT(DFT)/aVTZ+b

M1

M2

3.753 3.760

1.183 1.216

-4.298 -4.467

-3.492 -3.599

-2.527

-2.689

1.769

0.803

S8

-2.699 -2.7 ( 0.2

-2.795 -2.8 ( 0.1

4.271 4.223 4.215 -3.104 -3.252 -3.307 -1.534 -1.583 1.570 1.669 -1.682 -1.63 ( 0.10

-2.585 -2.742

-2.701 -2.771

-1.662 -1.850

a ∆CCSD(T) ) ECCSD(T) - EMP2. The acronyms aVXZ+b stand for int int the aug-cc-pVXZ basis set plus midbond functions described in the text. SAPT(DFT) results calculated in the MC+BS scheme. All calculations performed in the present work. b MP2/aVTZ+b plus ∆CCSD(T)/aVDZ+b.

V. Estimates of Accuracy When applying our PES and comparing with experimental results, one would like to know how accurate is the surface. Due to the size of the system, only very rough estimates of accuracy can be made. One has to take into account the uncertainties due to the truncation of the level of theory and due to the use of the finite basis set. The former factor can only be estimated by examining several levels of theory or by analyzing results from different theories. We will compare SAPT(DFT) results to those from the CCSD(T) method. The latter method is known to be, in general, the most accurate of practically applicable ab initio approaches. In many cases, one can evaluate the differences between two levels of theory simply by computing interaction energies in the same basis set. This would not work in comparisons of SAPT(DFT) with CCSD(T) because the two methods show distinctly different patterns of convergence in basis set size,37,38 with SAPT(DFT) converging faster. Thus, it is the best to compare the two theories at the CBS limit. The perusal of the SAPT(DFT) results presented in Table 2 shows that the differences between the interaction energies computed in bases aVDZ+b and aVTZ+b (where aVXZ stands for aug-cc-pVXZ and “+b” denotes midbond functions) range between 0.07 and 0.19 kcal/mol. We have also compared the basis set convergence of the MC+BS scheme used by us with the DC+BS one. The use of MC+BS, which contains 576 functions at the aVTZ level, leads to more than 3-fold savings of the computer time compared to the use of DC+BS, which gives 864 functions at the same level. These savings allowed us to compute the required number of points on the PES. To check the accuracy of the MC+BS scheme, we computed the DC+BS results for the stationary points. We found that the use of MC+BS leads to an underestimation of the magnitude of the interaction energies from 0.05 kcal/mol for the M2 structure, through 0.09 for S8, to 0.11 for M1, relatively small errors compared the computer time savings. We have also performed a 1/X3 extrapolation of the dispersion energy to the CBS limit using aVDZ+b and aVTZ+b DC+BS results (the remaining components are practically converged in aVTZ+b). The extrapolation deepens the results further by 0.05 kcal/mol for M2, 0.05 for S8, and 0.10 kcal/mol for M1. On the basis of the extrapolations and on the DC+BS differences with respect to

10350 J. Phys. Chem. A, Vol. 110, No. 34, 2006 the MC+BS results, we can assume that the basis set incompleteness error remaining in the basis aVTZ+b (MC+BS form) is, depending on configuration, about 0.1-0.2 kcal/mol in the well region. To investigate the convergence of the CCSD(T) method with respect to the basis set size, we performed the CCSD(T) calculations using several aVXZ+b bases. The results for the M1, M2, and S8 structures are presented in Table 2. The interaction energies for the remaining stationary points are included in the Supporting Information.55 The S8 sandwich structure has the highest symmetry, and therefore, we were able to perform in this case calculations in larger basis sets than for the other structures. In particular, the MP2 energies were computed in the aVQZ+b basis and CCSD(T) in aVTZ+b. These basis sets are larger than in any published work on the benzene dimer. The 1/X3 extrapolation was applied to the MP2 correlation contribution to the interaction energy, and the result was added to the SCF energy in the larger basis giving -3.29 and -3.34 kcal/mol at D-T and T-Q levels, respectively. An analogous extrapolation from D-T bases gives -1.58 kcal/mol at the CCSD(T) level. Adding to the latter result the difference between D-T and T-Q extrapolations for MP2, one obtains the CBS estimate of the S8 interaction energy as -1.63 ( 0.05 kcal/mol, where the error bar is equal to the difference used. For the remaining two configurations, we can perform only D-T extrapolations at the MP2 level, which give -4.54 and -3.66 kcal/mol for structures M1 and M2, respectively. If we were to add the values of ∆CCSD(T) computed in the basis aVDZ+b to these CBS limits, this would result in a significant overestimate of the CCSD(T) interaction energy, as seen for the S8 structure where this procedure gives -1.74 kcal/mol, 0.11 kcal/ mol from the CBS result obtained with larger basis sets. For S8, a better result, -1.68 kcal/mol, is obtained by taking the computed aVTZ+b result for MP2. Using this method, one obtains the best estimates of CCSD(T) interaction energy for the configurations M1 and M2 equal to -2.7 ( 0.1 and -2.8 ( 0.1 kcal/mol, respectively, with the error bars doubled compared to the S8 case. Our best estimates of the CCSD(T) interaction energies are in reasonable agreement with those obtained by Sinnokrot et al. in refs 23 and 24 (see the last two rows of Table 2 in ref 24). Part of the differences is due to the slightly different monomer and dimer geometries in the two sets of results. We have computed the MP2 interaction energy for the S8 configuration using the same geometry as in Table 2 of ref 24 and the aVTZ+b basis. This calculation gave the interaction energy 0.07 kcal/mol below the result in our Table 2. Thus, in our geometry, the two estimates by Sinnokrot et al.23,24 would be -1.63 and -1.73 kcal/mol. The former result is in excellent agreement with our estimate. The latter result is based on the MP2-R12 interaction energy, computed using an explicitly correlated basis set. In our S8 geometry, this interaction energy would amount to -3.57 kcal/mol, 0.23 kcal/mol below our T-Q extrapolated value. Similarly, there are large differences between orbital and MP2-R12 interaction energies for M1 configuration, but not for S3. We believe that the limit obtained by the extrapolation is more reliable than the MP2-R12 result. However, to partly accommodate the latter result, we have doubled the error bars of our limits for configurations S8 and M1, and these final estimates are listed in Table 2. The estimated limits of the CCSD(T) interaction energies allow us to judge the accuracy of the SAPT(DFT) interaction energies. Of course, the error of the CCSD(T) method with respect to the exact interaction energy of the benzene dimer is

Podeszwa et al. completely unknown. In fact, this error can be estimated only for a handful of systems. For example, it is about 1% near the minimum of the argon dimer.66 For the helium dimer, this error is known exactly and amounts to 3% (ref 67). Thus, if the relative errors were 3% also for the benzene dimer, for the S8 configuration the absolute error would amount to 0.05 kcal/ mol, smaller than the basis set truncation errors discussed above. On the other hand, some models of the benzene molecule indicate that effects beyond CCSD(T) might be fairly large for this system,68 which may also be reflected in larger uncertainties of CCSD(T) interaction energies. At the present time, we can only assume that these uncertainties will not change the error estimates given in Table 2 in a significant way. Comparing the SAPT(DFT) interaction energies in the basis aVTZ+b listed in Table 2 to the estimated CCSD(T) CBS limits, we see that these energies are virtually identical for configurations M1 and M2. For S8, the SAPT(DFT) result is 0.2 kcal/mol below the CCSD(T) limit. However, for this configuration the difference between the MP2 and CCSD(T) results is particularly large; therefore, the CCSD(T) theory truncation error may be relatively large there. In any case, we will assume that our potential is accurate to about 0.2 kcal/mol in the well region. The results of Table 1 show that the conclusions reached above based on configurations M1, M2, and S8 should extend to other geometries of the dimer. For the structures M3 and S1-S6, the agreement of SAPT(DFT) with the CCSD(T) results [computed as MP2/aVTZ+b + ∆CCSD(T)/aVDZ+b] is to within 0.06 kcal/mol. Only for the S7 sandwich structure is the difference larger, 0.18 kcal/mol, just like for the very similar S8 structure. Let us note that except for the S7 and S8 structures (which are almost isoenergetic in both approaches), SAPT(DFT) gives the same energy orderings as CCSD(T). This includes the ordering of the S3 structure (T-shape) and M2 (tilted T-shape). Even though the energetic difference between the two structures is within the error estimate of our PES, this ordering should not change at the complete basis set limit because the two structures are very similar and therefore the errors should be very close. The results of Table 1 show that the 0.2 kcal/mol estimate of the error of the SAPT(DFT) PES, based mainly on the results for the structure S8, may represent an upper bound for this error. Thus, structure S8 seems to represent the most difficult case for the convergence of both SAPT(DFT) and the supermolecular approach. A further confirmation of our estimates of the accuracy of our PES comes from comparisons of the computed and experimental virial coefficients; see section VII. Our results allow us to make some comment on the often used practice of adding the ∆CCSD(T) contribution computed in a small-size basis to the MP2 interaction energy computed in a much larger basis or possibly extrapolated to the CBS limit. Whereas this method usually works very well, for the benzene dimer it may result in some overestimation of the CCSD(T) interaction energy because the fairly large ∆CCSD(T) contribution is quite sensitive to the basis set size. This energy difference tends to increase with the increase of the basis set, and for the two basis sets considered in the Table 2 the increase is 0.10 kcal/mol, similar to the 0.15 kcal/mol change at the MP2 level. The CCSD(T) results obtained by adding ∆CCSD(T) computed in a small basis to MP2 computed in a moderate-size basis set, such as aVTZ+b, may actually give a more accurate result than the former approach. This is the case for the structure S8 discussed above. One should also note that due to the use of the midbond functions, our aug-cc-pVDZ and aug-cc-pVTZ results are more basis-set saturated than the literature data without midbond functions.

Potential Energy Surface for the Benzene Dimer

Figure 5. Cross section of the angular configuration of structure M1. “elst” ) electrostatic energy, “exch” ) the sum of first and second-order exchange energies, “ind” ) induction energy, “disp” ) dispersion energy.

Figure 6. Cross section of the angular configuration of structure S3. Labels as in Figure 5.

Figure 7. Cross section of the angular configuration of structure S8. Labels as in Figure 5.

VI. Potential Cross-sections and Interaction Energy Components Symmetry-adapted perturbation theory provides insight into the physical contributions to the interaction energy. In Figures 5-7, we present the interaction energy components and total energies for radial cross-sections through the potential energy surface for angular coordinates corresponding to the structures

J. Phys. Chem. A, Vol. 110, No. 34, 2006 10351 most extensively studied in the literature: M1, S3, and S8. Detailed numerical data for all the 491 grid points are presented in the Supporting Information.55 Both the M1 and S3 structures are electrostatically attractive near the minimum. At large distances, the electrostatic interaction remains attractive for S3 but becomes repulsive for M1, in the region where it becomes dominated by the quadrupolequadrupole interaction. It is interesting to note that the often used argument that the slipped-parallel configuration is the minimum due to the quadrupole-quadrupole interaction is actually not true. This interaction has a minimum for θ1 ) θ2 ≈ 41°, where the angles are defined in Figure 1. At θ1 ) θ2 ) 61°, the value for the M1 minimum, the quadrupole-quadrupole interaction is slightly positive because, as a function of θ1 ) θ2, it changes sign at θ1 ) θ2 ≈ 59.5°. The reasons the electrostatic energy of M1 becomes negative for smaller R are the negative hexadecapole-hexadecapole interaction and overlap effects. Asymptotically (not shown in the figure), the total interaction energy becomes positive for M1 as the O(R-5) decay of the quadrupole-quadrupole interaction is slower than the O(R-6) decay of the dispersion interaction. Compared to S3, in the M1 configuration all the contributions except for the electrostatic one are larger in magnitude at the minimum due to the shorter (on the average) interatomic distances. However, the larger magnitude of dispersion and induction energies are canceled by the larger exchange energy of M1 and the resulting total interaction energy in the minimum is almost equal to that in the S3 configuration. Because the induction component is almost completely quenched by its exchange counterpart, whereas the dispersion one is not (see below), it is in fact the dispersion interaction that stabilizes M1. Unlike M1 and S3, the S8 structure is electrostatically repulsive at a wide range of distances. The electrostatic energy becomes negative only at the short range due to the overlap effects, but at these distances the rapidly increasing exchange repulsion dominates the overall picture. Although the S8 structure is not as strongly bound as M1, its repulsive region rises more gradually and, at the shortest distance shown on the picture (3.25 Å), S8 is significantly less repulsive than M1. The stronger M1 repulsion is due to the larger overlap of electron densities because at the same distance between the centers of the mass of the monomers, the distance between the ring planes is shorter for M1. Comparing M2 and S3 structures, one sees that the tilted structure of M2 is a result of more favorable dispersion interactions. Although the electrostatic energy of M2 (-1.639 kcal/mol) is slightly less negative than the S3 one (-1.660 kcal/ mol), the dispersion energy is much more favorable (-4.773 vs -4.581 kcal/mol) with only slightly less favorable first-order exchange (3.448 vs 3.413 kcal/mol).55 It should be noted that the second-order exchange-induction energy quenches to a large extent the induction contribution and the sum of these two energies plays a minor role in the total energy of the dimer, as observed also in ref 24. All the investigated structures have a very large dispersion component (which is only slightly quenched by its exchange counterpart), and for all distances except the extreme long range, dispersion is the most significant attractive contribution. For the directionality of the interaction, it is the electrostatic component that depends mostly on the angular configuration, although the dispersion component plays also some role in the directionality. One example is the tilt of the M2 structure described above. Our current SAPT(DFT) results are generally similar to the SAPT2 results computed by Sinnokrot and Sherril24 for the

10352 J. Phys. Chem. A, Vol. 110, No. 34, 2006

Podeszwa et al.

sandwich, slipped-parallel, and T-shape structures, despite the fact that SAPT2 includes electron correlation at a level much lower than does SAPT(DFT). In particular, SAPT2 is asymptotically equivalent to the MP2 method, which for the benzene dimer gives results that differ significantly from the SAPT(DFT) or CCSD(T) results. The reason for the overall good performance of the SAPT2 results is cancellation of errors of the method with errors due to the use of a small basis set, as already noted in ref 24. There are a few significant discrepancies between the results from the two methods, though. For example, the parallel-displaced structure of Sinnokrot and Sherill is more electrostatically attractive than the T-shape structure, in contrast to our results. The simple model of Tsuzuki et al.,21 where the electrostatic component is estimated from monomer properties and dispersion is estimated from the difference between the CCSD(T) and SCF interaction energies, provides a physical picture of the interaction that only very roughly agrees with the accurate SAPT(DFT) decomposition. There are several important deficiencies, for example, the electrostatic interaction for the parallel-displaced structure is repulsive in ref 21, in contrast to our results. VII. Second Virial Coefficient The fitted potential energy surface has been used in calculations of the second virial coefficient, B(T), of benzene, defined by the virial expansion

B(T) pV )1+ + ... RT V

(6)

where p, V, R, and T are the pressure, the molar volume, the gas constant, and the temperature, respectively. The formula employed is given by eqs 22-24 of ref 69 and contains the lowest-order quantum correction (of the order of p2). It should be noted that due to a relatively large mass and large moments of inertia of the benzene molecule, the quantum contribution to B(T) is generally small and does not exceed 0.6% even for the lowest temperatures considered here. The calculation of B(T) involves a 6-dimensional integration over the relative coordinates of the dimer. The integral over the angular coordinates was computed using a simple Monte Carlo (MC) procedure sampling with 360 000 configurations. For each of these configurations, the radial integration was performed over the range from 0 to 70 Å by splitting this range into 15 intervals and then using a 31-point Stenger quadrature70 in each of these intervals. The radial integral computed in this way is practically exact compared to the uncertainty of the angular integration. The latter can be assessed from the MC procedure and amounts to about 0.1% for all temperatures considered. The second virial coefficients computed with the SAPT(DFT) potential over a range of temperatures are shown in Figure 8 and compared to two sets of experimental data. Also shown are the results obtained with the empirical OPLS potential of Jorgensen and Severance.62 It is seen that the agreement of the SAPT(DFT) results with experiment is generally very good over the whole range of temperatures, the discrepancies being typically between 2 and 6%. Both sets of experimental values are slightly above the theoretical ones, which may indicate that the SAPT(DFT) potential is a little too attractive on the average. The theoretical curve remains too low in the high-temperature region, which may indicate that the repulsive wall of the SAPT(DFT) potential is somewhat too soft. The better performance at low temperatures suggests higher accuracy of the minima wells. Let us note that although the B(T) values obtained by Cacelli et al.25 are also accurate, its agreement with experiment is the

Figure 8. Second virial coefficient of benzene. SAPT(DFT): this work. Experiment 1: ref 71. Experiment 2: ref 72. OPLS: calculated with the empirical potential of ref 62.

result of a fortuitous cancellation of errors because both the method and the basis set used are not adequate for benzene dimer. Moreover, their fit function, of the OPLS form, does not reproduce their ab initio data, as seen in Figure 2 of ref 25 for the bottom parts of all the wells and for the whole curve in the slipped-parallel configuration. The virial coefficient data obtained from the empirical OPLS potential62 show that this potential strongly underestimates the attractive interactions between benzene molecules, which results in values of B(T) significantly less negative than the experimental ones. The reason for this behavior is, of course, the OPLS potential has been fitted to represent bulk-phase properties of benzene, and as such, it must implicitly account for many-body interactions. Because empirical pair potentials account for such effects by a distortion of pair interaction energies, this typically leads to poor values of purely pairwise properties, such as B(T). The underbinding of OPLS is in contrast to the situation for polar liquids, like water, for which empirical pair potentials tend to be more attractive than the true pair potentials. We now know73,74 that the reason for this overbinding is that many-body effects in polar systems are dominated by the induction interactions, which increase the magnitude of the interaction energy. For benzene, the many-body forces are probably dominated by the third-order dispersion nonadditive effects and possibly also by the exchange effects, like in the argon trimer.75,76 The OPLS results indicate that the overall effect of these interactions is repulsive. Compared to the water case, where empirical pair potentials give B(T) that are often a factor of 2 or more larger in magnitude than the experimental values,69 the values of B(T) given by the empirical potential for benzene agree reasonably well with the experimental and ab initio ones. This suggests that the many-body effects in benzene are not as strong as in water, where the inclusion of three-body forces is absolutely necessary for the description of liquid. Recent results computed by Tauer and Sherill77 for a few benzene trimer configurations also show that the three-body effects are, in general, small for this system. Although accurate simulations of the condensed phases of benzene would require inclusion of the many-body effects, the two-body ab initio potential alone should lead to much smaller errors than in the case of water. VIII. Conclusions With the low cost of the SAPT(DFT) calculations, we were able to compute several hundred configurations of the benzene dimer. The accuracy of the calculations was estimated to be

Potential Energy Surface for the Benzene Dimer about 0.2 kcal/mol in the potential well region. The large number of points allowed us to obtain an accurate fit of the PES. With the fitted potential, it has been possible to determine all important stationary points, extending the information obtained in previous studies. We have found two almost isoenergetic global minima, the slipped parallel and the tilted T-shape ones. The latter structure, not reported before, achieves the stability due to a larger dispersion component compared to the C2V T-shape structure, which was found by us to be a saddle point. The relative stability of the structures has been confirmed by separate CCSD(T) calculations. The surface is very flat around these minima and the barriers between them are only about 0.1 kcal/mol. Another minimum of a higher energy has also been found. The complete characterization of the stationary points resolves the ongoing debate whether the slipped parallel, T-shape, and sandwich configurations are minima or saddle points.13,15,16,18,22,23 Unlike MP2, SAPT(DFT) energy results are very close to the CCSD(T) values. The good agreement of the computed second virial coefficient with experiment provides an independent verification of the quality of the potential. Our decomposition of the total energy confirms that the stabilization energy of benzene π-π interactions is dominated by the dispersion forces, although the exchange effects are of similar size in the minima. The directionality of the interaction is mainly due to the electrostatic forces with some effect from the balance between the dispersion and exchange terms. The induction component (strongly quenched by the exchange-induction term) plays a minor role in both stability and directionality. Because the standard implementations of the supermolecular DFT approach are unable to recover the dispersion component, numerous DFT studies of π-π interactions for larger systems cannot be expected to give reliable results. We believe the results for the benzene dimer shed some more light on understanding of the strength and directionality of the π-π interactions and may help to improve the effective potentials for larger π-π systems. Our potential provides accurate and complete information about the interaction energy of two benzene molecules and the correct asymptotic decomposition of the interaction energy components (the components of the fit cannot be related to the physical components in the region of the potential energy well, though). In particular, the analysis of the surface in the regions around the minima and along the tunneling paths may help to interpret the experimental results for this system. Acknowledgment. Funding for this work was provided by an ARO DEPSCOR grant and by NSF grants CHE0239611 and CHE0555979. Supporting Information Available: (a) The energy components, total energy, and geometric parameters of all points calculated in this study, (b) the MP2 and CCSD(T) results for the characteristic points, (c) the asymptotic parameters, (d) the parameters of the fit, (e) harmonic frequencies of the minima. This material is available free of charge via Internet at http:// pubs.acs.org. References and Notes (1) Kim, K. S.; Tarakeshwar, P.; Lee, J. Y. Chem. ReV. 2000, 100, 4145-4185. (2) Meyer, E. A.; Castellano, R. K.; Diederich, F. Angew. Chem., Int. Ed. 2003, 42, 1210-1250. (3) McGaughey, G. B.; Gagne´, M.; Rappe´, A. K. J. Biol. Chem. 1998, 273, 15458-15463. (4) Hobza, P.; Sponer, J. J. Am. Chem. Soc. 2002, 124, 11802-11808.

J. Phys. Chem. A, Vol. 110, No. 34, 2006 10353 (5) Remeta, D. P.; Mudd, C. P.; Berger, R. L.; Breslauer, K. J. Biochemistry 1993, 32, 5064-5073. (6) Janda, K. C.; Hemminger, J. C.; Winn, J. S.; Novic, S. E.; Harris, S. J.; Klemperer, W. J. Chem. Phys. 1975, 63, 1419-1421. (7) Grover, J. R.; Waiters, E. A.; Huit, E. T. J. Phys. Chem. 1987, 91, 3233-3237. (8) Arunan, E.; Gutowsky, H. S. J. J. Chem. Phys. 1993, 98, 42944296. (9) Venturo, V.; Felker, P. J. Phys. Chem. 1993, 97, 4882-4886. (10) Schaeffer, M.; Maxton, P. M.; Felker, P. M. Chem. Phys. Lett. 1994, 224, 544-550. (11) Krause, H.; Ernstberger, B.; Neusser, J. Chem. Phys. Lett. 1991, 184, 411. (12) Hobza, P.; Selzle, H. L.; Schlag, E. W. J. Chem. Phys 1990, 93, 5893-5897. (13) Hobza, P.; Selzle, H. L.; Schlag, E. W. J. Am. Chem. Soc. 1994, 116, 3500-3506. (14) Tsuzuki, S.; Uchimaru, T.; Tanabe, K.; Kuwajima, S. J. Phys. Chem. 1998, 98, 1830-1833. (15) Jaffe, R. L.; Smith, G. D. J. Chem. Phys. 1996, 105, 2780-2788. (16) Hobza, P.; Selzle, H. L.; Schlag, E. W. J. Phys. Chem. 1996, 100, 18790-18794. (17) Hobza, P.; Spirko, V.; Selzle, H. L.; Schlag, E. W. J. Phys. Chem. A 1998, 102, 2501-2504. (18) Sˇ pirko, V.; Engkvist, O.; Solda´n, P.; Selzle, H. L.; Schlag, E. W.; Hobza, P. J. Chem. Phys. 1999, 111, 572-582. (19) Tsuzuki, S.; Uchimaru, T.; Matsumura, K.; Mikami, M.; Tanabe, K. Chem. Phys. Lett. 2000, 319, 547-554. (20) Tran, F.; Weber, J.; Wesołowski, T. A. HelV. Chim. Acta 2001, 84, 1489-1503. (21) Tsuzuki, S.; Honda, K.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2002, 124, 104-112. (22) Tsuzuki, S.; Uchimaru, T.; Sugawara, K.; Mikami, M. J. Chem. Phys. 2002, 117, 11216-11221. (23) Sinnokrot, M. O.; Valeev, E. F.; Sherrill, C. D. J. Am. Chem. Soc. 2002, 124, 10887-10893. (24) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2004, 108, 10200-10207. (25) Cacelli, I.; Cinacchi, G.; Prampolini, G.; Tani, A. J. Am. Chem. Soc. 2004, 126, 14278-14286. (26) Zhikola, O. A.; Shishkin, O. V.; Lyssenko, K. A.; Leszczynski, J. J. Chem. Phys. 2005, 122, 144104. (27) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 4209-4212. (28) Sato, T.; Tsuneda, T.; Hirao, K. J. Chem. Phys. 2005, 123, 104307. (29) Diedrich, C.; Luchow, A.; Grimme, S. J. Chem. Phys 2005, 123, 184106. (30) Dobson, J. F.; McLennan, K.; Rubio, A.; Wang, J.; Gould, T.; Le, H. M.; Dinte, B. P. Aust. J. Chem. 2001, 54, 513-527. (31) Dion, M.; Rydberg, H.; Schro¨der, E.; Langreth, D. C.; Lundqvist, B. I. Phys. ReV. Lett. 2004, 92, 246401. (32) Misquitta, A. J.; Szalewicz, K. Chem. Phys. Lett. 2002, 357, 301306. (33) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Phys. ReV. Lett. 2003, 91, 033201. (34) Hesselmann, A.; Jansen G. Chem. Phys. Lett. 2002, 357, 464470; 2002, 362, 319-328; 2003, 367, 778-784. (35) Williams, H. L.; Chabalowski, C. F. J. Phys. Chem. A 2001, 105, 646-659. (36) Jeziorski, B.; Moszyn´ski, R.; Szalewicz, K. Chem. ReV. 1994, 94, 1887-1930. (37) Misquitta, A. J.; Szalewicz, K. J. Chem. Phys. 2005, 122, 214109. (38) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. J. Chem. Phys. 2005, 123, 214103. (39) Hesselmann, A.; Jansen, G.; Schu¨tz, M. J. Chem. Phys. 2005, 122, 014103. (40) Podeszwa, R.; Szalewicz, K. Chem. Phys. Lett. 2005, 412, 488493. (41) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 4993-4999. (42) Bukowski, R.; Podeszwa, R.; Szalewicz, K. Chem. Phys. Lett. 2005, 414, 111-116. (43) Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Chem. Theory Comput. 2006, 2, 400-412. (44) DALTON, a molecular electronic structure program, Release 2.0 2005, see http://www.kjemi.uio.no/software/dalton/dalton.html. (45) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865-3868. (46) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158-6170. (47) Tozer, D. J.; Handy, N. C. J. Chem. Phys. 1998, 109, 1018010189. (48) Lias, S. G. Ionization Energy Evaluation, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69 (http://webbook.nist.gov).

10354 J. Phys. Chem. A, Vol. 110, No. 34, 2006 (49) Tamagawa, K.; Iijima, T.; Kimura, M. J. Mol. Struct. 1976, 30, 243-253. (50) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796-6806. (51) Akin-Ojo, O.; Bukowski, R.; Szalewicz, K. J. Chem. Phys. 2003, 119, 8379-8396. (52) Williams, H. L.; Mas, E. M.; Szalewicz, K.; Jeziorski, B. J. Chem. Phys. 1995, 103, 7374. (53) Weigend, F.; Ko¨hn, A.; Ha¨ttig, C. J. Chem. Phys. 2002, 116, 31753183. (54) Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285-4291. (55) See the Supporting Information. (56) Available at http://www.physics.Udel.edu/∼szalewic/sapt/benzene. (57) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726-3741. (58) Wormer, P. E. S.; Hettema, H. POLCOR package, University of Nijmegen, 1992. (59) Kumar, A.; Meath, W. J. Mol. Phys. 1992, 75, 311-324. (60) Cerjan, C. J.; Miller, W. H. J. Chem. Phys. 1981, 75, 2800-2806. (61) Werner, H.-J.; et al. MOLPRO: a package of ab initio programs, version 2002.6. (62) Jorgensen, W.; Severance, D. J. Am. Chem. Soc. 1990, 112, 47684774. (63) Schmied, R.; Carcabal, P.; Dokter, A. M.; Lonij, V. P. A.; Lehmann, K. K.; Scoles, G. J. Chem. Phys. 2004, 121, 2701-2710. (64) Bacon, G. E.; Curry, N. A.; Wilson, S. A. Proc. R. Soc. London, Ser. A 1964, 279, 98-110.

Podeszwa et al. (65) Schweizer, W. B.; Dunitz, J. D. J. Chem. Theory Comput. 2006, 2, 288-291. (66) Patkowski, K.; Murdachaew, G.; Fou, C.-M.; Szalewicz, K. Mol. Phys. 2005, 103, 2031-2045. (67) Cencek, W.; Jeziorska, M.; Bukowski, R.; Jaszu_ski, M.; Jeziorski, B.; Szalewicz, K. J. Phys. Chem. A 2004, 108, 3211-3224. (68) Podeszwa, R.; Kucharski, S. A.; Stolarczyk, L. Z. J. Chem. Phys. 2002, 116, 480-493. (69) Mas, E. M.; Szalewicz, K.; Bukowski, R.; Jeziorski, B. J. Chem. Phys. 1997, 107, 4207-4218. (70) Stenger, F. SIAM ReV. 1981, 23, 165-224. (71) Smith, D. B.; Srivastava, R. Thermodynamic data for pure compounds: Part A Hydrocarbons and Ketones; Elsevier: Amsterdam, 1986. (72) Bich, E.; Hendl, H.; Lober, T.; Millat, J. Fluid Phase Equilib. 1992, 76, 199-211. (73) Mas, E. M.; Bukowski, R.; Szalewicz, K. J. Chem. Phys. 2003, 118, 4386-4403. (74) Mas, E. M.; Bukowski, R.; Szalewicz, K. J. Chem. Phys. 2003, 118, 4404-4413. (75) Lotrich, V.; Szalewicz, K. J. Chem. Phys. 1997, 106, 9668-9687. (76) Lotrich, V. F.; Szalewicz, K. Phys. ReV. Lett. 1997, 79, 13011304. (77) Tauer, T. P.; Sherrill, C. D. J. Phys. Chem. A 2005, 109, 1047510478.