Introduction

Two major factors for the structure and stability of DNA are the hydrogen bonds between two adjacent bases from opposite strands (forming a Watson-Crick pair), and the π-π stacking interactions between two bases within each of the two DNA strands. The hydrogen bonding interactions in DNA are well understood by high-level Density Functional Theory (DFT) calculations [1–3]. They were shown to be provided by nearly equally important contributions from electrostatic attraction and donor-acceptor orbital interactions (e.g., from an N lone pair to an N-H σ* orbital), and to be sensitive to effects from the molecular environment (counter-ions, solvent) [1, 2]. A similar understanding of the π-stacking interactions is currently lacking due to several reasons: (1) π-π stacking interactions have proven to be sensitive to the methodology used (see below); (2) π-π stacking interactions have proven to be very dependent on the size of the basis set used in (ab initio) calculations, resulting in a large basis-set superposition error (BSSE); (3) high-level ab initio benchmark studies (CCSD(T) with at least triple-zeta basis sets) have only recently become available; and (4) the notorious problems of DFT to correctly describe dispersion interactions, [4–16] which constitute one of the components involved in π-stacking interactions (see below). The latter problems of DFT resulted in potential energy surfaces (PES) that, erroneously, lack any equilibrium configuration, i.e., the stacked systems were incorrectly described as being purely repulsive.

Because of the technical difficulties involved in obtaining accurate benchmarks (CCSD(T) using large basis sets), Sponer, Hobza and co-workers have studied in a ground-breaking series of papers [17–22] the stacking energies of DNA bases and DNA base analogs using a variety of ab initio methods and basis sets. Earlier studies showed that reasonable interaction energies could be obtained when employing the MP2/6-31G*(0.25) method, which was confirmed by later and more accurate studies [17–22]; in the aforementioned 6-31G*(0.25) basis set, the standard d-polarization functions (with exponent of 0.80) were replaced by more diffuse ones (exponent of 0.25) with the aim to improve description of the dispersion attraction [22]. Thus, for the anti-parallel displaced pyrimidine dimer, the MP2/6-31G*(0.25) energy is only 1.1 kcal mol−1 lower than the CCSD(T)/6-31++G*(0.25,0.15) energy, for the anti-parallel cytosine dimer it is 0.6 kcal mol−1 higher than the CCSD(T)/6-31++G**(0.25,0.15) energy, and for the twisted uracil dimer the MP2/6-31G*(0.25) is equal to the CCSD(T)/6-31++G**(0.25,0.15) energy. Although there seems to be no systematic shift in energy when comparing MP2/6-31G*(0.25) with the more elaborate CCSD(T) stacking energies, the MP2 stacking energies are always larger in magnitude than the corresponding CCSD(T) ones when using the same basis set for both methods. Interestingly, in connection with the quest for reliable ab initio benchmarks, preliminary results [23–27] indicate that the recently developed spin-component scaled MP2 (SCS-MP2) yields systematically lower stacking energies than MP2. This suggests that SCS-MP2 is in good agreement with the CCSD(T) results at only a fraction of the computational cost.

Several approaches have been proposed to remedy the “dispersion problem” of DFT, ranging from Time-Dependent DFT frameworks [28] to more ad-hoc procedures such as including (C6) van der Waals terms in the interaction energy (DFT-D) [29] or fitting (reparameterizing) the form of the density functional to some reference energies and/or distances. A more elegant solution was recently proposed by Tozer and co-workers [30], who analyzed the DFT exchange-correlation (xc) potential in comparison with the near-exact multiplicative Zhang-Morrison-Parr (ZMP) potential from coupled cluster Brueckner doubles (BD) densities. Their resulting KT1, and the related KT2, functional [30] produces xc potentials that closely resemble those from ZMP. As a result, KT1 and KT2 produce NMR chemical shifts [4, 31] that in general show smaller deviations from experimental values than those from standard density functionals. Furthermore, the ZMP potential from Brueckner (BD(T)) coupled cluster densities gives dispersion forces in good agreement with near-exact dispersion forces [32]. It is therefore to be expected that the KT1 and KT2 functionals behave properly for dispersion forces as well. However, the application of these new functionals to π-stacking interactions has so far not appeared in the literature.

