Thermochemistry of Guanine Tautomers Re-Examined by Means of High-Level CCSD(T) Composite Ab Initio Methods*

All DNA bases exist in several tautomeric forms and DNA tautomerization is an important chemical process that can alter base pairing. Over the past two decades, the gas-phase energetic and spectroscopic properties of the guanine tautomers have been extensively studied both experimentally and theoretically.^[1–23] Nevertheless, as the most energetically stable guanine tautomers are separated by energy differences of only a few kilojoules per mole, it is of interest to use highly accurate ab initio methods for obtaining these energy separations. Hanus et al.^[5] found that as many as four guanine tautomers lie within a narrow energy difference of ~5.5 kJ mol⁻¹. In particular, they found that at the CCSD(T)/aug-cc-pVDZ level of theory, the energies of the keto and enol guanine tautomers are: 2.9 (2), 3.8 (3), and 5.5 (3′) kJ mol⁻¹ relative to tautomer 1 (structures are shown in Fig. 1). In contrast, Plekan et al.^[17] found that at the CCSD(T)/6–311+G** level of theory, the enol tautomer 3 is the second most stable tautomer after 1. In particular, they obtained energies of 0.9 (3), 1.2 (2), and 3.9 (3′) kJ mol⁻¹ relative to tautomer 1. They also showed that the tautomer order remains the same when using the following additive basis-set correction: CCSD(T)/6–311+G** + CCSD/cc-pVTZ – CCSD/6–311+G**. In particular, they obtained energies of: 1.7 (3), 1.9 (2), and 3.4 (3′) kJ mol⁻¹ relative to tautomer 1. Nevertheless, as these levels of theory are far from the CCSD(T)/CBS limit (coupled cluster with single, double, and quasiperturbative triple excitations at the complete basis set limit), it is of interest to obtain the energy separations closer to the CCSD(T)/CBS limit.

Fig. 1.  Approximated CCSD(T)/CBS relative enthalpies at 298 K (ΔH₂₉₈) for the guanine tautomers (in kJ mol⁻¹). The CCSD(T)/CBS energy is approximated by the formula HF/jul-cc-pV{D,T}Z + CCSD/jul-cc-pVTZ + (T)/jul-cc-pVDZ (see below). Atomic colour scheme: H, white; C, grey; N, blue; O, red.

In the present work, we use the high-level composite W1-F12 protocol to determine the energy separation and ordering of the low-lying guanine tautomers. W1-F12 theory approximates the all-electron, relativistic CCSD(T)/CBS energy and has been found to produce thermochemical properties such as heats of formation and reaction energies with mean-absolute deviations (MADs) from highly accurate experimental or theoretical data below the threshold of chemical accuracy (i.e. with MADs below 1 kcal mol⁻¹ (1 kcal mol⁻¹ = 4.184 kJ mol⁻¹)).^[24–29] It should be noted, however, that because the low-lying guanine tautomers are chemically similar to each other, W1-F12 theory should yield even better performance for their relative energies owing to a high degree of error cancellation between reactants and products.^[26,30–47] We find that at the all-electron, relativistic CCSD(T)/CBS limit, the four lowest-energy guanine tautomers are separated by less than 1 kcal mol⁻¹. This work demonstrates that in order to reliably predict the order and energy separation between the low-lying guanine tautomers, highly accurate composite ab initio methods have to be used and that lower-level composite approaches may not be sufficiently accurate. For example, the considered Gn composite procedures (G4, G4(MP2), G4(MP2)-6X, G3, G3B3, G3(MP2), and G3(MP2)B3) predict that the enol tautomer (3) is more stable than the keto tautomer (2).

