Modeling π-π Interactions with the EFP Method

Modeling π–π Interactions with the Effective Fragment Potential Method: The Benzene Dimer and Substituents This study compares the results of the general effective fragment potential (EFP2) method to the results of a previous combined coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] and symmetry-adapted perturbation theory (SAPT) study [Sinnokrot and Sherrill, J. Am. Chem. Soc., 2004, 126, 7690] on substituent effects in π-π interactions. EFP2 is found to accurately model the binding energies of the benzene-benzene, benzene-phenol, benzene-toluene, benzene-fluorobenzene, and benzene-benzonitrile dimers, as compared with high-level methods [Sinnokrot and Sherrill, J. Am. Chem. Soc., 2004, 126, 7690], but at a fraction of the computational cost of CCSD(T). In addition, an EFP-based Monte Carlo/simulated annealing study was undertaken to examine the potential energy surface of the substituted dimers. 1. Introduction Intermolecular π-π interactions are among the major noncovalent forces controlling structural organization and recognition processes in biomolecules.2 Interactions between aromatic rings are largely responsible for DNA base-pair stacking,3 host-guest complexation,4–7 and the tertiary structure of proteins.8,9 Certain drugs rely on π-π interactions for intercalation into DNA.10,11 While these interactions have been studied extensively,12–24 their relative weakness and shallow potential energy surface makes them challenging to describe by either experiment or theory.25–29 The binding energy of the gas-phase benzene dimer, for example, is 2–3 kcal/mol, and the dimer is stable only at low temperatures.30 Only the highest levels of electronic structure theory can accurately capture these weak interactions. The fundamental nature of π-π interactions, along with the difficulties that theory encounters when targeting the multifaceted patterns of intermolecular bonding, has made the benzene dimer and its substituents popular benchmark systems for the accuracy and feasibility of new computational techniques and approaches.1,30–39 Both experimental18,28,40 and theoretical30–35 studies suggest that the benzene dimer has two almost isoenergetic geometries: T-shaped and paralleldisplaced. At the coupled cluster with singles, doubles, and perturbative triples [CCSD(T)]41 level in an estimated complete basis set limit, the gas-phase binding energies De (Do) of these geometries were calculated to be 2.7 (2.4) and 2.8 (2.7) kcal/mol, respectively.30 A sandwich geometry (faceto-face stacking) was found to be 1 kcal/mol higher in energy and is a transition state.30 Adding a substituent (OH, CH3, F, or CN) to one of the benzene rings in the dimer can result in interesting changes to the π-π interaction energies.1,37 Sherrill and co-workers studied the sandwich geometry and two different T-shaped geometries for several substituted benzene dimers (Figure 1). An estimated CCSD(T)/aug-cc-pVTZ level of theory predicts that all substituted sandwich dimers have a larger binding energy than the unsubstituted benzene dimer. The sandwich configuration exhibits a larger stabilizing effect due to substituents than either of the T-shaped configurations does, so the energy difference between these configurations is smaller for the substituted dimers than it is for the unsubstituted benzene dimer. For the Tshaped(1) structures (referred to as “T-shaped” in ref 1), the presence of an electron-withdrawing substituent (F, CN) results in a stronger interaction compared to the unsubstituted benzene dimer, while electron-donating substituents (OH, CH3) interact more weakly than the unsubstituted dimer does. For the T-shaped(2) dimers, the opposite is true. It is noteworthy that, as discussed in ref 1, the binding patterns found for the substituted benzene dimers contradict the Hunter-Sanders model,42 which qualitatively predicts binding energies based on electrostatics alone. This indicates the importance of the nonelectrostatic terms in determining the binding energies of these dimers. The same set of substituted benzene dimers was recently studied using the Van der Waals density functional theory (vdW-DFT).36 Most DFT methods do not accurately describe dispersion interactions, and generalized gradient approximation (GGA) functionals often fail to predict any binding in molecular dimers.36 However, the vdW-DFT method,43 which incorporates a nonlocal correlation energy, has been shown36 to reproduce the energetic ordering of the substituents consistent with the CCSD(T) results of ref 1. But, the interaction energies of the sandwich dimers are overestimated by 0.6–0.7 kcal/mol, compared to CCSD(T), whereas the binding in the T-shaped dimers is underestimated by 0.2–0.4 kcal/mol. In this contribution, the interactions in mixed substituted benzene dimers are studied by means of the general effective fragment potential (EFP2) method.44,45 A benchmark study is presented, comparing the performance of EFP2 with the CCSD(T) and symmetry-adapted perturbation theory (SAPT)46 results of Sinnokrot and Sherrill.1 EFP2 is an ab initio-based model potential method for treating intermolecular interactions, wherein parameters for all major types of noncovalent forces;Coulomb, exchange-repulsion, polarization (induction), and dispersion;are obtained for each unique fragment from a single ab initio calculation. Contrary to many other model potential methods, EFP2 does not employ empirically fitted parameters. This method can then be used to model interactions with other EFP fragments or with fully ab initio molecules. EFP fragments have frozen internal geometries. Figure 1. Substituted benzene dimer geometries. The geometry that is called “T-shaped” in ref 1 has been renamed “T-shaped(1)”. The computational cost of EFP2 is orders of magnitude lower than that of ab initio techniques. For example, for the benzene dimer in the 6-311+G(3df,2p) basis set with 660 basis functions, a single-point energy calculation with second order Moller-Plesset (MP2) perturbation theory47 would require 142 min of CPU time on one IBM Power5 processor, whereas the analogous EFP calculation requires only 0.4 s.48 If EFP2 can be shown to produce results of sufficient accuracy, the advantage of using it in place of more computationally costly methods, such as CCSD(T) or MP2 theory, is clear. The effective fragment potential method has been shown to successfully model numerous hydrogen-bonded systems. For example, EFP1 (the EFP-type model designed specifically for water, with a fitted exchange-repulsion term) has been used to model chemical reactions in solution,49,50 solvent clusters,51,52 an SN2 reaction,53 and amino acid neutral/zwitterion equilibria.54,55 Styrene-styrene interactions56 and methanol-water mixing57 have been studied by the general EFP2 method. Recently, EFP2 was shown to successfully model the π-π interactions in the unsubstituted benzene dimer.48 An EFP2 study of benzene-water complexes is in progress.58 The main goal of the present work is to further investigate the performance of EFP2 for systems with π-π interactions, to confidently model various biological systems (e.g., interactions in DNA base pairs) in future studies. Therefore, the first part of this work analyzes the accuracy of EFP2 on a set of mixed benzene-substituted benzene dimers, comparing the EFP2 results with those of CCSD(T) and SAPT theory from ref 1. The second part of this paper undertakes an independent EFP2 study of the dimer potential energy surfaces and bonding patterns. This work is organized as follows. Section 2 describes the theoretical methods and computational details. Section 3 presents and discusses the results. The main conclusions are given in Section 4. 2. Theoretical Methods and Computational Details The EFP2 method has been implemented in the quantum chemistry program GAMESS,59 which was used for all calculations in this study. Dimer Geometries. To benchmark the performance of EFP2 theory versus CCSD(T) and SAPT theory, the monomer and dimer geometries were adopted from ref 1. The procedure used in ref 1 to find the dimer structures was as follows. The geometries of all monomers were optimized with MP2 theory and the aug-cc-pVTZ basis set.60 For benzene, a C-C bond length of 1.3942 Å and a C-H bond length of 1.0823 Å were obtained. The geometry of the toluene monomer was fixed at Cs symmetry, with one methyl hydrogen above the benzene ring and two below it, and phenol was chosen to be planar. In all dimer calculations in ref 1, the monomers were held rigid, having fixed internal coordinates. The dimer geometries were obtained by moving the rigid monomers relative to each other to find the optimal intermonomer distance (defined as the distance between ring centers) for the sandwich, T-shaped(1), and T-shaped(2) configurations (see Figure 1). These structures are not actual energy minima for the dimers in question, but are constrained model structures chosen in ref 1 for the purpose of assessing the influence of substituents on the π-π interactions. For example, to minimize direct interactions between the substituent on one ring and the other unsubstituted ring, the substituent in the T-shaped dimers is held in a para position, relative to the unsubstituted benzene. This serves to maintain focus on the effect of the substituent on the π-π interactions. In this study, the monomer geometries from ref 1 were used to generate EFP2 interaction parameters with the 6-311++G(3df,2p) basis set. The underlying 6-311G basis set is accurate for modeling electrostatic interactions, whereas polarization functions are important for modeling dispersion interactions57 and diffuse functions are required for accurate exchangerepulsion interactions.45 Distributed multipoles on atoms and bond midpoints were generated using a numerical integration scheme.48 The electrostatic charge-charge, charge-dipole, charge-quadrupole, and dipole-dipole energies were screened by charge-penetration damping functions, as described in ref 48. As in ref 1, the intermonomer distances in each of the chosen dimer configurations (Figure 1) were varied to find the optimal EFP2 dimer geometry. Performing the geometry optimizations in this way facilitates a direct comparison of the optimal intermonomer distances (R) and energies (Eint) with the corresponding CCSD(T) values from ref 1. Energy Component Comparison. To analyze the accuracy of each component (Coulomb, exchange-repulsion, polarization, and dispersion) of the EFP2 energy, separate EFP2 calculations were performed at the dimer geometries used in the SAPT/aug- TABLE 1: Optimized Intermonomer Distances (R) and Interaction Energies (Eint) for the Sandwich, T-Shaped(1), and T-Shaped(2) Structures of the Benzene-Substituted Dimersa Sandwich X H OH CH3 F CN a T-Shaped(1) T-Shaped(2) method R (Å) Eint (kcal/mol) R (Å) Eint (kcal/mol) R (Å) Eint (kcal/mol) MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ est. CCSD(T)/aug-cc-pVTZ EFP2 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ est. CCSD(T)/aug-cc-pVTZ EFP2 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ est. CCSD(T)/aug-cc-pVTZ EFP2 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ est. CCSD(T)/aug-cc-pVTZ EFP2 MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ est. CCSD(T)/aug-cc-pVTZ EFP2 3.80 3.70 3.90 3.95 3.70 3.60 3.80 3.90 3.70 3.65 3.80 3.90 3.70 3.70 3.80 3.90 3.70 3.60 3.80 3.85 -2.90 -3.26 -1.80 -2.17 -3.40 -3.75 -2.17 -2.72 -3.58 -3.96 -2.27 -2.78 -3.50 -3.81 -2.29 -3.02 -4.49 -4.86 -3.05 -3.91 5.01 4.89 4.99 5.15 5.00 4.90 5.00 5.15 5.00 4.90 5.00 5.20 4.95 4.90 5.00 5.15 4.90 4.80 4.90 5.15 -3.16 -3.46 -2.62 -2.42 -3.14 -3.42 -2.58 -2.54 -3.11 -3.39 -2.55 -2.47 -3.35 -3.61 -2.77 -2.79 -3.79 -4.11 -3.25 -3.20 5.01 4.89 4.99 5.15 4.95 4.90 5.00 5.15 4.90 4.80 5.00 5.15 5.00 4.90 5.00 5.20 5.00 4.90 5.00 5.15 -3.16 -3.46 -2.62 -2.42 -3.23 -3.52 -2.67 -2.45 -3.60 -3.89 -2.95 -2.95 -2.87 -3.17 -2.38 -2.30 -2.82 -3.08 -2.20 -2.23 MP2 and estimated CCSD(T) results taken from ref 1. cc-pVDZ calculations in ref 1. Sandwich dimers have an intermonomer separation of 3.70 Å , whereas a separation of 4.90 Å was used for the T-shaped(1) and T-shaped(2) dimers. Investigation of the EFP2 Potential Energy Surface. In addition to direct comparisons of the accuracy of EFP2 to that of CCSD(T), MP2, and SAPT, a Monte Carlo/simulated annealing (MC/SA)61 study was performed on the EFP2 dimers to investigate the potential energy surface of each dimer. To increase the conformational sampling, two temperature ranges were used in each case: 20000-100 K and 3000-100 K. A sandwich geometry was used as the starting structure for each dimer. Geometry optimizations were performed every 10 steps. 3. Results and Discussion Intermonomer Separations and Binding Energies. Table 1 compares the optimal EFP2 distances and binding energies of the substituted dimers with the distances and energies found with the MP2 and CCSD(T) methods, as given in ref 1 for the constrained sandwich, T-shaped(1), and T-shaped(2) structures. The optimal EFP2 intermonomer distances are consistently larger than the estimated CCSD(T)/aug-cc-pVTZ distances by 0.05-0.10 Å for the sandwich dimers, 0.15-0.25 Å for the T-shaped(1) dimers, and 0.15-0.21 Å for the T-shaped(2) dimers. Compared to CCSD(T), MP2 with the aug-cc-pVTZ basis set consistently underestimates the intermonomer distances by 0.1-0.2 Å for all dimer geometries. MP2/aug-cc-pVDZ results are in better agreement with those of CCSD(T), underestimating the intermonomer distances by, at most, 0.1 Å. This suggests that there is a fortuitous cancellation of errors for MP2 theory with the smaller basis set. For the sandwich dimers, EFP2 overestimates the magnitude of the binding energies by 0.36-0.86 kcal/mol, compared to CCSD(T). In all cases, this represents a significant improvement over MP2, which overestimates the binding energy by 1.1-1.8 kcal/mol, compared to the CCSD(T) values. EFP2 most accurately models the unsubstituted benzene sandwich dimer, giving an error of 0.37 kcal/mol in the binding energy, in comparison to the estimated CCSD(T) energy. The benzene– benzonitrile sandwich dimer binding energy has the greatest error, at 0.86 kcal/mol. Compared to CCSD(T), MP2 underestimates the intermonomer separation and overestimates the binding energies of the dimer, whereas EFP2 overestimates both the intermonomer distances and binding energies. The overestimation of binding energies occurs to a much smaller degree with EFP2 than with MP2. A disconcerting observation is that the MP2 error in the binding energy increases when the basis set is improved. EFP2 is very accurate in calculating the binding energies of the T-shaped(1) and T-shaped(2) dimers. The EFP2 binding energy of the T-shaped(1) unsubstituted benzene dimer (identical to the T-shaped(2) unsubstituted dimer) is higher (less strongly bound) than the corresponding CCSD(T) energy, by 0.20 kcal/ mol. Of the substituted T-shaped(2) dimers, benzene-phenol is the most in error, at 0.22 kcal/mol higher in energy. Discrepancies in the other EFP2 T-shaped dimers range from 0.02 kcal/mol lower in energy (more strongly bound) to 0.08 kcal/mol higher in energy, compared to CCSD(T). The EFP2 benzene-toluene energy is in perfect agreement with the estimated CCSD(T) energy for that dimer. Contrary to the excellent performance of EFP2, MP2 overestimates the binding energies of T-shaped(1) and T-shaped(2) dimers, becoming even less accurate when increasing the basis set from aug-cc-pVDZ to aug-cc-pVTZ. These EFP2 results also compare favorably with the results of a DFT study36 performed on the same set of substituted benzene dimers using a Van der Waals (vdW)-corrected density functional.43 This vdW-DFT method overestimates the magnitude of the sandwich dimer interaction energies by 0.56-0.68 kcal/mol compared with CCSD(T) and underestimates the binding in the T-shaped dimers, by 0.28-0.39 kcal/mol in the T-shaped(1) and by 0.19-0.34 kcal/mol in the T-shaped(2) dimers. In comparison, EFP2 overbinds the sandwich dimers by 0.36- 0.86 kcal/mol compared to CCSD(T) and differs from the CCSD(T) results for the T-shaped dimers by <0.1 kcal/ mol (except for the unsubstituted dimer and the T-shaped(2) benzene-phenol dimer, which are 0.2 kcal/mol less strongly bound). EFP2 requires considerably less computer time (on the order of seconds) after the MAKEFP potentials have been TABLE 2: Contributions to the Binding Energy (kcal/mol) for the Sandwich, T-Shaped(1), and T-Shaped(2) Dimersa,b Contribution to Binding Energy (kcal/mol) X H OH CH3 F CN Sandwich T-shaped(1) T-shaped(2) Energy Component SAPT EFP2 SAPT EFP2 SAPT EFP2 Coulomb ex.-repulsion polarization dispersion total energy Coulomb ex.-repulsion polarization dispersion total energy Coulomb ex.-repulsion polarization dispersion total energy Coulomb ex.-repulsion polarization dispersion total energy Coulomb ex.-repulsion polarization dispersion total energy -0.97 6.03 -0.33 -6.53 -1.80 -1.08 5.85 -0.34 -6.72 -2.29 -1.03 6.21 -0.40 -7.19 -2.41 -1.36 5.73 -0.29 -6.49 -2.40 -1.83 5.78 -0.29 -7.01 -3.36 -0.70 6.03 -0.41 -6.38 -1.46 -1.03 5.75 -0.39 -6.55 -2.22 -0.65 5.85 -0.42 -7.03 -2.26 -1.54 5.63 -0.31 -6.31 -2.53 -1.92 5.55 -0.29 -6.87 -3.54 -2.24 4.87 -0.67 -4.37 -2.42 -2.21 4.85 -0.66 -4.37 -2.39 -2.24 5.02 -0.67 -4.46 -2.34 -2.27 4.55 -0.68 -4.22 -2.63 -2.59 4.59 -0.83 -4.29 -3.12 -2.71 4.96 -0.31 -3.83 -1.89 -2.69 5.05 -0.33 -4.07 -2.04 -2.63 5.24 -0.34 -4.16 -1.89 -2.82 4.72 -0.37 -3.92 -2.39 -3.22 4.77 -0.53 -3.83 -2.81 -2.24 4.87 -0.67 -4.37 -2.42 -2.18 4.78 -0.66 -4.41 -2.47 -2.38 4.81 -0.70 -4.59 -2.85 -1.98 4.73 -0.57 -4.30 -2.12 -1.73 4.69 -0.51 -4.53 -2.09 -2.71 4.96 -0.31 -3.83 -1.89 -2.64 4.86 -0.27 -3.86 -1.91 -2.69 4.76 -0.35 -4.26 -2.54 -2.43 4.82 -0.20 -3.92 -1.73 -1.96 4.49 -0.10 -4.20 -1.77 a The intermonomer separations are 3.70 Å in the sandwich dimers and 4.90 Å in the T-shaped(1) and T-shaped(2) dimers. b SAPT results are taken from ref 1. generated. In contrast, the GGA portion of the vdW-DFT calculations requires 1.3 h on a single Opteron processor.36 Energy Component Comparison. Table 2 and Figures 2, 3, and 4 compare the SAPT (from ref 1) and EFP2 Coulomb, exchange-repulsion, polarization (or induction), and dispersion energies for the three types of constrained geometries. The energy component comparison was performed at fixed intermonomer distances (3.70 Å for the sandwich structures, 4.90 Å for both types of T-shaped structures) that were smaller than the EFP2 equilibrium intermonomer distances but close to the estimated CCSD(T) equilibrium distances (see Table 1). As described in ref 1, SAPT predicts that dispersion is the single greatest attractive contribution to the overall binding energy of the dimers, especially for the sandwich structures. The sandwich dimers with electron-withdrawing substituents (CN and F) have the most favorable Coulomb interaction energy coupled with the lowest exchange-repulsion energy, giving them the strongest sandwich binding energies. The order of the total binding energies for the other sandwich dimers corresponds to the order of their dispersion energies. Because of the greater intermonomer separation in the T-shaped dimers, the exchangerepulsion and dispersion components of their binding energies are smaller than those of the sandwich structures. The Coulomb attractions in the T-shaped(1) and T-shaped(2) dimers are stronger, compared to the sandwich structures, because the negatively charged π-cloud of one monomer interacts with the positively charged hydrogens of the other. Overall, good agreement is found between the EFP and SAPT energy terms (see Table 2 and Figures 2, 3, and 4). The EFP2 dispersion interaction is consistently lower in magnitude than that predicted by SAPT, but always by <0.5 kcal/mol; dispersion remains the dominant attractive contribution to the binding energy for both the sandwich and T-shaped dimers. The EFP2 polarization terms for the sandwich dimers agree to within ∼0.1 kcal/mol with those predicted by SAPT. For the T-shaped(1) Figure 2. (a) Electrostatic, exchange-repulsion, polarization (induction), and dispersion components and (b) total interaction energies for the sandwich dimers by EFP2 and SAPT. SAPT data are taken from ref 1. and T-shaped(2) dimers, the polarization terms are smaller in magnitude by 0.3–0.4 kcal/mol, compared with the SAPT results. The EFP2 Coulomb term is slightly more negative than Figure 3. (a) Electrostatic, exchange-repulsion, polarization (induction), and dispersion components and (b) total interaction energies for the T-shaped(1) dimers by EFP2 and SAPT. SAPT data are taken from ref 1. the corresponding SAPT term in the sandwich dimers with electron-withdrawing substituents (CN, F), but less negative in the benzene-benzene and benzene-toluene dimers. The EFP2 Coulomb interaction for the T-shaped structures is ∼0.5 kcal/ mol lower in energy, on average, compared to SAPT. The EFP2 exchange-repulsion generally differs from that obtained from SAPT by <0.1 kcal/mol. For the sandwich structures (Figure 2), the SAPT binding energy increases in magnitude in the following order: H < OH < F ≈ CH3 < CN. The order predicted by both EFP2 and CCSD(T) is a very similar: H < OH < CH3 < F < CN (see Table 1). Overall, the trends predicted by EFP2 for the sandwich dimers agree with those found with SAPT in ref 1. In the first set of T-shaped dimer structures (see panels a and b in Figure 3), the magnitudes of the dispersion and the exchange-repulsion terms are significantly reduced, compared to their values in the sandwich dimer structures. This is reflected in both the SAPT and the EFP2 results. The order of increasing strength of the total interaction energy is CH3 ≈ OH ≈ H < F < CN and CH3 ≈ H < OH < F < CN for SAPT and EFP2, respectively. Taking into account the very small (<0.1 kcal/ mol) binding energy differences between CH3-, OH-, and H-substituted dimers, the agreement between SAPT and EFP2 is very reasonable. The overall EFP2 binding energies are 0.2-0.5 kcal/mol lower in magnitude than SAPT binding energies. This reflects the fact that the chosen intermonomer separations are shorter than the optimal EFP2 geometries. SAPT predicts the order of increasing binding energies for the T-shaped(2) structures to be CN ≈ F < H ≈ OH < CH3, whereas EFP2 predicts F < CN (see panels a and b in Figure Figure 4. (a) Electrostatic, exchange-repulsion, polarization (induction), and dispersion components and (b) total interaction energies for the T-shaped(2) dimers by EFP2 and SAPT. SAPT data are taken from ref 1. 4). However, these energy differences are very small. Generally, the EFP2 exchange-repulsion energy in T-shaped(2) dimers closely resembles the SAPT exchange-repulsion. The EFP Coulomb terms are more attractive by 0.2-0.5 kcal/mol, compared to SAPT. As with the T-shaped(1) dimers, the magnitude of the T-shaped(2) EFP2 polarization and dispersion energies are slightly lower than those for SAPT; this results in some error cancellation with the Coulomb term. To summarize, trends that are due to substituent effects in each interaction energy component, as well as in the total binding energies, are well-reproduced by EFP2 (see Figures 2, 3, and 4). This is encouraging for future EFP2 studies of more complex species. The Potential Energy Surface. The conformational space of each dimer was explored using Monte Carlo/simulated annealing with the EFP2 method. The lowest energy structures for each dimer are shown in Figures 5, 6, 7, 8, and 9. All geometries shown have positive definite Hessians. The three minima found by EFP2 on the benzene dimer potential energy surface are parallel-displaced (Figure 5b) structures, tilted T-shaped structures (Figure 5a), and edge-toedge structures (Figure 5c), in agreement with recent SAPT(DFT) studies by Podeszwa et al.33 The parallel-displaced and tilted T-shaped configurations are predicted to be almost isoenergetic, with CCSD(T) favoring the latter structure by 0.1 kcal/mol.30,33 EFP2 is in agreement with CCSD(T), predicting the tilted T-shaped structure to be the global minimum, with the parallel-displaced configuration being ∼0.4 kcal/mol higher in energy. The edge-to-edge configuration is ∼1 kcal/mol higher Figure 5. Lowest-energy benzene-benzene structures found with EFP2 Monte Carlo/simulated annealing. Energies are given in units of kcal/mol. Figure 6. Lowest energy benzene-phenol structures, using a planar phenol monomer, found with EFP2 Monte Carlo/simulated annealing. Energies are given in units of kcal/mol. in energy than the parallel-displaced one. The EFP2 paralleldisplaced configuration is slightly less displaced48 than is predicted by MP2,30 CCSD(T),31 or SAPT(DFT).33 As discussed in ref 48, the discrepancy in the R2 distance (the “displacement coordinate” defined in ref 30) between MP2 or CCSD(T) and EFP2 is 0.4 Å. That is, the EFP2 parallel-displaced structure is less relaxed (closer to the sandwich structure) than are the MP2 or CCSD(T) isomers. This results in a slightly higher EFP2 energy for the parallel-displaced structure and, consequently, a slightly larger (0.4 kcal/mol) energy difference between the tilted T-shaped and parallel-displaced configurations, compared to 0.1 kcal/mol for both CCSD(T) and SAPT(DFT). For the substituted dimers, unconstrained geometry optimizations lead to minima that are similar to those of the benzene dimer; that is, configurations resembling the parallel-displaced, T-shaped, or edge-to-edge structures. However, the predicted energetic ordering of these conformations depends on the substituent. The origins of these differences are analyzed in the following discussion. 5292 J. Phys. Chem. A, Vol. 112, No. 23, 2008 Smith et al. Figure 7. Lowest energy benzene-toluene structures found with EFP2 Monte Carlo/simulated annealing. Energies are given in units of kcal/mol. Figure 8. Lowest-energy benzene-fluorobenzene structures found with EFP2 Monte Carlo/simulated annealing. Energies are given in units of kcal/mol. The parallel-displaced structure is common to all dimers. As in the benzene dimer, binding in the substituted paralleldisplaced dimers is dominated by dispersion. As discussed previously for the constrained sandwich dimers, substituents significantly increase the binding energies in the paralleldisplaced dimers (by 1.3–2.8 kcal/mol). The CN-substituted dimer (Figure 9a) is the most strongly bound, followed by OH (Figure 6b), F (Figure 8a), and CH3 (Figure 7a). The binding energies of the latter three structures are very similar and 1.3–1.4 kcal/mol larger than that in the parallel-displaced benzene dimer. Binding energies in the unconstrained substituted paralleldisplaced dimers are 0.6–1.1 kcal/mol stronger than those in the corresponding constrained sandwich dimers. The origin of the stronger binding is partially the electron-donating/electronwithdrawing effects of the substituents and partially the interaction between the substituent and the benzene ring in the unconstrained structures. For example, stabilization in the OH, F, and CN substituted dimers is partly due to favorable Coulomb interactions between the electronegative substituent and the positively charged hydrogens of the unsubstituted benzene. This is reflected in both the significantly larger Coulomb energies, and the slightly tilted dimer structures in which the substituent is angled closer to the benzene. In contrast, the CH3 substituted dimer is perfectly parallel, and the additional stabilization in Figure 9. Lowest-energy benzene-benzonitrile structures found with EFP2 Monte Carlo/simulated annealing. Energies are given in units of kcal/ mol. this species is mainly due to the dispersion forces between the CH3 groupandthebenzenering.Consequently,thebenzene-toluene dimer (Figure 7a) has the smallest binding energy of the substituted parallel-displaced species. Binding in the T-shaped dimers is governed by a balance between the dispersion and Coulomb interactions. Among the unconstrained EFP2 T-shaped structures, the CN-substituted dimer (Figure 9b) has the largest binding energy of 4.7 kcal/ mol, closely followed by the OH dimer (Figure 6a). Two CH3substituted dimers (panels b and c in Figure 7) and the F-substituted dimer (Figure 8b) have binding energies of 3.0–3.3 kcal/mol. As in the parallel-displaced structures, the benzene– benzonitrile and benzene–fluorobenzene dimers are stabilized through Coulomb interactions between the electronegative substituents and the unsubstituted benzene ring. The T-shaped benzene-phenol structure is different from other T-shaped dimers, because its binding is dominated by the interaction between the partially positive hydrogen of the hydroxyl group and the negative π-cloud of the unsubstituted benzene. The interaction in benzene-phenol is thus similar to that in the benzene-water dimer.62 The important role of the hydroxyl group in the binding of benzene-phenol in the T-shaped isomer is also reflected in the modestly larger polarization term, as compared to other dimers. Both conformations of the T-shapedlike benzene-toluene dimer are slightly stabilized by dispersion. All substituted dimer edge-to-edge structures are also bound more strongly than is the analogous benzene dimer isomer. Additional stabilization in OH-, F-, and CN-substituted dimers occurs due to Coulomb interactions between the positively charged hydrogens of the benzene ring and the partial negative charge on the substituent group. Indeed, in these dimers, a weak edge-to-edge interaction between the two benzene rings is augmented by a stronger edge-to-substituent interaction. As for the T-shaped dimers, the strength of the interaction in the edgeto-edge-type dimers decreases in the following order: CN > OH > F. Now, consider how the substituents affect the relative energies of the different conformations of the dimers. Based on the preceding analyses, it is expected that the EFP2 interactions are within 0.5 kcal/mol or less of the CCSD(T) values, so the trends are generally qualitatively and quantitatively captured by the EFP2 method. In all but the benzene-benzene and phenol-benzene dimers, the most strongly bound EFP2 structure is parallel-displaced rather than T-shaped. A rationale for this is that, in the offset structures, the optimal benzene-ring orientation remains almost unchanged, but additional stabilization occurs because of favorable Coulomb and dispersion interactions of the substituent group with the unsubstituted benzene. In the substituted T-shaped dimers the interaction between the benzene rings is sterically less favorable than in the unsubstituted T-shaped benzene dimer. The exception to this trend is the benzene-phenol dimer, in which the T-shapedlike configuration is lower in energy than the parallel-displaced configuration, by 0.8 kcal/mol. This is because the phenol hydroxyl group strongly interacts with the π-cloud of the benzene. 4. Conclusion This study explores the applicability and accuracy of the general effective fragment potential (EFP2) method, the firstprinciples-based model potential, to systems with aromatic π-π interactions (i.e., benzene–benzene, benzene–toluene, benzene– phenol, benzene–fluorobenzene, and benzene–benzonitrile dimers). The first part of the paper compares the EFP binding energies for the constrained geometries based on the work of Sinnokrot and Sherrill.1 The EFP2 binding energies are in excellent agreement with the CCSD(T) results in the constrained T-shaped dimers and overestimate those in the sandwich dimers by <0.9 kcal/mol. For these aromatic complexes, EFP2 is more accurate than MP2. The latter, compared to CCSD(T), overestimates binding in both the T-shaped and sandwich structures by up to 1.0 and 1.8 kcal/mol, respectively. These results are excellent in view of the computational cost of each method. For example, a single-point energy calculation of the benzene dimer in the 6-311++G(3df,2p) basis with MP2 would require 142 min of CPU time on one IBM Power5 processor, whereas the analogous EFP2 calculation requires only 0.4 s. Some of the discrepancy between EFP2 and higher levels of theory such as CCSD(T) likely occurs because of the omission of some higher-order terms in the EFP2 expansions. For instance, the inclusion of induced quadrupoles in the polarization energies and higher-order terms in the dispersion expansion can be important. These terms will be targeted in future developments of the EFP2 method. It is very encouraging that EFP2 correctly reflects the energy changes due to different substituents. EFP2 also provides excellent agreement with symmetry-adapted perturbation theory (SAPT) for the different components of the total binding energy;Coulomb, exchange-repulsion, polarization, and dispersion;demonstrating the utility of EFP2 as an accurate and computationally inexpensive tool for the analysis of binding patterns in molecular complexes and (potentially) liquids. Finally, an independent EFP2 study was performed of the potential energy surface of each dimer by employing a Monte Carlo/simulated annealing technique. The CN, F, and CH3 substituents stabilize and favor the parallel-displaced configurations. Binding in the benzene-phenol dimer resembles that of the water-benzene complex. Acknowledgment. This work was supported in part by a Scientific Discovery through Advanced Computing (SciDAC) Department of Energy grant and in part by a NIRT grant from the National Science Foundation. The authors are grateful to Professor Brett Bode for providing the MacMolPlt program [J. Molec. Graphics 1999, 16, 133] for visualization. Supporting Information Available: The coordinates of the structures found with Monte Carlo/simulated annealing, as shown in Figures 5–9, are available as a text file. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Sinnokrot, M. O.; Sherrill, C. D. J. Am. Chem. Soc. 2004, 126, 7690. (2) Meyer, E. A.; Castellano, R. K.; Diederich, R. Angew. Chem., Int. Ed. Engl. 2003, 42, 1210. (3) Saenger, W. Principles of Nucleic Acid Structure; Springer-Verlag: New York, 1984. (4) Askew, B.; Ballester, P.; Buhr, C.; Jeong, K. S.; Jones, S.; Parris, K.; Williams, K.; Rebek, J., Jr J. Am. Chem. Soc. 1989, 111, 1082. (5) Smithrud, D. B.; Diederich, F. J. Am. Chem. Soc. 1990, 112, 339. (6) Hunter, C. A. Chem. Soc. ReV. 1994, 23, 101. (7) Rebek, J., Jr Chem. Soc. ReV. 1996, 25, 255. (8) Burley, S. K.; Petsko, G. A. Science 1985, 229, 23. (9) Hunter, C. A.; Singh, J.; Thornton, J. M. J. Mol. Biol. 1991, 218, 837. (10) Lerman, L. S. J. Mol. Biol. 1961, 3, 18. (11) Brana, M. F.; Cacho, M.; Gradillas, A.; Pascual-Teresa, B.; Ramos, A. Curr. Pharm. Des. 2001, 7, 1745. (12) Tsuzuki, S.; Uchimaru, T.; Tanabe, K. J. Mol. Struct. (THEOCHEM) 1994, 307, 107. (13) Hobza, P.; Selzle, H. L.; Schlag, E. W. J. Phys. Chem. 1996, 100, 18–790. (14) Tsuzuki, S.; Luthli, H. P. J. Chem. Phys. 2001, 114, 3949. (15) Steed, J. M.; Dixon, T. A.; Klemperer, W. J. Chem. Phys. 1979, 70, 4940. (16) Arunan, E.; Gutowsky, H. S. J. Chem. Phys. 1993, 98, 4294. (17) Law, K. S.; Schauer, M.; Bernstein, E. R. J. Chem. Phys. 1984, 81, 4871. (18) Felker, P. M.; Maxton, P. M.; Schaeffer, M.W. Chem. ReV. 1994, 94, 1787. (19) Venturo, V. A.; Felker, P. M. J. Chem. Phys. 1993, 98, 4294. (20) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2002, 124, 104. (21) Hobza, P.; Jurecka, P. J. Am. Chem. Soc. 2003, 125, 15–608. (22) Claessens, C. G.; Stoddart, J. F. J. Phys. Org. Chem. 1997, 10, 254. (23) Fyfe, M. C. T.; Stoddart, J. F. Acc. Chem. Res. 1997, 30, 393. (24) Dahl, T. Acta Chem. Scand. 1994, 48, 95. (25) Rappé, A. K.; Bernstein, E. R. J. Phys. Chem. A 2000, 104, 6117. (26) Kim, K. S.; Tarakeshwar, P.; Lee, J. Y. Chem. ReV. 2000, 100, 4145. (27) Dunning, T. H. J. Phys. Chem. A 2000, 104, 9062. (28) Müller-Dethlefs, K.; Hobza, P. Chem. ReV. 2000, 100, 143. (29) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2006, 110, 10656. (30) Sinnokrot, M. O.; Valeev, E. F.; Sherrill, C. D. J. Am. Chem. Soc. 2002, 124, 10887. (31) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2004, 108, 10200. (32) Jaffe, R. L.; Smith, G. D. J. Chem. Phys. 1996, 105, 2780. (33) Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Phys. Chem. A 2006, 110, 10345. (34) Puzder, A.; Dion, M.; Langreth, D. C. J. Chem. Phys. 2006, 124, 164105. (35) Hill, J. G.; Platts, J. A.; Werner, H. J. Phys. Chem. Chem. Phys. 2006, 8, 4072. (36) Thonhauser, T.; Puzder, A.; Langreth, D.C. J. Chem. Phys. 2006, 124, 164106. (37) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2003, 107, 8377. (38) DiStasio, R. A., Jr.; von Helden, G.; Steele, R. P.; Head-Gordon, M. Chem. Phys. Lett. 2007, 437, 277. (39) Tauer, T. P.; Derrick, M. E.; Sherrill, C. D. J. Phys. Chem. A 2005, 109, 191. (40) Sun, S.; Bernstein, E. R. J. Phys. Chem 1996, 100, 13348. (Also see ref 28.). (41) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (42) Hunter, C. A.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112, 5525. (43) Dion, M; Rydberg, H.; Schrder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. ReV. Lett. 2004, 92, 246401. (44) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. J. Phys. Chem. A 2001, 105, 293. (45) Jensen, J. H.; Gordon, M. S. Mol. Phys. 1996, 89, 1313. (46) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. ReV 1994, 94, 1887. (47) Möller, C; Plesset, M. S. Phys. ReV. 1934, 46, 618. (48) Slipchenko, L. V.; Gordon, M. S. J. Comput. Chem. 2007, 28, 276. (49) Chen, W.; Gordon, M. S. J. Chem. Phys. 1996, 105, 11081. (50) Webb, S. P.; Gordon, M. S. J. Phys. Chem. A 1999, 103, 1265. (51) Merrill, G.N.; Gordon, M. S. J. Phys. Chem. A 1998, 102, 2650. (52) Day, P. N.; Pachter, R.; Gordon, M. S. J. Chem. Phys. 2000, 112, 2063. (53) Adamovic, I.; Gordon, M. S. J. Phys. Chem. A 2005, 109, 1629. (54) Bandyopadhyay, P.; Gordon, M. S. J. Chem. Phys. 2000, 113, 1104. (55) Bandyopadhyay, P.; Gordon, M. S; Mennucci, B.; Tomasi, J. J. Chem. Phys. 2002, 116, 5023. (56) Adamovic, I.; Li, H.; Lamm, M. H.; Gordon, M. S. J. Phys. Chem. A 2006, 110, 519. (57) Adamovic, I.; Gordon, M. S. J. Phys. Chem. 2006, 110, 10267. (58) Slipchenko, L. V.; Gordon, M. S. Manuscript in preparation. (59) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr J. Comput. Chem. 1993, 14, 1347. (60) Woon, D. E.; Dunning, T. H., Jr J. Chem. Phys. 1993, 98, 1358. (61) Day, P. N.; Pachter, R.; Gordon, M. S.; Merrill, G. N. J. Chem. Phys. 2000, 112, 2063. (62) Suzuki, S.; Green, P. G.; Bumgarner, R. E.; Dasgupta, S.; Goddard, W. A.; Blake, G. A. Science 1992, 257, 942. JP800107Z