2Department of Chemistry

Oral Roberts University

Tulsa, OK 74171

This paper reports the analysis of adenine spectra using both harmonic and anharmonic approximations to the vibrational frequencies reported in matrix isolation studies. The harmonic approximation procedure consists of the application of a scaled ab initio calculated harmonic force field to predict the frequencies, and infrared intensities, of adenine. Theoretical calculations were made using Hartree-Fock - density functional theory (DFT) B3-LYP/6-31G* and GGA/DNP computational methods. The equilibrium calculated force constants were scaled according to the method of Pulay and compared with the experimentally determined frequencies, and intensities, to assess the accuracy and fit of the theoretical calculation. Good agreement is found except for the in-plane X-H bending or stretching and the out-of-plane X-H bending or wagging modes (X=C,N) which exhibit cubic and quartic anharmonicity respectively. The NH₂ puckering mode, is an out-of plane mode that is poorly modeled by both DFT methods, which we show is most probably due to its quartic-quadratic anharmonicity. In this work, we document the anharmonicity and show that intermode coupling can be estimated using harmonic shift analysis. The goal of this study is to compare various DFT approaches with the aim of determining their limits within the harmonic approximation. We present the method of harmonic shift analysis as a tool for the estimation of mode anharmonicity and for the determination of intermode coupling in the DFT calculation of adenine.

The search for accurate vibrational analysis of nucleobases has engaged the attention of many scientists for several decades. Although, reasonable results can be achieved for heterocyclic molecules, nucleobases have presented many anomalies that have prevented accurate normal coordinate analysis. Uracil presents the possible exception, as there are a number of theoretical studies that have demonstrated reasonable agreement with vibrational data. One can hypothesize that the difficulty in assigning nucleobase vibrational spectra are not simply due to the complexity of the infrared and Raman spectra, but to an aspect of the computational methods that does not account for the spectral features exhibited in nucleobases. These molecules contain interacting low frequency modes with significant anharmonicity, extensive hydrogen bonding, and physical properties that preclude easy spectral measurement in the liquid or gas state. This is particularly the case for exocyclic amino groups, which present a greater difficulty for harmonic analysis than carbonyl groups. The success in modeling uracil spectra may be attributable to the absence of exocyclic amino groups in this molecule. Adenine provides an interesting test case of a molecule that contains no carbonyl groups and one exocyclic amino group. By applying both harmonic and anharmonic analysis to adenine, insight can be gained into the role that these properties play in the vibrational analysis of polynucleotides.

The present study focuses on matrix isolation data for adenine that has been obtained previously. The analysis of purines has received far less attention than that of pyrimidines. Adenine has been studied intensely, by infrared and Raman spectroscopy, in solution, in the solid state, in the gas phase, in an Argon matrix and by neutron inelastic scattering. In spite of numerous vibrational analyses, the harmonic force field and its corresponding normal mode analysis has been reliably assigned only recently using matrix isolated isotopomers. This recent study used a single scale factor, for scaling the Hessian matrix from a B3LYP/6-31G** density functional calculation to fit the matrix isolation data. They found reasonable agreement for some modes, but noted discrepancies in low frequency modes. Baker used a multiple scale factor scheme that is similar to the method presented here. Their multiple scale factor method fits the matrix isolation significantly better, but used qualitative intensity comparisons to assess the accuracy of the intensity predictions. The goal of developing a vibrational model that can be used to interpret the vibrational spectra of purine nucleobases has been stymied, by the discrepancies in modes below 1000 cm^-1. The disagreement in the low frequency modes and the systematic error in the harmonic approximation applied to high frequency modes can be addressed by a model that includes anharmonic coupling. Such an approach is key to further advancement of the vibrational analysis of nucleobases.