In order to obtain accurate energies for the guanine tautomers, calculations were carried out using the high-level ab initio W1-F12 thermochemical protocol. W1-F12 theory combines explicitly correlated F12 techniques^[48] with basis-set extrapolations in order to approximate the all-electron, relativistic CCSD(T)/CBS energy. The computational details of W1-F12 theory have been specified elsewhere.^[24] In short, the Hartree–Fock (HF) and valence CCSD-F12 correlation components are extrapolated from the VDZ-F12 and VTZ-F12 basis sets, where VnZ-F12 denotes the cc-pVnZ-F12 basis sets of Peterson et al.^[49] The complementary auxiliary basis set (CABS) singles correction is included in the self-consistent field (SCF) energy.^[50–52] An optimal extrapolation exponent of 3.38 was used in the CCSD-F12/cc-pV{D,T}Z-F12 extrapolation as recommended in reference [24]. The (T) valence correlation energy is obtained from the original W1 theory,^[53] i.e. it is extrapolated from the jul-cc-pVDZ and jul-cc-pVTZ basis sets.^[54–56] The diagonal, fixed-amplitude 3C(FIX) ansatz,^[51,57–59] and the CCSD-F12b approximation are used in all of the explicitly correlated coupled cluster calculations.^[52,60] The CCSD inner-shell contribution is calculated with the core-valence weighted correlation-consistent cc-pwCVTZ basis set of Peterson and Dunning,^[61] while the (T) inner-shell contribution is calculated with the cc-pwCVTZ(no f) basis set (where cc-pwCVTZ(no f) indicates the cc-pwCVTZ basis set without the f functions). The scalar relativistic contribution (in the second-order Douglas–Kroll–Hess approximation)^[62,63] is obtained as the difference between non-relativistic CCSD(T)/jul-cc-pVDZ and relativistic CCSD(T)/jul-cc-pVDZ-DK calculations.^[64] The diagonal Born–Oppenheimer correction (DBOC) is calculated at the HF/cc-pVDZ level of theory using the CFOUR program suite.^[65] The energy calculations involved in W1-F12 theory are carried out with the Molpro 2016.1 program.^[66,67] The relativistic, all-electron, DBOC-inclusive W1-F12 CCSD(T)/CBS tautomerization energies are converted to tautomerization enthalpies at 298 K using zero-point vibrational energies (ZPVEs) and enthalpic corrections calculated at the B3LYP/jul-cc-pVTZ level of theory.^[68–70] The harmonic ZPVEs were scaled as recommended in references [25] and [71].

In addition, the performance of more approximate Gaussian-n^[72] and CBS-type^[73] composite thermochemical procedures is assessed.^[74,75] We consider the following composite procedures: G4,^[76] G4(MP2),^[77] G4(MP2)-6X,^[78] G3,^[79] G3(MP2),^[80] G3B3,^[81] G3(MP2)B3,^[81] CBS-QB3,^[82] and CBS-APNO.^[83] The calculations for these composite methods were performed using the Gaussian 16 program suite.^[84]

Fig. 1 gives the relative enthalpies at 298 K for the guanine tautomers. The energies are calculated at the CCSD(T)/CBS level of theory, using the following formula HF/jul-cc-pV{D,T}Z + CCSD/jul-cc-pVTZ + (T)/jul-cc-pVDZ, which has been found to give excellent performance relative to W1-F12 theory (see below). For the enol tautomers, we have considered two possible conformers, one in which the O–H bond is pointing towards the pyrimidine ring (labelled n, Fig. 1) and the other in which the O–H bond is pointing towards the imidazole ring (labelled n′, Fig. 1). Inspection of the relative tautomerization energies in Fig. 1 reveals that the energy separations between the low-lying tautomers are spread over a narrow range; namely, the four most stable tautomers are found within a small energy window of less than 1.0 kcal mol⁻¹. It is therefore clear that in order to reliably and accurately predict the order and energy separation between these tautomers, we have to use a high-level composite ab initio method such as W1-F12 theory. However, it should be noted that because tautomers 2 and 3 are separated by merely 0.3 kJ mol⁻¹, the scaled ZPVE component of W1-F12 theory^{[24,25,37,71,85]} may not be sufficiently accurate to unambiguously determine the energy separation between these two tautomers.

