1 Introduction

Besides the well-known Watson–Crick pairing of nucleic acid bases in the helical DNA molecule, these bases are known in telomeres to form cyclic hydrogen-bonded quartets [1–11] (see Scheme 1), which, in turn, are able to stack together through π-stacking interactions. These quartet structures of DNA bases are increasingly recognized for their biological importance [1–11].

Scheme 1
scheme 1

Guanine and adenine quartets

Full size image

By far the longest-known is the guanine quartet (G4), which has an exceptional stability due to its capacity to form two hydrogen bonds between neighboring guanines, i.e., the N1–H···O6 and N2–H···N7 hydrogen bonds (see Scheme 1). Furthermore, the stabilization of guanine quartets can be enhanced by cations [2]. This high stability of G4 facilitates through stacking the formation of other quartets, for instance of A4 [4–7]. However, the experimental evidence for the existence of adenine quartets is rare. To our knowledge, adenine quartets have only been observed experimentally sandwiched between two guanine quartets in [4–7]. The adenine quartets are formed with one hydrogen bond between adjacent bases (see Scheme 1). The three quartets have as proton-donor the N6–H bond of the amino-group and as proton-acceptor the N1 atom (A4-N1) [4, 5], the N3 atom (A4-N3) [6], and the N7 atom (A4-N7) [7], respectively. The guanine quartet, as well as the three adenine quartets, has been the subject of computational studies [12–19].

In this work, we analyze the stability of these four quartets in the gas phase and in aqueous solution using dispersion-corrected density functional theory (DFT-D). First, we explore the performance of different DFT-D variants for the hydrogen-bonded (see Scheme 2) and stacked DNA base pairs (see Fig. 1, later on) by comparing them with high-level ab initio results [20–22]. In our previous work, we showed that the hydrogen bonds between adjacent bases from opposite strands in DNA are described adequately by a number of GGA, meta-GGA and hybrid DFT approaches (for instance BP86) [23, 24]. However, most of the current DFT approaches, ranging from GGA via meta-GGA to hybrid DFT, fail in describing π-stacking interactions between two bases within each of the two DNA strands [25, 26].

Scheme 2
scheme 2

Watson–Crick base pairs

Full size image
Fig. 1
figure 1

Stacked AT and GC dimers in the gas phase shown from two side views (BLYP-D/TZ2P)

Full size image

After having established the proper dispersion-corrected functional, we investigate for the four quartets three different configurations which are of C4h, C4 and S4 point group symmetry and which all have the cyclic hydrogen bonds as shown in Scheme 1. The adenine quartets have been observed experimentally in a planar configuration [4–7], and we examine whether these quartets remain in a planar configuration when the stabilization by the guanine quartets is not present. Finally, we test the applicability of the widely used B3LYP functional for these quartets in the different symmetries and we show that the neglect of the dispersion correction will even lead to erroneous conclusions.

2 Computational methods

All calculations were performed using the Amsterdam density functional (ADF) program developed by Baerends et al. [27–38], and the QUantum-regions Interconnected by Local Descriptions (QUILD) program by Swart and Bickelhaupt [39, 40]. The QUILD program is a wrapper around ADF (and other programs) and is used for its superior geometry optimizer which is based on adapted delocalized coordinates [39]. The numerical integration was performed using the procedure developed by te Velde et al. [34, 35].

The MOs were expanded in a large uncontracted set of Slater type orbitals (STOs) containing diffuse functions: TZ2P (no Gaussian functions are involved) [36]. The basis set is of triple-ζ quality for all atoms and has been augmented with two sets of polarization functions, i.e., 3d and 4f on C, N, O and 2p, 3d on H. The 1s core shells of carbon, nitrogen and oxygen were treated by the frozen-core approximation [30]. 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 self-consistent field cycle [37].

Calculations were done using dispersion-corrected DFT-D as developed by Grimme [41–47] for a correct treatment of the stacking interactions between the DNA bases. The density functionals are augmented with an empirical correction for long-range dispersion effects, described by a sum of damped interatomic potentials of the form C6R−6 added to the usual DFT energy [41–47]. Equilibrium structures were optimized using analytical gradient techniques [38]. Geometries and energies were calculated with the generalized gradient approximation (GGA) at BLYP [48, 49], BP86 [48, 50] and PBE [51] as well as with the dispersion-corrected variants BLYP-D, BP86-D and PBE-D. Furthermore, the hybrid functional M06-2X [52] has been used for the computation of geometries and energies of AT and GC (both hydrogen-bonded and stacked dimers) and for the systems with stacking interactions involved. For the quartets the B3LYP [49, 53, 54] has been used to compute energies in a single-point fashion using the BLYP-D geometries.