This work examines the role of anharmonic analysis via harmonic shift calculation, redundant coordinate Pulay type scaling and refinement (SQM scaling) and matrix isolation measured infrared intensities in application to matrix isolation data. Comparison of two different methods of calculation shows that some of the calculated vibrational modes are outside of the expected error for quantum chemical calculations. It is shown that the discrepancy arises due to a breakdown in the harmonic approximation. If these anharmonic bands are too numerous, then using normal coordinates, as a basis for modeling the molecular vibrations would be difficult if not impossible. Thus, it is important to establish the number of vibrational bands that deviate from harmonic behavior. Here we show that the number is reasonably small and that it is possible to design anharmonic corrections that present a possible route for complete assignment of nucleobases.

In this study, we compare the B3LYP density functional using Gaussian98 and the GGA density functional using DMol3. Both functionals give comparable results and both fail to fit high frequency modes and certain low frequency out-of-plane modes. Using eigenvector projections, we show that these modes are anharmonic. Out of 39 modes, we find 11 modes that have significant anharmonicity. We establish that modes having C-H or N-H stretching character have large cubic and quartic anharmonic terms and that modes with (C, N)-H and (C-N or N-C)-H out-of-plane bending or wagging character have large quartic anharmonic terms. In addition to exhibiting significant mode anharmonicity, the stretching or wagging modes are shown to have intermode coupling which involves a hydrogen atom common to both modes. The conclusion of this type of harmonic shift analysis points in a direction for the development of a comprehensive approach to the interpretation of vibrational spectra of not only nucleobases, but also all other complex molecular systems.

Experimental frequencies and relative intensities were taken from matrix isolation experiments and are summarized in columns 10 & 11 of Table 1. The data span the range from 200 – 3600 cm^-1. A molecular structure of adenine was constructed using crystallographic coordinates obtained from the Cambridge Structural Database (Cambridge Crystallographic Data Centre, Cambridge, UK), for adenine trihydrate (Figure 1.). The optimized ground state geometries were obtained using both the generalized gradient approximation (GGA) of Perdew and Wang as implemented in DMol3 using a double-numeric basis with polarization functions (DNP) and using the B3LYP density functional as implemented in Gaussian98 with a 6-31G* basis set. All calculations were carried out on the SGI/Cray Origin 2000 or the IBM RS/6000 SP supercomputers at the North Carolina Supercomputer Center (NCSC). The geometry optimizations with no symmetry constraints were carried out until the energy difference was less than 10^-6 a.u. on subsequent iterations. Following geometry optimization the Hessian matrix was constructed by finite difference of the analytic gradient for DMol3 and the Hessian was found by analytic second derivative methods as implemented in Gaussian98. Frequency calculation results and calculated infrared intensities from both DFT methods are summarized in columns 2 & 4 of Tables 1 & 2 respectively. Conversion from Cartesian to internal coordinates, automatic generation of the redundant internal coordinates, SQM (Pulay) scaling, least square refinement of scale factors, and decomposition of the potential energy distribution (PED) were carried out using the program FCART01, a major modification of previous software.

Normal coordinate analysis was conducted on both computational models using FCART01 as described earlier. The basic procedure is outlined as follows: The final output and coordinate files were read into FCART01. The information in these output files consists of Cartesian coordinates (including atom type and number sequencing), atomic forces, Hessian matrix, infrared dipole and Raman polarizability intensity derivatives, and vibrational frequencies calculated within the harmonic approximation. The input data were then used by FCART01 to generate scaled frequencies and calculated intensities, including the normal coordinate potential energy distribution (PED). The calculated frequencies can then be further refined by nonlinear least squares on the scale factors to fit entered experimental frequencies and intensities if necessary. The scaled frequencies and relative intensities for both models are presented in columns 5 & 7 of Tables 1 & 2 respectively. The potential energy distribution for the normal coordinates can be found in the Supporting Information along with a detailed description of scaling procedure used and normal coordinate vector images for both methods. All of the experimental frequencies were used in the FCART01 treatment except the 591 and 242 cm^-1 modes. In the 591 and 583 cm^-1 modes appear as a single intensity and hence only one was used. The 242 cm^-1 mode does not match the calculated frequencies. The anharmonicity of low frequency modes is likely to be reasonable as demonstrated in the following.

