Toward True DNA Base-Stacking Energies: MP2, CCSD(T), and Complete Basis Set Calculations Pavel Hobza* and Jir�� �� �� Sponer �� Contribution from the J. HeyroVsky �� Institute of Physical Chemistry, Academy of Sciences of the Czech Republic and Center for Complex Molecular Systems and Biomolecules, DolejskoVa �� 3, 182 23 Prague 8, Czech Republic Received May 2, 2002 Abstract: Stacking energies in low-energy geometries of pyrimidine, uracil, cytosine, and guanine homodimers were determined by the MP2 and CCSD(T) calculations utilizing a wide range of split-valence, correlation-consistent, and bond-functions basis sets. Complete basis set MP2 (CBS MP2) stacking energies extrapolated using aug-cc-pVXZ (X ) D, T, and for pyrimidine dimer Q) basis sets equal to -5.3, -12.3, and -11.2 kcal/mol for the first three dimers, respectively. Higher-order correlation corrections estimated as the difference between MP2 and CCSD(T) stacking energies amount to 2.0, 0.7, and 0.9 kcal/mol and lead to final estimates of the genuine stacking energies for the three dimers of -3.4, -11.6, and -10.4 kcal/mol. The CBS MP2 stacking-energy estimate for guanine dimer (-14.8 kcal/mol) was based on the 6-31G*(0.25) and aug-cc-pVDZ calculations. This simplified extrapolation can be routinely used with a meaningful accuracy around 1 kcal/mol for large aromatic stacking clusters. The final estimate of the guanine stacking energy after the CCSD(T) correction amounts to -12.9 kcal/mol. The MP2/6-31G*(0.25) method previously used as the standard level to calculate aromatic stacking in hundreds of geometries of nucleobase dimers systematically underestimates the base stacking by ca. 1.0-2.5 kcal/mol per stacked dimer, covering 75-90% of the intermolecular correlation stabilization. We suggest that this correction is to be considered in calibration of force fields and other cheaper computational methods. The quality of the MP2/6-31G*- (0.25) predictions is nevertheless considerably better than suggested on the basis of monomer polarizability calculations. Fast and very accurate estimates of the MP2 aromatic stacking energies can be achieved using the RI-MP2 method. The CBS MP2 calculations and the CCSD(T) correction, when taken together, bring only marginal changes to the relative stability of H-bonded and stacked base pairs, with a slight shift of ca. 1 kcal/mol in favor of H-bonding. We suggest that the present values are very close to ultimate predictions of the strength of aromatic base stacking of DNA and RNA bases. Introduction Stacking of aromatic systems plays an important role in nature and is responsible for the structure and dynamics of many complexes. The best-known example represents stacking of nucleic acid (NA) bases which fundamentally contributes to the stability and conformational variability of nucleic acids.1 The physicochemical origin of stacking differs considerably from nucleobase H-bonding. Whereas the H-bonding is mainly of electrostatic origin, the stacking interaction is due to the London dispersion energy.2-4 The electrostatic interactions are correctly described already at the Hartree-Fock (HF) level of quantum chemical description with rather small basis sets and that is why H-bonding was studied earlier and more extensively. The very popular density functional theory (DFT), including some por- tions of electron correlation energy, yields accurate data on structure, stabilization energy, and vibration frequencies of H-bonded complexes while the semiempirical methods are less satisfactory. In contrast, for base stacking HF, DFT and semiempirical methods fail completely.2 We nevertheless wish to underline that there have been considerable recent efforts to improve the DFT predictions for stacked molecular clusters though, to our opinion, ultimate success for stacking (quality comparable to conventional electron correlation methods and available to all stacked NA base pairs) has not yet been achieved.4b-g The only exception represents to our best knowl- edge our recent attempt4h where we combined the self- consistent-charge, density-functional tight-binding method with * Corresponding author. E-mail: hobza@indy.jh-inst.cas.cz. (1) (a) Mathews, D. H. Sabina, J. Zuker, M. Turner, D. H. J. Mol. Biol. 1999, 288, 911-940. (b) Bommarito, S. Peyret, N. SantaLucia, J. Nucleic Acids Res. 2000, 28, 1929-1934. (c) Suzuki, M. Amano, N. Kakinuma, J. Tateno, M. J. Mol. Biol. 1997, 274, 421-435. (2) (a) Hobza, P. Sponer, �� J. Chem. ReV. 1999, 99, 3247-3276. (b) Sponer, �� J. Leszczynski, J. Hobza, P. Biopolymers 2001, 61, 3-31. (3) (a) Sponer, �� J. Leszczynski, J. Hobza, P. J. Phys. Chem. 1996, 100, 5590- 5596. (b) Sponer, �� J. Gabb, H. A. Leszczynski, J. Hobza, P. Biophys. J. 1997, 73, 76-87. (c) Sponer, �� J. Leszczynski, J. Hobza, P. J. Phys. Chem. A 1997, 101, 9489-9495. (4) (a) Sponer, �� J. Leszczynski, J. Hobza, P. J. Comput. Chem. 1996, 17, 841- 850. (b) Adamo, C. Barone, V. J. Chem. Phys. 1998, 108, 664-675. (c) Wesolowski, T. A. Parisel, O. Ellinger, Y. Weber, J. J. Phys. Chem. A 1997, 101, 7818-7825. (d) Williams, H. L. Chabalowski, C. F. J. Phys. Chem. A 2001, 105, 646-659. (e) Kurita, N. Kobayashi, K. Comput. Chem. 2000, 24, 351-357. (f) Kurita, N. Araki, M. Nakao, K. Kobayashi, K. Chem. Phys. Lett. 1999, 313, 693-700. (g) Wesolowski, T. A. Morgantini, P. Y. Weber, J. J. Chem. Phys. 2002, 116, 6411-6421. (h) Elstner, M. Hobza, P. Frauenheim, T. Suhai, S. Kaxiras, E. J. Chem. Phys. 2001, 114, 5149-5155. Published on Web 09/05/2002 11802 9 J. AM. CHEM. SOC. 2002, 124, 11802-11808 10.1021/ja026759n CCC: $22.00 �� 2002 American Chemical Society