At the BLYP-D level of theory, all energy minima of hydrogen-bonded AT and GC pairs, stacked AT and GC dimers, and the global energy minima of DNA-base quartets have been verified in the gas phase and in water to be equilibrium structures through vibrational analysis [55–57]. The lowest energy minima were found to have zero imaginary frequencies (see also Electronic Supplementary Material). For the dispersion-corrected functionals the basis set superposition error (BSSE) on the bond energy was not calculated because the dispersion correction [42] has been developed such that the small BSSE effects are absorbed into the empirical potential. For the bond energy calculated with the uncorrected functionals the BSSE has been computed and corrected through the counterpoise method [58], using the individual DNA bases as fragments.

Solvent effects in water have been estimated using the conductor-like screening model (COSMO) [59, 60], as implemented in the ADF program [61]. For settings see [62]. The continuum solvent model performs adequately for the determination of geometries as has been done in this work. However, for the calculation of spectra of molecules in a solvent, it is essential to take into account the first shell of solvent molecules explicitly [63–69]. For instance, Nicu et al. [63] showed that the VCD spectrum of 2-benzoic acid is influenced by the hydrogen bonding with the solvent due to the donor–acceptor interaction (which of course is not taken into account by continuum solvent models). According to the work by Riley et al. [70] the dispersion correction does not need to be modified for the solvated systems.

3 Results and discussion

3.1 AT and GC pairs

The hydrogen-bond distances and energies of the Watson–Crick base pairs AT and GC in the gas-phase and in aqueous solution are given in Tables 1 and 2, respectively, together with the best ab initio estimate by Sponer et al. [20]. The base pairs have been optimized in Cs symmetry and the bases in C1 in accordance with our previous work [23, 24]. The hydrogen bond energy for AT and GC is defined as the difference in energy between the optimized pair and the fully optimized bases. From Tables 1 and 2, we see that the hydrogen bond energies for AT and GC obtained with the dispersion-uncorrected density functionals, including the M06-2X functional, are all a few kilocalorie per mole more weakly bound than the CCSD(T) reference value with a clear underestimation of the hydrogen-bond energies by the B3LYP functional. The hydrogen-bond energies acquired with the dispersion-corrected density functionals are all slightly more strongly bound than the CCSD(T) reference value. Although the BP86-D and the PBE-D results already agree well with the ab initio results, the best agreement with CCSD(T) is obtained with BLYP-D. At BLYP-D/TZ2P, the hydrogen-bond energy for the Watson–Crick AT and GC pairs amounts to −16.7 and −30.1 kcal/mol, which has to be compared with −15.4 and −28.8 kcal/mol, respectively, at CCSD(T)/aug-cc-pVQZ//RI-MP2/cc-pVTZ [20–22]. This close agreement between dispersion-corrected functionals and ab intio results for the hydrogen-bonded and stacked AT and GC pair was also found by Grimme and co-workers [43] with Gaussian-type triple and quadruple-zeta basis sets. Furthermore, the inclusion of water as a solvent has been investigated. The hydrogen bonds between the AT and GC in water are half as strong as the same hydrogen bonds in the gas phase.

Table 1 Hydrogen bond distances (Å) and bond energies (kcal/mol) for AT computed at various levels of theory
Full size table
Table 2 Hydrogen bond distances (Å) and bond energies (kcal/mol) for GC computed at various levels of theory
Full size table

The π-stacking interaction in the AT and GC dimers in the gas phase and in aqueous solution are given in Table 3, together with the best ab initio estimate by Sponer et al. [20]. For the stacked AT and GC pairs, the dispersion corrected functionals are in close agreement with the ab initio results. The M06-2X energies are about 2 kcal/mol too weakly bound, which has also been found by Kabelac et al. [71, 72] for the M05-2X functional. The BLYP-D energies of −11.7 and −16.9 kcal/mol, respectively, are in excellent agreement with the CCSD(T)/aug-cc-pVQZ//RI-MP2/TZVPP ones of −11.6 and −16.9 kcal/mol, respectively [20–22].