The geometry optimized coordinates were used as the starting point for the eigenvector projection procedure. The distorted Cartesian coordinate geometries were written as input files, in the appropriate format for calculation of frequencies or energies by either DMol3 or Gaussian98. The procedure for eigenvector projections uses mass-weighted Cartesian eigenvectors, h_j calculated using GGA / DNP and B3LYP / 6-31G* methods, respectively, divided by the square root of the atomic mass, m_j.

(3)

The Cartesian normal coordinate Q_j is used for the eigenvector projections. In practice, these eigenvectors consist of the three columns of the similarity transform matrix; obtained by matrix diagonalization of the mass-weighted force constant matrix.

S^-1FS = L (4)

The similarity transform matrix, S diagonalizes the mass-weighted force constant matrix, F. The eigenvalue matrix L is a diagonal matrix consisting of 3N eigenvalues l_j = 4p²c²n_j² columns of eigenvectors for each of the N atoms in the molecule. Each successive three columns of the S matrix correspond to one eigenvector, the displacement along x, y and z for each of the N atoms in the molecule. In this study, 21 uniformly displaced geometry files for each of the eigenvectors were generated over the range ±0.5 in the normal coordinate units defined by Eqn. 3. The frequencies or energies are then extracted and plotted against the displaced normal coordinate.

The SCF energy was calculated for each of the 21 distorted geometries. Potential energy surface (PES) plots of each eigenvector were then fitted to a fourth order polynomial. In principle, the linear term (column 2 of Table 3) should be zero. The cubic and quartic anharmonic terms are reported in columns 4 and 5 of Table 6. Column 6 of Table 3 shows the calculated harmonic frequency based on the quadratic term using the formula

The quadratic parameter A has units of cm^-1/a.m.u./Ų. (5)

appropriate to the coordinate system as used in Eqn. 3 for a potential energy surface plotted in cm^-1. The frequencies determined by Eqn. 5 agree quite well with the original calculated values obtained from the GGA/DNP as implemented by Dmol3. For those potential energy surfaces that show significant cubic and quartic terms, the eigenvalues and eigenfunctions can be determined using the Numerov-Cooley method for numerical calculation of the wave function. The computational code used for the Numerov-Cooley procedure was adapted from program #407, titled "Systems for the Numerical Solution of the Radial Schroedinger Equation" which was obtained from Quantum Chemistry Program Exchange (QCPE). For PESs that are nearly harmonic the Numerov-Cooley eigenvalues are in good agreement ±2% with the density functional determined frequencies obtained from the DMol3 calculation. The numeric generation of eigenfunctions and eigenvalues for an arbitrary potential energy surface provides an estimate for anharmonic corrections to harmonic eigenvalues. This procedure even gave a well-defined value for the eigenfunctions of mode 1, which has a negative eigenvalue (i.e. an imaginary frequency).

Harmonic Shift Analysis (HSA) is a procedure that consists of calculation of frequencies at the same distorted geometries used in single-mode anharmonicity analysis. In the harmonic approximation, eigenvectors are an orthonormal basis and therefore projection along any eigenvector by definition should not affect the eigenvalues or frequencies of any other harmonic mode. The harmonic assumption can be tested, by performing a vibrational frequency calculation on a geometry distortion along a particular eigenvector. Any mode that has a frequency shift is not orthogonal to the distorted eigenvector and is therefore anharmonically coupled to it. This simple procedure provides a test for the relative strength of intermode coupling or multimode anharmonicity in large or small molecules. This procedure was used to determine the anharmonic multimode coupling and consisted of a vibrational frequency calculation for a series of distorted geometries that corresponded to an eigenvector projection along mode 1 of adenine. The frequencies where then extracted and plotted against the displaced normal coordinate.

There are a number of possible methods for interpreting the HSA output. In this study, we determine an average frequency shift by averaging over the nuclear probability in the distorted mode using the wave function obtained from the Numerov-Cooley method, above. The shifted frequency is fit to a polynomial. The harmonic frequency shift is then calculated by integrating the product of the frequency function and the square of the nuclear wave function in the first vibrational state (v=1) minus that of the ground vibrational state (v=0). Only the even powers of the shifted harmonic frequency will contribute to the integral. For simplicity, we use harmonic wave functions that are fit to the exact wave functions as determined by Numerov-Cooley procedure:

