Computational chemistry applied to vibrational spectroscopy: A tool for characterization of nucleic acid bases and some of their 5-substituted derivatives

5.2.1 Scaling procedures

The calculated wavenumbers generally have a large error for many reasons [1, 8, 9], such as anharmonicity, errors in the computed geometry, Fermi resonance, solvent effects, etc., and it is necessary to correct them. For this purpose, the first step after the calculation of the vibrational wavenumbers is the scaling, which helps to carry out an accurate assignment. To improve these computed wavenumbers, two scaling procedures can be used [8, 27].

(a) The first procedure uses a single overall scale factor for the calculated wavenumbers, ν^exp/ν^calc. This is the easiest way and thus is the procedure generally used in the bibliography to scale wavenumbers. To correct the overestimation of the calculated wavenumbers, several authors have reported scale factors for different levels. The most complete set of values has been determined by Scott and Radom [31], with the particularity of using two scale factors, one for the high- and medium-wavenumber vibrations, and another for low wavenumbers.

(b) A remarkable improvement in the accuracy of the scaled wavenumbers is obtained if a linear relationship is established between the calculated and experimental wavenumbers. This procedure called linear scaling equation (LSE) uses a scaling equation to correct the computed wavenumbers of a molecule at a specific level of theory. This scaling equation, obtained previously for a specific molecule (or group of molecules) with experimental gas-phase wavenumbers available (or Ar matrix values), can be used in related molecules and in their derivatives due to the good transferability of their parameters, which enables accurate assignment of the experimental values for many compounds. This procedure, developed by us, represents a compromise between accuracy and simplicity, and the results obtained are accurate enough for the standard today in large molecules [1, 8, 9].

Following this last procedure, the scaling equations for Cy and U molecules collected in Table 2 were obtained at different levels of computation.

5.2.2 Assignment of all the bands in the cytosine molecule

The vibrational bands computed in tautomer C1 of the Cy molecule with B3LYP and two different basis sets are examined in Table 3. The first column refers to the numbers assigned to the calculated vibrations, and they are given in increasing order of wavenumbers. These normal modes of Cy structure appear plotted in Figure 4. The second and the third column list the calculated wavenumbers with the 6-31G(d,p) and 6-311+G(2d,p) basis sets, respectively. The relative IR intensities of the fourth column were obtained by normalizing the computed values to the intensity of the strongest line, no. 28. Small differences appear between the calculated values with both basis sets. The value of the IR intensity helps to carry out the match of theoretical to experimental wavenumbers. Although wavenumber calculations on Cy have been reported with other quantum chemical methods [22, 32], the values shown in Table 3 represent the most accurate today.

The characterization established by B3LYP for each calculated wavenumber is shown in the fifth column. The percent contribution of the different modes to a computed wavenumber appears in parentheses. Contributions lower than 10 % were not considered. The sixth and seventh columns collect the experimental wavenumbers reported in neon and argon [33, 34] matrices. The most detailed study corresponds to an Ar matrix [34], and thus their values were selected as reference.

The scaled wavenumbers of the eight to tenth columns were obtained following the scaling procedures mentioned above. These scaled wavenumbers refer to an isolated molecule, and thus they can be directly compared to the experimental wavenumbers in an Ar matrix. From this match, the experimental bands can be assigned. This scaling was carried out in the calculated wavenumbers of tautomers C1, C2a, C2b, C3a and C3b of Cy. When the scaled wavenumbers of these tautomers are compared with the experimental values reported in an Ar matrix, it was not possible to identify the characteristic bands corresponding to the enol (C2a, C3b) or imino (C3a, C3b) tautomers. Thus, only the keto form C1 appears clearly in gas phase. However, in other studies, both C1 and C2 tautomers have been detected with a small amount of C3 [35].