Table 1 gives an overview of the basis set convergence of the HF component of the guanine tautomerization energies relative to the HF/CBS limit from W1-F12 theory. As expected, the jul-cc-pVDZ basis set results in very large deviations from the CBS values, with deviations of up to 1.8 kJ mol⁻¹. The jul-cc-pVTZ basis set performs well with deviations on the order of 0.1 kJ mol⁻¹ and a maximum deviation of 0.2 kJ mol⁻¹. The jul-cc-pV{D,T}Z basis set extrapolation results in deviations on the order of 0.05 kJ mol⁻¹ and a maximum deviation of 0.08 kJ mol⁻¹. The jul-cc-pVQZ basis set gives tautomerization energies that are practically identical to those at the HF/CBS values from W1-F12 theory (Table 1).

Table 1.  Overview of basis set convergence of the HF and valence CCSD and (T) correlation contributions to the guanine tautomerization energies (energies are relative to tautomer 1, kJ mol⁻¹)

Similarly to the HF energy, calculating the CCSD correlation energy with the jul-cc-pVDZ basis set results in very large deviations of up to 1.8 kJ mol⁻¹ from W1-F12 CBS limit values. The jul-cc-pVTZ basis set performs much better, with deviations on the order of 0.1 kJ mol⁻¹. However, the jul-cc-pV{D,T}Z basis set extrapolation results in significantly larger deviations on the order of 0.8 kJ mol⁻¹. Thus, the jul-cc-pVDZ basis set should not be used as a stand-alone basis set or in basis set extrapolations for the calculation of the the CCSD correlation component of the guanine tautomerization energies.

Calculating the (T) correlation energy with the jul-cc-pVDZ basis set results in fairly good performance, with deviations of 0.0–0.2 kJ mol⁻¹ from W1-F12 CBS limit values.

Table 2 gives the component breakdown of the W1-F12 tautomerization energies for the five energetically most stable guanine tautomers. The HF/CBS level of theory predicts the wrong order of the three most stable tautomers with relative energies of −0.34 (3), −0.20 (2), 0.00 (1), 2.77 (3′), and 16.31 (4) kJ mol⁻¹. The CCSD/CBS level of theory gives the correct order of all tautomers with relative energies of 0.00 (1), 1.28 (2), 2.56 (3), 4.13 (3′), and 16.40 (4) kJ mol⁻¹. The (T) correlation component affects the relative energies by significant amounts of up to 1.52 kJ mol⁻¹. At the valence CCSD(T)/CBS level of theory, we obtain relative energies of 0.00 (1), 2.80 (2), 3.21 (3), 4.43 (3′), and 15.18 (4) kJ mol⁻¹. The CCSD(T) core-valence correction still affects the relative energies by up to 0.17 kJ mol⁻¹, whereas the scalar relativistic and DBOC corrections do not affect the tautomerization energies by a chemically significant amount; namely, they are ≤0.03 kJ mol⁻¹.

Table 2.  Component breakdown of the W1-F12 guanine tautomerization energies (relative to tautomer 1, kJ mol⁻¹)

We proceed to evaluate the performance of computationally economical procedures for the relative energies of the guanine tautomers. Table 3 gives the root-mean-square deviations (RMSDs), MADs, and mean-signed deviations (MSDs) for several standard and composite ab initio procedures for the five most stable guanine tautomers relative to W1-F12 reference values. We note that as the four low-lying tautomers are separated by energy differences as small as 0.3 kJ mol⁻¹, it is desirable to identify computationally economical methods with RMSDs below this threshold.

Table 3.  Overall error statistics for standard and composite ab initio procedures for the tautomerization energies of the five most stable guanine tautomers (relative to reference values from W1-F12 theory, in kJ mol⁻¹)