(6)

If the form of the shift is

(7)

The formulae used for the v=0 and v=1 vibrational states are:

(8)

and

, (9)

respectively. The limit l is determined by the extent of the polynomial fit to the frequency shift or potential energy surface used. The harmonic shift terms are then given by:

(10)

The above integrals provide analytic solutions, for the general case of an arbitrary l. These can be integrated analytically to obtain

(11)

and

(12)

For l = infinity the integrals are:

(13)

The calculated harmonic shift depends on the magnitude of l. It is not practical to numerically extend the procedure to l = ¥, since normal modes have a finite possible extent of displacement. Thus, the procedure consists of carrying out the calculation at various values of l to ascertain that the shift from the harmonic frequency is consistent. Ideally, the harmonic shift will converge as calculations are carried out with larger and larger values of l.

The frequencies obtained from the harmonic approximation are presented in Tables 1 and 2 for both GGA / DNP (DMol3) and B3LYP / 6-31G* (Gaussian98) calculations. In column 2 of Tables 1 and 2 the harmonic frequencies are reported for the DMol3 and Gaussian98 calculations, respectively. Column 5 in each table lists respective scaled frequency obtained using FCART01. The GGA / DNP calculated frequencies fit the experimental data remarkably well even before scaling is applied, the average percent difference was found to be 1.54%, with only six modes having a percent difference greater than 2.5%. Scaling the GGA / DNP frequencies did not improve the agreement. The average percent difference increased to 2.19% due to a poorer correspondence between the lower frequency modes (mode 14 and below) and the experimental data. For the B3LYP / 6-31G* calculation, the unscaled frequencies deviated by 2.48% but the scaled frequencies differed by only 0.77%. The B3LYP / 6-31G* unscaled frequencies have a greater deviation in the higher modes, than those of the GGA / DNP method, with almost all modes, above mode 21 having a percent difference greater than 2.5%. Scaling the B3LYP / 6-31G* frequencies improves the fit to the lower modes, with only three modes having a percent difference greater than 2.0%. It should be noted that the scaling factors used on the B3LYP / 6-31G* method had been refined using a larger number of molecules, than those used on the GGA / DNP method. Overall the agreement for both methods is encouraging, however, specific mode classes e.g. (high frequency C-H and N-H stretching modes) all show significantly poorer agreement than the average, suggests that the origin of discrepancies arises from a breakdown in the harmonic approximation.

There are several trends seen with both DFT methods. Generally, the agreement between experimental frequencies and the unscaled frequencies is poor for the C-H and N-H stretching region (3100 - 3700 cm^-1), but is correctable with scaling. Modes involving C-C and C-N stretching (600 - 1700 cm^-1) generally fit well before scaling and only improve slightly upon scaling. However, all of the out-of-plane C-H and N-H wags in the 500 – 1000 cm^-1 region show strong anharmonicity (see Single Mode Anharmonicity). The lower frequency modes (< 600 cm^-1), are fitted poorly in both models with five modes out of ten not even being observed. The DFT assignments of modes above 800 cm^-1 are reasonably reliable, whereas for modes below 800 cm^-1 the correspondence between calculated and observed frequencies is less reliable. Thus, in the absence of isotopomer studies, DFT assignment of modes below 800 cm^-1 should be done with caution.