In addition, with the help of computational methods, it is possible to modify the assignment reported by other authors. Thus, in Cy, the calculated C2=O stretch (scaled at 1,712 cm^–1, ninth column) is predicted as the strongest IR band in accordance with that observed in an Ar matrix spectrum at 1,718 cm^–1 [34]. However, Radchenko et al. [36] pointed out that the most intense IR band in an Ar matrix at 1,620 cm^–1 corresponds to the ν(C=C) stretch. We have corrected all errors in Table 3. For the discussion of the remaining bands, Table 3 is self-explanatory.

5.2.3 Assignment of all the bands in the uracil molecule

Table 4 shows the results obtained in the U molecule with the corresponding scaling equations of Table 2. The first column refers to the notation used for the ring normal modes [8, 27]. The main characterization of these modes appears in the second column. The third and fourth columns list the calculated wavenumbers with the B3LYP method and the 6-31G(d,p) and 6-311+G(2d,p) basis sets, respectively. The column with relative intensities was obtained by dividing the computed values by the intensity of the strongest line. Although wavenumber calculations on uracil have been determined by other authors using MP2 [37] and B3LYP DFT methods [38, 39], the values shown here represent the most accurate today.

Unfortunately, few studies have reported the gas-phase vibrational spectrum of U. The experimental values selected in the seventh to ninth columns were those reported in an argon matrix [40, 41] and in the gas phase [40, 42]. For the determination of the scaling equations, the values reported in the gas phase by Colarusso et al. [42] were selected as reference, since the assignments given there correspond most closely to our own. These wavenumbers appear underlined in Table 4. In the lack of these data, the values obtained in an Ar matrix were used. The scaling equations obtained in this way and listed in Table 2 have a good transferability to uracil derivatives [27]. With this correction (scaling) of the calculated wavenumbers, in general, they are remarkably close to the experimental wavenumbers, and thus they can be used for the assignments.

In the assignment of the uracil modes, Table 4 is self-explanatory, and only a few comments need be made: The N-H stretches appear as pure modes and with strong intensity, and thus they were clearly characterized and assigned. The N-H in-plane bendings are more complex because they appear mixed with other modes. These N-H modes have noticeable contributions to seven calculated vibrations. The description of mode 22 is complex. Mainly, it was characterized as δ(C6-H)+ δ(N-H) with some contribution of ν(C-N), in agreement with the assignment reported by Harsányi et al. [40] and Colarusso et al. [42]. Mode 21 is also complex and it has been characterized as ν(C-N)+ δ(N3-H) in agreement with Harsányi et al. [40] but in contrast to Colarusso et al. [42] and Lés et al. [43].

In the gas phase, the γ(N1-H) out-of-plane bending mode is clearly related to the band at 556 cm^–1 [40] or at 545 cm^–1 [42], in disagreement with Aamouche et al. [37], who related this band to the two calculated wavenumbers corresponding to δ(C2=O) and γ(N1-H) modes. The assignment of the very weak gas-phase band at 672 cm^–1 [40] or at 692 cm^–1 [42] should be changed and related to the γ(N3-H) mode.

The calculated C2=O and C4=O stretches are the strongest IR bands, and they are predicted with similar intensity. The scaled wavenumbers of both modes are at 1,739.1 and 1,706.7 cm^–1 (fourteenth column in Table 4), modes C2=O and C4=O, respectively, in good agreement with the gas-phase results. The C=O out-of-plane bending appears in general mixed with γ(C-H) modes. The description of mode 13 as γ(C5-H)+ γ(C4=O) is in agreement with Colarusso et al. [42] and other authors [38, 43, 44]. The band reported in the gas phase [42] at 717.4 cm^–1 and in an argon matrix at 707.4 and 719 cm^–1 by Ivanov et al. [41] and Szczepaniak et al. [38], respectively, were related here to mode 10, defined as γ(C4=O) + γ(C5-H).

As we show here, with accurate scaling, it is possible to carry out a good assignment and in some cases to correct that reported by other authors.