Here, we report the results of a systematic investigation of π-stacking interaction energies, potential energy surfaces (PES) and stacking distances for dimers of, among others, the five DNA/RNA bases. A long list of standard as well as recent density functionals (including the KT1 and KT2 functionals) has been considered in this study, ranging from the Local Density Approximation (LDA), via the Generalized Gradient Approximation (GGA), and meta-GGA, to hybrid functionals. We have also examined BHandH, which was recently [12] shown to give reasonable agreement with ab initio data. The DFT results are compared to benchmark data that were taken from the literature, which in some cases contained energy profiles as a function of the torsion associated with the rotation of one base with respect to the other.

The present investigation is part of a larger undertaking. We aim at developing an extended QM/QM scheme, designated QUILD [33], which allows a (complex) molecular system to be decomposed into arbitrary, interpenetrating domains such that each type of interaction can be described with a different method, i.e., the density functional that performs best for that type of interaction or process, for example, electron-pair as well as hydrogen bonding (BP86: see [2, 3]), π-π stacking (LDA, KT1, KT2: this work) or reaction barriers (OLYP, OPBE: see [34, 35]). This is a pragmatic (of course not a fundamental) approach to cope with the present-day shortcomings of DFT, which features density functionals that perform satisfactorily for a number but not all types of interactions and phenomena. The present study focuses on identifying which DFT approaches serve best to describe the π-π stacking interactions in the above-mentioned QM/QM scheme.

Furthermore, we have analyzed the π-stacking interactions in more detail by decomposing the interaction energy into three physically meaningful components, namely, electrostatic energy, Pauli repulsion (i.e., the destabilizing interactions between occupied orbitals), and bonding orbital interactions [36]. Our analyses reveal, among others, that the electrostatic attraction between the stacked molecules is the most important component of the π-π stacking interaction (see also [37]) which determines the shape and depth of the PES. We have also analyzed the π-π interactions between two stacked Watson–Crick base pairs with the purpose to differentiate between inter- and intra-strand stacking interactions.

Methods

Computational details

All calculations were performed with the Amsterdam Density Functional (ADF) program. [38–49] Molecular orbitals (MOs) were expanded using a large, uncontracted set of Slater-type orbitals: TZ2P [40]. The TZ2P basis is an all-electron basis of triple-ζ quality, augmented by two sets of polarization functions (2p and 3d on H; d and f on heavy atoms). An auxiliary set of s, p, d, f, and g STOs was used to fit the molecular density and to represent the Coulomb and exchange potentials accurately in each SCF cycle.

The equilibrium structure of individual benzene analogs, DNA bases and Watson-Crick base pairs were optimized at the BP86 level of DFT, i.e., using the local density approximation (LDA; Slater exchange and VWN [42] correlation) with non-local corrections due to Becke [43] (exchange) and Perdew [44] (correlation) added self-consistently, which was previously shown to be one of the best DFT methods for the accuracy of geometries [41, 45] and hydrogen-bonding interactions [1–3]. The π-π stacking interaction energies for all density functionals were evaluated post-SCF using LDA as SCF functional. This procedure was recently shown to introduce an error of only a few tenths of a kcal mol−1 [46].

Stacking distances (vertical separation) and orientations (twist angle) were explored with the various density functionals through scans of the potential energy surface (PES) in which the BP86 geometries of the monomers (e.g., DNA bases or Watson–Crick base pairs) were kept frozen. PES scans as a function of the twist angle (see below) were done using steps of 30° in case of homo-base stacks, 60° in the case of hetero-base stacks and 36° for stacks of Watson–Crick base pairs. PES scans as a function of the vertical separation (see below) were done in steps of 0.1 Å.