Table 3 Stacking energies (kcal/mol) and distances between the bases (Å) for AT and GC computed at various levels of theory
Full size table

The B3LYP functional appears to be completely inadequate for describing stacked DNA bases, similar to what we have found previously for regular GGA functionals such as BP86 [25]. This follows from numerical experiments in which we evaluate the stacking energies in the gas phase and in water with B3LYP/TZ2P at the BLYP-D/TZ2P equilibrium structures, i.e., B3LYP/TZ2P//BLYP-D/TZ2P. In the gas phase, this yields for the stacked AT a slightly repulsive result of 2.1 kcal/mol. For GC the stacking energy in the gas phase amounts to −7.5 kcal/mol. The same numerical experiment has been done for PBE, which gives the AT pair bonded by only 1.4 kcal/mol and the GC pair by 10.1 kcal/mol at the PBE/TZ2P//BLYP-D/TZ2P level. Note that the reason for the attractive B3LYP and PBE energies of stacked GC is the presence of two hydrogen-bonding like interactions in the stacked GC system in which G and C are, in fact, not exactly parallel but under a slight angle. This is because of the aforementioned “partial” hydrogen bonds which B3LYP or PBE are able to describe. Stacked GC differs in this respect from stacked AT which is bound by a more pure π-stacking interaction which B3LYP or PBE is not able to describe properly (see Fig. 1). The stacked GC pair is therefore not a good benchmark to evaluate the adequacy of functionals for stacking interactions. However, starting from a stacked AT or GC pair the optimization with the B3LYP or PBE functional leads eventually to a hydrogen bonded base pair because of the erroneous description of the stacking interactions by these functionals. In water, both the AT and GC stacks would, according to B3LYP, not exist, which is evidently incorrect.

3.2 Adenine and guanine quartets

The results of our BLYP-D/TZ2P study on the formation of the G4, A4-N1, A4-N3 and A4-N7 quartets in the point group symmetries C4h, C4 and S4 in the gas phase and in water are summarized in Table 4 (energies), Figs 2, 3, 4, 5 and Tables S1–S4 of the ESM (geometries). The geometries are in both phases almost similar, except for the A4-N3 minima in the S4 symmetry. For the C4h symmetry the hydrogen bond lengths are given for the gas phase and the condensed phase in Fig. 2. In the C4 symmetry we find for all quartets bowl-shaped structures (see Fig. 3). The S4 structures in water are presented in Fig. 4. For G4 and A4-N3 quartets the systems are slightly distorted from planarity and for A4-N1 and A4-N7 we find stacked dimers with cyclic hydrogen bonds in between (for A4-N1 with four N6(H)···N1 and for A4-N7 with four N6(H)···N7 hydrogen bonds as presented in Scheme 1). The number of imaginary frequencies from the vibrational analysis of the structures in different symmetries can be found in Table S6 of the ESM. The hydrogen bond energies for B3LYP have been calculated with the BLYP-D geometries to give a qualitative impression of how B3LYP is unable to reproduce energies for stacked complexes (see Table 4).

Table 4 Hydrogen bond energies (kcal/mol) for DNA quartets in the gas phase and in water
Full size table
Fig. 2
figure 2

The C4h structures in water of G4, A4-N1, A4-N3 and A4-N7 at the BLYP-D/TZ2P level of theory. Hydrogen-bond distances (Å) are given for aqueous solution (gas-phase values in parentheses)

Full size image
Fig. 3
figure 3

Side view and top view of the C4 structures in water of A4-N1, A4-N3 and A4-N7 at the BLYP-D/TZ2P level of theory

Full size image
Fig. 4
figure 4

Structures of the S4 global minima in water of G4, A4-N1, A4-N3 and A4-N7 at the BLYP-D/TZ2P level of theory (for G4 and A4-N3 a side and a top view are shown)

Full size image
Fig. 5
figure 5

Side view of the local S4 minimum in the gas phase of A4-N1 at the BLYP-D/TZ2P level of theory

Full size image