5.2.4 Scaling in Cy and U molecules: accuracy of the different methods

The accuracy of the two scaling procedures mentioned above on the Cy molecule and the two theoretical methods used is shown in Table 3. The absolute error obtained in the wavenumbers is shown in the last five columns. They were calculated with the experimental values in an Ar matrix. In the last row, the root mean square (rms) error of each level and procedure is determined. With the LSE procedure, the errors obtained in the predicted wavenumbers were very small; the mean deviation was 13 cm^–1 (1 %). The values of these rms errors were very close to those obtained in other related compounds studied by us [9, 45, 46]. This good match of scaled to experimental wavenumbers helps in the assignment and analysis of the experimental fundamental modes.

Similar to the Cy molecule, the absolute errors obtained in the calculated and scaled wavenumbers of the U molecule are shown in Table 4. The largest values correspond to the calculated wavenumbers and they are marked in bold type. However, with the scaling equation they are remarkably reduced as shown in the last column of the Table. The bottom of the table shows the rms error obtained for the calculated and scaled wavenumbers at several computational levels.

Calculations at other levels were also carried out on the U and Cy molecules. The errors obtained in the U molecule and by the two scaling procedures are shown in Table 5. As can be seen, remarkable differences appeared between the HF and B3LYP methods. In both cases, the large rms error in the calculated wavenumbers is remarkably reduced with the scaling, especially by HF. However, the error obtained after scaling by HF continues to be higher than by B3LYP. That is the HF method should not be used for assignment of the experimental wavenumbers.

With the use of a scaling equation, a remarkable reduction of the error is obtained by B3LYP and thus of the risk to a mistake. The best for this purpose is the 6-311+G(2d,p) basis set shown in the last three columns of Table 5. Larger basis sets represent an excessive increase in the computational cost for a very slight improvement. The largest absolute error corresponds to vibration no. 12, a ring bending. Large errors also appear in the vibrations nº 21 and 25, corresponding to C-N and C-C stretches, respectively.

Finally, it can be concluded that the LSE procedure leads to very low errors, and lower than using an overall scaling factor. With B3LYP, the wavenumbers are significantly close to the experimental and much better than by HF. This B3LYP method is also the best among the DFT methods used, and it is the recommended method for carrying out an assignment of the vibrational bands.

5.3.1 Tautomerism in nucleic acid bases

The nucleic acid bases can undergo keto-enol tautomerism. Much of the interest of the tautomerism is due to the fact that tautomers induce alterations in the normal base pairing, leading to the possibility of spontaneous mutations in the DNA or RNA helices. On the basis of geometry, ionization potential and dipole moment, it is virtually impossible to decide which arrangement is actually the most stable one. However, computational methods applied to the vibrational spectra may help in deciding the most stable arrangement (tautomer). The absorption bands due to C=O, NH and OH groups give the most straightforward information about the tautomeric forms present, because they correspond to the characteristic, well-localized vibrations of the functional groups directly involved in tautomeric changes. Thus, special attention was paid on these modes.

All the tautomers in the Cy molecule appear with relative energies much lower than their counterparts in the U molecule. Two features can explain it: (i) the negative charge on N1 atom is lower in Cy than in U, with longer N1-H bond length. (ii) The negative charge on the O2 atom is higher in Cy than in U, with longer C=O bond length. Both features favor tautomerism in Cy molecule. Thus, there is a special interest to study this molecule and its derivatives [15, 16, 47].

The Cy molecule has been studied intensively by IR spectroscopy in the free monomeric form (in low-temperature inert matrices [34, 36]) and in the crystalline phase at room [48] or at low temperatures [49]. The different aggregate states can contain two tautomeric forms [33, 50]: the gas phase with both the enol and the keto forms, whereas in the crystalline state and polar solvents only the keto form has been observed [34]. For its interpretation, the theoretical vibrational spectrum of Cy has been predicted at different levels of approximation [34, 51].