G3 theory, in which the geometries and ZPVEs are obtained at the MP2 and HF levels respectively, results in a large RMSD of 0.93 kJ mol⁻¹. Using G3B3, in which the geometries and ZPVEs are obtained with B3LYP, results in a slightly lower RMSD of 0.73 kJ mol⁻¹. In contrast, the more recent G4 theory results in better performance at a similar computational cost, with an RMSD of 0.54 kJ mol⁻¹. Both G3(MP2) and G3(MP2)B3 result in relatively poor performance with RMSDs >1 kJ mol⁻¹, whereas G4(MP2) and G4(MP2)-6X give essentially the same performance, with an RMSD of 0.7 kJ mol⁻¹. We find that all the considered Gn procedures predict that the enol tautomer (3) is energetically more stable than the keto tautomer (2) by amounts ranging from 0.36 (G4) to 1.28 (G3(MP2)) kJ mol⁻¹ (see Table S1 in the Supplementary Material for the tautomerization energies of the Gn-type methods).

Let us move to the CBS-type composite methods. CBS-QB3 results in excellent performance with an RMSD of 0.31 kJ mol⁻¹. However, the older CBS-APNO method results in a significant deterioration of performance (RMSD = 0.81 kJ mol⁻¹). Both CBS-QB3 and CBS-APNO predict the correct energy order of the low-lying tautomers.

A cost-effective approach for approximating the CCSD(T)/CBS energy is using an additivity-based scheme in which this energy is estimated from the CCSD(T)/jul-cc-pVDZ energy and an MP2-based basis-set correction (ΔMP2 = MP2/jul-cc-pV{T,Q}Z – MP2/jul-cc-pVDZ, where MP2/jul-cc-pV{T,Q}Z energy is extrapolated to the basis-set limit with an extrapolation exponent of 3).^[86] This additivity scheme has been found to be an efficient way for approximating non-covalent interactions (e.g. hydrogen-bonding and dispersion interactions),^[87–93] as well as reaction energies^[45,94,95] at the CCSD(T)/CBS level. Nevertheless, this simple additivity scheme does not perform particularly well for the guanine tautomers, with an overall RMSD of 0.69 kJ mol⁻¹. Furthermore, it should be noted that the CCSD(T)/jul-cc-pVDZ level of theory alone results in much better performance, with an RMSD of 0.40 kJ mol⁻¹.

The results presented in the previous section indicate that a cost-effective approach to improve on the CCSD(T)/jul-cc-pVDZ result is to extrapolate the SCF energy from the jul-cc-pV{D,T}Z basis set pair and calculate the CCSD and (T) correlation energies with the jul-cc-pVTZ and jul-cc-pVDZ basis sets respectively. Using an extrapolation exponent of 5.0 in the SCF extrapolation results in an RMSD of 0.23 kJ mol⁻¹. Thus, this approach outperforms all the composite and standard methods considered so far, including the CBS-QB3 method (Table 3). We note that using an optimal extrapolation exponent of α = 3.5 in the SCF extrapolation cuts this RMSD by ~50 % and gives an RMSD of 0.11 kJ mol⁻¹. In the next section, we use this method to obtain reference values for evaluating the performance of approximated theoretical procedures for all the 14 guanine tautomers considered in this work (these reference values are given in Fig. 1).

Table 4 gives the error statistics for several standard and composite ab initio procedures for the 14 guanine tautomers relative to approximated CCSD(T)/CBS reference values. The reference HF/jul-cc-pV{D,T}Z + CCSD/jul-cc-pVTZ + (T)/jul-cc-pVDZ values are given in Fig. 1 and the tautomerization energies for all methods are given in Table S1 (Supplementary Material). Before proceeding to a detailed discussion of the performance of the various methods, we note that with the composite methods that do not use B3LYP geometries, i.e. G3, G3(MP2), and CBS-APNO, we were unable to optimize the geometry of the highly energetic tautomer 8 (Fig. 1). Thus, in these cases, this tautomer is not included in the error statistics. Having said that, exclusion of one tautomer is not expected to drastically affect the overall error statistics in these cases.