The comparison shown in Tables 1 and 2 indicate that both DFT calculations have similar shortcomings. They are accurate in the mid-frquency region and show similar frequency and intensity discrepancies with the experimental data. The large IAWPE reported for the GGA / DNP calculation results from a number of modes (e.g. Modes 35, 31, 22, 15 & 11) that are predicted to have significant IR intensity but are hardly observed in the spectrum. The B3LYP / 6-31G* IAWPE value of 39.4 compares favorably with the average planar IAWPE of 22 and nonplanar IAWPE of 72 reported for gas phase IR spectra of heterobicyclics. It can be concluded that both the frequencies and IR intensities of matrix isolation spectra of adenine are reasonably well modeled by the gas phase calculations carried out here. The modes involved are similar in both calculations. The FCART01 scale factors used for scaling the GGA / DNP and the B3LYP / 6-31G* calculations are quite different from each other for C-C/N-H bending and C-C-C-H waging coordinates which suggests that there is a systematic error in the calculation (see supplementary material for the scale factors used in each method). In the following, we show that the discrepancies are due to a failure of the harmonic approximation rather than a failure of the either of DFT methods. The consideration of anharmonic contributions is divided into a single-mode and a intermode coupling or multimode contribution.

Adenine has thirty-nine normal modes of vibration, nominally twenty-seven are in the molecular plane and twelve are out of the molecular plane. Single mode anharmonicity is distinguishable by large cubic and quartic terms, obtained from fitting each modes potential energy surface as described the methods, the results of which are reported in columns 4 and 5 of Table 3. The N-H and C-H stretching vibrations (modes 35 – 39) exhibit significant cubic anharmonicities (column 4 of Table 3). This is important since the unscaled frequencies deviate substantially, for example mode 39, is calculated as 108 cm^-1 and 178.5 cm^-1, higher than the experimental frequencies, for the GGA / DNP and B3LYP / 6-31G* calculations, respectively. Modes 19, 17, 9, 6 and 1 have the largest quartic anharmonicity values. These modes all involve out-of-plane wagging of hydrogen atoms in C-H or N-H bonds. Mode 1 involves N10 - H14 & H15 puckering, mode 6 involves N9 - H12 wagging, mode 9 involves N10 - H14 & H15 twisting, mode 17 involves C8 - H13 wagging and mode 19 involves C2 - H11 wagging. Modes 34, 31, 24, 21 and 20, also exhibit quartic anharmonicities. These modes all involve in-plane bending or stretching motions, of C-C-H or C-N-H internal coordinates. Mode 20 involves N10 - H14 & H15 rocking, Mode 21 involves N9 – C8 - H13 & C8- N9 - H12 bending, Mode 24 involves C8 - H13 bending, Mode 31 involves N7 – C8 - H13 stretching and Mode 34 involves C6 – N10 stretching and N10- H14&H15 scissoring. All of these modes (34, 31, 24, 21, 20, 19, 17, 9, 6 and 1) consist of relatively isolated movements. An analysis of the preceding results indicates that there is a correlation between specific classes of normal modes and the magnitude of there anharmonic terms. This is significant since it suggests that there are a limited number and type of normal modes that need to be examined in order to provide a better understanding of where anharmonicity is important in nucleic acids and other important biological molecules.

Mode 1 is of particular interest for it is the C6 – N10 - H14 & H15 puckering mode that involves concerted out-of-plane motion of H14 and H15 and the smaller motion of N10 in the opposite sense relative to the plane of the adenine ring. We refer to this mode as the NH₂ puckering mode in the following. The NH₂ puckering mode is a classic example of a quartic-quadratic potential surface (Figures 2A and 2B). Both the GGA / DNP and the B3LYP / 6-31G* calculations result in a negative frequency (i.e. imaginary eigenvalue) for this mode. This corresponds to the fact that the mode oscillates about a local maximum rather than a local minimum as shown in Figure 2B. Figure 2A shows the PES of the NH₂ puckering mode with the first five Numerov-Cooley wave functions. The displacements represent the eigenvalues determined by the procedure. The exocyclic NH₂ puckering mode frequency can be estimated from the 0à 1 transition as calculated by the Numerov-Cooley procedure and was found to be 420.0 cm^-1. The amino group of adenine (N10 - H14 & H15) is directly involved in six of the fifteen single mode anharmonicities (1, 9, 20, 34, 37 & 39) and all ten of the multimode anharmonicities (4, 7, 11, 20, 22, 23, 37 & 39). In summary, there is strong correlation between the type of motion for a given normal mode and the magnitude of its anharmonicity. This correlation affects five high frequency stretching modes, five low frequency out-of-plane wagging modes, and five low frequency in-plane bending modes, which have cubic, large quartic and moderate quartic anharmonicities, respectively.