5.3.2 Tautomerism in 5-bromocytosine

Similar to the Cy molecule, 5-BrCy can exist in various tautomeric forms. However, this molecule has been much less studied than Cy, and thus it is discussed here. The computational methods have been applied to vibrational spectroscopy to simulate (scaled) the IR spectrum of the different tautomers of 5-BrCy. In the isolated state, the enol form C2b of 5-BrCy is the most stable one. The next most stable tautomer is the enol form C2a, 3.18 kJ/mol above C2b in 5-BrCy, and 12.96 kJ/mol above C2b in the Cy molecule with the MP2 method. The bromine atom in position 5 favors the tautomerism with much lower relative energies (about 66 %) in 5-BrCy than in Cy.

In the isolated state, the harmonic vibrational bands computed in the 5-BrCy ring are shown in Table 6. The second column lists the calculated wavenumbers with the 6-31G(d,p) basis set in tautomer C1. The third column collects their relative IR intensities (A) in percent. They were obtained by normalizing the computed value to the intensity of the strongest band. The assignment with the calculated percent potential energy distribution (PED) of the different modes for each vibration appears in the fourth column. Contributions lower than 10 % were not considered. The scaled wavenumbers in the fifth column correspond to tautomer C1 and those in the ninth column to tautomer C2b. With these data obtained theoretically, the experimental spectrum in an Ar matrix reported by Jaworski et al. [52] can be now analyzed in detail. For this task, the scaled values were directly compared to the experimental values, and from this comparison it was possible to separate those corresponding to tautomer C1 (sixth column), from those to tautomer C2b (tenth column). In Table 6, only the experimental bands with high intensity are shown, and those selected in the comparison with the scaled values are underlined.

Through the comparison of these spectra, we have identified, separated and assigned the bands corresponding to the different tautomers, which in this case are only the amino-oxo form (tautomer C1) and the amino-hydroxy form C2b which is the most stable form in accordance to the computations, Figure 5. Due to many bands between both conformers appear very close in frequency, the IR intensity criterion was followed in these cases. Thus, the assignment reported previously by Jaworski et al. [52] has been improved by using scaling, and all bands accurately assigned to both tautomers. Bands corresponding to tautomer C3b, simulated theoretically, were not found in the experimental spectrum.

Figure 5:

Assignment to tautomers C1 and C2b of the experimental IR spectrum in argon matrix of 5-bromocytosine.

5.4.1 5-Chlorocytosine

The crystal unit cell of 5-ClCy was simulated at the B3LYP/6-31G(d,p) level by a tetramer form, Figure 7. Its simulated spectrum appears closer to the experimental spectrum in the solid state than the simulated spectrum of the monomer form. Figure 8 only shows the comparison between the experimental (IR) and that scaled with the tetramer form. Similar agreement is observed with the Raman spectrum.

Figure 7:

Simulation of the crystal unit cell of 5-chlorocytosine by a tetramer form at the B3LYP/6-31G(d,p) level.

Figure 8:

Experimental and theoretical IR spectra of 5-chlorocytosine in the 4,000–200 cm^–1 range. The predicted scaled spectrum was obtained using the LSE procedure in the tetramer form at the B3LYP/6-31G(d,p) level.

The scaled and experimental wavenumbers in the monomer and tetramer forms are collected in Table 7, together with the corresponding assignment. They are the most accurate today for this molecule. The scaling equation ν^scaled = 4.8 + 0.9671·ν^calculated deduced from the Cy molecule (Table 3) at the B3LYP/6-311+G(2d,p) level was used for the monomer, while ν^scal = 34.6 + 0.9447·ν^cal at the B3LYP/6-31G(d,p) level was used for the tetramer. Four wavenumbers appear in the tetramer form for each characterized vibration corresponding to its four molecules. Of these four wavenumbers, those determined in the molecules II and III (see Figure 7) are underlined. These wavenumbers appear affected by the internal H-bonds of the crystal, and they represent better the solid state. Among these wavenumbers, that with the highest IR intensity appears in bold type, while that with the highest Raman intensity is shown in italic type. The main characterization carried out, and corresponding to the tetramer form, is shown in the last column.