In accordance with previous work of Meyer et al. [18] at the B3LYP level, we find in the gas phase for the G4 quartet that the difference between the C4h and the S4 structures is small, only 0.6 kcal/mol (see Table 4; Fig. 4). Also for the A4-N3 quartet, we find geometrically equivalent structures as in previous work [17] at the B3LYP level. However, we find that all structures are energetically very similar (within 0.4 kcal/mol equal), whereas Meyer et al. [17] found that at the B3LYP level the S4 structure is 8.6 kcal/mol higher in energy than the C4h structure.

For the A4-N1 and the A4-N7 quartets, we find that the S4 structure is more than 16 and 12 kcal/mol, respectively, lower in energy than the other symmetric structures. This is in sharp contrast with previous work [18] at the B3LYP level, where the S4 structure of A4-N1 is only a few kilocalorie per mole lower in energy than the C4 and C4h structures, and for the A4-N7 quartet all symmetric structures are within 0.5 kcal/mol energetically equal. The difference can be ascribed to the difference in the geometries. Our S4 minima are stacked dimers with hydrogen bonds in between, whereas the B3LYP minima [18] look more like a distorted C4h structure. For the A4-N1 quartet we were able to find a local minimum in the S4 symmetry at the BLYP-D/TZ2P level (see Fig. 5), which is very similar to the B3LYP structure. The bond energy of the local S4 minimum amounts only to −31.2 kcal/mol, while the bond energy of the global S4 minimum is −46.0 kcal/mol. To analyze if the B3LYP functional could reproduce this result, we calculated the B3LYP/TZ2P bonding energies with our BLYP-D/TZ2P structures: the bond energy of the stacked dimers with hydrogen bonds in between (Fig. 4) amounts then to −9.7 kcal/mol and of the “distorted C4h” structure (Fig. 5) to −14.4 kcal/mol. Thus, in addition to a general underestimation of the stability of our stacked model systems, the B3LYP functional also yields the wrong energetic ordering of the minima (see also Table 4). This demonstrates once more how the B3LYP functional leads to qualitatively wrong chemical conclusions.

3.3 To stack or not to stack

In aqueous solution, our results show again that for the G4 and the A4-N3 all three symmetries are very close in energy and that it costs only 0.2 kcal/mol in both cases to go from the S4 geometry (as presented in Fig. 4) to the flat C4h symmetric geometry. Such planarization is necessary in order to form stacked systems.

The bonding energy (−33.6 kcal/mol) of the G4 quartet in the C4h symmetry is twice as large as the bonding energy (−16.8 kcal/mol) of the A4-N3 quartet and therefore G4 is, as expected, more stable. The A4-N1 and the A4-N7 quartets bind in the S4 geometry with −33.2 and −27.2 kcal/mol, respectively, which is of the same magnitude as the bonding energy of G4. To go from the nonplanar stacked system in the S4 geometry (see Fig. 4) to the planar C4h geometry, which is needed to form stacks of quartets, it costs 22.3 kcal/mol for the A4-N1 quartet and 11.0 kcal/mol for the A4-N7 in aqueous solution.

The relatively low stability of planar A4-N1 and A4-N7 quartets and their strong tendency to form S4-symmetric nonplanar arrangements might be the decisive factor behind the fact that these A4 quartets are not so often observed experimentally. This also explains why these quartets are found sandwiched between two G4 quartets, which supplies the required stabilization of the planar configuration through stacking interactions.

4 Conclusions

We have shown that the dispersion-corrected density functionals, especially BLYP-D, are very well suited for describing both hydrogen bonded AT and GC pairs and stacked AT and GC dimers, with bond energies and distances that are in excellent agreement with the best ab initio CCSD(T) benchmark data in the literature.

The B3LYP functional fails completely to describe the stacked AT and GC dimers. Furthermore, B3LYP fails to reproduce the strong tendency of A4 quartets to adopt nonplanar structures. Our findings imply that B3LYP should not be used for DNA systems that involve π-stacking interactions.

Our BLYP-D results show that in water G4 is the most strongly bound quartet and easily adopts a planar geometry. The planar geometries of all A4 quartets studied are two to three times less stable than planar G4. In addition, by far the most stable A4 structure is nonplanar. In other words, A4 has, unlike G4, a strong tendency to adopt a geometry that is unsuitable for stacking. This explains the fact that A4 quartets have been experimentally observed in stacks only in between G4 quartets: they are less stable and nonplanar and therefore need the stabilization of planar G4 in order to form larger stacks.