The interaction energy (see below) was studied as a function of the vertical separation between the two monomers in the π-stacked dimer, and as a function of the twist angle between these monomers. The latter is defined as the right-handed rotation of the upper base around the axis that runs through the center of mass of both the upper and lower base, and perpendicular to the plane of the base (see Scheme 1).

Scheme 1
scheme 1

Definition of the twist angle (TW). The black dot is the center of mass

Full size image

In the stacked DNA systems, the twist angle of 0° is defined as that twist angle at which the respective “glycosidic” N-H bonds (more precisely, the N-H bonds that replace the glycosidic N-C bonds to the sugar in the backbone) are oriented in parallel.

Bond-energy decomposition

The overall π-stacking energy ΔE is made up of two major components [36] (Equation 1):

$$\Delta E = \Delta E_{{{\text{prep}}}} + \Delta E_{{\operatorname{int} }} $$
(1)

In this formula, the preparation energy ΔE prep is the energy needed to deform the separate molecular fragments from their equilibrium structure to the geometry that they attain in the overall molecular system. The interaction energy ΔE int is the energy released when the prepared fragments are brought together into the position they have in the overall molecule. It is analyzed for our model systems in the framework of the KS-MO model using a Morokuma-type decomposition [47] into electrostatic interaction, Pauli repulsion (or exchange repulsion), and (attractive) orbital interactions (Equation 2).

$$\Delta E_{{\operatorname{int} }} = \Delta V_{{{\text{elstat}}}} + \Delta E_{{{\text{Pauli}}}} + \Delta E_{{{\text{oi}}}} $$
(2)

The term ΔV elstat corresponds to the classical electrostatic interaction between the unperturbed charge distributions of the prepared (i.e. deformed) fragments and is usually attractive. The Pauli-repulsion, ΔE Pauli, comprises the destabilizing interactions between occupied orbitals and is responsible for the steric repulsion. The orbital interaction ΔE oi in any MO model, and therefore also in Kohn-Sham theory, accounts for electron-pair bonding, charge transfer (i.e., donor-acceptor interactions between occupied orbitals on one fragment with unoccupied orbitals of the other, including the HOMO-LUMO interactions), and polarization (empty-occupied orbital mixing on one fragment due to the presence of another fragment). The orbital interaction energy can be further decomposed into the contributions from each irreducible representation ℘ of the interacting system (Equation 3) using the extended transition state (ETS) scheme developed by Ziegler and Rauk [48, 49].

$$\Delta E_{{{\text{oi}}}} = \sum _{\Gamma } {\text{ }}\Delta E_{\Gamma } $$
(3)

Results and discussion

We have analyzed π-π stacking interactions using different systems that were taken mainly from reference data present in the literature. We start with the prototypical system for studying π-π stacking, namely, the benzene dimer. The second set of reference data was taken from the supporting information of a paper by Mignon and co-workers [50], who studied hydrogen bonding and stacking in complexes of substituted benzene C6H5X with pyridine at MP2/6-31G*(0.25) with substituents X = H, F, NH2, Cl, CH3, OH, CN, COOH, CHO and NO2. All substitutions of H by another X in benzene were done at the para-position with respect to the nitrogen of pyridine (see Scheme 2). The third set of reference data was taken from a paper by Jurecka and co-workers [21], who investigated hydrogen bonding and stacking interactions in the cytosine dimer with MP2 and CCSD(T) methods, using several orientations (see Scheme 3) to investigate the potential energy surfaces. The fourth set of reference data was taken from a paper by Wu and Yang [14] (based on previous work by Hobza, Sponer and co-workers [17, 19, 22]) who studied the dependence of the stacking energy of DNA bases on the twist angle (see Scheme 4) and vertical separation using MP2/6-31G*(0.25) calculations.

Scheme 2
scheme 2

Stack of substituted benzene and pyridine

Full size image
Scheme 3
scheme 3

Orientations of the stacked cytosine dimer from Jurecka and co-workers ([21]). The upper base is drawn in bold, the center of mass is indicated by a black dot, and the carbonylic oxygen and the hydrogen replacing the glycosidic bond are indicated for clarity