The absence of a ν(OH) band in the 3,500–3,700 cm^–1 range and the appearance of ν(C=O) modes as strong bands indicate that in the solid state this compound is only in the keto form. Another feature is that the majority of the fundamental modes absorb in relatively narrow wavenumber ranges. We have changed the assignments of several bands reported by other authors [20] according to our scaled values. A good match between the scaled wavenumbers in the tetramer to the experimental ones was found, and much better than that with the monomer form. That is the simulation of the crystal unit cell is appropriate for an accurate assignment of the experimental bands in the solid state.

Figure 9:

Simulation by a tetramer form of two possible crystal unit cells of 5-bromocytosine with three tautomeric forms.

5.4.2 5-Bromocytosine

X-ray data of 5-BrCy has not yet been reported, and thus the tautomer present in the solid state is not known, but it is possible to resolve this problem by simulation of the crystal unit cell of all the possible tautomers using as reference the structure reported in related molecules, which in the present case is 5-ClCy. The crystal data of this unit cell (the tetramer form) was the starting point used for the simulation of the solid state (Figure 9).

In the isolated,state the molecule of 5-BrCy is calculated to be full planar, and the small deviations from planarity observed in crystals of Cy [53] are due to the intermolecular H-bonds. Our simulations in the tetramer form also show these deviations, in general with values slightly higher than in the crystal. Due to the orientation of the bromine atom in the tetramer form, and its almost zero interaction with other atoms, its bond length is almost the same as that in the isolated state. The simulations were carried out at the B3LYP/6-31G(d,p) level [54]. The tetramer form with tautomer C1 is 30.5 kJ/mol more stable than with tautomers C2b-C3b. The simulated spectra of these tetramer forms appear closer to the experimental spectrum in the solid state than the simulated spectrum of the monomer form. As an example of this concordance, the experimental (IR) and the scaled spectrum of this molecule are compared in Figure 10, for simplicity only in the dimer form and in the 3,700–2,700 cm^–1 range. Similar agreement is observed in the other regions, as well as in the Raman spectrum.

Figure 10:

Experimental and theoretical scaled Raman spectra in the dimer form of 5-bromocytosine in the 3,650–2,700 cm^–1 range and in the C1 and C2b tautomers using the scale equation ν^scaled = 34.6 + 0.9447 ν^calculated

The scaled values of the ring vibrations in two of the tetramer forms are collected in Table 8. Four wavenumbers appear for each vibration corresponding to the four 5-BrCy molecules of the tetramer. Of these four wavenumbers, that with the highest IR intensity is shown in bold type, while that with the highest Raman intensity is shown in italic type. These scaled tetramer values can be directly compared to the experimental IR and Raman data reported in the solid state [20], [54]. We have underlined the experimental bands selected for this comparison. The absolute error (Δ) between the observed and the scaled wavenumbers is lower in tetramer C1 than in C2b-C3b. I.e. The spectrum of crystalline solid 5-BrCy shows that this molecule exists in the crystal in only the amino-oxo form (C1).

The preference for only one tautomeric form in the crystal than in the isolated state is a clear indication of the importance of intermolecular interactions, in particular H-bonding, to determine the structure of the condensed phase. Because molecules in the crystal are involved in hydrogen bonding and in other intermolecular interactions, it is expected that these H-bonds significantly change the spectrum obtained in the solid state compared to that in the gas phase (or in Ar/Ne matrices). The bands corresponding to the stretching and bending vibrations of the N1-H, C=O and NH₂ groups change more drastically than those related to the ring or to the CH group. Thus, the calculations of the tetramer form help to interpret the spectra in the crystalline solid.