The single-mode anharmonicities given in Table 3 do not report on intermode anharmonic coupling. To ascertain whether there is strong intermode coupling, we implemented the harmonic shift analysis (HSA) procedure as described in the methods section. Within the harmonic approximation, all normal coordinates are orthogonal and therefore, in theory, a distortion along one normal coordinate does not affect frequencies of the other normal coordinates. The observation of a frequency shift for a projected geometry is evidence that a normal mode is coupled anharmonically to other modes. The determination of intermode anharmonic coupling for all of the normal modes in a molecule would require that an eigenvector displacement and a subsequent frequency calculation be performed each mode, but this would be computationally expensive process. Therefore, we choose to perform an eigenvector displacement only on the NH₂ puckering mode, for this mode results in a negative frequency (i.e. imaginary eigenvalue). Moreover, Table 3 indicates that only certain classes of modes are likely to have strong intermode coupling, and furthermore that these motions are relatively isolated. Thus, the HSA method can be applied to a subsystem of a complicated molecule.

The HSA procedure was implemented for the NH₂ pucker mode of adenine. The results are shown in Figures 3A-E for modes 2-39. In Figures 3A-E those modes that are flat show no shift for displacements along the NH₂ puckering mode eigenvector and behave as expected for the harmonic approximation. On the other hand, those modes In Figure 3 that have curvature show a harmonic shift. HSA reveals that there are 11 modes out of a total of 38 that should be considered in intermode anharmonic coupling that involves the -NH₂ puckering mode. The remaining 27 modes have negligible coupling to mode 1. It is not surprising that all of the 11 modes have motions that include the exocyclic -NH₂ group. Three modes stand out as strongly coupled to the NH₂ puckering mode. Mode 9, Figure 3D, which has a very unusual pattern, that being a triple minimum, and modes 37 and 39 Figure 3A, that have very large negative shifts. Modes 4, 7, 11, 20, 22, 23, 32 and 34 are coupled weakly and they exhibit moderate shifts over the range of nuclear coordinate space sampled.

There are two methods for evaluating the relative magnitude of the anharmonicity. The first method is the approach described above in Eqns. 6-12, which involves averaging over the nuclear displacement of mode 1 using the empirically determined Numerov-Cooley wave function as the probability distribution function for the nuclei. The results of this method are presented in Column 4 of Table 4 labeled HSA. Since this HSA correction involves the contribution of only a single mode it is not expected to completely fit the data. However, there is an improvement in the agreement of all of the modes in adenine that are anharmonically coupled to the NH₂ puckering mode. The second method is a fixed geometry calculation, Column 5 of Table 4, which compares the frequencies as calculated at an eigenvector projection corresponding to one of the minima in the PES, Figure 2B. When the -NH₂ group was forced into this configuration there were no longer any negative eigenvalues. The significance of this calculation is that a modest change in environment (e.g. hydrogen bonding to either a nucleotide or water) can produce a similar geometry distortion of the -NH₂ group out of the plane. For example, the X-ray structure of adenine trihydrate indicates that there is an out-of-plane geometry distortion of the -NH₂ group. The frequencies and percent differences from the experimental data are reported in Table 4 for the GGA / DNP, HSA, Fixed Distorted and Scaled GGA / DNP models. Considering even a single mode, the NH₂ puckering mode, HSA provides an overall improvement in the agreement with experiment (column 7, Table 4). The fixed distorted geometry model (column 8, Table 4) also improves the agreement with experiment in cases where there is more than one stable equilibrium conformation as is the case for the NH₂ pucker mode (Figure 2).