Full size image
Scheme 4
scheme 4

Geometries of stacked base dimers with a twist angle of 0 and 120 degrees. The upper base is drawn in bold, the center of mass is indicated by a black dot, and the hydrogens replacing the glycosidic bond are indicated for clarity. Note that for the homo dimers at zero twist angle the lower base is hidden behind the upper base

Full size image

As mentioned already, the prototypical system for studying π-π stacking interactions is posed by the (parallel sandwich) benzene dimer. Because of its high symmetry (D6h), CCSD(T) energies using large basis sets (aug-cc-pVQZ) are available [51] to compare with (see Table 1). The minimum is found at 3.9 Å for this system with a stacking energy of −1.7 kcal mol−1. The standard density functionals are unable to give a proper description of this system, i.e., they either predict a purely repulsive energy surface (e.g., BP86; see Table 1), or give a shallow well (e.g., OLYP). Interestingly, the KT1 functional predicts a minimum that is very close (3.8–3.9 Å) to that of the CCSD(T) method, and has a well depth that is equally close (−1.6 kcal mol−1; see Table 1). The related KT2 and, surprisingly, the LDA functionals have a minimum at the same distance as CCSD(T) and KT1, but with a somewhat smaller well depth (−1.3 kcal mol−1). Also, the BHandH functional performs reasonably well with a similar equilibrium distance, but a reduced well depth (−0.9 kcal mol−1; see Table 1).

Table 1 Stacking energy (in kcal mol−1) of the benzene-benzene complex as a function of the vertical separation Rvert (in Å), computed with various density functionalsa and CCSD(T)/aug-cc-pVQZ*b
Full size table

Reference data from Mignon and co-workers

Scheme 2 shows the structure of the π-stacked complexes of pyridine with (substituted) benzene, i.e. with benzene, fluorobenzene, aminobenzene, chlorobenzene, toluene, phenol, cyanobenzene, benzoic acid, benzaldehyde and nitrobenzene. For all of these complexes the coordinates were taken from the supporting information of the original paper by Mignon and co-workers [50]. The geometry of each dimer had been fully optimized and therefore the vertical separation is different for each complex; it is found in the range of 3.2–3.4 Å.

The interaction energies for several typical density functionals are compared to the original MP2 data in Table 2. Consistent with previous reports, the mean absolute deviations (MAD) with respect to the reference MP2 data are much larger for density functionals containing the Becke88 [43] exchange functional (e.g., BP86, MAD value 7.68 kcal mol−1) than for functionals based on PW91 [52] or PBE [53] exchange (e.g., PW91, MAD value 5.33 kcal mol−1). However, these standard density functionals all show repulsive interactions, in contrast to the MP2 results that indicate bound systems (see Table 2). Much better agreement with the MP2 data is observed for the KT1, KT2 and, surprisingly, LDA functionals, with MAD values of 0.25 kcal mol−1 (see Table 2). The BHandH functional, recently [12] shown to give good agreement with MP2/CCSD(T) energies, shows a larger MAD value of 0.4 kcal mol−1 (see Table 2). Because of the uncertainty of ca. ± 0.5 kcal mol−1 in the reference MP2 data with respect to the high-level CCSD(T) data (see above), it is uncertain which of the four density functionals (LDA, KT1 , KT2 or BHandH) performs best. However, the improvement over standard density functionals is obvious and promising.

Table 2 Stacking energy (in kcal mol−1) of substituted benzene-pyridine complexes (see Scheme 2), computed with various density functionalsa and MP2/6-31G*(0.25)b
Full size table

Reference data from Jurecka and co-workers