the empirical expression describing the London dispersion energy. This technique successfully described stacking in all 10 NA base pairs. While stabilization energies of stacked nucleobase dimers are typically around -10 kcal/mol, the HF and most DFT approaches give only repulsive potential.2,4a What makes the theoretical analysis of base-stacking inter- actions especially difficult is the fact that higher-order electron correlation effects substantially affect the strength of aromatic stacking.5 The higher-order contributions are not included in the second-order Moeller-Plesset (MP2) perturbational theory conventionally used to study molecular interactions.5 This is not the case of H-bonding where the higher-order electron correlation terms cancel each other and, therefore, the MP2 method is usually very accurate.5 Further, evaluation of reason- able stacking energies requires the use of diffuse functions and, especially, diffuse polarization functions. In preceding studies on NA base stacking, we used the MP2 method combined with the 6-31G*(0.25) basis set.2-4a Here, the standard d-polarization functions (with exponent of 0.8) were replaced by more diffuse ones (exponent of 0.25)6 with the aim to improve description of the dispersion attraction.2-4a This basis set is obviously not fully balanced and cannot be used, for example, for geometry optimizations. Nevertheless, the utilization of this basis set was crucial in early electron correlation studies of base stacking as the standard 6-31G* basis set would provide a highly distorted picture of stacking.2 Stacking energies of NA bases are sensitive to the quality of the basis set and very extended basis sets increase (in absolute value) the calculated stacking energies significantly. The final answer to this problem can be obtained by evaluating the complete basis set (CBS) stacking energies, and the first attempt was made recently when Nielsen et al.7 studied two gradiently optimized geometries of stacked uracil dimer. The authors obtained CBS MP2 binding energies from basis set extrapola- tions using cc-pVXZ and aug-cc-pVXZ basis sets. Higher-order correlation contributions were obtained from the difference of MP2 and CCSD(T) stacking energies determined with the 6-31G*(0.25) basis set. CCSD(T) stands for the coupled cluster method with noniterative evaluation of triple excitations. The final estimate of the stacking energy was by 1.2 and 1.7 kcal/ mol higher (in absolute values) compared with MP2/6-31G*- (0.25) calculations.7 In the present paper, we reevaluate the stacking interaction of smallest NA base pairs (uracil and cytosine homodimer) and a related model system (pyrimidine dimer) using the MP2 procedure combined with extended basis sets. Higher-order correlation contributions are obtained with the CCSD(T) procedure using basis sets larger than the 6-31G*(0.25) one. We are aware of the fact that the higher-level correlation treatment like CCSD(T) requires the use of extended basis sets containing higher polarization functions such a calculation for the present NA base pairs is, however, clearly impractical. Using a small (symmetrical) model system, we demonstrate that reasonable values of the difference between MP2 and CCSD- (T) interaction energies can be obtained already using smaller basis sets. Further, an efficient procedure for estimation of genuine stacking energies of larger NA dimers (guanine dimer) will be suggested. We will also demonstrate the ability of resolution of identity MP2 (RI-MP2) procedure8a-c to correctly describe stacking interactions. It was suggested in the past that the RI-MP2 correlation interaction energy might be smaller than that in the accurate MP2 method. Our recent calculations using the RI-MP2 method for NA base pairs indicated an excellent performance of the RI-MP2 method and it is confirmed here using a more systematic comparison.8d Accurate characterization of nucleobase stacking allows proper calibration of the balance between stacking and hydrogen bonding of nucleic acid bases. This information is important for studies of a wide variety of systems ranging from advanced gas-phase physicochemical experiments up to the condensed- phase experiments and simulations.9 Advanced ab initio calcula- tions of base stacking are vital for parametrization of molecular mechanics (empirical) potentials and other computational techniques since gas-phase experiments on the energetics of nucleobase association are presently not available. We suggest that the present values of base stacking are very close to ultimate predictions of the strength of aromatic base stacking of DNA and RNA bases. Methods Geometries. Four aromatic stacked structures were investigated (Figure 1). For pyrimidine and cytosine homodimers, we used anti- parallel undisplaced face-to-back dimers with a vertical separation between the coplanar monomers of 3.3 ��, assuming rigid monomers.3a,5 (5) Sponer, �� J. Hobza, P. Chem. Phys. Lett. 1997, 267, 263-270. (6) Kroon-Batenburg, L. M. J. van Duijneveldt, F. B. J. Mol. Struct. THEOCHEM 1985, 121, 185-199. (7) Leininger, M. L. Nielsen, I. M. B. Colvin, M. E. Janssen, C. L. J. Phys. Chem. A 2002, 106, 3850-3854. (8) (a) Feyereisen, M. Fitzgerald, G. Komornicki, A. Chem. Phys. Lett. 1993, 208, 359-363. (b) Vahtras, O. Almlo ��f, J. Feyereisen, M. W. Chem. Phys. Lett. 1993, 213, 514-518. (c) Bernholdt, D. E. Harrison, R. J. Chem. Phys. Lett. 1996, 250, 477-484. (d) Jurec ��ka, P. Nachtigall, P. Hobza, P. Phys. Chem. Chem. Phys. 2001, 3, 4578-4582. (9) (a) Nir, E. Janzen, C. Imhof, P. Kleinermanns, K. de Vries, M. S. Phys. Chem. Chem. Phys. 2002, 4, 732-739. (b) Nir, E. Janzen, C. Imhof, P. Kleinermanns, K. de Vries, M. S. Phys. Chem. Chem. Phys. 2002, 4, 740- 750. (c) Desfrancois, C. Abdoul-Carime, H. Schultz, L. P. Schermann, J. P. Science 1995, 269, 1707-1709. (d) Schnier, P. D. Klassen, J. S. Strittmatter, E. E. Williams, E. R. J. Am. Chem. Soc. 1998, 120, 9605- 9613. (e) Dunger, A. Limbach, H.-H. Weisz, K. J. Am. Chem. Soc. 2000, 122, 10109-10114. Figure 1. Structures of pyrimidine dimer (1), cytosine dimer (2), uracil dimer (3), and guanine dimer (4). Toward True DNA Base-Stacking Energies A R T I C L E S J. AM. CHEM. SOC. 9 VOL. 124, NO. 39, 2002 11803