The application of the HSA method provides insight into the exocyclic amino group of adenine. The large discrepancy between the calculated and experimental values of the N-H symmetric (mode 37) and asymmetric (mode 39) modes is due in large part to the intermode coupling with the lowest frequency NH₂ puckering mode. The NH₂ scissoring mode (mode 34) has a higher experimental frequency 1633 cm^-1 than that calculated by the GGA/DNP method within the harmonic approximation (1620 cm^-1) consistent with an increase in frequency of +13 cm^-1 due to anharmonic intermode coupling to mode 1 (Table 4). Anharmonic intermode coupling to mode 9 (the amino twist mode) provides a substantial increase in the frequency compared to the harmonic frequency moving this mode closer to the observed pair of modes at 583/592 cm^-1. The frequency of the NH₂ pucker mode itself is estimated to be ~420 cm^-1 based on the single-mode anharmonicity. However, HSA suggests that the frequency is significantly lower and may correspond to the mode at 242 cm^-1 observed in matrix isolation data. The results suggest that the exocyclic amino group of adenine can be treated as a separate subsystem. The method described here provides a conceptually simple framework for approaching the complexity of anharmonic coupling in large molecules.

The HSA approach provides a new means to enhance the computational accuracy of DFT modeling of vibrational spectra. Matrix isolation data present a good starting point since the complexity introduced by hydrogen bonding can be ignored. However, there are number of challenges that remain, even for the isolated adenine molecule. Given that normal modes of vibration often tend to occur in clusters, an experimental approach to the separation of harmonic and anharmonic modes is needed. The combination of infrared and Raman spectroscopy, including Raman polarization and depolarization ratios, provides an experimental handle to aid in assignment of all modes since it will help to determine the true mode ordering in a congested region. Temperature and isotopic labeling are two other important experimental factors that can be varied in order to obtain information on these modes. The Pulay matrix scaling method provides a starting point for fitting of the experimental data. However, an approach that includes anharmonic multimode coupling reveals that low frequency modes are not expected to be well-represented by the harmonic approximation. The scaled quantum mechanical (SQM) approach will never satisfactorily calculate these modes since they are not harmonic. The important result from the present study is that modes can be treated as subsystems that are anharmonically coupled. For example, the asymmetric and symmetric N-H stretching modes of the exocyclic amino group are coupled to the NH₂-pucker mode. Multimode coupling accounts for a significant amount of the discrepancy observed between the calculated values and experimental data for N-H stretching modes. The low frequency pucker mode is not determined at all within the harmonic approximation. By contrast, the approach taken here using the Numerov-Cooley method provides an estimate for its frequency for the first time. However, multimode coupling will be important for this mode and will probably lower the calculated frequency. The future goal is to combine the HSA method with matrix scaling by using the fact that only certain modes are anharmonically coupled while other modes obey the harmonic approximation. The key application for this approach is interpretation of vibrational spectra in DNA, RNA and complexes of polynucleotides with drugs and proteins.

Acknowledgments: SF acknowledges support by North Carolina State University startup funds. The computational work carried out by WC was supported by a Research Opportunity Award supplement to NSF grant MCB-9874895. We gratefully acknowledge a Faculty Award from the North Carolina Supercomputer Center.

Supporting Information Available: The potential energy distribution for the normal coordinates, normal coordinate vector images and scaling factors for both the GGA/DNP and B3LYP/6-31G* methods are provided in the supplemental for this publication. This material is available free of charge via the internet at http://pubs.acs.org

Table 1. Calculated frequencies and intensities of adenine using the GGA/DNP approach. Experimental data are from the matrix isolation experiments of Nowak et al.. Calculated intensities are in units of km/mol.

Table 2. Calculated frequencies and intensities of adenine using the B3LYP/6-31G* approach. Experimental data are from the matrix isolation experiments of Nowak et al.. Calculated intensities are in units of km/mol.

Table 3. Based on GGA / DNP DFT calculation. Wavenumber (n) values in the table were obtained from the quadratic term (column 3) of fitting to the potential energy surface to a quartic polynomial using Eqn. 3. No frequency, (n) is listed in column 6 for mode 1 since it is not possible to calculate a real harmonic frequency with A in Eqn. 4 is negative.