Scheme 3 shows the fourteen orientations of the stacked cytosine dimer that have been considered in the paper by Jurecka and co-workers [21] (with a vertical separation of 3.3 Å for the first orientation, and 3.4 Å for all others). The corresponding reference energies and our DFT interaction energies are collected in Table 3. The latter also contains the MAD values with respect to two sets of reference energies, i.e. the usual MP2/6-31G*(0.25) results (MAD1) as well as the MP2 energies obtained after complete-basis-set (CBS) extrapolation corrected for the difference between MP2 and CCSD(T) with the 6-31G*(0.25) basis set (MAD2). Similar to what was found for the benzene-pyridine systems (see above), the MAD values for standard density functionals are found between 5 and 14 kcal mol−1, resulting in many cases in erroneous repulsive stacking interactions between the cytosine fragments (see Table 3). Improved results are again obtained with LDA, KT1, KT2, and to a lesser extent BHandH, which have MAD1 values (i.e. deviations with respect to MP2/6-31G*(0.25)) of 0.9, 0.8, 0.6 and 1.5 kcal mol−1 respectively, and MAD2 values of 0.4, 0.5, 0.7 and 0.5 kcal mol−1 respectively with the CCSD(T) data (see Table 3).

Table 3 Stacking energies (in kcal mol−1) for several orientations of the cytosine dimer (see Scheme 3), computed with various density functionals and MP2a,b
Full size table

Reference data from Wu and Yang

In the paper by Wu and Yang, [14] the potential energy surface (PES) was calculated as a function of the twist angle between two stacked bases (see Scheme 1 and methodological section for the definition of the twist angle) using MP2/6-31G*(0.25), which was chosen as a reliable reference method based on the paper by Hobza and co-workers [22]. Although thymine was not included by Wu and Yang, we report here the energy profiles for all combinations of the five RNA/DNA bases, i.e., thymine, adenine, guanine, cytosine, uracil (see Scheme 4 for the geometries). The vertical separation is 3.4 Å for the C-C, G-C and G-G stacked base dimers, and 3.3 Å for all other systems. This difference in vertical separations was chosen based on earlier papers by Sponer and co-workers [19].

From the tests with the benzene dimer (Table 1), the Mignon (Table 2) and Jurecka reference data (Table 3), it is already evident that standard density functionals do not give reliable interaction energies, which is now shown to be true also for the energy profiles as function of the twist angle. Figure 1 shows the PES for two stacked bases, namely C-C and G-U, which are representative for all possible base pair combinations (the complete figure for all combinations of stacked bases can be found in the Supporting Information). Interestingly, although the standard density functionals typically underestimate the interaction energy, the shape of the PES is in all cases highly similar to that based on the more accurate MP2 data. This suggests that, at least trends and qualitative features in the rotational energy profile can be correctly reproduced with DFT which otherwise shows a (functional and basis-set dependent) constant shift with respect to the more accurate MP2 rotational profile.

Fig. 1
figure 1

Stacking interaction (in kcal mol−1) as function of the twist angle (in degree) for π-π stacks of two bases, C-C and G-U, computed with various density functionals. The thick lines with filled circles, triangles, squares and diamonds show the results obtained with LDA, PW91, BLYP, and BP86, respectively. The dashed lines with open squares, triangles and diamonds show the results obtained with OLYP, B3LYP and X3LYP, respectively. The bold line with crosses represents the MP2/6-31G*(0.25) reference data

Full size image

It is, however, more instructive to look at functionals that by themselves already give accurate interaction energies, such as LDA, KT1, KT2 and BHandH. Figure 2 shows the corresponding rotational-energy profiles, again for the stacked bases C-C and G-U. From Fig. 2 it is evident that the LDA, KT1, KT2, and to a lesser extent BHandH are visually (and virtually) indistinguishable from the reference MP2 data.

Fig. 2
figure 2

Stacking interaction (in kcal mol−1) as function of the twist angle (in degree) for π-π stacks of two bases, C-C and G-U, computed with various density functionals. The red, blue and green lines show the results obtained with BHandH, KT1 and KT2, respectively. The bold line with crosses represents the MP2/6-31G*(0.25) reference data

Full size image