Table 4.  Overall error statistics for standard and composite ab initio procedures for the tautomerization energies of the 14 guanine tautomers in Fig. 1 (relative to approximated CCSD(T)/CBS reference values, in kJ mol⁻¹)

Inspection of the RMSDs obtained for the 14 guanine tautomers (Table 4) reveals that they are fairly similar to those obtained for the five most stable tautomers (Table 3). In particular, the RMSDs for G3, G3B3, and G4(MP2) remain practically unchanged, and the RMSDs for G4, G4(MP2)-6X, and CCSD(T)+ΔMP2 differ by less than 0.2 kJ mol⁻¹ from each other. The largest changes between the two sets of RMSDs are observed for the CCSD(T)/jul-cc-pVDZ and CBS-QB3 methods, where the performance for the larger set deteriorates by 0.65 and 0.93 kJ mol⁻¹, respectively. The methods that emerge as best performers for the entire set of tautomers are G4(MP2)-6X and CCSD(T)+ΔMP2 with nearly identical RMSDs of 0.51 and 0.52 kJ mol⁻¹, respectively. Nevertheless, we note that the poorest performers (CCSD(T)/jul-cc-pVDZ, CBS-QB3, and CBS-APNO) still attain RMSDs of 1.05–1.38 kJ mol⁻¹.

We obtained the gas-phase tautomerization energies at the W1-F12 level of theory for the five most stable guanine tautomers. We find that as many as four tautomers populate a small energetic window of less than 1.0 kcal mol⁻¹. The W1-F12 relative enthalpies at 298 K for the five lowest-lying guanine tautomers are: 0.00 (1), 2.37 (2), 2.63 (3), 4.03 (3′), and 14.31 (4) kJ mol⁻¹. We use these highly accurate W1-F12 tautomerization energies to benchmark the performance of a wide range of lower-level standard and composite ab initio procedures. The best performing procedures with RMSDs ≤0.5 kJ mol⁻¹ are: 0.54 (G4), 0.40 (CCSD(T)/jul-cc-pVDZ), 0.31 (CBS-QB3), and 0.23 (approximately CCSD(T)/CBS) kJ mol⁻¹. Here, approximated CCSD(T)/CBS is calculated via the additive formula HF/jul-cc-pV{D,T}Z + CCSD/jul-cc-pVTZ + (T)/jul-cc-pVDZ, where the HF energy is extrapolated with an exponent of 5.0. Using instead an effective extrapolation exponent of α = 3.5 in the SCF extrapolation results in an RMSD of merely 0.11 kJ mol⁻¹. We use this approximate CCSD(T)/CBS method for calculating the tautomerization energies for all 14 tautomers. The resulting tautomerization energies are: 0.00 (1), 2.20 (2), 2.51 (3), 4.06 (3′), 14.30 (4), 25.65 (5), 43.78 (4′), 53.50 (6′), 61.58 (6), 77.37 (7), 82.52 (8′), 86.02 (9), 100.70 (10), and 121.01 (8, Fig. 1) kJ mol⁻¹. Finally, we evaluate the performance of the standard and composite methods relative to these reference values. For the entire set of tautomers, the best performing procedures emerge as (RMSDs are given in parenthesis): G4(MP2)-6X (0.51), CCSD(T)+ΔMP2/CBS (0.52), and G4(MP2) (0.64 kJ mol⁻¹).

Tautomerization energies for the 14 guanine tautomers calculated with the standard and composite CCSD(T) methods considered in this work (Table S1) and B3LYP/jul-cc-pVTZ optimized geometries for all guanine tautomers (Table S2) are available on the Journal’s website.

The author declares no conflicts of interest.