Figure 3 shows the interaction energy as a function of the vertical separation for both the standard (upper diagrams) and promising (lower diagrams) density functionals, together with MP2 reference curve in each of the four diagrams. Not surprisingly, the density functionals that gave larger errors for the rotational energy profile, also give larger errors for the vertical separation profiles: almost all standard density functionals show potential energy surfaces that are repulsive, or at best weakly attractive with a very shallow minimum at elongated distances. In contrast, the promising functionals (i.e., KT1, KT2, LDA and BHandH) do provide well-defined minima with equilibrium distances (3.2–3.3 Å) close to the reference MP2 data. Similar to the uncertainty in the reference energy data of some 0.5–1.0 kcal mol−1 (see above), we may expect a similar uncertainty for the vertical separation, which we very roughly estimate at 0.1–0.2 Å. The differences between the equilibrium distance of the promising density functionals and that of the reference MP2 falls well within this estimated uncertainty. Thus, KT1, KT2, LDA and BHandH perform well not only for rotational energy profiles but also for the PES as a function of the vertical separation.

Fig. 3
figure 3

Stacking interaction (in kcal mol−1) as function of the vertical separation (in Å) for π-π stacks of two bases, A-A and G-G, computed with various density functionals. The thick lines with filled triangles, squares and diamonds in the upper diagrams show the results obtained with PW91, BLYP, and BP86, respectively. The dashed lines with open squares, triangles and diamonds in the upper diagrams show the results obtained with OLYP, B3LYP and X3LYP, respectively. The red, blue, green and orange lines in the lower diagrams that truly show well-defined minima are obtained with BHandH, KT1, KT2 and LDA respectively. The bold black line with crosses represents the MP2/6-31G*(0.25) reference data

Full size image

Decomposition of interaction energy

Now that we have established the reliability of the KT1, KT2, LDA and BHandH functionals for the interaction energy of stacked DNA bases, it is interesting to study the actual interaction into greater detail. This is done through a quantitative decomposition of the interaction energy ΔE int into electrostatic attraction ΔV elstat, Pauli repulsion ΔE Pauli and orbital interactions ΔE oi (see Methods). Table 4 provides the results of the energy decomposition at KT1/TZ2P for the five π-stacked homo dimers, for each one at a twist angle of 0° and 180° (which is the lowest-energy conformation) and at two different vertical separations. Note that the lowest-energy conformation is always achieved at a twist angle of 180° but that the optimal stacking distance (vertical separation) varies by 0.1 – 0.2 Å between the different π-stacked homo dimers. The energy changes by −6.8 (A-A), −11.8 (C-C), −10.2 (G–G), −11.2 (T-T) and −7.9 (U-U) kcal mol−1, respectively, when going from a twist angle of 0° to 180° (the lowest-energy conformation). This change in energy is almost completely resulting from a change in electrostatics that change by −5.4 (A-A), −13.1 (C-C), −10.3 (G-G), −6.6 (T-T) and −8.3 (U-U) kcal mol−1 respectively. The sum of Pauli repulsion and orbital interactions changes much less, by ca. 1.4 kcal mol−1 or less, except for T-T where the sum of Pauli repulsion and orbital interactions changes by 4.6 kcal mol−1 due to the influence of the peripheral methyl groups. But also the changes in each term individually, i.e., either ΔE Pauli or ΔE oi, are much smaller (typically 2 to 5 times, but in some cases even more) than the changes in ΔV elstat (again, only in the case of T-T, the change in ΔE Pauli is nearly as large as that in ΔV elstat). Here we note that the same energy decomposition analyses for the five π-stacked homo dimers at LDA/TZ2P instead of KT1/TZ2P yield values that differ by only 1 kcal mol−1 or less (not shown in Table 4).

Table 4 Decomposition of interaction energiesa (in kcal mol−1) between stacked DNA bases in different geometries, computed at KT1/TZ2P
Full size table

Additivity of interaction energies of stacked base pair dimers

With the good performance of KT1 for stacking interactions kept in mind, we went one step further, and looked also at the interactions within stacked dimers of hydrogen-bonded base pairs (see Scheme 5) [54]. The geometry was obtained by taking the geometry of one hydrogen-bonded base pair and putting another one at 3.4 Å on top of the former. The base pair dimer is stabilized by hydrogen bonding between the bases within one layer, and stacking interactions between the bases in different layers of the stack. In this case, because there are two bases per layer, there are also interactions present between bases that are, in a sense diagonally, on opposite sides and in different layers [55], the so-called cross terms. For example for the AT-AT system as shown in Scheme 5, this refers to the interaction between the thymine of the top layer with the adenine of the bottom layer.

Scheme 5
scheme 5

Structure of stacked DNA base pair dimer AT-AT with twist angle of 36°. The dot represents the B-DNA helical axis of rotation.[54]

Full size image

In order to investigate in more detail the importance of the various energy terms, in particular, the cross-terms, we approximate the total stacking interaction between the base pairs in the stacked base pair dimer by the sum of four pairwise interactions between two bases, i.e., we approximate the total (stacking) interaction energy ΔE int between the layers AB and CD in Scheme 6 by the sum ΔE add of interactions A-C, B-D, A-D and B-C as shown in Equations 4a,b:

$$\Delta E_{{{\text{add}}}} = \Delta E_{{{\text{AC}}}} + \Delta E_{{{\text{BD}}}} + \Delta E_{{{\text{AD}}}} + \Delta E_{{{\text{BD}}}} $$
(4a)
$$\Delta E_{{\operatorname{int} }} = \Delta E_{{{\text{add}}}} + \Delta E_{{{\text{coop}}}} $$
(4b)
Scheme 6
scheme 6

Additivity approximation for the π-π interaction between two stacked Watson-Crick base pairs in terms of pairwise interactions between individual bases

Full size image

The remaining term ΔE coop that, added to the approximate ΔE add yields again the real stacking interaction ΔE int, consists of cooperativity effects between the different interactions.

In Table 5, we report the energy values resulting from the additivity scheme for the stacked AT-AT, GC-GC and AU-AU base pair dimers at three different values of the twist angle (0, 36 and 180°). Importantly, the additive approximation ΔE add of the real stacking interaction ΔE int largely accounts for the latter. ΔE coop has in most cases a positive sign and is small compared to either ΔE add or ΔE int. Furthermore, the dependence of the ΔE coop term on the twist angle is much smaller than that of the interaction energy ΔE int. For example, ΔE coop varies from 1.0 to 1.6 kcal mol−1 whereas ΔE int varies from −1.4 to −10.7 kcal mol−1 for the stacked AT-AT base pair dimer (see Table 5). Larger cooperativity effects are found in the case of the stacked GC-GC base pair dimer, where ΔE coop values are found between 2.0 and 4.4 kcal mol−1. Note however that the changes in ΔE coop are again small compared to the changes in ΔE int.

Table 5 Interaction energy (kcal mol−1) between stacked DNA base pairs in terms of the additivity scheme, computed at KT1/TZ2Pa,b
Full size table

The cross-terms in the additivity scheme are as important as the non-cross terms, especially when rotating the base pair dimer from a twist angle of 0 degrees to finally a twist angle of 180 degrees. At the end point of 180°, the A and B base have switched, and the cross terms are no longer AD and BC, but AC and BD; this switch occurs after a twist angle of 90°. For instance for the AT-AT base pair dimer at 0°, the cross terms are stabilizing (−2.5 kcal mol−1) and dominate the ΔEadd of −2.2 kcal mol−1. At 180°, the non-cross terms dominate (−10.9 kcal mol−1) with still a substantial contribution from the cross terms (−1.4 kcal mol−1). Note also the dominance of the cross terms (−10.2 kcal mol−1) over the non-cross terms (−3.0 kcal mol−1) for the GC-GC base pair dimer at 180°. This is mainly due to favorable hydrogen-bonding interactions between the cross G-G pair.

When rotating base pair dimers from a twist angle of 0° to 36°, the cross terms in the additivity scheme reduce slightly. However, this is dominated by the increase of favorable interactions between the stacked bases. For instance for AT-AT, the cross terms drop from −2.5 to −2.0 kcal mol−1, but the non-cross terms increase from +0.0 to −9.2 kcal mol−1. The same trend is observed for the AU-AU and GC-GC base pair dimers.

Again we note that the same energy decomposition analyses for stacked DNA base pairs at LDA/TZ2P instead of KT1/TZ2P yield values that differ by only 1 kcal mol−1 or less (not shown in Table 5).

Significance of the observed trends for the structure of DNA

The stabilizing interactions for the structure of DNA mainly originate from hydrogen-bonding interactions between bases in one plane and the π-π stacking interactions between stacked bases. However, we have shown here that in addition to these terms there is also an important contribution of the cross terms that arise between bases that are neither hydrogen-bonded within one layer, nor stacked above each other. All these interaction energies can now be understood and quantified with density functional theory using appropriate functionals. For the hydrogen-bonding interactions, the BP86 and PW91 functionals seem to perform best, [3] while for π-stacking interactions KT1, KT2 and (surprisingly) LDA work well.

The experimentally observed value for the twist angle (36°) is retrieved in our study on base pair dimers as the orientation with optimal stacking interactions. We have shown that this is predominantly determined by the classical electrostatic component of the stacking energy. The rotation to a twist angle of 36° destroys part of the favorable cross terms (see above), which is however completely overcome by the increase in favorable stacking interactions between the stacked bases. The twist angle and vertical separation that are experimentally observed are of course influenced by both the sugar-phosphate backbone and the presence of solvent and/or counter ions. These factors will be tackled in forthcoming work.

Conclusions

We have analyzed π-π stacking interactions between two benzenes or benzene analogs, between two DNA bases, and between two Watson-Crick base pairs using Density Functional Theory (DFT) in combination with large basis sets. The interaction energies for a large number of density functionals have been compared with ab initio reference data. In line with previous studies, most standard density functionals recover, at best, only part of the favorable stacking interactions.

An exception is the new KT1 functional, which has been constructed to closely resemble the near-exact Zhao–Morrison–Parr (ZMP) exchange-correlation potential. Thus, not unexpectedly, KT1 gives a good description for π-π stacking energies and potential energy surfaces (PES) resulting in bound systems. Surprisingly, however, a similarly good performance is achieved with the Local Density Approximation (LDA). In view of the accurate prediction of PESes and the low computational cost, we recommend the use of either KT1 or LDA for treating the π-π stacking interactions in larger systems, using a QM/QM approach (e.g., our QUILD approach which is under development [33]) that allows this type of interaction to be separated from other types of interaction, such as hydrogen bonding but also regular covalent bonds.

To gain insight into the origin of π-π stacking interactions, we have decomposed the interaction energies into the classical electrostatic attraction, Pauli repulsion and orbital interactions. Interestingly, the electrostatic interactions appear to be the most important factor that determines the shape and depth of the PES. In the case of two stacked Watson-Crick base pairs, this classical electrostatic attraction causes a minimum to occur along the energy profile at a twist angle of 36°. Furthermore, the stabilizing contributions to the stacking interaction between two Watson-Crick base pairs is shown to originate from the inter-strand stacking terms, that is, from the interaction between two bases that are in different Watson-Crick pairs and also not directly stacked on top of each other.

The above-mentioned insight that electrostatic attraction plays an essential role in π-π stacking interactions nicely consolidates and extends the recent finding by Krapp, Bickelhaupt and Frenking [56] that classical electrostatic attraction is essential for understanding trends in bond strength of diatomic molecules. It is interesting to note that qualitative models of both “regular” chemical bonds and those of π-π stacking interactions often completely ignore classical electrostatic attraction. These qualitative models attempt to rationalize trends in bonding entirely in terms of orbital interactions (in the case of regular bonds) and dispersion interactions (in the case of π-π stacking).